Smoothed Particle Hydrodynamics (SPH) : An Overview and Recent Developments

Arch Comput Methods Eng (2010) 17: 2576
DOI 10.1007/s11831-010-9040-7
Smoothed Particle Hydrodynamics (SPH): an Overview

and Recent Developments
M.B. Liu G.R. Liu
Received: 1 August 2009 / Accepted: 1 August 2009 / Published online: 13 February 2010
CIMNE, Barcelona, Spain 2010
Abstract Smoothed particle hydrodynamics (SPH) is a

meshfree particle method based on Lagrangian formulation,
and has been widely applied to different areas in engineering and science. This paper presents an overview on the
SPH method and its recent developments, including (1) the
need for meshfree particle methods, and advantages of SPH,
(2) approximation schemes of the conventional SPH method
and numerical techniques for deriving SPH formulations
for partial differential equations such as the Navier-Stokes
(N-S) equations, (3) the role of the smoothing kernel functions and a general approach to construct smoothing kernel
functions, (4) kernel and particle consistency for the SPH
method, and approaches for restoring particle consistency,
(5) several important numerical aspects, and (6) some recent
applications of SPH. The paper ends with some concluding
remarks.
1 Introduction
1.1 Traditional Grid Based Numerical Methods
Computer simulation has increasingly become a more and
more important tool for solving practical and complicated
G.R. Liu is SMA Fellow, Singapore-MIT Alliance.

M.B. Liu ()
Key Laboratory for Hydrodynamics and Ocean Engineering,
Institute of Mechanics, Chinese Academy of Sciences,
#15, Bei Si Huan Xi Road, Beijing 100190, China
e-mail: liumoubin@imech.ac.cn
G.R. Liu
Centre for Advanced Computations in Engineering Science
(ACES), Department of Mechanical Engineering, National
University of Singapore, 10 Kent Ridge Crescent, Singapore
119260, Singapore
problems in engineering and science. It plays a valuable role

in providing tests and examinations for theories, offering insights to complex physics, and assisting in the interpretation
and even the discovery of new phenomena. Grid or mesh
based numerical methods such as the finite difference methods (FDM), finite volume methods (FVM) and the finite element methods (FEM) have been widely applied to various
areas of computational fluid dynamics (CFD) and computational solid mechanics (CSM). These methods are very
useful to solve differential or partial differential equations
(PDEs) that govern the concerned physical phenomena. For
centuries, the FDM has been used as a major tool for solving partial differential equations defined in problem domains
with simple geometries. For decades, the FVM dominates in
solving fluid flow problems and FEM plays an essential role
for solid mechanics problems with complex geometry [13].
One notable feature of the grid based numerical models is to
divide a continuum domain into discrete small subdomains,
via a process termed as discretization or meshing. The individual grid points (or nodes) are connected together in a
pre-defined manner by a topological map, which is termed
as a mesh (or grid). The meshing results in elements in FEM,
cells in FVM, and grids in FDM. A mesh or grid system
consisting of nodes, and cells or elements must be defined
to provide the relationship between the nodes before the approximation process for the differential or partial differential
equations. Based on a properly pre-defined mesh, the governing equations can be converted to a set of algebraic equations with nodal unknowns for the field variables. So far the
grid based numerical models have achieved remarkably, and
they are currently the dominant methods in numerical simulations for solving practical problems in engineering and
science [15].
Despite the great success, grid based numerical methods
suffer from difficulties in some aspects, which limit their
26
applications in many types of complicated problems. The

major difficulties are resulted from the use of mesh, which
should always ensure that the numerical compatibility condition is the same as the physical compatibility condition for
a continuum. Hence, the use of grid/mesh can lead to various difficulties in dealing with problems with free surface,
deformable boundary, moving interface, and extremely large
deformation and crack propagation. Moreover, for problems
with complicated geometry, the generation of a quality mesh
has become a difficult, time-consuming and costly process.
In grid based numerical methods, mesh generation for the
problem domain is a prerequisite for the numerical simulations. For Eulerian grid methods such as FDM constructing
a regular grid for irregular or complex geometry has never
been an easy task, and usually requires additional complicated mathematical transformation that can be even more
expensive than solving the problem itself. Determining the
precise locations of the inhomogeneities, free surfaces, deformable boundaries and moving interfaces within the frame
of the fixed Eulerian grid is also a formidable task. The
Eulerian methods are also not well suited to problems that
need monitoring the material properties in fixed volumes,
e.g. particulate flows [1, 2, 6]. For the Lagrangian grid methods like FEM, mesh generation is necessary for the solids
and structures, and usually occupies a significant portion of
the computational effort. Treatment of extremely large deformation is an important issue in a Lagrangian grid based
method. It usually requires special techniques like rezoning.
Mesh rezoning, however, is tedious and time-consuming,
and may introduce additional inaccuracy into the solution
[3, 7].
The difficulties and limitations of the grid based methods
are especially evident when simulating hydrodynamic phenomena such as explosion and high velocity impact (HVI).
In the whole process of an explosion, there exist special
features such as large deformations, large inhomogeneities,
moving material interfaces, deformable boundaries, and free
surfaces [8]. These special features pose great challenges to
numerical simulations using the grid based methods. High
velocity impact problems involve shock waves propagating
through the colliding or impacting bodies that can behave
like fluids. Analytically, the equations of motion and a highpressure equation of state are the key descriptors of material behavior. In HVI phenomena, there also exist large deformations, moving material interfaces, deformable boundaries, and free surfaces, which are, again, very difficult for
grid based numerical methods to cope with [9].
The grid based numerical methods are also not suitable
for situations where the main concern of the object is a set
of discrete physical particles rather than a continuum. Typical examples include the interaction of stars in astrophysics,
movement of millions of atoms in an equilibrium or nonequilibrium state, dynamic behavior of protein molecules,
M.B. Liu, G.R. Liu
and etc. Simulation of such discrete systems using the continuum grid based methods is often not a good choice [10,
11].
1.2 Meshfree Methods
Over the past years, meshfree methods have been a major
research focus, towards the next generation of more effective computational methods for more complicated problems
[4, 6]. The key idea of the meshfree methods is to provide
accurate and stable numerical solutions for integral equations or PDEs with all kinds of possible boundary conditions using a set of arbitrarily distributed nodes or particles
[4, 6]. The history, development, theory and applications of
the major existing meshfree methods have been addressed in
some monographs and review articles [4, 6, 7, 1215]. The
readers may refer to these literatures for more details of the
meshfree methods. To avoid too much detour from our central topic, this article will not further discuss these meshfree
methods and techniques, except for mentioning briefly some
of the latest advancements.
For solid mechanics problems, instead of weak formulations used in FEM, we now have a much more powerful
weakened weak (W2) formulation for general settings of
FEM and meshfree methods [1619]. The W2 formulation
can create various models with special properties, such upper bound property [2023], ultra-accurate and supper convergent solutions [22, 2429], and even nearly exact solutions [30, 31]. These W2 formulations have a theoretical
foundation on the novel G space theory [1619]. All most
all these W2 models work very well with triangular mesh,
and applied for adaptive analyses for complicated geometry.
For fluid dynamics problems, the gradient smoothing
method (GSM) has been recently formulated using carefully
designed gradient smoothing domains [3235]. The GSM
works very well with unstructured triangular mesh, and can
be used effective for adaptive analysis [34]. The GSM is an
excellent alternative to the FVM for CFD problems.
One distinct meshfree method is smoothed particle hydrodynamics or SPH. The SPH is a very powerful method
for CFD problems governed by the Navier-Stokes equations.
1.3 Smoothed Particle Hydrodynamics
Smoothed particle hydrodynamics is a truly meshfree,
particle method originally used for continuum scale applications, and may be regarded as the oldest modern meshfree
particle method. It was first invented to solve astrophysical
problems in three-dimensional open space [36, 37], since
the collective movement of those particles is similar to the
movement of a liquid or gas flow, and it can be modeled by
Smoothed Particle Hydrodynamics (SPH): an Overview and Recent Developments
the governing equations of the classical Newtonian hydrodynamics.

In SPH, the state of a system is represented by a set
of particles, which possess material properties and interact
with each other within the range controlled by a weight function or smoothing function [6, 38, 39]. The discretization of
the governing equations is based on these discrete particles,
and a variety of particle-based formulations have been used
to calculate the local density, velocity and acceleration of the
fluid. The fluid pressure is calculated from the density using
an equation of state, the particle acceleration is then calculated from the pressure gradient and the density. For viscous
flows, the effects of physical viscosity on the particle accelerations can also be included. As a Lagrangian particle
method, SPH conserves mass exactly. In SPH, there is no
explicit interface tracking for multiphase flowsthe motion
of the fluid is represented by the motion of the particles, and
fluid surfaces or fluid-fluid interfaces move with particles
representing their phase defined at the initial stage.
SPH has some special advantages over the traditional grid
based numerical methods.
1. SPH is a particle method of Lagrangian nature, and the
algorithm is Galilean invariant. It can obtain the time history of the material particles. The advection and transport
of the system can thus be calculated.
2. By properly deploying particles at specific positions at
the initial stage before the analysis, the free surfaces, material interfaces, and moving boundaries can all be traced
naturally in the process of simulation regardless the complicity of the movement of the particles, which have been
very challenging to many Eulerian methods. Therefore,
SPH is an ideal choice for modeling free surface and interfacial flow problems.
3. SPH is a particle method without using a grid/mesh.
This distinct meshfree feature of the SPH method allows a straightforward handling of very large deformations, since the connectivity between particles are generated as part of the computation and can change with
time. Typical examples include the SPH applications in
high energy phenomena such as explosion, underwater
explosion, high velocity impact, and penetrations.
4. In SPH method, a particle represents a finite volume in
continuum scale. This is quite similar to the classic molecular dynamics (MD) method [11, 40] that uses a particle
to represent an atom or a molecule in nano-scale, and the
dissipative particle dynamics (DPD) method [41, 42] that
uses a particle to represent a small cluster of molecules
in meso-cale. Thus, it is natural to generalize or extend
SPH to smaller scales, or to couple SPH with molecular dynamics and dissipative particle dynamics for multiple scale applications, especially in biophysics, and biochemistry.
27
5. SPH is suitable for problems where the object under

consideration is not a continuum. This is especially true
in bio- and nano- engineering at micro and nano scale,
and astrophysics at astronomic scale. For such problems,
SPH can be a natural choice for numerical simulations.
6. SPH is comparatively easier in numerical implementation, and it is more natural to develop three-dimensional
numerical models than grid based methods.
The early SPH algorithms were derived from the probability theory, and statistical mechanics are extensively used
for numerical estimation. These algorithms did not conserve
linear and angular momentum. However, they can give reasonably good results for many astrophysical phenomena.
For the simulations of fluid and solid mechanics problems,
there are challenges to reproduce faithfully the partial differential equations governing the corresponding fluid and solid
dynamics. These challenges involve accuracy and stability
of the numerical schemes in implementing the SPH methods.
With the development of the SPH method, and the extensive applications to a wide range of problems, more attractive features have been showcased while some inherent drawbacks have also been identified. Different variants
or modifications have been proposed to improve the original SPH method. For example, Gingold and Monaghan
found the non-conservation of linear and angular momentum of the original SPH algorithm, and then introduced
an SPH algorithm that conserves both linear and angular
momentum [43]. Hu and Adams also invented an angularmomentum conservative SPH algorithm for incompressible
viscous flows [44].
Many researchers have conducted investigations on the
SPH method on the numerical aspects in accuracy, stability, convergence and efficiency. Swegle et al. identified the
tensile instability problem that can be important for materials with strength [45]. Morris noted the particle inconsistency problem that can lead to poor accuracy in the SPH
solution [46]. Over the past years, different modifications or
corrections have been tried to restore the consistency and
to improve the accuracy of the SPH method. Monaghan
proposed symmetrization formulations that were reported
to have better effects [4749]. Johnson and his co-workers
gave an axis-symmetry normalization formulation so that,
for velocity fields that yield constant values of normal velocity strains, the normal velocity strains can be exactly
reproduced [50, 51]. Randles and Libersky derived a normalization formulation for the density approximation and a
normalization for the divergence of the stress tensor [52].
Chen et al. proposed a corrective smoothed particle method
(CSPM) which improves the simulation accuracy both inside the problem domain and around the boundary area [53,
54]. The CSPM has been improved by Liu et al. in resolving problems with discontinuity such as shock waves in a
28
discontinuous SPH (DSPH) [55]. Liu et al. also proposed

a finite particle method (FPM), which uses a set of basis
function to approximate field variables at a set of arbitrarily
distributed particles [56, 57]. FPM can be regarded as an improved version of SPH and CSPM with better performance
in particle consistency. Batra et al. concurrently developed a
similar idea to FPM, and it is named modified SPH (MSPH)
[58] with applications mainly in solid mechanics. Fang et al.
further improved this idea for simulating free surface flows
[59], and they later developed a regularized Lagrangian finite point method for the simulation of incompressible viscous flows [60, 61]. A stress point method was invented to
improve the tensile instability and zero energy mode problems [6265]. Other notable modifications or corrections
of the SPH method include the moving least square particle hydrodynamics (MLSPH) [66, 67], the integration kernel correction [68], the reproducing kernel particle method
(RKPM) [69, 70], the correction for stable particle method
[71, 72], and several other particle consistency restoring approaches [6, 56, 73]. Belytschko and his co-workers have
conducted a series of stability and convergence analyses on
meshfree particle methods, and some of the numerical techniques and analyses can also be applicable to SPH [13, 72,
74].
This article is organized as follows. In Sect. 1, the background of meshfree particle methods is first addressed with
highlights on overcoming the limitations of the grid based
numerical models. The invention, features and developments of the SPH method are then briefly introduced. In
Sect. 2, the approximation schemes of the SPH method
are discussed, while some numerical techniques for developing SPH formulations are presented. The SPH formulations for the N-S equation, which governs the general fluid
flow problems, are also given. Section 3 presents a review
on the smoothing kernel function. Conditions for constructing smoothing functions are developed with examples of
smoothing functions constructed. Section 4 introduces consistency concept of SPH including kernel consistency and
particle consistency, and also provides an in-depth review
on the existing approaches for restoring particle consistency.
Several new approaches are also presented, which include a
discontinuous SPH for simulating problems with discontinuity, a general approach to restore particle inconsistency,
and a finite particle method. In Sect. 5, some important numerical topics in SPH are discussed. These special topics
include (1) solid boundary treatment, (2) representation of
solid obstacles, (3) material interface treatment, and (4) tensile instability. Different applications of the SPH method
have been reviewed in Sect. 6. Some concluding remarks
are given in Sect. 7.
M.B. Liu, G.R. Liu
2 SPH Approximation Techniques

The conventional SPH method was originally developed for
hydrodynamics problems in which the governing equations
are in strong form of partial differential equations of field
variables such as density, velocity, energy, and etc. There are
basically two steps in obtaining an SPH formulation. The
first step is to represent a function and/or its derivatives in
continuous form as integral representation, and this step is
usually termed as kernel approximation. In this kernel approximation step, the approximation of a function and its
derivatives are based on the evaluation of the smoothing kernel function and its derivatives. The second step is usually
referred to as particle approximation. In this step, the computational domain is first discretized by representing the domain with a set of initial distribution of particles representing the initial settings of the problem. After discretization,
field variables on a particle are approximated by a summation of the values over the nearest neighbor particles.
2.1 Kernel Approximation of a Function
The kernel approximation in the SPH method involves representation of a function and its derivatives using a smoothing function. The smoothing function should satisfy some
basic requirements, and it has been called kernel, smoothing
kernel, smoothing kernel function, or sometimes even weight
function in some SPH literature [38, 46, 49, 75]. A detailed
discussion on smoothing function, basic requirements and
constructing conditions will be given in Sect. 3.
The kernel approximation of a function f (x) used in the
SPH method starts from the following identity

f (x) =
f (x )(x x )dx ,
(1)
where f is a function of the position vector x, and (x x )

is the Dirac delta function given by

1, x = x ,

(x x ) =
(2)
0, x = x .
In (1), is the volume of the integral that contains x. Equation (1) implies that a function can be represented in an integral form. Since the Dirac delta function is used, the integral
representation in (2) is exact and rigorous, as long as f (x)
is defined and continuous in .
The Delta function (x x ) with only a point support, and hence (1) cannot be used for establishing discrete
numerical models. If replacing the Delta function (x x )
by a smoothing function W (x x , h) with a finite spatial
dimension h, the kernel approximation of f (x), f (x), becomes

.
f (x) =
f (x )W (x x , h)dx ,
(3)
where h is the smoothing length defining the influence or

support area of the smoothing function W . Note that as long
as W is not the Dirac delta function, the integral representation shown in (3) can only be an approximation, except for
special cases. Therefore (3) can be written as
29
where r stands for the residual. Note that W is an even function with respect to x, and (x x)W (x x , h) should be
an odd function. Hence we should have

(x x)W (x x , h)dx = 0.
(9)
f (x )W (x x , h)dx .
f (x) =
(4)
A smoothing function W is usually chosen to be an even

function for reasons given later in Sect. 3. It should also satisfy a number of conditions. The first one is the normalization condition that states

W (x x , h)dx = 1.
(5)
This condition is also termed as unity condition since the

integration of the smoothing function produces the unity.
The second condition is the Delta function property that
is observed when the smoothing length approaches zero
lim W (x x , h) = (x x ).
(6)
h0
The third condition is the compact condition

W (x x , h) = 0 when |x x | > h,
f (x) =
(7)
[f (x) + f (x)(x x) + r((x x)2 )]
W (x x , h)dx

= f (x) W (x x , h)dx + f (x) (x x)

W (x x , h)dx + r(h ),
2
f (x) = f (x) + r(h2 ).
(10)
It is clear that SPH kernel approximation of an arbitrary field

function is of second order accuracy.
2.2 Kernel Approximation of Derivatives
The approximation for the spatial derivative f (x) is obtained simply by substituting f (x) with f (x) in (4),
which gives

(11)
f (x) = [ f (x )]W (x x , h)dx ,
where the divergence in the integral is operated with respect

to the primed coordinate. Considering
[ f (x )]W (x x , h)
where is a constant related to the smoothing function for a

particle at x, and h defines the effective (non-zero) area of
the smoothing function. This effective area is usually called
as the support domain of the smoothing function for a point
at x (or the support domain of that point). Using this compact condition, integration over the entire problem domain
is localized as integration over the support domain of the
smoothing function. Therefore, the integration domain
can be the same as the support domain.
In the SPH literatures, the kernel approximation is often
said to have h2 accuracy or second order accuracy [46, 47,
49, 7577]. The observation can be obtained easily using
Taylor series expansion on (4). Note from (7) that the support domain of the smoothing function is |x x| h, the
errors in the SPH integral representation can be roughly estimated by using the Taylor series expansion of f (x ) around
x in (4). If f (x) is differentiable, we can get
Using (5) and (9), (8) becomes
(8)
= [f (x )W (x x , h)]
f (x ) [W (x x , h)],
(12)
the following equation is obtained

[f (x )W (x x , h)]dx
f (x) =
f (x ) [W (x x , h)]dx .
(13)
The first integral on the right hand side (RHS) of (13) can
be converted using the divergence theorem into an integral
over the surface S of the domain of the integration, , as

dS
f (x) = f (x )W (x x , h) n
S
f (x ) W (x x , h)dx ,
(14)
is the unit vector normal to the surface S. Since the

where n
smoothing function W is usually defined to have compact
support (see (7)), the value of W on the surface of the integral in (14) is zero in SPH. Therefore, the surface integral
on the right hand side of (14) is also zero. Hence, the kernel
approximation of the derivatives can be written from (14) as

(15)
f (x) = f (x ) W (x x , h)dx .
30
It is clear that the differential operation on a function is

transformed into a differential operation on the smoothing
function. In other words, the SPH kernel approximation of
the derivative of a field function allows the spatial gradient
to be determined from the values of the function and the
derivatives of the smoothing function W , rather than from
the derivatives of the function itself.
Kernel approximation of higher order derivatives can be
obtained in a similar way by substituting f (x) with the corresponding derivatives in (4), using integration by parts, divergence theorem and some trivial transformations. Another
approach is to repeatedly use (15) to obtain the kernel approximation of the higher order derivatives, since any higher
order derivative can always be regarded as the first order
derivative of its next lower order derivative.
Following similar analyses based on Taylor series expansion, it is easy to show that the kernel approximation of the
derivative is also of second order accuracy. Since the SPH
kernel approximations for a field function and its derivatives
are of second order accuracy, that is why the SPH method
has usually been referred as a method of second order accuracy. However, (10) is not always true because (5) and (9)
are sometimes not satisfied. For example, in a 1D problem
space, if the support domain is within the problem domain
under consideration, the integration of the smoothing function is unity (see (5)), and the integration of the first moment
of the smoothing function (see (9)) is zero. Also the surface integral in (14) is zero. Hence the SPH kernel approximations are of second order accuracy, and this is shown in
Fig. 1.
However, there are scenarios in which the support domain intersects with the problem domain boundary, as
shown in Fig. 2. Therefore, the smoothing function W
is truncated by the boundary, and the integration of the
smoothing function is no longer unity. The integration of the
first moment term of the smoothing function and the surface
integral in (14) are also no longer zero. At such scenarios,
the SPH kernel approximations are not of second order accuracy.
M.B. Liu, G.R. Liu
Fig. 1 Schematic illustration of the scenarios in which the support

domain is located within the problem domain. For such scenarios, the
SPH kernel approximations are of second order accuracy
Fig. 2 Schematic illustration of the scenarios in which the support

domain intersects with the problem domain. For such scenarios, the
SPH kernel approximations are not exactly of second order accuracy
2.3 Particle Approximation

The second step of SPH method is the particle approximation, which involves representing the problem domain using
a set of particles, and then estimating field variables on this
set of particles. Considering a problem domain filled with
a set of particles (usually arbitrarily distributed, see Fig. 3
for illustration in a two-dimensional domain). These particles can either be centered particles initially generated using existing mesh generation tools or concentrated particles
initially generated using some kind of space discretization
model such as the particle-fill model in AUTODYN [78].
The state of the system is represented by these particles,
Fig. 3 SPH particle approximations in a two-dimensional problem domain with a surface S. W is the smoothing function that is used to
approximate the field variables at particle i using averaged summations over particles j within the support domain with a cut-off distance
of hi
each associated with field properties. These particles can be

used not only for integration, interpolation or differencing,
but also for representing the material. The volume of a subsection is lumped on the corresponding particle. Therefore
one particle i is associated with a fixed lumped volume Vi

without fixed shape. If the particle mass and density are concerned, the lumped volume can also be replaced with the
corresponding mass to density ratio mi /i . These particles
can be fixed in an Eulerian frame or move in a Lagrangian
frame.
After representing the computational domain with a finite
number of particles, the continuous form of kernel approximation expressed in (4) can be written in discretized form
of a summation of the neighboring particles as follows
f (x) =
N

mj
j =1
f (x j )W (x x j , h),
(16)
where N is the total number of particles within the influence

area of the particle at x. It is the total number of particles
that are within the support domain which has a cut-off distance, characterized by the smoothing length, h, multiplied
by a scalar constant . This procedure of summation over the
neighboring particles is referred to as particle approximation, which states that the value of a function at a particle can
be approximated by using the average of the values of the
function at all the particles in the support domain weighted
by the smoothing function. Following the same procedure,
the particle approximation of a derivative can be obtained as
f (x) =
N

mj
j =1
f (x j ) W (x x j , h),
(17)
where the gradient W in the above equation is evaluated at

particle j . Equation (17) states that the value of the gradient
of a function at a particle located at x can be approximated
by using the summation of those values of the function at
all the particles in the support domain weighted by the gradient of the smoothing function. The particle approximation
in (16) and (17) converts the continuous form of kernel approximation of a field function and its derivatives to the discrete summations over a set of particles. The use of particle summations to approximate the integral is, in fact, a key
approximation that makes the SPH method simple without
using a background mesh for numerical integration, and it is
also the key factor influencing the solution accuracy of the
SPH method.
One important aspect is that the particle approximation
in the SPH method introduces the mass and density of the
particle into the equations. This can be conveniently applied to hydrodynamic problems in which the density is a
key field variable. This is probably one of the major reasons
for the SPH method being particularly popular for dynamic
fluid flow problems. If the SPH particle approximation is
applied to solid mechanics problems, special treatments are
required. One of the ways is to use the SPH approximation
to create shape functions, and to establish the discrete system equations [4].
31
The particle approximation is, however, related to some

numerical problems inherent in the SPH method, such as
the particle inconsistency and the tensile instability, as will
be addressed in the following sections. One basic reason
is that the discrete summation is only taken over the particles themselves (collocation). In general, in meshfree methods, to achieve stability and accuracy, the number of sampling points for integration should be more than the field
nodes (particles). This is especially true for meshfree methods based on weak forms for solid mechanics problems [4].
Otherwise, it may (not always) lead to some kind of instability problems.
2.4 Techniques for Deriving SPH Formulations
By using the above-described procedure of kernel approximation and particle approximation, SPH formulations for
partial differential equations can always be derived. There
are in fact a number of ways to derive SPH formulation of
PDEs. Benz used one approach to derive the SPH equations
for PDEs that is to multiply each term in the PDEs with the
smoothing function, and integrate over the volume with the
use of integration by parts and Taylor expansions [79]. Monaghan employed a straightforward approach of directly using (16) and (17) [49]. In that approach, the following two
identities are employed to place the density inside the gradient operator,
1
[ (f (x)) f (x) ],

f (x)
f (x)
+ 2 .
f (x) =
f (x) =
(18)
(19)
The above two identities may be substituted into the integral

in (11). The same procedure of the particle approximation
to obtain (17) is applied to each gradient term on the right
hand side of (18) and (19). Note that each expression at the
outside of every gradient term is evaluated at the particle
itself, the results from (18) and (19) for the divergence of
f (x) at particle i are obtained as
N
1
f (x i ) =
mj [f (x j ) f (x i )] i Wij , (20)
i
j =1
and
f (x i )
N
f (x j ) f (x i )
= i
mj
+
i Wij .
j2
i2
j =1
(21)
One of the good features for the above two equations is

that the field function f (x) appears pairwisely and involves
32
M.B. Liu, G.R. Liu
asymmetric and symmetric SPH formulations. These asymmetric and symmetric formulations can help to improve the
numerical accuracy in SPH simulations [6, 49, 56].
Besides the above-mentioned two identities, some other
rules of operation can be convenient in deriving the SPH
formulations for complex system equations [6]. For example, for two arbitrary functions of field variables f1 and f2 ,
the following rules exist.
f1 f2 = f1 f2 ,
(22)
f1 f2 = f1 f2 .
(23)
Hence, an SPH approximation of the sum of functions

equals to the sum of the SPH approximations of the individual function, and an SPH approximation of a product of
functions equals to the product of the SPH approximations
of the individual functions.
If f1 is a constant denoted by c, we should have
cf2 = cf2 .
(24)
It is clear that the SPH approximation operator is a linear

operator. It is also easy to show that the SPH approximation
operator is commutative, i.e.,
f1 + f2 = f2 + f1 ,
(25)
and
f1 f2 = f2 f1 .
(26)
For convenience, the SPH approximation operator is

omitted in later sections.
2.5 SPH Formulations for Navier-Stokes (N-S) Equations
Using the afore-mentioned kernel and particle approximation techniques with necessary numerical tricks, it is possible to derive SPH formulations for partial differential equations governing the physics of fluid flows. For example, for
Navier-Stokes equations controlling the general fluid dynamic problems, we have
D
v
Dt = x ,

Dv
1
(27)
Dt = x + F,
De = v ,
Dt
x
where the Greek superscripts and are used to denote
the coordinate directions, the summation in the equations is
taken over repeated indices, and the total time derivatives are
taken in the moving Lagrangian frame. The scalar density
, and internal energy e, the velocity component v , and
the total stress tensor are the dependent variables. F is
the external forces such as gravity. The spatial coordinates
x and time t are the independent variables. The total stress

tensor is made up of two parts, one part of isotropic
pressure p and the other part of viscous stress , i.e., =
p + . For Newtonian fluids, the viscous shear stress
should be proportional to the shear strain rate denoted by
through the dynamic viscosity , i.e., = , where
=
v
v
2
+ ( v) .
x
x
3
Substituting the SPH approximations for a function and

its derivative (as shown in (16) and (17)) to the N-S equations, the SPH equations of motion for the N-S equations
can be written as
N
Wij
Di
j =1 mj v ij x ,
Dt =

N
Dv i
i
Wij
j
(28)
+ Fi ,
j =1 mj 2 + 2
Dt =
x i
i
j

pj
Wij
pi
i
1 N
i
De
+ 2i i i ,
2 + 2 v ij
j =1 mj
Dt = 2
i
x i
where v ij = v i v j . Equation (28) is a set of commonly

used SPH equations for the N-S equations. It should be
noted that by using different numerical tricks, it is possible to get other different forms of SPH equations for the
same partial differential equations. The obtained SPH formulations may have special features and advantages suitable
for different applications [6]. One typical example is the approximation of density. If the field function is the density,
(16) can be re-written as
i =
N
mj Wij .
(29)
j =1
This is another approach to obtain density directly from the

SPH summation of the mass of all particles in the support
domain of a given particle, rather than from the continuum
equation. Compared to the SPH formulations on density
change in (28), this summation density approach conserves
mass exactly, but suffers from serious boundaries deficiency
due to the particle inconsistency. A frequently used way to
remediate the boundaries deficiency is the following normalization form by the summation of the smoothing function
itself [52, 53]
N
j =1 mj Wij
.
(30)
i = N mj
j =1 ( j )Wij
3 SPH Smoothing Function
3.1 Review on Commonly Used Smoothing Functions
One of the central issues for meshfree methods is how to effectively perform function approximation based on a set of
nodes scattered in an arbitrary manner without using a predefined mesh or grid that provides the connectivity of the
nodes. In the SPH method, the smoothing function is used
for kernel and particle approximations. It is of utmost importance in the SPH method as it determines the pattern to
interpolate, and defines the cut-off distance of the influencing area of a particle.
Many researchers have investigated the smoothing kernel, hoping to improve the performance of the SPH method,
and/or to generalize the requirements for constructing the
smoothing kernel function. Fulk numerically investigated a
number of smoothing kernel functions in one-dimensional
space, and the obtained results are basically valid for regularly distributed particles [39, 75]. Swegle et al. revealed the
tensile instability, which is closely related to the smoothing kernel function [45]. Morris studied the performances
of several different smoothing functions, and found that by
properly selecting the smoothing function, the accuracy and
stability property of the SPH simulation can be improved
[46, 80]. Omang provided investigations on alternative kernel functions for SPH in cylindrical symmetry [81]. Jin and
Ding investigated the criterions for smoothed particle hydrodynamics kernels in stable field [82]. Capuzzo-Dolcetta
gave a criterion for the choice of the interpolation kernel in SPH [83]. Cabezon and his co-workers proposed a
one-parameter family of interpolating kernels for SPH studies [84].
Different smoothing functions have been used in the SPH
method as shown in published literatures. Various requirements or properties for the smoothing functions have been
discussed. Major properties or requirements are now summarized and described in the following discussion.
1. The smoothing function must be normalized over its support domain (Unity)

W (x x , h)dx = 1.
(31)
This normalization property ensures that the integral of

the smoothing function over the support domain to be
unity. It can be shown in the next section that it also ensures the zero-th order consistency (C 0 ) of the integral
representation of a continuum function.
2. The smoothing function should be compactly supported
(Compact support), i.e.,
W (x x ) = 0,
for |x x | > h.
(32)
The dimension of the compact support is defined by the

smoothing length h and a scaling factor , where h is
the smoothing length, and determines the spread of the
specified smoothing function. |x x | h defines the
support domain of the particle at point x. This compact
supportness property transforms an SPH approximation
33
from a global operation to a local operation. This will

later lead to a set of sparse discretized system matrices,
and therefore is very important as far as the computational efforts are concerned.
3. W (x x ) 0 for any point at x within the support domain of the particle at point x (Positivity). This property
states that the smoothing function should be non-negative
in the support domain. It is not mathematically necessary
as a convergent condition, but it is important to ensure
a meaningful (or stable) representation of some physical phenomena. A few smoothing functions used in some
literatures are negative in parts of the support domain.
However in hydrodynamic simulations, negative value of
the smoothing function can have serious consequences
that may result in some unphysical parameters such as
negative density and energy.
4. The smoothing function value for a particle should be
monotonically decreasing with the increase of the distance away from the particle (Decay). This property is
based on the physical consideration in that a nearer particle should have a bigger influence on the particle under
consideration. In other words, with the increase of the
distance of two interacting particles, the interaction force
decreases.
5. The smoothing function should satisfy the Dirac delta
function condition as the smoothing length approaches
zero (Delta function property)
lim W (x x , h) = (x x ).
h0
(33)
This property makes sure that as the smoothing length

tends to be zero, the approximation value approaches the
function value, i.e. f (x) = f (x).
6. The smoothing function should be an even function
(Symmetric property). This means that particles from
same distance but different positions should have equal
effect on a given particle. This is not a very rigid condition, and it is sometimes violated in some meshfree particle methods that provide higher consistency.
7. The smoothing function should be sufficiently smooth
(Smoothness). This property aims to obtain better approximation accuracy. For the approximations of a function and its derivatives, a smoothing function needs to be
sufficiently continuous to obtain good results. A smoothing function with smoother value of the function and
derivatives would usually yield better results and better
performance in numerical stability. This is because the
smoothing function will not be sensitive to particle disorder, and the errors in approximating the integral interpolants are small, provided that the particle disorder is
not extreme [6, 49, 75].
Any function having the above properties may be employed as SPH smoothing functions, and many kinds of
34
M.B. Liu, G.R. Liu
smoothing functions have been used. Lucy in the original

SPH paper [36] used a bell-shaped function
W (x x , h) = W (R, h)

(1 + 3R)(1 R)3 ,
= d
0,
R 1,
R > 1,
(34)
where d is 5/4h, 5/h2 and 105/16h3 in one-, two- and

three-dimensional space, respectively, so that the condition
of unity can be satisfied for all the three dimensions. R is
the relative distance between two points (particles) at points
|
x and x , R = hr = |xx
h , where r is the distance between
the two points.
Gingold and Monaghan in their original paper [37] selected the following Gaussian kernel to simulate the nonspherical stars
W (R, h) = d eR ,
2
(35)
where d is 1/ 1/2 h, 1/h2 and 1/ 3/2 h3 , respectively,

in one-, two- and three-dimensional space, for the unity
requirement. The Gaussian kernel is sufficiently smooth
even for high orders of derivatives, and is regarded as a
golden selection since it is very stable and accurate especially for disordered particles. It is, however, not really
compact, as it never goes to zero theoretically, unless R approaches infinity. Because it approaches zero numerically
very fast, it is practically compact. Note that it is computationally more expensive since it can take a longer distance
for the kernel to approach zero. This can result in a large
support domain with more particles for particle approximations.
The most frequently used smoothing function may be the
cubic B-spline function, which was originally used by Monaghan and Lattanzio [85]
2
1 3
2
3 R + 2R ,
W (R, h) = d 1 (2 R)3 ,
6
0,
0 R < 1,
1 R < 2,
(36)
R 2.
In one-, two- and three-dimensional space, d = 1/ h,

15/7h2 and 3/2h3 , respectively. The cubic spline function has been the most widely used smoothing function in
the emerged SPH literatures since it closely resembles a
Gaussian function while having a narrower compact support. However, the second derivative of the cubic spline is
a piecewise linear function, and accordingly, the stability
properties can be inferior to those of smoother kernels.
Morris has introduced higher order (quartic and quintic)
splines that are more closely approximating the Gaussian
and more stable [46, 80]. The quartic spline is
W (R, h)
(R + 2.5)4 5(R + 1.5)4 + 10(R + 0.5)4 ,
0 R < 0.5,
(2.5 R)4 5(1.5 R)4 ,

(37)
= d
0.5 R < 1.5,
(2.5 R)4 , 1.5 R < 2.5,
0, R > 2.5,
where d is 1/24h in one-dimensional space. The quintic

spline is
(3 R)5 6(2 R)5 + 15(1 R)5 ,
0 R < 1,
(3 R)5 6(2 R)5 ,

W (R, h) = d
(38)
1 R < 2,
(3 R) , 2 R < 3,
0, R > 3,
where d is 120/ h, 7/478h2 and 3/359h3 in one-, twoand three-dimensional space, respectively.
Johnson et al. used the following quadratic smoothing
function to simulate high velocity impact problems [86]

3 2 3
3
R R+
, 0 R 2,
(39)
W (R, h) = d
16
4
4
where in one-, two- and three-dimensional space, d = 1/ h,
2/h2 and 5/4h3 , respectively. Unlike other smoothing
functions, the derivative of this quadratic smoothing function always increases as the particles move closer, and always decreases as they move apart. This was regarded by the
authors as an important improvement over the cubic spline
function, and it was reported to relieve the problem of compressive instability.
Some higher order smoothing functions that are devised
from lower order forms have been constructed, such as the
super-Gaussian kernel [85]

3
2
2
R eR , 0 R 2,
(40)
W (R, h) = d
2
where d is 1/ in one-dimensional space. One disadvantage of the high order smoothing function is that the kernel
is negative in some region of its support domain. This may
lead to unphysical results for hydrodynamic problems [75].
The smoothing function has been studied mathematically
in detail by Liu and his co-workers. They proposed a systematical way to construct a smoothing function that may meet
different needs [38]. A new quartic smoothing function has
been constructed to demonstrate the effectiveness of the approach for constructing a smoothing function as follows.
2 9 2 19 3
5 4
R ), 0 R 2,
( R + 24 R 32
W (R, h) = d 3 8
0,
R > 2,
(41)
where d is 1/ h, 15/7h2 and 315/208h3 in one-, twoand three-dimensional space, respectively. Note that the centre peak value of this quartic smoothing function is defined
as 2/3. The quartic smoothing function behaves very much
like the widely used cubic B-spline function given in (36),
but has only one piece, and hence is much more convenient and efficient to use. More discussions on this quartic
smoothing function will give in the next section.
3.2 Generalizing Constructing Conditions
Major requirements of an SPH smoothing function have
been addressed in Sect. 3.1. Some of these requirements can
be derived by conducting Taylor series analysis. This analysis is carried out at the stage of the SPH kernel approximation for a function and its derivatives. It shows that, to
exactly approximate a function and its derivatives, certain
conditions need to be satisfied. These conditions can then be
used to construct the smoothing functions.
Considering the SPH kernel approximation for a field
function f (x) as show in (4), if f (x) is sufficiently smooth,
applying Taylor series expansion of f (x ) in the vicinity of
x yields
1
f (x ) = f (x) + f (x)(x x) + f (x)(x x)2 +
2

n
k
k
(k)

(1) h f (x) x x k
=
k!
h
k=0

x x
,
(42)
+ rn
h
where rn is the remainder of the Taylor series expansion.
Substituting (42) into (4) leads to

n

x x
,
(43)
Ak f (k) (x) + rn
f (x) =
h
k=0
where
Ak =
(1)k hk
k!

x x
h
k
W (x x , h)dx .
(44)
Comparing the LHS with the RHS of (43), in order for

f (x) to be approximated to n-th order, the coefficients Ak
must equal to the counterparts for f (k) (x) at the LHS of
(43). Therefore, after trivil transformation, the following
conditions for the smoothing function W can be obtained
as follows

M0 = W (x x , h)dx = 1

M1 = (x x )W (x x , h)dx = 0

2

(45)
M2 = (x x ) W (x x , h)dx = 0 ,
..
.
Mn = (x x )n W (x x , h)dx = 0
35
where Mk is the k-th moments of the smoothing function.

Note that the first equation in (45) is, in fact, the unity condition expressed in (31), and the second equation in (45)
stands for the symmetric property. Satisfaction of these two
conditions ensures the first order consistency for the SPH
kernel approximation for a function.
Also performing Taylor series analysis for the SPH kernel approximation of the derivatives of a field function f (x),
using the concept of integration by parts, and divergence theorem with some trivial transformation, the following equations
W (x x , h)|S = 0
(46)
and

M1 = (x x )W (x x , h)dx = 1

2

M2 = (x x ) W (x x , h)dx = 0
..
.

n

Mn = (x x ) W (x x , h)dx = 0
M0 =
W
(x
x , h)dx = 0
(47)
can be obtained. Equation (46) actually specifies that the

smoothing function vanishes on the surface of the support
domain. This is compatible to the compactness condition of
the smoothing function. Equation (47) defines the conditions
with which the derivatives of the smoothing function should
be satisfied. Note that (45) and (47) are actually compatible
considering integration by parts, divergence theorem and the
boundary value vanishing effects (see (46)) of the smoothing
function.
Performing Taylor series analysis on the SPH kernel approximation for the second derivatives, similar equations
can be obtained. Except for the requirements on the second derivatives of the momentums, the first derivative of
the smoothing function also needs to vanish on the surface,
which is
W (x x , h)|S = 0.
(48)
Equations (45)(48) can be used to construct smoothing

functions. It can be seen that the conditions of smoothing
functions can be classified into two groups. The first group
shows the ability of a smoothing function to reproduce polynomials. Satisfying the first group, the function can be approximated to n-th order accuracy. The second group defines
the surface values of a smoothing function as well as its first
derivatives, and is the requirements of the property of compact support for the smoothing function and its first derivative. Satisfying these conditions, the first two derivatives of
the function can be exactly approximated to the n-th order.
36
M.B. Liu, G.R. Liu
3.3 Constructing SPH Smoothing Functions

By using above-mentioned conditions, it is possible to have
a systematic way to construct the SPH smoothing functions.
If the smoothing function is assumed to be a polynomial
dependent only on the relative distance of the concerned
points, it can be assumed to have the following form in the
support domain with an influence width of h.
W (x x , h) = W (R) = a0 + a1 R + a2 R 2 + + an R n .
(49)
It is clear that a smoothing function in the abovementioned form is a distance function since it depends on the
relative distance. It is easy to show that for the second derivative of the smoothing function to exist, a1 should vanish.
Substituting this polynomial form smoothing function into
the conditions (see (45)(48)), the parameters a0 , a2 , . . . , an
can be calculated from the resultant linear equations, and
then the smoothing function can be determined.
There are several issues that need further consideration.
Firstly, a smoothing function derived from this set of conditions (see (45)) will not necessarily be positive in the entire
support domain, especially when high order reproducibility is required. Such a negative smoothing function may result in unphysical solutions, for example, negative density
(mass) and negative energy. For this reason, smoothing functions used in SPH literatures are generally non-negative for
CFD problems. On the other hand, for the even moments
(k = 2, 4, 6, . . .) to be zero, a smoothing function has to
be negative in some parts of the region. This implies that
one cannot have both non-negativity and high-order reproducibility at the same time.
Secondly, in constructing a smoothing function, the center peak value is a factor that needs to be considered. The
center peak value of a smoothing function is very important
since it determines how much the particle itself will contribute to the approximation. Revisiting (45), if a positive
smoothing function is used, the highest order of accuracy
for the function approximation
is second order. Therefore,
the second momentum (M2 = (x x )2 W (x x , h)dx )
can be used as a rough indicator to measure the accuracy
of the kernel approximation. The smaller the second moment M2 is, the more accurate the kernel approximation is.
The center peak value of a smoothing function is closely
related to M2 . A positive smoothing function with a large
center peak value will have a smaller second moment M2 .
This implies that a smoothing function is closer to the Delta
function, and therefore is more accurate in terms of kernel
approximations.
Thirdly, in some circumstances, a piecewise smoothing function is preferable since the shape of the piecewise
smoothing function is easier to be controlled by changing
the number of the pieces and the locations of the connection

points. For example, consider the general form of a smoothing function with two pieces,
W1 (R), 0 R < R1 ,
(50)
W (R) = W2 (R), R1 R < R2 ,
0
R2 R.
The function itself and the first two derivatives at the
connection points should be continuous, i.e., W1 (R1 ) =
W2 (R1 ), W1 (R1 ) = W2 (R1 ) and W1 (R1 ) = W2 (R1 ). Considering the requirements at these points as well as the compact support property, one possible form of the smoothing
function is
n
n
b1 (R1 R) + b2 (R2 R) , 0 R < R1 ,
W (R) = d b2 (R2 R)n , R1 R < R2 ,
(51)
0, R2 R.
It is also feasible to construct smoothing function with more
pieces using similar expressions.
To show the effectiveness of this approach to constructing general SPH smoothing functions Liu et al. [38] derived
a new quartic smoothing function using the following conditions
the unity condition,
compact support of the smoothing function,
compact support of the first derivative of the smoothing
function,
centre peak value.
The constructed smoothing function is given as W (R, h) =
5 4
3
d ( 23 98 R 2 + 19
24 R 32 R ), for 0 R 2, where d is
2
1/ h, 15/7h and 315/208h3 in one-, two- and threedimensional space, respectively. Note that the centre peak
value of this quartic smoothing function is defined as 2/3.
As defined, this quartic function satisfies the normalization condition, while the function itself and its first derivative have compact support. It is very close to the most commonly used cubic spline (see (36)) with a same center peak
value of 2/3, and monotonically decreases with the increase
of the distance as show in Fig. 4. However, this quartic function produces a smaller second momentum than the cubic
spline function, and therefore can produce better accuracy
for kernel approximation. Also this quartic smoothing function has a smoother second derivative than the cubic spline
smoothing function, thus the stability properties should be
superior to those of the cubic spline function, as reported by
many researchers that a smoother second derivative can lead
to less instability in SPH simulation [80, 87].
4 Consistency
In early studies, the SPH method has usually been reported
to have second order accuracy. This is reasonable, as the ker-
37
4.1 Consistency in Kernel Approximation (Kernel

Consistency)
For a constant (0th order polynomial) function f (x) = c
(where c is a constant) to be exactly reproduced by the SPH
kernel approximation, following (1), we require

(52)
f (x) = cW (x x , h)dx = c,
or
Fig. 4 The quartic smoothing function constructed by Liu et al. by using the smoothing function constructing conditions [38]. The shapes of
the quartic function and its first derivative are very close to the shapes
of the cubic spline function and its first derivative. However, this one
piece smoothing function is expected to produce better accuracy as it
has smaller second momentum. It is also expected to be more stable
since it has a continuous second derivative
nel approximation of a field function has a remainder of second order (see (8)). However, the kernel approximation of
a function and/or its derivatives may not necessarily be of
second order accuracy as the smoothing function may not
be an even function since it is truncated when approaching
the boundary region. Moreover, high order accuracy of the
kernel approximation does not necessarily mean high order
accuracy of the SPH simulations, as it is the particle approximation rather than the kernel approximation that eventually
determines the accuracy of the SPH simulations.
Following the Lax-Richtmyer equivalence theorem, we
know that if a numerical model is stable, the convergence of
the solution to a well-posed problem will then be determined
by the consistency of the function approximation. Therefore
the consistency of SPH approximation is crucial. In finite
element method, to ensure the accuracy and convergence of
an FEM approximation, the FEM shape function must satisfy a certain degree of consistency. The degree of consistency can be characterized by the order of the polynomial
that can be exactly reproduced by the approximation using
the shape function. In general, if an approximation can reproduce a polynomial up to k-th order exactly, the approximation is said to have k-th order consistency or C k consistency. The concept of consistency from traditional finite element methods can also be used for meshfree particle methods such as smoothed particle hydrodynamics [3, 6, 13, 72,
74, 88]. Therefore it is similarly feasible to investigate the
consistency of SPH kernel and particle approximations in
reproducing a polynomial by using the smoothing function.
W (x x , h)dx = 1.
(53)
Equation (53) is exactly the normalization condition described previously.

Further, for a linear function f (x) = c0 + c1 x (where c0
and c1 are constants) to be exactly reproduced, we must have

f (x) = (c0 + c1 x )W (x x , h)dx = c0 + c1 x.
(54)
Using (52), (54) can be simplified as

x W (x x , h)dx = x.
(55)
Multiplying x to both side of (53), we have the following

identity

xW (x x , h)dx = x.
(56)
Subtracting (55) from the above identity yields

(x x )W (x x , h)dx = 0.
(57)
Equation (57) is the symmetric condition described in

Sect. 3.1.
More generally, by performing Taylor series analyses
of a function f (x) (=
on the kernel approximation
f (x )W (x x , h)dx ) in a one-dimensional space, we
have already obtained a set of requirements of the smoothing function in Sect. 3.2, as expressed in (45). Equations (53)
(normalization condition) and (57) (symmetric condition)
are actually components in (45), which describes the 0th
and 1st moments.
For the integrations expressed in (45), the integration domain is assumed to be a fully continuous support domain
that is not truncated by the boundaries. Equation (45) states
the requirements on the moments of a smoothing function to
reproduce certain order of polynomials.
Equation (45) can thus be used as an approximation accuracy indicator. If a smoothing function satisfies (45), a
function can be approximated to n-th order accuracy. Furthermore the 0th moment in (45) states the normalization
38
M.B. Liu, G.R. Liu
condition, and the 1st moment states the symmetry property

of the smoothing function.
Similar to the consistency concept in the traditional FEM,
if an SPH approximation can reproduce a polynomial of up
to nth order exactly, the SPH approximation is said to have
nth order or C n consistency. If the consistency of an SPH
kernel approximation in continuous form is termed as kernel consistency, the kernel consistency of an SPH kernel approximation is of nth order when the smoothing function
satisfies (45). Therefore the expressions in (45) are also the
kernel consistency conditions of the smoothing function for
an SPH kernel approximation.
Note that if the SPH kernel approximations are carried
out for regions truncated by boundaries, constant and linear functions can not be reproduced exactly since (53) and
(57) are not satisfied for these regions. Therefore we can
conclude that, since a conventional smoothing function satisfies the normalization and symmetric conditions, the conventional SPH method has up to C 1 consistency for the interior regions. However, for the boundary regions, it even does
not have C 0 kernel consistency.
4.2 Consistency in Particle Approximation (Particle
Consistency)
Satisfying the consistency conditions at the kernel approximation stage does not necessarily mean that the discretized
SPH model will have such a consistency. This is because
such a consistency can be distorted by the particle approximation process in discrete SPH model. Therefore, the consistency analysis should be conducted for the discrete SPH
model in the particle approximation process, and this consistency can be termed as particle consistency.
The discrete counterparts of the constant and linear consistency conditions as expressed in (53) and (57) are
N
W (x x j , h)vj = 1,
(58)
j =1
and
N

(x x j )W (x x j , h)vj = 0.
(59)
j =1
These discretized consistency conditions are not satisfied

in general. One obvious and simple example is the particle
approximations at the boundary particles (Fig. 5a). Even for
uniform particle distribution, due to the unbalanced particles
contributing to the discretized summation, the LHS of (58)
is smaller than 1 and the LHS of (59) will not vanish, due
to the truncation of the smoothing function by the boundary. For cases with irregularly distributed particles (Fig. 5b),
it is also easy to verify that even for the interior particles
Fig. 5 SPH particle approximations in one-dimensional cases. (a) Particle approximation for a particle whose support domain is truncated by
the boundary. (b) Particle approximation for a particle with irregular
particle distribution in its support domain
whose support domains are not truncated, the constant and

linear consistency conditions in discretized forms may not
be exactly satisfied. Therefore the original SPH method does
not even have C 0 consistency in the particle approximation.
It is clear that the inconsistency caused by the particle approximation is closely related to the corresponding kernel
approximation and particles involved in the approximation.
Such an inconsistency problem results in directly the solution inaccuracy in the original SPH method.
Besides the particle approximation features associated
with boundary particles or irregular distributed particles, the
choice of the smoothing length is also important in the particle approximation process. In a one-dimensional domain
with the cubic spline smoothing function, it is easy to verify
that for uniformly distributed interior particles, the original
SPH method has C 0 particle consistency if the smoothing
length is taken exactly as the particle spacing (h = x) since
(58) is satisfied. However, varying the smoothing length can
results in a dissatisfaction of (58), leading to poor accuracy
in the original SPH method. This is a reason why we often
need to examine the influence of the smoothing length on
the SPH approximation results.
In summary, the original SPH models, in general, do not
have even C 0 consistency. Such an inconsistency originates
from the discrepancy between the SPH kernel and particle approximations. Boundary particles, irregular distributed
particles, and variable smoothing length can usually produce
inconsistency in the particle approximation process. In the
next section, we discuss ways to restore the consistency in
SPH models.
4.3 Review on Approaches for Restoring Consistency
It has been shown that the original SPH method even does
not have 0th particle consistency. Different approaches have
been proposed to improve the particle inconsistency and
hence the SPH approximation accuracy. Some of them involve reconstruction of a new smoothing function so as
to satisfy the discretized consistency conditions. However,
these approaches are usually not preferred for hydrodynamic
simulations because the reconstructed smoothing function
can be partially negative, non-symmetric, and not monotonically decreasing. Approaches which improve the particle
consistency without changing the conventional smoothing
function are usually more preferable in simulating hydrodynamics.
One early approach [49, 52] is based on the antisymmetric assumption of the derivative of a smoothing function
N
Wi, vj = 0,
N

(fj fi )Wi, vj ,
(61)
j =1
or
fi, =
f (x) = fi + (x xi )fi,x +
N

(fj + fi )Wi, vj .
(62)
j =1
It should also be noted that (60) is not necessarily

valid, even

if its corresponding continuous counterpart Wi, dx = 0 is
valid (for interior regions). This is also a manifestation of
the particle inconsistency. Therefore (61) and (62) actually
use the particle inconsistency in approximating the derivative of the smoothing function to offset or balance the particle inconsistency in approximating the derivatives of a field
function, with a hope to improve the accuracy of the approximations.
Randles and Libersky [52] derived a normalization formulation for the density approximation
N
j =1 j Wij vj
i = N
,
(63)
j =1 Wij vj
and a normalization for the divergence of the stress tensor
N
j =1 (j i ) i Wij vj
( )i = N
,
(64)
j =1 (x j x i ) i Wij vj
where is the tensor product. Again, (63) and (64) also use
the inconsistency in approximating the smoothing function
and its derivatives to offset the inconsistency in approximating a field function and its derivatives, also with an aim to
improve the accuracy of the approximations.
Based on Taylor series expansion on the SPH approximation of a function, Chen et al. [54] suggested a corrective smoothed particle method (CSPM). In one-dimensional
space, the process of CSPM can be briefed as follows.
(x xi )2
fi,xx + .
2!
(65)
Multiplying both sides of (65) by the smoothing function W

and integrating over the entire computational domain yield

f (x)Wi (x)dx

= fi
where Wi, = Wi (x)/x , in which is the dimension

index repeated from 1 to d (d is the number of dimensions). Therefore when approximating the derivative of a
function f , the particle approximation can be rewritten as
fi, =
Performing Taylor series expansion at a nearby point xi ,

a sufficiently smooth function f (x) can be expressed as
(60)
j =1
39

Wi (x)dx + fi,x
fi,xx
2
(x xi )Wi (x)dx

(x xi )2 Wi (x)dx + .
(66)
If the terms involving derivatives in this equation are neglected, a corrective kernel approximation for function f (x)
at particle i is obtained as

f (x)Wi (x)dx
fi =
.
(67)
Wi (x)dx
For a conventional smoothing function (non-negative and
symmetric), the second term at the RHS of (66) is zero for
interior region and not zero for boundary region. Therefore
the corrective kernel approximation expressed in (67) is also
of 2nd order accuracy for interior region and 1st order accuracy for boundary region. Comparing (67) with (16), it is
found that for the interior regions, the kernel approximations
in the original SPH and CSPM are actually the same due
to the satisfaction of the normalization condition (in continuous form). For the boundary regions, since the integral
of the smoothing function is truncated by the boundary, the
normalization condition cannot be satisfied. By retaining the
non-unity denominator, CSPM restores the C 0 kernel consistency.
The corresponding particle approximation for function
f (x) at particle i can be obtained using summation over
nearest particles for each term in (66) and again neglecting
the terms related to derivatives
N
j =1 fj Wij vj
fi = N
.
(68)
j =1 Wij vj
It is noted that the particle approximation of the second term
at the RHS of (66) is not necessarily zero even for the interior particles due to the irregularity of the particles. Therefore strictly speaking, the particle approximation expressed
in (68) is of 1st order accuracy for both the interior and
boundary particles. Only if the particles are uniformly distributed can the particle approximation of the second term
at the RHS of (66) be zero. In this case, the particle approximation expressed in (68) is of 2nd order accuracy for the
uniformly distributed interior particles.
40
M.B. Liu, G.R. Liu
If replacing Wi (x) in (66) with Wi,x and neglecting the

second and higher derivatives, a corrective kernel approximation for the first derivative is generated as

[f (x) f (xi )]Wi,x (x)dx

.
(69)
fi,x =
(x xi )Wi,x (x)dx
half of the support domain, i.e., x i < d b, the integration

of the multiplication of f (x) and the smoothing function W
over the entire support domain can be divided into two parts
b
d
f (x)Wi (x)dx =
f (x)Wi (x)dx
a
The particle approximations corresponding to (69) is

N
j =1 (fj
fi )Wi,x vj
j =1 (x j
x i )Wi,x vj
fi,x = N
(71)
f (x)Wi (x)dx.
d
(70)
Similarly, The CSPM kernel approximations for the

derivatives are also of second order accuracy (or 1st order consistency) for interior regions, but 1st order accuracy (or 0th order consistency) for boundary regions. Except for cases with uniformly distributed interior particles,
the CSPM particle approximations for the derivatives are of
1st order accuracy (or 0th order consistency) for both the
interior and boundary particles.
Expanding f (x) in the first integral on the right hand side

around point x i , and around another arbitrary point x k in the
second integral, where d x k b gives
b
f (x)Wi (x)dx
a

= f (x i )
+ f (x i )
Wi (x)dx + f (x k )
Wi (x)dx
d
(x x i )Wi (x)dx
+ f (x k )
4.4 An SPH Formulation for Discontinuity
(x x k )Wi (x)dx + r(h2 ).
(72)
It is clear that CSPM can have better accuracy than the

conventional SPH method as it improves the boundary deficiency problem. It was reported that the CSPM can also
reduce the so-called tensile instability inherent in the traditional SPH method [89, 90].
One notable point in CSPM is that it is based on Taylor series analysis on the kernel approximation of a field
function and its derivatives. To perform Taylor series analysis, the function under consideration should be sufficiently
smooth. Therefore, CSPM is not applicable to problems
with discontinuities such as hydrodynamic problems that
generate shock waves.
Liu and his co-workers proposed a further improvement
on SPH in resolving discontinuity problems, and they call
this as a discontinuous SPH, or DSPH [55]. The idea of
DSPH one-dimensional space is briefed here as follows.
Examine the kernel approximation for any function f (x)
in the support domain of x i . The support domain is bounded
by a and b with a dimension of 2h, shown in Fig. 6. Assuming that function f (x) has an integrable discontinuity
at d in the support domain and that it is located in the right
Rearranging by combining some similar terms with some

transformations yields
b
f (x)Wi (x)dx
a

= f (x i )
a
Wi (x)dx
d
+ f (x i )
Wi (x)dx + [f (x k ) f (x i )]
(x x i )Wi (x)dx
a
b
[(x x k )f (x k ) (x x i )f (x i )]Wi (x)dx
+ r(h2 ).
(73)
Since the kernel W is assumed to be even, normalized, and

has a compact support, we have |x x k | h. In the above
equation, we assume that f (x) must exist and be bounded
in [a, d) (d, b]. In other words, the derivative of f (x) exists, and when x approaches d from both sides, it is bounded
within a finite limit. Hence, the last two terms (excluding
the residual term) on the RHS of the above equation can be
bounded by terms of order of h respectively. Therefore, the
above equation can be rewritten as:
b
f (x)Wi (x)dx
a
Fig. 6 Kernel approximations for a one-dimensional function with a

discontinuity at point d. The support domain of x i is bounded by a and
b with a dimension of 2h

= f (x i )
a
+ r(h).
Wi (x)dx + [f (x k ) f (x i )]
Wi (x)dx
d
(74)
41
This equation can be re-written as

b
f (x i ) =
f (x)Wi (x)dx
b
a Wi (x)dx
b

[f (x k ) f (x i )] d Wi (x)dx
+ r(h), (75)
b
a Wi (x)dx
a
which is the kernel approximation of the field function with

a discontinuity.
Similarly, the kernel approximation of the derivative is
obtained as
f (x i )
b
=
[f (x) f (x i )]Wix (x)dx

b
a (x x i )Wix (x)dx
b

[f (x k ) f (x i )] d Wix (x)dx
b
a (x x i )Wix (x)dx
b

[(x x k )f (x k ) (x x i )f (x i )]Wix (x)dx
+ d
b
a (x x i )Wix (x)dx
a
+ r(h).
(76)
The above two equations are the kernel approximations

of a field function (and its derivative) with a discontinuity
in the support domain. The kernel approximations at the
RHS of (75) and (76) consist of two parts. The first parts
are the same as those at the RHS of (67) and (69), which describe the kernel approximations of a field function and its
derivatives. It is the second parts in the big brackets that describe the behavior of the discontinuity. Neglecting the second parts of (75) and (76) should result in numerical errors
since the resultant kernel approximations in the presence of
a discontinuity are inconsistent. If the second terms are retained, the resultant kernel approximations are consistent up
to the first order.
The particle approximations for the field function and its
derivatives at particle i are a little bit different from the particle approximations in the conventional SPH and CSPM due
to the presence of the discontinuity. When the domain is discretized by particles, since two particles cannot be located
at the same position, the discontinuity should always be located between two particles. In deriving (75) and (76), since
the point x k is arbitrarily selected, it can be taken as the
particle that is nearest to and on the right hand side of the
discontinuity (Fig. 7).
The particle approximation of the discontinuous function
is
N mj
j =1 ( j )fj Wij
fi = N m j
j =1 ( j )Wij
Fig. 7 Particle approximations for a function with a discontinuity at

point d. In the process from the kernel approximation to particle approximation, an arbitrary point x k is associated with a particle k that
is the nearest particle on the right hand side of the discontinuity. The
total number of particles in the support domain of [a, b] is N
[f (x ) f (x )] N ( mj )W
k
i
ij
j =k j
.
N mj
j =1 ( j )Wij
(77)
And similarly, the particle approximation of the derivative

of the discontinuous function is
N mj
j =1 ( j )(f (x j ) f (x i ))i Wij
fi =
N mj
j =1 ( j )(x j x i )i Wij
[f (x ) f (x )] N ( mj ) W
k
i
i ij
j =k j
N mj
N
mj
j =k ( j
)[(x j x k )f (x k ) (x j x i )f (x i )]i Wij

.
N mj
(78)
Also, the particle approximations consist of two parts, the
first primary parts similar to those in the CSPM approximation and the second additional parts in the big brackets developed for the discontinuity treatment. It should be
noted that the summations of the numerators of the additional parts are only carried out for particles in the right portion of the support domain [d, b] shaded in Fig. 7. The coordinates of these particles are x j (j = k, k + 1, . . . , N ) that
satisfy (x k x j )(x k x i ) 0. Particles at the same side as
i do not contribute to the summation in the additional parts.
Otherwise, the additional part would become zero, and the
modification term would disappear.
Equations (75) and (76), (77) and (78) are the essential DSPH formulations for approximating a discontinuous
function and its first derivative. It is possible to extend the
DSPH to higher order derivatives, and multi-dimensional
space.
The different performances of the SPH, CSPM and
DSPH in resolving discontinuity have been demonstrated
by a number of numerical tests. One of them is the sod
tube problem, which has been studied by many researchers
for validation purposes [47, 91]. The shock-tube is a long
straight tube filled with gas, which is separated by a membrane into two parts of different pressures and densities. The
42
M.B. Liu, G.R. Liu
Fig. 8 Density profiles for the shock tube problem obtained using different versions of SPH formulation. DSPH gives the best result. CSPM
failed to capture the major shock physics (t = 0.2 s)
Fig. 9 Pressure profiles for the shock tube problem obtained using
different versions of SPH formulation
gas in each part is initially in an equilibrium state of constant

pressure, density and temperature. When the membrane is
taken away suddenly, a shock wave, a rarefaction wave and
a contact discontinuity will be produced. The shock wave
moves into the region with lower density gas; the rarefaction
wave travels into the region with higher density gas; while
the contact discontinuity forms near the center and travels
into the low-density region behind the shock. Exact solution
is available for comparison for this one-dimensional problem.
In this example, the cubic spline function is used as the
smoothing function to simulate this shock tube problem. The
initial conditions are same as those introduced by Sod [92]
x 0,
= 1,
p = 1,
x = 0.001875,
x > 0,
= 0.25,
p = 0.1795,
v = 0,
v = 0,
e = 2.5,
e = 1.795,
x = 0.0075
where , p, e, and v are the density, pressure, internal energy, and velocity of the gas, respectively. x is the particle
spacing. A total of 400 particles are deployed in the onedimensional problem domain. All particles have the same
mass of mi = 0.001875. 320 particles are evenly distributed
in the high-density region [0.6, 0.0], and 80 particles are
evenly distributed in the low-density region [0, 0.6]. The initial particle distribution is to obtain the required discontinuous density profile along the tube. The equation of state for
the ideal gas p = ( 1)e is employed in the simulation
with = 1.4.
In the simulation, the time step is 0.005 s and the simulation is carried out for 40 steps. In resolving the shock,
the Monaghan type artificial viscosity is used, which also
Fig. 10 Velocity profiles for the shock tube problem obtained using
serves to prevent unphysical penetration [49]. Figures 8, 9,

10 and 11 show, respectively, the density, pressure, velocity
and internal energy profiles.
It can be seen that the obtained results from the DSPH
agree well with the exact solution in the region [0.4, 0.4].
The shock is observed at around x = 0.3 within several
smoothing lengths. The rarefaction wave is located between
x = 0.3 and x = 0. The contact discontinuity is between
x = 0.1 and x = 0.2. Just as reported by Hernquist and Katz
[91] and Monaghan [47], the traditional SPH with artificial viscosity yields comparable results in capturing shock
physics. The DSPH gives slightly better results than the traditional SPH method though only 40 steps passed. If more
dynamic steps are involved, the better performance of the
DSPH is more evident [55].
43
By substituting the above smoothing function into (45), and

after some trivial transformation, the discretized form of
(45) can be written as
N xx j I
k
I =0 bI (x, h)
j =1 ( h ) x j = 1
k
N xx j I +1
b
(x,
h)
(
)
x
=
0
I
j
I =0
j =1
h
.
(80)
..
N xx j I +k
k
x j = 0
I =0 bI (x, h)
j =1 ( h )
The k + 1 coefficients bI (x, h) can then be determined by
solving the following matrix equation
m0 (x, h)
m (x, h)
Fig. 11 Energy profiles for the shock tube problem obtained using
..
mk (x, h)

m1 (x, h)
m2 (x, h)
..
.
mk+1 (x, h)

..
.

b0 (x, h)
mk (x, h)

m1+k (x, h)
b1 (x, h)
= .. ,
..
..
.
.

0
bk (x, h)
mk+k (x, h)

(81)
CSPM failed to captures the major shock physics. In this
example, CSPM performs poorer even than the traditional
SPH method. The poor performance of the CSPM in simulating shock waves should arise from the denominator that
in fact acts as a normalization factor. Revisiting (67), it is
found that the summation of the smoothing function around
the discontinuity region is far from the unity. When it is
used to normalize the summation of the density, the local
features of shock wave can be smoothed out, and hence the
shock physics can be concealed. In contrast, the traditional
SPH does not suffer from this smoothing effect, and the
DSPH benefits from the addition part in resolving the shock
physics. Their performance, therefore, should be better than
that of the CSPM. It is clear that the different performances
of CSPM and DSPH originate from the additional corrective
parts in DSPH that approximate the discontinuity.
4.5 A General Approach to Restore Particle Inconsistency
Liu et al. gave a general approach to restore particle inconsistency through reconstructing the smoothing function [6].
In general, a smoothing function can be written in the following form

x xj
W (x x j , h) = b0 (x, h) + b1 (x, h)
h
2

x xj
+
+ b2 (x, h)
h

k

x xj I
=
bI (x, h)
.
h
I =0
(79)
or
Mb = I ,
(82)
where
mk (x, h) =

N

x xj k
j =1
x j ,
(83)
M is a moment matrix, b is a vector of coefficients, I is a

vector of given constants.
After determining the coefficients bI (x, h), the smoothing function expressed in (79) can be calculated. The procedure ensures particle consistency to kth order. Therefore,
the particle consistency restoring process actually gives an
approach to construct some kind of smoothing function for
the SPH methods.
Comparing with the traditional smoothing function,
which is only dependent on the particle distance and applicable for all the particles, the consistency restored smoothing
function is particle-wise. It therefore depends on both the
distance and positions of the interacting particles. The costeffectiveness for this approach in constructing particle-wise
smoothing functions needs to be considered since it will
require additional CPU time to solve the particle-wise equation (79) for all the particles. Moreover, since all particles
are moving, the particle location is changing as well. Hence,
the particle-wise smoothing functions need to be computed
for every time step. Another problem is that, to solve (81),
the moment matrix M is required to be non-singular. Therefore, the particle distribution must satisfy certain conditions
to avoid singular momentum matrix. This implies that when
we enforce on consistency, we will face the stability problem shown as the bad-conditioned moment matrix in the
44
M.B. Liu, G.R. Liu
SPH settings. In contrast, in the original SPH method, particles can be arbitrary distributed, though the obtained results
may be less accurate.
As far as the approximation is concerned, restoring particle consistency is an improvement on the accuracy of the
particle approximation, provided that the moment matrix M
is not singular. However, it is noted that restoring the consistency in discrete form leads to some problems in simulating hydrodynamic problems. Firstly, the resultant smoothing function is negative in some parts of the region. Negative value of smoothing function can leads to unphysical
representation of some field variables, such as negative density, negative energy that can lead to a breakdown of the
entire computation. Secondly, the resultant smoothing function may not be monotonically decreasing with the increase
of the particle (node) distance. Moreover, the constructed
smoothing function may not be symmetric and using this
non-symmetric smoothing function violates the equal mutual interaction in physics.
4.6 Finite Particle Method
Considering the disadvantages of the above-mentioned particle inconsistency restoring approach in constructing a
point-wise smoothing function, Liu et al. devised another
particle consistency restoring approach, which retains the
conventional non-negative smoothing function instead of
reconstructing a new smoothing function [56, 57]. This approach has been termed as Finite Particle Method (FPM), in
which a set of basis functions can be used in the numerical
approximation.
Performing Taylor series expansion at a nearby point
x i = {xi , yi , zi } and retaining the second order derivatives,
a sufficiently smooth function f (x) at point x = {x, y, z}
can be expressed as follows
following equation

f (x)1 (x x i )dx

1 (x x i )dx + fi,
= fi
fi,
+
2

(x x i )1 (x x i )dx
(x x i )(x x i )1 (x x i )dx
+ r((x x i )3 ).
(88)
It is seen that the integration is carried out over the entire

problem space, and can be quite time-consuming. One usual
assumption is that a field variable at point x i is only strongly
influenced by the field variables at nearby points and that the
influence of the field variables at points far away from point
x i is very weak and hence can be neglected. Therefore, the
global integration can be converted into a local integration
by defining a support domain for point x i in which the field
variables at point x i can be determined. The shape of the
support domain can be conveniently taken as a circle (in 2D)
or a sphere (in 3D) with a radius of h, in which is a
constant scalar factor, and h is a length characterizing the
support domain. The function 1 (x x i ) is also limited to
the local support domain, and can be rewritten as 1 (x
x i , h).
Since the points distributed in the problem space are actually particles, each occupying individual lumped volume,
(88) can be numerically approximated by summation over
the particles surrounding point xi as follows
N
f (x j )1 (x j x i , h)Vj
j =1
= fi
N
1 (x j x i , h)Vj
j =1
f (x) = fi + (x x i )fi, +
(x x i )(x x i )
fi,
2!
+ r((x x i )3 ),
+ fi,
j =1
(84)
where , are the dimension indices repeated from 1 to 3

(or from x to z). r((x x i )3 ) is the remainder of the expansion. fi , fi, and fi, are defined as
fi = f (x i ),
(85)
fi, = f (x i ) = (f/x )i ,
(86)
fi, = f (x i ) = ( 2 f/x x )i .
(87)
Multiplying both sides of (84) with a function 1 (x

x i ) and integrating over the problem space can yield the
N

(x j x i )1 (x j x i , h)Vj
N
fi,
(x j x i )(x j x i )1 (x j x i , h)Vj ,
2
j =1
(89)
where N is the number of particles within the support domain of particle i. The remainder term r((x x i )3 ) in (88)
is omitted in (89) for the sake of conciseness. The summation over the particles is illustrated in Fig. 3.
Equation (89) can be further simplified as the following
equation at point x i
B1i = A1ki Fki ,
(90)
where
Fki = [ fi fi, fi, ]T ,

B1i =
45
N
f (x j )1 (x j x i , h)Vj ,
(91)
(92)
A1ki =
N
1 (x j x i , h)Vj
N
(x j x i )1 (x j x i , h)Vj
j =1
N
1
(x j x i )(x j x i )1 (x j x i , h)Vj .
2
j =1
(93)
Corresponding to 1, 2, and 3 dimensional cases, there are
one function value, 1, 2 and 3 first derivatives, and 1, 3 and 6
second derivatives that will be approximated. It is clear that
k in (91), (93) and (90) is 3, 6, and 10 respectively corresponding to 1, 2, and 3 dimensional cases. To calculate the
function value, the first and the second derivatives at x i , 2, 5,
and 9 other equations similar to (90) are required. Therefore
in 1, 2, and 3 dimensional cases, totally 3, 6, and 10 functions (M (x x i , h), M = 3, 6, or 10) are required in order
to approximate the function value, the first and second derivative. These functions form a set of basis functions used for
approximating the function value, its first and second derivatives. A conventional SPH smoothing function and its first
and second derivatives can form a set of basis function in the
FPM. For example, in 2D space, a smoothing function W ,
its two first order derivatives, W and W , and its three second order derivatives, W , W and W forms a set of 6
basis functions.
In summary, multiplying a set of basis functions on both
sides of (84), integrating over the problem domain, summing
over the nearest particles within the local support domain of
particle i, a set of matrix equation can be produced to approximate the function value as well as the first and second
derivatives at particle i. The matrix equations at particle i at
can be written as
BMi = AMki Fki
or B = AF
(94)
f (x j )M (x j x i , h)Vj ,
(95)
where
N

j =1
AMki =
N
1
(x j x i )(x j x i )M (x j x i , h)Vj .
2
(96)
j =1
BMi =
j =1
j =1
j =1
N

(x j x i )M (x j x i , h)Vj
N

j =1
M (x j x i , h)Vj
Equation (94) is the basis of the finite particle method and

can be used to approximate a function value and its derivatives for a field variable. It is seen that only if the coefficient
matrix A is not singular, can these M equations determine
a unique set of solutions at particle i for the M unknowns
in vector F . Solving the above pointwise matrix equations,
the function value as well as the first and second derivatives
at every particle can be simultaneously approximated. Note
that the conditioning of A matrix reflects the stability of the
FPM model.
Since the governing equations in CFD only involve the
first and second derivatives, only the derivatives up to the
second order are retained in (84). For problems in other
areas such as computational solid mechanics, high order
derivatives may be involved. If third or higher order derivatives are to be approximated, in expanding f (x) at x i , the
interested derivatives need to be retained in (84). To obtain the increased number of unknowns, more functions like
M (x x i , h) are necessary to complete the matrix equation (94). Therefore, except for the increased number of unknowns, increased number of basis functions, and therefore
more computational efforts, the solution procedure for the
interested unknowns is the same.
Comparing conventional SPH and FPM, it is clear that
both FPM and SPH are meshfree particle methods in which
particles with lumped volumes are used to represent the state
of a system. The particles form a frame for interpolation,
differencing or integration in a certain approximation. Both
FPM and SPH can be used as Lagrangian methods if allowing the particles to move in the problem space. However, the
difference between FPM and SPH is obvious.
1. FPM uses a set of basis function to approximate the function value and its derivatives, whereas SPH employs a
smoothing function and its derivatives to approximate
a function value and the corresponding derivatives. The
smoothing functions in SPH should have some special
properties as described in Sect. 3. However, the basis
functions in FPM are more general. Any set of functions
which do not lead to a singular coefficient matrix A can
be used as basis functions. Therefore the smoothing function and its specific derivatives actually can be one possible choice as a suitable set of basis functions.
2. SPH can be regarded as a special case of FPM, whereas
FPM is a generalized version of SPH with modifications.
46
3.
4.
5.
6.
M.B. Liu, G.R. Liu
In (84), if all the terms related to derivatives are neglected, multiplying both sides of (84) with the smoothing function W , and integrating over the problem space
can lead to the approximation in SPH. Summation over
the nearest particles within the support domain of a particle further produces the particle approximation of the
field variable at that particle (see (89)).
FPM should have better accuracy than SPH. Since no
derivative term is retained in (84), the SPH method actually is of first order accuracy. If a symmetric smoothing function is used, the terms related to the first order
derivatives are actually zero for the interior particles in
the problem domain. Therefore SPH is of second order
accuracy in interior parts. In contrast, since up to second
order derivatives are retained in the expansion process,
the accuracy of FPM is of third order. Moreover, if higher
derivatives are retained, better accuracy can be achieved.
FPM should have a better accuracy than SPH both for the
interior particles and boundary particles.
The accuracy of FPM is not sensitive by the selection of
smoothing length, and extremely irregular particle distribution. This has been demonstrated by Liu and his coworkers in testing one- and two-dimensional cases.
As the solution is to be obtained from solving the matrix
(94), a good matrix inversion algorithm is necessary to
prevent the co-efficient matrix to be negative.
The matrix (94) is to be solved at every particle, and
every time step. It can be more computationally expensive than the conventional SPH method.
One simpler example of FPM is to only consider the first

derivative in (84). Using the smoothing function and its first
derivatives as the basis functions, the following equations
can be obtained

f (x)Wi (x)dx

= fi

Wi (x)dx + fi,
(x x i )Wi (x)dx,
(97)
(x x i )Wi, (x)dx.
(98)
and

f (x)Wi, dx

= fi

Wi, dx + fi,
Again is the dimension index repeated from 1 to d. The

corresponding discrete forms for (97) and (98) are
N
fj Wij vj
j =1
= fi
N

j =1
N

Wij vj + fi,
(x j x i )Wij vj ,
j =1
(99)
Table 1 Kernel consistency of SPH, CSPM, and FPM

Interior domain
Boundary area
SPH
1st order
less than 0th order
CSPM
1st order
0th order
FPM
1st order
1st order
and
N
fj Wij, vj
j =1
= fi
N

j =1
Wij, vj + fi,
N

(x j x i )Wij, vj .
j =1
(100)
There are d + 1 equations for d + 1 unknowns (fi and fi, ).
Equations (99) and (100) are therefore complete for solving
with respect to fi and fi, , and the solutions for fi and fi,
are

fi
fi,

1
N
N
j =1 Wij vj
j =1 (x j x i )Wij vj
=
N
N
x )W
W
v
(x
v
j
ij,
j
j =1 ij,
j =1 j
i

N
j =1 fj Wij vj
.

(101)
N
f
W
v
j
ij,
j
j =1
In (99) and (100), the terms related to the function and the
first order derivatives are all retained, only the terms related
to second or high order derivatives are neglected. Therefore
the resultant particle approximations for a function and its
derivatives (see (101)) are able to exactly reproduce a constant and a linear function (C 0 and C 1 consistency). Hence
the algorithm shown in (101) actually restores the particle
consistency that conventional SPH method does not have.
Again, this particle consistency restoring approach is independent of the particle distribution (either regular or irregular), and the choices of the smoothing kernel and smoothing length. Another advantage is that this particle consistency restoring approach does not change the conventional
smoothing function and should be preferable in simulating
hydrodynamics.
In summary, the consistency of the conventional SPH
method, CSPM, and FPM (if using (101), rather than (94),
which is associated with higher order particle consistency
than (101)) are described in Tables 1 and 2.
Figure 12 shows the approximation results for a linear
function f (x, y) = x + y using discrete models of SPH
and FPM. The particles are uniformly distributed. It is seen

Table 2 Particle consistency of
SPH, CSPM, and FPM
47
Interior domain
Boundary area
Regular distribution
Irregular distribution
SPH
1st order
less than 0th order
less than 0th order
CSPM
1st order
0th order
0th order
FPM
1st order
1st order
1st order
Fig. 12 Approximation results for a linear function f (x, y) = x + y using (a) SPH and (b) FPM in a domain of [0, 1, 0, 1]
that the SPH results are associated with oscillation especially near the boundary, whereas the particle consistency
restoring approach can obtain much better results, which are
in very good agreement with the analytical solution. It is
worth noting that if the particles are randomly or irregularly
distributed, FPM can also obtain accurate results, whereas,
SPH results can be far away from the analytical solutions.
It is also worth noting that in this simple numerical test, the
approximation is carried out by obtaining numerical values
at current step from the existing numerical values at the particles at the previous step. There is no physical dynamics
involved in these two numerical performance studies. Only
the accuracy of the numerical approximations of the given
function is studied using different approximation methods. It
can be expected that numerical errors existing in the present
approximated results can propagate to the next approximation steps and can even be magnified. However, when these
methods are used to simulate a well-posed physical problem, the physics is governed by the conservation equations.
If the discrete model is stable, the error in the current step is
controlled by the consistency. Therefore, the errors should
be more or less within some level, and will not be magnified
from one step to the next step.
4.7 Consistency vs. Stability
We have seen in the two previous sections that a typical
dilemma exists for many numerical methods: consistency or
stability. For a given setting of a numerical model such

as a particle method, we may choose one over another, but
probably a balanced one can be difficult to choose. The
original SPH has clearly chosen the stability (and also flexibility) over the consistency, which gives the SPH a distinct feature of working well for many complicated problems with good efficiency, but less accuracy. It always tries
to deliver some reasonably good results for the price paid.
This seems to be a very practical choice for many practical
engineering problems, as should be regarded as an advantage of the SPH method. Attempts to improve the accuracy
of SPH via restoring the consistency can be helpful, provided that the stability and efficiency is not too much compromised.
The consistency restoring approaches, such as the FPM,
put more emphases on the consistency (hence hopefully accuracy), but the stability (the conditioning of the M or A matrix) can be in question for some types of problems. Hence
proper measures are needed to establish the stability ensuring the accuracy of the solution.
The question is that can we have both the consistency
and the stability at the same time? The answer is yes, provided we are willing to change the setting and pay the
price. The recently proposed GSM method [3235] is a typical example that guarantees both excellent stability for arbitrary grids and 2nd order accuracy. It uses also the gradient
smoothed technique, but in a very carefully designed fashion. However, the GSM is not a particle method any more;
48
M.B. Liu, G.R. Liu
the simplicity and efficiency features of the particle methods are lost. The GSM requires precise evaluation of the
integrals for carefully chosen types of smoothing domains,
and it works more like FVM. Therefore, the final question
depends on what we want and at what cost: choosing a numerical method should be closely related to the nature of the
problem, the requirement for the solutions and the resources
we have.
5 Special Topics
Fig. 13 SPH kernel and particle approximations for interior and

boundary particles. The boundary particles suffer from particle insufficiency, and no contribution comes from outside since there are no
particles beyond the boundary
5.1 Solid Boundary Treatment

Particle methods such as MD, DPD, and SPH are very different from the grid based numerical methods in treating solid
boundary, and care must taken in this regards. For finite element methods, proper boundary conditions can be well defined and properly imposed without affecting the stability of
the discrete model. In the finite difference method, because
the domain is relatively regular, proper ways to handling the
boundary conditions have already been developed. In either
FEM or FDM, the implementations of boundary conditions
(either Neumann boundary conditions or Dirichlet boundary conditions, or mixed boundary conditions) are generally
straightforward [2, 3, 5, 93]. In contrast, in the MD, DPD
and SPH, the implementation of the boundary conditions is
not as straightforward as in the grid based numerical models. This has been regarded historically as one weak point of
the particle methods [12, 15, 77, 79, 9499]. For molecular
dynamics and dissipative particle dynamics, it is common to
fill the solid obstacle areas with frozen particles. These fixed
particles can prevent the mobile particles from penetrating
the solid walls, and interact with the mobile particles with
proper interaction models. Different complex solid matrix
and solid-fluid interface physics can thus be modeled by appropriately deploying fixed particles plus a suitable interaction model with mobile particles. Preventing fluid particles
from penetrating solid walls by exerting a repulsive force
simulates the boundary condition with zero normal velocity. Different models of reflection or mirror of fluid particles
or interaction models of solid particles with fluid particles
can also be applied to exert slip or non-slip boundary condition [11, 40, 97, 100103]. In principal, the boundary treatment techniques in MD and DPD can also be used for SPH.
However, SPH is a continuum scale particle method, and
the field variables on boundaries such as pressure need to be
directly calculated. Therefore solid boundary treatment in
SPH is more difficult and numerical techniques for implementing solid wall boundary conditions are also more diversified. Since the invention of the SPH method, a lot of works
have been published addressing the treatments of solid wall
boundary conditions [94, 97, 104111].
As discussed in Sect. 4, the conventional SPH model suffers from particle inconsistency that is closely related to the
selection of smoothing function, smoothing length, and most
importantly the distribution of particles. An extreme case is
the particle distribution near the boundary area. Only particles inside the boundary contribute to the summation of
the particle interaction, and no contribution comes from the
outside since there are no particles beyond the boundary.
This one-sided contribution does not give correct solutions,
because on the solid surface, although the velocity is zero,
other field variables such as the density are not necessarily
zero.
Some improvements have been proposed to treat the
boundary condition. Monaghan used a line of ghost or virtual particles located right on the solid boundary to produce
a highly repulsive force to the particles near the boundary, and thus to prevent these particles from unphysically
penetrating through the boundary [112]. Campbell treated
the boundary conditions by including the residue boundary
terms in the integration by parts when estimating the original kernel integral involving gradients [113], as discussed
in Sect. 3 in deriving the constructing conditions for the
smoothing kernel function. Libersky and his co-workers introduced ghost or image (or imaginary) particles to reflect
a symmetrical surface boundary condition with opposite velocity on the reflecting image particles [114]. Later, Randles and Libersky proposed a more general treatment of the
boundary condition by assigning the same boundary value
of a field variable to all the ghost particles, and then interpolating smoothly the specified boundary ghost particle value
and the calculated values of the interior particles [52].
In general, it is feasible to use virtual particles to implement the solid boundary conditions. These virtual particles can be allocated on and outside the boundary. These
virtual particles can be classified into two categories, virtual particles that are located right on the solid boundary
(Type I) and virtual particles that are outside the boundary
(Type II) (see Fig. 14). Type I virtual particles are similar to
what Monaghan used [112], and they are used to exert a repulsive boundary force on approaching fluid (real) particles
49
cles N N (i) that are within its influencing area of hi can be

categorized into three subsets (see Fig. 14)
(a) I (i): all the interior particles that are the neighbors of i
(real particles);
(b) B(i): all the boundary particles that are the neighbors
of i (virtual particles of type I);
(c) E(i): all the exterior particles that are the neighbors of i
(virtual particles of type II).
Fig. 14 Schematic illustration of flow region, solid boundary, real
fluid particles and virtual particles
to prevent the interior particles from penetrating the boundary. The penalty force is calculated using a similar approach
for calculating the molecular force of Lennard-Jones form
[115]. If a Type I virtual particle is the neighboring particle
of a real particle that is approaching the boundary, a force is
applied pairwisely along the centerline of these two particles
as follows

x
D[( rrij0 )n1 ( rrij0 )n2 ] 2ij , ( rrij0 ) 1,
rij
PBij =
(102)
0,
( rrij0 ) > 1.
where the parameters n1 and n2 are empirical parameters,
and are usually taken as 12 and 4 respectively. D is a problem dependent parameter, and should be chosen to be in the
same scale as the square of the largest velocity. The cutoff
distance r0 is important in the simulation. If it is too large,
some particles may feel the repulsive force from the virtual particles in the initial distribution, thus leads to initial
disturbance and even blowup of particle positions. If it is
too small, the real particles may have already penetrated the
boundary before feeling the influence of the repulsive force.
In most practices, r0 is usually selected approximately close
to the initial particle spacing. To more effectively prevent
fluid particles from penetrating the solid wall, Type I virtual
particles are usually more densely distributed at the solid
boundary than the fluid particles (Fig. 14).
Type II virtual particles are usually obtained through reflecting the real fluid particles along the solid boundary.
Scalar properties of the virtual particles can be taken as the
same as the reflected real particles, while velocity of the virtual particles are usually taken as opposite of the reflected
counterparts. As such the non-slip boundary condition is implemented. Due to the compactness of the smoothing kernel function, only the real particles that are within h to
the solid boundary are necessary to be reflected to form
Type II virtual particles. Hence there are only several layers of Type II virtual particles outside the solid boundary.
For example, if the cubic spline smoothing function is used,
there are two layers of Type II virtual particles (if the real
particles are regularly distributed). Therefore, for a fluid particle i close to the solid boundary, all the neighboring parti-
These two types of virtual particles are specially marked

for contribution in the later summation on the real particles.
Type I virtual particles exert repulsive force on an approaching particle, but they usually do not contribute in the particle approximation. Type II virtual particles take part in the
kernel and the particle approximations of the real particles.
It is clear that by using type II virtual particles, the onesided particle distribution problem in the conventional SPH
method can be somewhat remedied. In different implementations, the position and physical variables for type II virtual
particles can be fixed (determined once a set of virtual particles are obtained), or can adapt in the simulation process
with the evolution of the corresponding counterpart real particles. In general, adapting type II virtual particles can result in better results. However, the adapting process may be
time-consuming, and may be complicated especially when
the solid wall is not able to be defined in a simple curve.
In many circumstances, the computational domain can
be complicated, and can contain solid obstacles of arbitrary
shapes. Generating type II virtual particles may not be an
easy task. For such cases, a feasible way is to deploy particles into the computational domain either by placing all the
particles in some kind of lattice, or injecting particles into
the computational domain and then running SPH simulation
for the entire system to reach equilibrium, just as in MD and
DPD simulation. The particles outside of the boundary (or
within the solid obstacles) can be used as type II virtual particles. This approach has some advantages since the system
starts from a quiet or equilibrium state of the same number
density without unwanted noise. For implementing the nonslip boundary condition, these virtual particles can move
with the velocity as the obstacles to which they are attached.
For better accuracy, the velocity of the virtual particles can
be obtained from extrapolating the velocity of a neighboring
fluid particle, as illustrated in Fig. 15. Assuming the normal
distance of a virtual particle to the solid boundary is dv , and
the normal distance of a fluid particle to the solid boundary is df , with non-slip boundary condition (velocity on the
boundary is zero), the velocity of the virtual particle under
consideration should be v v = (dv /df )v f . Hence the velocity difference between the fluid and virtual particles is
v f v = (1 + ddfv )v f , which can be used in calculating viscous forces as addressed in Sect. 5. Note that the positions
of the virtual particles do not need to evolve. Also if the
50
Fig. 15 A general approach to construct type II virtual particles and

assign velocities. Only one virtual particle is plotted in the obstacle
region for illustration
obstacle is a moving object, the velocity of the fluid particle should be taken as the relative velocity to the obstacle.
To prevent numerical singularities caused by a fluid particle
closely approaching the solid boundary (df 0), a safety
parameter can be used in calculating the velocity difference between the fluid and virtual particles as v f v = v f ,
where = min(max , 1 + ddfv ). max is an empirical parameter, it is reported that taking 1.5 can result in good results
[111].
5.2 Representation of Solid Grains
To model the geometries of a complicated solid matrix,
the entire computational domain can be discretized using
a shadow grid. The grid cells are labeled 0 for regions
occupied by pore spaces and 1 for solid-filled regions
(Fig. 16a). This simple identification of fluid and solid cells
can be used to represent any arbitrary complex geometry.
The unit vectors normal to the solid-fluid interfaces define
the local orientation of the fluid-solid interface and can be
obtained by simply calculating the surface gradient from the
Fig. 16 Schematic illustration
of the treatment of solid
obstacles. (a) The cells in the
entire computational domain are
first labeled, 0 for fluid (void)
cells and 1 for solid (obstacle)
cells. (b) The SPH particles in
the obstacle cells are frozen.
(c) Only the frozen particles that
are close to the fluid cells
(within 1 cut-off distance) are
retained as boundary particles
M.B. Liu, G.R. Liu
indicator numbers (0 for liquid regions and 1 for solid

regions). As discussed in Sect. 6.1 in constructing type II virtual particles, we can place particles in the entire computational domain regularly in some kind of lattice according to a
designated number density. It is also applicable to inject particles into the computational domain, and run SPH simulation until the system arrives at an equilibrium. At the beginning of each simulation, the particles are initialized and positioned randomly within the entire computational domain
until a pre-defined particle number density is reached, and
the system is then run to equilibrium. The particles within
the solid cells (marked as 1) are then frozen (their positions are fixed) to represent the solid grains (Fig. 16b). The
solid grains in obstacles can occupy a considerable fraction
of the entire computational domain, and hence the number
of frozen particles representing the solids can be very large.
Most of the frozen particles inside the solid grains are more
than 1 cut-off distance away from the adjacent fluid cells.
These particles do not contribute to the solid-fluid interactions and consequently they have no influence on the movement of the flow particles within the fluid cells. Therefore,
only the frozen particles that are within 1 cut-off distance
from the solid-fluid interface are retained as boundary particles (Fig. 16c), and the rest of the particles further inside the
solid grains are removed from the model domain.
Figures 17 and 18 show two example of solid obstacle
representation. In Fig. 17, the flow region is a porous media system, in which the circular obstacles occupy a great
amount of the computational domain. Figure 17a shows the
original porous media system where the circular obstacles
are filled with (solid) particles. Figure 17b shows the porous
media system for later simulation where only the (solid) particles within 1 cut-off distance away from the adjacent fluid
cells are retained. In Fig. 18, the flow region is a fractured
network system and the obstacles also occupy most of the
computational domain (Fig. 18a). After treatment, only the
fracture particles within 1 cut-off distance away from the adjacent fluid cells are retained. For both cases, it is convenient
and flexible to treat the complex geometries by identifying
the flow and obstacle regions separately with the indicator
51
5.3 Material Interface Treatment
Fig. 17 An example of solid obstacle representation in a porous media

flow region with circular obstacles. (a) Original porous media system.
(b) Solid particle distribution for later simulation
Fig. 18 An example of solid obstacle representation in a fracture network in which most of computational domain is full of obstacles.
(a) Original fractured network system. (b) Solid particle distribution
for later simulation
numbers. Also since only the solid particles that are within
1 cut-off distance away from the adjacent fluid cells are retained for later simulation, the number of particles can be
greatly reduced, and therefore the computational efficiency
can be significantly improved.
The above approach of deploying frozen particles to represent solid grains and interact with flow particles is applicable to particle methods such as MD, DPD, and SPH,
and it is especially effective for complex boundaries [99,
116118]. For problems with comparatively simpler geometries, some other approaches are commonly used. These approaches may also involve particles fixed on the solid boundaries with or without equilibrium process, and may involve
non-fixed boundary particles updated with the neighboring
flow particles. Although fixed boundary particles are more
frequently used in SPH literature, non-fixed boundary particles updated with the neighboring flow particles may produce better results [6, 110, 111, 119].
For modeling systems with multiple materials, an effective

algorithm is required for efficiently describing the material
interface dynamics. It is somewhat different from the grid
based numerical methods in correctly detecting contact and
effectively treating it in particle methods like SPH. Traditionally SPH uses particles to represent the state of a system under consideration. One point is that the SPH particles are interpolation points on which mass of each particle has been lumped. An SPH particle has a finite mass and
density, and hence can have a finite volume. However, in
contrast to FEM, FDM and FVM, in which a cell or an element has a definite shape, an SPH particle has never assumed as a definite shape. On the other hand, the shape of
a particle can also be associated with the influencing area
(or support domain, that is kh) of the particle. Generally an
SPH particle is assumed to be circular (or spherical) when
using a smoothing kernel function with a scalar smoothing
length (Fig. 19a). In some cases, an SPH particle can be elliptical (or ellipsoidal) when a tensor smoothing length is
used [120124] (Fig. 19b). In few simulations, it even can
be taken as a square or a cube [125] (Fig. 19c). In the SPH
simulation, the mass of a particle is usually assumed to be
constant, and updating the density means change of the volume. If the particle takes a definite shape such as a circular,
then the size of the particle also changes. Therefore different shapes and sizes of the SPH particles increase the difficulty in implementing the contact algorithm. During the
course of SPH development, different contact algorithms
have emerged [126142]. For simplicity, only SPH models with circular or spherical shaped particles are discussed
here. The SPH particles are assumed to be hard-spheres
rather than soft balls, as models with soft balls also increase
the difficulty in implementing the contact algorithm.
In the SPH particle models, if particles from different
materials are not bonded together, particles from one material can influence, and also be influenced by the particles
from the other material, according to conventional SPH algorithm. Figure 20 shows an example of the material interface at the initial stage, and at a later stage in the evolution.
Initially the SPH particles are distributed regularly, and there
is no intersection and deformation of particles. With the SPH
evolution, particles from different materials influence each
other, and the conventional SPH algorithm can introduce errors since particles from material A influence the particles
from material B in calculating the strain and strain rate, and
vice versa. This interaction between particles from different materials can introduce shear and tensile stress, which
prevent sliding and separation of different materials. Special
algorithms are required if an SPH model is used for simulation of problems with sliding and separation.
52
M.B. Liu, G.R. Liu
Fig. 19 Possible shapes of an

SPH particle. (a) Circular shape
corresponding to a scalar
smoothing length, (b) elliptical
shape corresponding to a tensor
smoothing length at different
directions, (c) square shape
whereas the smoothing length
can also be a tensor
Fig. 20 Illustration of material interface with SPH particle model.

(a) Initial SPH particle distribution, and (b) SPH particle evolution. Interaction between particles from different materials can introduce shear
and tensile stress and can prohibit sliding and separation of different
materials
Fig. 21 An illustration of particle and particle contact
pairwisely on the two approaching particles along the centerline of the two particles

Major steps to be considered in a contact algorithm in the

SPH simulation include (1) identifying the boundary, (2) detecting the contact and (3) applying repulsive contact forces.
Boundary (or interface) has usually been considered to be a
surface one half of the local smoothing length away from the
boundary SPH particles [52, 127]. A convenient way to detect particle-particle contact is to check the proximity of two
approaching particles from different materials, and compare
with the summation of the smoothing lengths of the particles. As shown in Fig. 21, the contact and penetration between two particles can be identified by
pe =
khi + khj
1,
rij
(103)
where rij is the distance between particle i and particle j .

After detection of the contact and penetration, a restoring
contact force needs to be applied along the centerline of the
two particles (position vector of particle i and particle j ).
There are different choices for the restoring contact force,
which is usually a function of penetration. One possible
choice is the penalty force in Lennard-Jones form. Liu et al.
had used this penalty force for simulating underwater explosion when dealing with different material interface [143]. In
their work, the Lennard-Jones form penalty force is applied
f ij =
n1 pen2 )
p(pe
0,
x ij
rij2
, pe 1,
(104)
pe < 1,
where the parameters p,

n1 , n2 are usually taken as 105 ,
6 and 4 respectively. In fact, these parameters can be adjusted to suit the needs of different problems, and the involved driven force in underwater explosions can vary in a
very big range [144]. The application of the penalty force in
combination with the summation between the interface particles has been shown to well prevent the unphysical penetration in the simulation of underwater explosion, though
numerical oscillations near the vicinity of the interface can
still exist.
It is reasonable that when the penetration of two contacting particle increases (or particle distance decreases),
the magnitude of the restoring force should also increase.
Monaghan [145] used the following form of repulsive force,
which is controlled by the kernel and the global distribution
of the particles
fij =
Wij
,
W (p)
(105)
where p denotes the average particle spacing in the neighborhood of particle i (proportional to smoothing length h).
53
Vignjevic et al. [146] used the following repulsive force

for treating possible contact and penetration between solid
particles
f i = Kf
N

mj Wij
,
i j x
j =1
i
(106)
where Kf is a user defined empirical parameter. In their

work, the force was applied to boundary particles that are
within the influencing area kh of each other. This repulsive force introduces the derivatives of the smoothing kernel, which in general reaches the maximum at a certain location between the origin and the boundary of the influencing
area, and vanishes at the origin and the boundary point (see
Fig. 4). Considering that decreasing the inter-particle distance means bigger repulsive force, using the kernel itself
rather than its derivative physically would be more reasonable. For example, one possible choice can be
f i = Kf
N

mj
xij
Wij 2 .
i j
rij
Fig. 22 Schematic illustration of the tensile instability in a

two-dimensional space. (a) Initial particle distribution; (b) Particle distribution in a certain stage under the effect of tensile instability
(107)
j =1
5.4 Tensile Instability

When using the SPH method for hydrodynamics with material strength, one numerical problem called tensile instability [45, 62, 63, 74, 145, 147152] may arise. The tensile
instability is the situation that when particles are under tensile stress state, the motion of the particle becomes unstable. It could result in particle clumping or even complete
blowup in the computation. Swegle et al. demonstrated the
tensile instability effects in a two-dimensional space as illustrated in Fig. 22 [45]. It is a simple test in which the particles representing the object initially are distributed regularly and under uniform initial stress. The boundary particles are fixed (two layers) and the interior particles are free
to move due to some kind of external perturbation. If applying a very small velocity perturbation on a center particle,
theoretically, it will take a long time for a particle under consideration to move over single particle spacing. However, in
an SPH simulation with a common cubic spline smoothing
kernel function, if the stress is tensile, after a number of time
steps (while the number of steps is not sufficiently to move
the particles even a single particle spacing), the interior particles can clumped together unphysically, and form a void.
It shows that the particle system exhibits a numerical instability from tensile stress.
According to Swegle, the tensile instability depends neither on the artificial viscosity, no on the time integration
scheme. It is closely related to the selection of smoothing
kernel function. In a one-dimensional von Neumann stability analysis, Swegle et al. gave a criterion for being stable or
instable in terms of the stress state and the second derivative
Fig. 23 Schematic illustration of tensile instability with the cubic

spline function and its first and second derivatives
of the smoothing function, i.e., a sufficient condition for the

unstable growth is
W > 0,
(108)
where W is the second derivative of the smoothing function.

In the SPH method, the cubic spline smoothing function
(illustrated in Fig. 4) is mostly commonly used. The initial smoothing length is usually set to be equal to the particle spacing. Under such circumstances, the first nearest
neighbor particles are located at r/ h = 1; and the next nearest neighbor particles are at r/ h = 2. As can be seen from
Fig. 23, the second derivatives of the cubic spline function
from r/ h = 1 to r/ h = 2 are always positive. Therefore it is
expected that, according to (108), the SPH method with the
cubic spline function would be stable in a compressed state
but could be unstable in a tensile state in this region.
54
M.B. Liu, G.R. Liu
Fig. 24 Two possible setups of

velocity and stress particles
Several remedies have been proposed to improve or

avoid such tensile instability. Morris suggested using special
smoothing functions since the tensile instability is closely
related to the second order derivative of the smoothing function [46]. Though successful in some cases, they do not always yield satisfactory results generally. Chen and his coworkers proposed the corrective smoothed particle method
(CSPM), which was reported to improve the tensile instability [89]. Recently Monaghan and his colleagues proposed an
artificial force to stabilize the computation [145, 153].
The basic reason of tensile instability is that the SPH
method is essentially a collocation method, in which the particle approximations are conducted ONLY over the particles
that represent the entire system. This leads to insufficient
sampling points for establishing equations, and can result
in numerical instability problem [4]. The situation is very
much similar to the so-called node integration in the implementation of the element free Galerkin method (EFG)
[154]. In the EFG method, the instability is restored by
adding stabilization terms in the Galerkin weak form. Based
on this analysis, we need somehow bring in more information from other points, in addition to these particles.
One of such a method is to make use of the information at additional points in the support domain, rather than
use only these collocation particles. Dyka et al. first introduced additional stress points other than the normal particles
in a one-dimensional algorithm aimed at removing the tensile instability in SPH [62, 63]. The stress points were also
shown to be stable in tension and contributed considerably
to the accuracy in wave propagation problems. Later, this
approach has been further extended to multi-dimensional
space by staggering the SPH particles using stress points
so that there are essentially an equal number in each set of
points [64, 108, 151]. Basically in this approach, two sets of
particles are used. One set of SPH particles carry velocity,
and are referred to as velocity particles. The other companion set of particles carry all required field variables except for the velocity, and are referred to as stress particles.
Figure 24 shows two possible setups of velocity and stress
particles.
Randles and Libersky pointed out that, the tensile instability for problems involving material strength generally is
latent. The growth rate of damages in solid continuum models is often much faster than the grow rate of the tensile instability [64].
Except for problems with material strength which can experience tensile instability, fluid mechanics problems sometimes can also meet tensile instability. Melean et al. had
showed the tensile instability in a formation of viscous drop
[155], and the instability can be removed by using the artificial stress proposed by Monaghan [145, 153].
6 Applications
The original applications of the SPH method is in astrophysical phenomena, such as the simulations of binary stars and
stellar collisions [49, 156, 157], supernova [158, 159], collapse as well as the formation of galaxies [160, 161], coalescence of black holes with neutron stars [162, 163], single
and multiple detonation of white dwarfs [164], and even the
evolution of the universe [165]. It also has been extended to
a vast range of problems in both fluid and solid mechanics
because of the strong ability to incorporate complex physics
into the SPH formulations [6]. The applications of SPH to
many other engineering applications include
multi-phase flows [119, 166174],
coastal hydrodynamics including water wave impact, dam
break, sloshing and overtopping [112, 175200],
environmental and geophysical flows including flood and
river dynamics, landslide, flow in fractures and porous
media, seepage, soil mechanics and mudflow [117, 201
221],
heat and/or mass conduction [53, 222227],
ice and cohesive grains [228234],
microfluidics and/or micro drop dynamics [123, 152, 169,
235246],
high explosive detonation and explosion [55, 125, 143,
247256],
underwater explosions and water mitigation [125, 143,

247249, 255],
elastic and/or plastic flow [52, 114, 257259],
fracture of brittle solids [260],
metal forming and high pressure die casting [68, 261
270],
magneto-hydrodynamics and magnetic field simulation
[271282],
problems with fluid-solid interactions [176, 196, 210,
283288], and
many other problems like blood flow [289292], traffic
flow [293].
In this review, we will discuss applications of SPH to the
following areas:
high strain hydrodynamics with material strength,
high explosive detonation and explosions, and underwater
explosion,
microfluidics and micro drop dynamics,
coast hydrodynamics and offshore engineering, and
environmental and geophysical flows.
6.1 High Strain Hydrodynamics with Material Strength
High strain hydrodynamics is generally characterized by
the presence of shock waves, intense localized materials
response and impulsive loadings. Numerical simulation of
high strain hydrodynamics with material strength such as
high velocity impact (HVI) and penetrations is one of the
formidable but attractive tasks in computational solid mechanics. Most of the wave propagation hydro-codes use traditional grid based methods such as finite difference methods and finite element methods to simulate high strain hydrodynamics. Some of them are associated with advanced
features which attempt to combine the best advantages
of FDM and FEM. Examples include arbitrary LagrangeEulerian (ALE) coupling and coupling Eulerian-Lagrangian
(CEL). Though many successful achievements have been
made using these methods, some numerical difficulties still
exist. These numerical difficulties generally arise from large
deformations, large inhomogeneities, and moving interfaces, free or movable boundaries [256, 294297].
SPH method has been intensively used for simulating
such phenomena of high strain hydrodynamics with material strength, due to its special features of meshfree, Lagrangian and particle nature. Many researchers have investigated impact (either high velocity or not) and penetration
using SPH method [50, 51, 86, 129, 132, 134, 139, 140,
298335]. The applications are mainly in defense industry,
and typical numerical examples range from Taylor bar impacting on a solid wall and spherical projectile/bumper collision/penetration [114, 330], projectile impact and penetration and an aircraft fuselage nose impacting a rigid target
55
[300], to bird strike impacting on an aircraft Wing [335],

and even to a high-quality high-fidelity visualization of the
September 11 attack on the world trade center [328].
Libersky and his co-workers have carried out the pioneering work of applying the SPH method to high strain hydrodynamic problems including hyper velocity impact (HVI),
fracture and fragmentation [52, 114, 258]. The group from
the Applied Physics Division of Los Alamos National Laboratory (LANL) has modeled high velocity impacts ranging
from the very small size (femtogram-scale projectiles that
have been accelerated by a Van de Graaf machine to create craters in various types of targets) to the very large size
(Shoemaker-Levy comet impact with the planet Jupiter).
Johnson et al. [51, 314] proposed a normalized smoothing function (NSF) for axisymmetric problems based on the
condition of uniform strain rate, and have made outstanding contributions in the application of SPH to impact problems. Attaway, Pampliton and Swegle et al. have worked in
coupling the SPH processor with a transient-dynamics FEM
code, PRONTO, in which high-strain areas that typically
tangle or break conventional finite element meshes are resolved using the SPH method [126, 336].
Zhou and Liu have revisited the Taylor-Bar-Impact with
the focus on the variation of results corresponding to the
different model parameters which represent varied SPH
implementation in a series of three-dimensional computational simulations, and have provided informative data on
appropriate SPH implementation options [259]. They also
have investigated normal and oblique hypervelocity impacts
of a sphere on the thin plate, producing large deformation of structures. Figure 25 shows the snapshots of threedimensional simulation of normal hypervelocity impact of
a sphere on the thin plate at t = 0, 10, 15, and 20 s. Figure 26 shows the top view (a) and three dimensional view
(b) of oblique hypervelocity impact of a sphere on the thin
plate at t = 20 s with a striking angle of 30o . They observed
that for oblique hypervelocity impacts of a sphere on the thin
plate, the shape of the crater is no longer in circle, and the
shape of the debris cloud changes with different striking angles while the ratio of the debris cloud (length over width)
remains approximately 1.3.
The standard SPH method uses an isotropic smoothing kernel which is characterized by a scalar smoothing
length. One of the problems associated with the standard
SPH is that the isotropic kernel of SPH can be seriously
mismatched to the anisotropic volume change that generally
occur in many problems. To closely match the anisotropic
volume changes, an anisotropic smoothing kernel which can
be characterized by a matrix (in two-dimensional space) or
tensor (in three-dimensional space) smoothing length can be
efficacious. This leads to the development of the adaptive
smoothed particle hydrodynamics (adaptive SPH or ASPH)
in which the smoothing length can be adapted with the volume changes or other dimension-dependent features. The
56
M.B. Liu, G.R. Liu
Fig. 25 Snapshots of
three-dimensional simulation of
normal hypervelocity impact of
a sphere on the thin plate at
t = 0, 10, 15, and 20 s (from
[259])
Fig. 26 Top view (a) and three

dimensional view (b) of oblique
hypervelocity impact of a sphere
on the thin plate at t = 20 s
with a striking angle of 30
(from [259])
idea of using anisotropic kernel with SPH dates back to

Bicknell and Gingold [337]. Fulbright et al. also presented
three-dimensional SPH designed to model systems dominated by deformation along a preferential axis using spheroidal kernels [338]. Later on Shapiro et al. systematically
introduced anisotropic kernels, tensor smoothing and shocktracking to SPH to create ASPH [121]. Owen et al. presented an alternative formulation of the ASPH algorithm
for evolving anisotropic smoothing kernels [120]. Except for

problems with anisotropic deformations, the concept of elliptical kernel has also been applied to channel flows with
very large length width ratio for saving computational efforts [123]. The numerical results presented in the references
further demonstrated that ASPH has significantly better performance than the standard SPH in terms of resolving ability
for a wide range of problems.
Fig. 27 Initial particle distribution in the vicinity of the contact area

for both the SPH and ASPH simulation of a cylinder impacting on a
plate (from [124])
Liu and his co-workers have developed an adaptive

smoothed particle hydrodynamics method (ASPH) for high
strain Lagrangian hydrodynamics with material strength
[124, 339]. In ASPH, the isotropic kernel in the standard
SPH is replaced with an anisotropic kernel whose axes
evolve automatically to follow the mean particle spacing
as it varies in time, space, and direction around each particle. Except for the features inherited from the standard
SPH, ASPH can capture dimension-dependent features such
as anisotropic deformations with a more generalized elliptical or ellipsoidal influence domain. A series of comparative
studies show that ASPH has better accuracy than the standard SPH when being used for high strain hydrodynamic
problems with inherent anisotropic deformations.
A typical example is the simulation of an aluminum (Al)
cylinder on an Al thin plate, as shown in Fig. 27. In the simulation, the cylinder is of 1.0 cm diameter. The plate is 0.4 cm
thick. The plate length is of 10 cm. The particles are all initialized as squares of 0.02 cm in side dimensions. The particles in the cylinder, which is an infinite cylinder in plane
symmetry, are arranged in circumferential rings as this gives
a most realistic representation of the geometry. The particles in the plate are arranged in a rectangular Cartesian array. There are 500 particles along the length and 20 particles
along the thickness of the plate. There are 1956 particles in
the cylinder and 10000 particles in the plate, for a total of
11956 particles. The cylinder is initially in contact with the
center of the plate. The problems are run with the plate free
of constraints. The impact speed of the cylinder is 6180 m/s.
The problems are run to 20 s after the impact. The environmental and initial temperature of the cylinder and plate
are set to 0C. Figure 28 shows the particle distributions obtained by using SPH at 0, 5, 15, and 20 s respectively. Fig-
57
Fig. 28 Particle distributions obtained using the SPH method for simulating an Al cylinder penetrating an Al plate at 0, 5, 15, and 20 s
respectively (from [124])
Fig. 29 Particle distributions obtained using ASPH for simulating an

Al cylinder penetrating an Al plate at 0, 5, 15, and 20 s respectively
(from [124])
ure 29 shows the particle distributions obtained using ASPH

at the same stages. In general, both SPH and ASPH can get
agreeable results with the experimental observations [340].
The symmetry of the problem is well preserved. Figure 30
shows the close-up view of the particle distributions near
the penetrated edge of the plate obtained using SPH (a) and
ASPH (b) at 20 s. In comparison with the SPH simulation,
the orientation and anisotropy of the deformation of the particles can be clearly seen in the ASPH simulation.
6.2 Detonation, Explosion and Underwater Explosion
The explosion of a high explosive (HE) charge can rapidly
convert the original explosive charge into gaseous products
58
M.B. Liu, G.R. Liu
Fig. 30 Close-up view of the

particle distributions near the
penetrated edge of the plate
obtained using SPH (a) and
ASPH (b) at 20 s (from [124])
with extremely high pressure through a chemical reaction

process. The high pressure can lead to damages to nearby
personnel and structures. A typical HE explosion consists of
the detonation process through the HE at a constant detonation velocity and the later expansion process of the gaseous
products to the surrounding medium (Fig. 31). The detonation process is accompanied by the propagation of the
reactive wave that advances through the explosive with a
constant velocity related to the particular type of explosive concerned. In a steady state detonation process, the reaction rate is essentially infinite and the chemical equilibrium is attained. After the completion of the detonation, the
detonation-produced explosive gas expands outwards. This
gas expansion generally involves moving material interfaces
if the explosive is surrounded by outside medium or free surfaces if the explosion occurs in the vacuum.
When the detonation and explosion happen underwater,
the explosion physics can be more complicated. The underwater explosion (UNDEX) [341, 342] produced by the detonation of a submerged high explosive poses a serious threat
to the integrity of nearby structures. Issues related to the underwater explosion include (1) the detonation process of the
high explosive (HE) charge, (2) the expansion process of the
detonation-produced explosive gas into the surrounding water, and (3) the interaction of underwater shocks with the
nearby structures. In the detonation process, the high explosive is converted into gaseous products at very high temperature and pressure through a violent chemical reaction, occurring with extreme rapidity and releasing a great deal of heat.
The propagation of the detonation wave through an explosive is so rapid that the gaseous products directly behind the
wave front are not in pressure equilibrium with the gas further behind the wave front. As the detonation wave reaches
the interface between the explosive and the surrounding water, a high pressure shock wave of step exponential type is
transmitted to propagate through the water, followed by a
series of bubble pulsation associated with the repeating expansion and contraction of the bubble of the explosive gas.
In the entire process of the underwater explosions, some
special features such as large deformations, large inhomo-
Fig. 31 Detonation of a 1D high explosive. The reaction end plane is

an interface of the pressurized high explosive charge and the explosive
gas produced in the detonation process
geneities, moving material interfaces, deformable boundaries, and free surfaces usually exist.
Theoretical solutions to the detonation and explosion of
high explosive, and underwater explosion are only limited
to some simple cases. Experimental studies need to resort to
dangerous and expensive firing trials, and sometimes certain physical phenomena related to the explosions cannot
be scaled in a practical experimental setup. Recently, more
and more analyses of detonation, explosion and underwater explosions are based on numerical simulations with the
advancement of the computer hardware and computational
techniques [8, 343, 344]. However, numerical simulations of
the high energy phenomena are generally very difficult for
the conventional grid based numerical methods. First, during the detonation process in the explosion, a very thin reaction zone divides the domain into two inhomogeneous parts
and produces large deformations. Second, in the expansion
process, there are free surfaces and moving interfaces involved. It is even more difficult to simulate underwater explosion, as underwater explosion phenomena are subject to a
number of physical laws and properties, including the physical conditions at the interface of the explosive gas and the
surrounding water. The surrounding water has such properties as large density that is approximately 1000 times of
the air density, low compressibility, and large sound speed.
Owing to the dynamic properties of the water (especially in

the regions surrounding the explosive gas), the pressures are
generally very high and the wave velocities are dependent
on the magnitude of the pressure and the displacement of
the water as it progresses. These complications for waves of
finite amplitude are expressed in much more difficult mathematical statements than those which suffice to explain the
propagation of small amplitude waves whose velocities are
practically independent on the magnitude of the pressure.
Traditional Lagrangian techniques such as the finite element methods are capable of capturing the history of the
detonation and explosion events associated with each material. It is, however, difficult to apply practically, since the
severely distorted mesh may result in very inefficient small
time step, and may even lead to the breakdown of the computation. Traditional Eulerian techniques, such as the finite
difference methods or finite volume methods, can well resolve the problem due to the large deformations in the global
motions, but it is very difficult to analyze the details of the
flow because of the lack of history and the smearing of
information as the mass moves through the fixed-in-space
Eulerian mesh [256]. Considering the difficulties of grid
based numerical models, the meshfree particle methods can
be a good alternate for simulating detonation, explosion and
underwater explosions [125, 136, 143, 247, 248, 252, 255,
256, 345350].
Liu and his co-workers have conducted a series of original work in extending the SPH method for simulating explosion phenomena including high explosive detonation, explosion [248, 351], shaped charge [247], and underwater explosions [143, 255]. Since the detonation and expansion speed
are extremely high, the gaseous products can be assumed to
be inviscid and the explosion process is adiabatic. The Euler
equations can be used to model the explosion process together with a suitable equation of state. The SPH equations
of motion derived in Sect. 2.5 can thus be used to model high
explosive explosion phenomena, with the viscosity terms in
the momentum and energy equations ignored. For the explosive gas, the standard Jones-Wilkins-Lee (JWL) [352] equation of state has been employed. These works not only show
the feasibility of applying the SPH method to explosion phenomena, but also identify some important and unresolved
issues. A Gruneisen form of equation of state for water has
been used, which is a polynomial form either in compressed
or expanded state [353].
To show the effectiveness of the SPH method in dealing with detonation and explosion, a benchmark onedimensional TNT slab denotation and explosion is provided
here, in which a 0.1 m long TNT slab detonates at one end
of the TNT slab.
If the solid wall boundary condition is used to prevent
material transport from everywhere, a symmetric setup can
be employed to deploy the particles, and thus makes the detonation of the 0.1 m long slab from one end to the other
59
Fig. 32 Pressure profiles along the TNT slab during the detonation
process (from [248])
end equivalent to the detonation of a 0.2 m long slab from

the middle point to both ends. Before detonation, particles
are evenly distributed along the slab. The initial smoothing
length is one and a half times the particle separation. After
ignition, a plane detonation wave is produced. According
to the detonation velocity, it takes around 14.4 s to complete the detonation to the end of the slab. Figure 32 show
the pressure profile along the slab at 2 s interval from 2 to
14 s by using 4000 particles. The dashed line in Fig. 32 represents the experimentally determined C-J detonation pressure, which is, according to the Chapman and Jouguets hypothesis, the pressure at the tangential point of the Hugoniot curve and the Rayleigh line, and represents the pressure at the equilibrium plane at the trailing edge of the very
thin chemical reaction zone. For this one-dimensional TNT
slab detonation problem, the experimental C-J pressure is
2.1 1010 N/m2 . It can be seen from Fig. 32 that, with the
process of the detonation, the detonation pressure converges
to the C-J pressure.
If the free boundary applies to the end of the detonation, the explosive gas behind the C-J plane disperses outwards with the forward propagating detonation wave. A theoretical solution for this one-dimensional TNT detonationdispersion problem exists. Figure 33 shows the comparisons
of velocity profiles between theoretical values and the presented SPH results at 1 and 2 s. The presented SPH results
are in close agreement with the theoretical values.
Figure 34 shows the pressure evolution in an underwater explosion process in a confined square chamber [143].
The initial outward propagating shock wave, reflection wave
from the solid wall, explosive gas expansion and later compression can all be seen from the figure. Right after the detonation, a shock wave is generated in the water and propagates outwards. At the same time, the rarefaction wave is
60
Fig. 33 Velocity transients at 1 and 2 s for the TNT slab detonation-dispersion process (from [248])
also produced within the explosive gas and advances inward.

With the advancement of the shock wave through the water, the gas bubble expands within the surrounding water.
The shock wave reaches the solid wall around the instant of
200 s, and then reflects from the solid wall. This reflection
wave propagates inward and tends to compress the expanding gas. At around 400 s, the gas bubble reaches to the
maximum and then reduces. With the continued process of
shock reflection and explosive gas contraction, the gas bubble then will reach a minimum size, and then expands afterwards. This repeated process of expansion and contraction
maintains for many circles, and will finally reach equilibrium at larger time step.
6.3 Microfluidics and Micro Drop Dynamics
By integrating mechanical elements, sensors, actuators, and
electronic components using micro-fabrication technology,
micro-electro-mechanical systems (MEMS) are fast in response, capable of achieving high spatial resolution, and
cost-effective due to the batch micromachining techniques.
Characterization of fluid flows in microfluidic devices has
increasingly becoming a very important topic since the fluidic behavior in MEMS is very different from what observed in daily life. Flows in microfluidic devices usually
involve small or ignorable inertial force, but dominant viscous, electro-kinetic and surface effects especially when the
surface-to-volume ratio increases [354]. Analytical or semianalytical solutions for microfluidics are generally limited
to a very few simple cases, whereas experimental studies
are usually expensive. Numerical simulation of flows in microfluidic devices, as an effective alternate, has been attracting more and more researchers. However, simulation of microfluidic devices is not easy due to the involved complex
features including movable boundaries (free surfaces and
M.B. Liu, G.R. Liu
moving interfaces), large surface-to-volume ratio, and phenomena due to micro scale physics. Numerical studies with
reliable models are needed to develop a better understanding
of the temporal and spatial dynamics of multiphase flows in
microfluidic devices.
On the other hand, drop formation and break-up in micro/nano scales are fundamentally important to diverse practical engineering applications such as ink-jet printing, DNA
and protein micro-/nano-arraying, and fabrication of particles and capsules for controlled release of medicines. Numerical studies provide an effective tool to improve better
understanding of the inherent physical dynamics of drop formation and breakup. Computational models for drop formation and breakup in micro/nano scales must be able to handle
movable boundaries such as free surfaces and moving interfaces, large density ratios, and large viscosity ratios. These
requirements together with micro scale phenomena and possible complex boundaries (fluid-fluid-solid contact line dynamics and fluid-fluid interface dynamics) in microfluidic
devices present severe challenges to conventional Euleriangrid based numerical methods such as FDM and FVM which
require special algorithms to treat and track the interfaces.
Algorithms based on Lagrangian-grid based methods such
as FEM have been shown to agree quantitatively with experimental measurements, but they are only capable of modeling the dynamics of formation of a single drop or the dynamics until the occurrence of the first singularity.
A number of meso-scale methods have been developed
for simulating micro- and nano-fluidics such as the direct
simulation Monte Carlo (DSMC) for rarefied gas flows
[355359], and dissipative particle dynamics for complex
fluid flows [41, 42, 118, 360]. Hoover and his co-workers
have first noted the similarity between the particle methods of molecular dynamics and smoothed particle hydrodynamics, and described the inherent links between them
[361363]. Espanol and his colleagues, when studying the
dissipative particle dynamics method for meso-scale applications, proposed a fluid particle dynamics model, which
is actually a synthesis of dissipative particle dynamics and
smoothed particle dynamics [96, 364]. Later on, they invented a smoothed dissipative particle dynamics for microor meso-scale applications, by introducing a fluctuation
term into the conventional smoothed particle hydrodynamics [365].
There have been a lot of literatures addressing the applications of SPH method to simulating microfluidics and
micro drop dynamics. Nugent and Posch described an approach to modeling liquid drops and surface tension for
a van der Waals fluid [236]. The cohesive pressure in the
equation of state for the van der Waals fluid actually acts
an attractive force between SPH particles. Melean and his
co-workers investigated the formation of micro drops using
61
Fig. 34 Pressure distributions

in a confined underwater
explosion (from [143])
SPH method [155]. It was reported that the tensile instability also exists when using SPH for simulating viscous liquid drops and using an artificial stress proposed by Monaghan et al. [145] can help greatly to remove the tensile
instability. Later on the group extended their work to the
simulation of coalescence of colliding van der Waals liquid
drops [240] and oscillation of viscous drops [239]. Liu et al.

applied the SPH method to multiphase fluid flow in micro
channels with applications to flip-chip underfill encapsulation process with both isotropic and anisotropic smoothing
kernels [123, 366]. Li et al. proposed an SPH model for
simulating droplet collision and coalescence [367]. Zhang
62
M.B. Liu, G.R. Liu
et al. developed an SPH model for free surface and solidification problems, which is also applicable to microfluidics
and micro drop dynamics such as droplet spreading, splashing and solidification, substrate melting and deformation
[243245]. A revised surface tension model was developed
for macro-scale particle methods including SPH by Zhou
et al. [246]. Fang et al. also developed an improved SPH
model for the simulations of droplet spreading and solidification [238]. Tartakovsky, Meakin and their co-workers
have conducted a series of excellent work in applying the
SPH method to the modeling of surface tension and contact angle [235], miscible flow in three-dimensional fractures and the two-dimensional Rayleigh Taylor instability
[368], unsaturated flow in complex fractures [117], reactive
transport and precipitation [369], mixing-induced precipitation [370] and non-aqueous phase liquid flow and dissolution [371].
There are basically two approaches in modeling the multiphase fluid flow in micro fluidic devices, and multiphase
drop dynamics. The first approach is to use the continuum surface force (CSF) model proposed by Brackbill et al.
[372], and introduced the surface tension force into the momentum equation as presented in Sect. 5. This approach is
straightforward, and real physical parameters are used in the
simulation. It is important to note that the surface tension
parameters in this approach are user-input parameters, and
surface curvature needs to be calculated to get the surface
tension force. Some researchers have used this approach in
SPH method to model multiphase fluid flow [123, 373].
With this CSF model, the surface tension is added to the
momentum equation as an external source force. The surface
tension force F can be formulated as follows
F = 2 k(x)V V ,
(109)
where is the surface tension coefficient, k is the surface

curvature. V is the volume of a fluid cell. Let n be the surface
normal which can be computed from the gradient of the fluid
cell as n = V , the surface curvature k is defined as

n
1
|n| ( n) ,
(110)
k = n =
|n| |n|
where n is the unit surface normal defined as
n =
n
.
|n|
(111)
For fluid cells near a solid wall, a wall adhesion model is

applied to adjust the surface normal through using the contact angle that a fluid is assumed to make with the wall.
The dynamic contact angle can be different for different contacting points. However, in micro channel flows with a small
buck velocity, the dynamic contact angle can be assumed to
be constant along the flow direction as the starting contact
Fig. 35 Snapshots of the moving flow leading edge in an SPH simulation of a multiphase fluid flow in a micro channel using the CSF model
(from [123])
angle. The unit normal n is adjusted as follows for fluid cells

close to the wall
n = n w cos + n t sin
(112)
where n w and n t are the unit vectors of the surface normal

(n = V ) normal and tangential to the wall, respectively.
Thereby the unit surface normal n for fluid cells near the
wall is adjusted using the contact angle. Whereas the unit
surface normal n for fluid cells one cell away from the wall
is normally calculated. After determination of the unit surface normal, the local surface curvature and therefore the
surface tension can be calculated.
Figure 35 gives an example of CSF model in simulating micro channel fluid flow. It is an underfill encapsulation
process. The encapsulant is used to fill the space between the
solder joints under the chip and the encapsulation process is
either driven by a capillary action or by a pressurized injection [123].
Another approach is to introduce an inter-particle interaction force (IIF) implicitly in the SPH equations, and can be
termed as IIF model. This attractive force between every pair
of SPH particles contribute to the surface tension. Considering the SPH equations presented in Sect. 2.5, an equation of
state can be used to close the equation system. An equation
of state describes the relationship of the pressure p, density
and the internal energy per unit mass, e, respectively. For
example, the van der Waals (vdW) equation of state can be
used to model the behavior of the fluid under consideration.
The van der Waals equation of state was derived from statistical mechanics as the mean-field limit for the free energy
density of a system of hard particles with a superimposed
long-range, attractive pair potential. It is realistic to display
a gas-to-liquid phase transition similar to that of a real fluid.
The van der Waals equation of state can written as
p=
kT
a 2 ,
1 b
(113)
63
and
a.
e = kT
(114)
In the above two equations, k = kB /T , where kB is the

Boltzmanns constant. T is the system temperature. a =
a/m2 , and b = b/m, where a and b are the parameter describing a van der Waals fluid. a controls the strength of the
attractive force, and b relates to the size of the particle.
The second part in (113) describes the cohesion between
particles. Using SPH particle approximation for pressure
(see (28)), and considering this cohesive pressure part separately, we can get

Dv i
Wij
mj
.
= 2a
Dt
x i
Fig. 36 SPH simulation of the formation of a liquid drop using the IIF
model. (a) initial stage, (b) final stage of liquid drop
(115)
j =1
Similarly substituting the cohesive pressure into the energy (see (28)), it is possible to get its analogous contribution
to energy

Dei
Wij
= a
mj v ij
.
Dt
x
N
j =1
(116)
It is clear that the attractive, long-range, inter-atomic

vdW force can be transformed into similarly attractive
forces between SPH particles. Equation (115) describes a
volume force on continuum scale to account for the formation of co-existing liquid-gas phases. Comparing with the
above-mentioned approach in using CSF model in calculating surface force, this approach does not need to locate
the surface, and then to calculate the local surface curvature. Also, the surface tension is not user-input parameters,
while it is implicitly obtained from inter-particle interactions. If the particles are distributed regularly (please refer
to Sect. 4.2 on particle consistency), the force obtained from
(115) vanishes for interior particles in both liquid and gas
phases. While for boundary particles, e.g., particles near
the gas-liquid interface, since the inter-particle interaction
force between different phases are generally different, a surface force is produced, which is basically perpendicular to
the surface, pointing towards the dense phase. Nugent and
Posch identified that to obtain acceptable results, a large influencing area (e.g., two times the smoothing length for approximating other field variables and the first part in (113))
is necessary to conduct the particle approximations in (115)
and (116).
Figure 36 shows an example of IIF model for the formation of a liquid drop. The initial square-shaped liquid
(Fig. 36a) gradually becomes rounded at the corner due to
the inter-particle interaction force, and finally forms a circular liquid drop.
Except for the inter-particle interaction force implicitly
calculated from the cohesive pressure in the van der Waals
Fig. 37 The shape (cos(1.5/kh)) of the inter-particle interaction

force expressed in (117) with k = 2
equation of state, Tartakovsky et al. added another interparticle interaction force, which can be written as
F ij = s ij cos(1.5/kh)r ij ,
(117)
where s ij is an interaction coefficient. Figure 37 shows the

shape of the inter-particle interaction force expressed in the
above equation (cos(1.5/kh) only, where k = 2). It is clear
that this inter-particle force is repulsive at short range, and
attractive at long-distance. Similar to inter-particle force calculated from the cohesive pressure in the van der Waals
equation of state, this inter-particle also describes a volume
force on continuum scale and can account for the formation
of co-existing liquid-gas phases. The force obtained from
(117) also vanishes for interior particles in both liquid and
gas phases, with small fluctuation around the overall direction of the macroscopic density gradient. For boundary particles, e.g., particles near the gas-liquid interface, a surface
force is produced, pointing towards the denser phase. The
IIF approach is comparatively simple and straightforward
since it does not need to calculate the surface curvature,
64
which is not an easy task for particle methods like SPH. One
problem is that as the IIF model implicitly calculates the surface tension force with parameters from atomistic level, it
needs parameter calibration, which relates the physical parameters from atomistic level to continuum level.
6.4 Ocean and Coastal Hydrodynamics and Offshore
Engineering
Flow phenomena in ocean and coastal hydrodynamics and
offshore engineering are significantly important as they can
greatly influence the nearly personnel and structures. The
flow phenomena include
wave dynamics (waver generation, wave breaking, and
wave interaction with other structures),
dam breaking,
water filling and water discharge (to and from a water tank
or reservoir),
shallow water flows,
entry of water, sloshing phenomena with fluid-solid interaction, and
different other problems.
These phenomena involve special features, which make it
difficult for numerical simulation. For example, water waves
can propagate shoreward where they undergo changes induced by the near-shore topography and increase in height.
Upon reaching the shoreline, they can break into pieces, and
travel inland for large distances with potential damage of
property and loss of life. Experimental setups for fluid flow
in coast hydrodynamics and offshore engineering are expensive and only limited to laboratory applications. Numerical
simulation has thus become a great tool to predicting fluid
flow in ocean and coast hydrodynamics and offshore engineering.
However, numerical simulation of fluid flow in these related areas is a formidable task as it involves not only complex geometries and free surfaces, but also fluid-solid interaction as well as other complex physics in a comparably very large scale. In many circumstances, violent fluid
structure interactions lead to air entrapment and multi-phase
flows, where the dynamics of the entrapped air at the impact may play a dominant role during the process and contribute to the high pressure maxima and pressure oscillations. Though conventional grid based methods like FDM,
FVM and FEM have achieved greatly in simulating fluid
flow in coast hydrodynamics and offshore engineering, there
is still a long way to go for practical engineering applications.
Smoothed particle hydrodynamics, due to its meshfree,
Lagrangian and particle nature, has been attracting more and
more researchers in coast hydrodynamics and offshore engineering. From the very early simulation of a simple dam
M.B. Liu, G.R. Liu
break problem [112], there have been a lot of literatures addressing the applications of SPH method in related areas.
Shao and his colleagues simulated near-shore solitary
wave mechanics by an incompressible SPH method [374].
They later extended the work to simulation of wave breaking and overtopping with turbulence modeling [375], plunging waves using a 2-D sub-particle scale (SPS) turbulence
model [376], solitary wave interaction with a curtain-type
breakwater [377], and water entry of a free-falling object
[190].
Frank and Reich addressed the conservation properties
of SPH applied to shallow water equation [378]. Ata and
Soulaimani proposed a stabilized SPH method for inviscid
shallow water flows [379]. Rodriguez-Paz and Bonet also
developed a corrected smooth particle hydrodynamics formulation of the shallow-water equations [380].
Investigation of water wave including water wave generation, water-structure interaction, and water wave breaking
has been a very attractive application of the SPH method.
Gomez-Gesteira and Dalrymple had investigated wave impact on a tall structure using a three-dimensional SPH
method [180]. The research group also investigated water
waves and waves breaking using the SPH method [175,
189]. Idelsohn and his co-workers applied the particle finite element method to solve incompressible flows with
free-surfaces and breaking waves [381]. Gotoh et al. developed an SPH-LES model for numerical investigation of
wave interaction with partially immersed breakwater [178].
They also simulated the coupled motion of progressive
wave and floating curtain wall by using the developed SPHLES model [177]. Gotoh and Sakai also addressed some
key issues in the particle method for computation of wave
breaking [183]. Recently, the research team developed a
corrected incompressible SPH method for accurate watersurface tracking in breaking waves [185]. Crespo et al. presented an example of 3D SPH simulation of large waves mitigation with a dike [176]. Qiu studied the water waves generated by landslide [188]. Yim investigated the water wave
generation by a vertical plunger using RANS and SPH models [193]. Cleary and Prakash discussed the feasibility of using SPH method for modeling tsunami effects [205, 266].
There are some references addressing the entry of water,
e.g., air cushion effects in a wedge water entry [382], waves
produced by a falling mass into a reservoir [383], water entry
of a free-falling object [190], and wedge water entries [384].
Simulation of sloshing problems using SPH is a promising research direction. Iglesias et al. simulated the anti-roll
tanks and sloshing type problems [196]. Rhee and Engineer
studied liquid tank sloshing with Reynolds-averaged NavierStokes [199]. Souto-Iglesias et al. assessed the liquid moment amplitude in sloshing type problems with smooth particle hydrodynamics [200]. Anghileri investigated the fluidstructure interaction of water filled tanks during the impact
Fig. 38 SPH simulation of the dam-break flow and impact against a

vertical wall at (a) single phase (water only) with free surface, and
(b) two phases (water and air) with moving interface (from [387])
with the ground [285]. Delorme et al. simulated the sloshing

loads in LNG tankers with SPH [385].
Figure 38 gives an example of the SPH simulation of the
dam-break flow and impact again a vertical wall at (a) single phase (water only) with free surface, and (b) two phases
(water and air) with moving interface. The flow pattern of
flow along the surge front, the impact of fluid against the
right hand side vertical wall and the entrapment of air bubble in water are in good agreement with the results from
other sources [386].
6.5 Environmental and Geophysical Flows
Environmental and geophysical flows are fundamentally important for human life and economic development. In large
scale, there are many typical environmental flow problems
including
flood due to dam break, or overflow of river and reservoir,
landslide, mudslide, mud-rock flow triggered by heavy
rain infiltration or even earth quake,
sand/dust or dust storm,
polluted liquid (e.g. oil) and air transport, and
many others.
Unwanted environmental flows can be hazardous to nearby
people and structures. In general, these large scale environmental flows are difficult to be accurately simulated and precisely predicted. On one hand, these large scale environmental flows involve complex geometries and arbitrarily moving interfaces. It is a formidable task to identify the origin
of the flows. Therefore, it is hard to determine the initial
and boundary condition for numerical simulations. On the
other hand, large scale environmental flows are associated
with complicated multiple fluid phases and multiple physics
with solid-liquid-air phase at different scales. Theses also
add difficulties to the numerical simulations for grid based
numerical models.
It is very appealing to use SPH to model such largescale environmental flow problems due to its meshfree, Lagrangian and particle nature. Actually in Sect. 6.4, some
large-scale environmental flow problems such as dam breaking and tsunami effects have been discussed. They are also
other references addressing the SPH applications in largescale environmental flow problems. Ghazali and Kamsin
65
presented a real time simulation and modeling of flood hazard [206]. Shen et al. conducted SPH simulation of river ice
dynamics [233]. Kipfer and Westermann also investigated
the realistic and interactive flows in river hydrodynamics
[207]. Xu and Shen studied the fluid-structure interaction of
hydrodynamic damper during the rush into the water channel [388].
Landslide and mudslide as well as mud-rock flows can be
simulated using the SPH method. Cleary and his co-workers
discussed the feasibility of using SPH method for modeling dam break, tsunami, landslide and volcano flows [205,
266]. The SPH method is coupled with discrete element
method (DEM) for modeling solid-fluid interaction. AtaieAshtiani and Shobeyri simulated landslide impulsive waves
by using incompressible smoothed particle hydrodynamics [389]. Qiu presented a two-dimensional SPH simulations of landslide-generated water waves [188]. McDougall
and Hungr developed a numerical model for the analysis
of rapid landslide motion across three-dimensional terrain
[216]. Pastor et al. proposed a depth-integrated, coupled
SPH model for flow-like landslides and related phenomena
[220].
There are more references with applications to other
large-scale environmental flow problems. Bui et al. developed an SPH model for large deformation and failure flows
of geo-material using elastic-plastic soil constitutive model
[209]. They also presented a numerical simulation of soilwater interaction using the SPH method [210]. Laigle and
his colleagues developed an SPH-based numerical investigation of mudflow and other complex fluid flow interactions
with structures [214].
Small scale environmental and geophysical flows are also
very important, but are usually difficult to simulate because
of the associated multiple fluid phases and multiple physics,
as well as the existence of complex geometries and arbitrarily moving interfaces. For example, fluid motion in the vadose zone is very important for groundwater recharge, fluid
motion and contaminant transport. Flow through fractures
and fractured porous media can lead to exceptionally rapid
movement of liquids and associated contaminants [99, 390,
391]. The physics of fluid flows in unsaturated fractures and
porous media is still poorly understood due to the complexity of multiple phase flow dynamics. Experimental studies
of fluid flow in fractures and fractured porous media are
limited, and in computer simulations it is usually difficult
to take into account the fracture surface properties and microscopic roughness. Predictive numerical models can be
divided into two general classes: volume-averaged continuum models (such as those based on Richards equation)
and discrete mechanistic models. Knowledge of the physical properties of the fluids and the geometry of the fracture
apertures is required in both classes. Volume-averaged continuum models are more suitable for large-scale systems,
66
M.B. Liu, G.R. Liu
Fig. 39 Snapshots of SPH

simulation of multiphase fluid
motion in a fracture at 4
representative stages
and they usually involve the representation of fractures as

porous media with porosity and permeability parameters adjusted to mimic flow within fractures. However, volumeaveraged continuum models are unable to describe the details of flow dynamics in fractures, they do not reproduce
the spatio-temporal complexity of multiphase fluid flow in
fractures, and they often fail to predict the rapid fluid motion
and contaminant transport observed in the fractured vadose
zone. Small-scale studies with discrete mechanistic models
are needed to develop a better understanding of the temporal
and spatial dynamics of fracture flows. However, the complexity of fracture flow dynamics makes it difficult to develop successful numerical models for fluid flows in fracture
networks. A broadly applicable model must be able to simulate a variety of phenomena including film flow with free
surfaces, stable rivulets, snapping rivulets, fluid fragmentation and coalescence (including coalescence/fragmentation
cascades), droplet migration and the formation of isolated
single-phase islands trapped due to aperture variability.
Realistic mechanistic models for multiphase fluid flows
in fracture and fractured porous media must be able to handle moving interfaces, large density ratios (e.g., 1000 : 1
for water and air), and large viscosity ratios (e.g., 100 : 1
for water and air). These requirements combined with the
complex geometries of natural fractures present severe challenges to mechanistic models. Grid based numerical methods such as finite difference methods, finite volume methods and Eulerian finite element methods require special algorithms to treat and track the interface between different
phases. However, continuum grid based numerical models
usually do not take account of the detailed void and obstacle
geometries, fluid-fluid interface dynamics within pores and
complex fluid-fluid-solid contact line dynamics. They rely
on constitutive equations that describe the coarse-grained
behavior and can, at least in principle, be derived from the
results of pore scale simulations or experiments. Therefore,
small-scale simulations with mechanistic models are needed
to develop a better understanding of the temporal and spatial
dynamics of multiphase flow through porous media. Porescale flows have been studied extensively using grid based
methods including finite difference method [392], finite volume method [393, 394], and finite element method [395].
However, due to the difficulties associated with geometrically complex boundaries, fluid-fluid-solid contact line dynamics, and fluid-fluid interface dynamics, it is difficult to
apply conventional grid based multiphase simulation methods based on computational fluid dynamics (CFD) coupled
with interface tracking algorithms [372, 396398] to porescale multiphase flow modeling.
The SPH method has been recently modified for much
smaller scale, typically low Reynolds number, applications
The performance of the SPH method was demonstrated for
two-dimensional single-phase flows through idealized, porescale porous media composed of spatially periodic square
and hexagonal arrays of cylinders [202, 399], and two-phase
(miscible and immiscible) flows through pore-scale fractures and fracture junctions [117, 368]. One obvious advantage of SPH over conventional grid based methods is
that SPH does not require explicit and complicated interface tracking algorithms, and thus there is no need to explicitly track the material interfaces, and processes such as fluid
fragmentation and coalescence can be handled without difficulty. SPH also does not require contact angle models since
contact angles can be naturally calculated from the shape of
the moving particle distributions, and can vary spatially and
temporally, depending on the dynamic balance of viscous,
capillary and gravitational forces. Figure 39 illustrates the
snapshots of SPH simulation of multiphase fluid motion in
a fracture at 4 representative stages. Water was injected into
the top entrance of the fracture using a syringe pump and
drained out through the bottom of the fracture.
7 Summary
Smoothed particle hydrodynamics (SPH) is a meshfree particle method, which not only uses particles as the computa-
tional frame for interpolation or differencing, but also uses

the particles to carry the material properties. In this survey, the smoothed particle hydrodynamics method, and its
recent development in numerical algorithms and applications have been reviewed. In methodology and numerical
approximation, the basic concepts of kernel and particle approximations as well as different techniques for developing
SPH formulations have been addressed. Different smoothing kernel functions have been reviewed with constructing
conditions provided. The emphasis is on the consistency
problem, which traditionally limits the accuracy of the SPH
method. A number of consistency restoring approaches have
been surveyed. An SPH formulation for discontinuity, a general approach for restoring particle consistency, and a finite
particle method are described. Several important numerical
topics have been investigated, including (1) solid boundary
treatment, (2) representation of solid obstacles, (3) material interface treatment and (4) tensile instability. In applications, different applications of the SPH mEthod have been
reviewed, with emphasis on (1) high strain hydrodynamics
with material strength, (2) high explosive detonation and explosions, and underwater explosion, (3) microfluidics and
micro drop dynamics, (4) coast hydrodynamics and offshore
engineering, and (5) environmental and geophysical flows.
The last decades have witnessed the great success of SPH
developments in methodology and applications, and there
are still many tasks and challenges remaining. To achieve
a reliable solution, the computational accuracy, consistency,
efficiency, stability and convergence need to be incorporated
into good SPH algorithms.
In general, SPH has great potential in many problems in
engineering and science. It has salient advantages over traditional grid based numerical models in treating large deformation, tracking free surfaces, moving interfaces, and deformable boundaries, resolving moving discontinuities such
as cracks and shock waves. The attraction of the SPH
method has been showcased in diversified applications, as
reviewed in Sect. 7.
One obvious advantage is that the SPH method provides
a feasible physical model for non-continuum, as it has some
same features as the classic molecular dynamics method
and the dissipative particle dynamics method. Therefore it
would be very attractive to apply the SPH method to simulating problems where the main concern of the object is a set
of discrete physical particles rather than a continuum, e.g.,
the interaction of stars in astrophysics, movement of millions of atoms in an equilibrium or non-equilibrium state,
dynamic behavior of protein molecules, and environmental flows with solid and fluid particles such as landslide and
mudslide. There is no need for discretization to begin with
for such situations. Also since SPH shares many similarities with other particle methods such as MD and DPD, it is
natural to couple these three methods, from molecular dynamics (at nano- and micro-scales), to dissipative particle
67
dynamics (at meso-scales), and then to smoothed particle

hydrodynamics (at macro-scales) for simulations with multiple scale physics. It is attractive to develop reliable models
of length scale coupling for problems with multiple physics
and multiple scales, such as ink-jet printing, DNA and protein micro-/nano-arraying, and fabrication of particles and
capsules for controlled release of medicines.
Acknowledgements The first author, M.B. Liu, has been supported
by the National Natural Science Foundation of China (NSFC Grant No.
10942004).
References
1. Chung TJ (2002) Computational fluid dynamics. Cambridge
University Press, Cambridge
2. Anderson JD (2002) Computational fluid dynamics: the basics
with applications. McGraw Hill, New York
3. Zienkiewicz OC, Taylor RL (2000) The finite element method.
Butterworth-Heinemann, Stonham
4. Liu GR (2002) Meshfree methods: moving beyond the finite element method. CRC Press, Boca Raton
5. Hirsch C (1988) Numerical computation of internal & external flows: fundamentals of numerical discretization. Wiley, New
York
6. Liu GR, Liu MB (2003) Smoothed particle hydrodynamics: a
meshfree particle method. World Scientific, Singapore
7. Liu GR, Gu YT (2005) An Introduction to meshfree methods and
their programming. Springer, Dordrecht, p 479
8. Zhang SZ (1976) Detonation and its applications. Press of National Defense Industry, Beijing
9. Zukas JA (1990) High velocity impact dynamics. Wiley, New
York
10. Hockney RW, Eastwood JW (1988) Computer simulation using
particles. Institute of Physics Publishing, Bristol
11. Allen MP, Tildesley DJ (1987) Computer simulation of liquids.
Oxford University Press, Oxford
12. Li S, Liu WK (2002) Meshfree and particle methods and their
applications. Appl Mech Rev 55(1):134
13. Belytschko T, Krongauz Y, Organ D, Fleming M, Krysl P (1996)
Meshless methods: an overview and recent developments. Comput Methods Appl Mech Eng 139(14):347
14. Idelsohn SR, Onate E (2006) To mesh or not to mesh? That is the
question. Comput Methods Appl Mech Eng 195(3740):4681
4696
15. Nguyen VP, Rabczuk T, Bordas S, Duflot M (2008) Meshless
methods: a review and computer implementation aspects. Math
Comput Simul 79(3):763813
16. Liu GR (2008) A generalized Gradient smoothing technique
and the smoothed bilinear form for Galerkin formulation of a
wide class of computational methods. Int J Comput Methods
5(2):199236
17. Liu GR (2009) A G space theory and weakened weak (W2) form
for a unified formulation of compatible and incompatible methods, part I: theory. Int J Numer Methods Eng. doi:10.1002/nme.
2719
18. Liu GR (2009) A G space theory and weakened weak (W2) form
for a unified formulation of compatible and incompatible methods, part II: applications to solid mechanics problems. Int J Numer Methods Eng. doi:10.1002/nme.2720
19. Liu GR (2009) On the G space theory. Int J Comput Methods
6(2):257289
68
20. Liu GR, Nguyen-Thoi T, H. N-X, Lam KY (2008) A node-based
smoothed finite element method (NS-FEM) for upper bound solutions to solid mechanics problems. Comput Struct 87:1426
21. Liu GR, Zhang GY (2008) Upper bound solution to elasticity problems: a unique property of the linearly conforming
point interpolation method (LC-PIM). Int J Numer Methods Eng
74:11281161
22. Liu GR, Nguyen-Thoi T, Lam KY (2009) An edge-based
smoothed finite element method (ES-FEM) for static, free and
forced vibration analyses in solids. J Sound Vib 320:11001130
23. Zhang GY, Liu GR, Nguyen TT, Song CX, Han X, Zhong ZH,
Li GY (2007) The upper bound property for solid mechanics of
the linearly conforming radial point interpolation method (LCRPIM). Int J Comput Methods 4(3):521541
24. Liu GR, Zhang GY (2009) A normed G space and weakened
weak (W2) formulation of a cell-based smoothed point interpolation method. Int J Comput Methods 6(1):147179
25. Liu GR, Zhang GY (2008) Edge-based Smoothed Point Interpolation Methods. Int J Comput Methods 5(4):621646
26. Liu GY, Zhang GY (2009) A novel scheme of strain-constructed
point interpolation method for static and dynamic mechanics
problems. Int J Appl Mech 1(1):233258
27. Liu GR, Xu X, Zhang GY, Nguyen-Thoi T (2009) A superconvergent point interpolation method (SC-PIM) with piecewise linear strain field using triangular mesh. Int J Numer Methods Eng
77:14391467
28. Liu GR, Xu X, Zhang GY, Gu YT (2009) An extended Galerkin
weak form and a point interpolation method with continuous
strain field and superconvergence using triangular mesh. Comput Mech 43:651673
29. Nguyen-Thoi T, Liu GR, Lam KY, Zhang GY (2009) A facebased smoothed finite element method (FS-FEM) for 3D linear
and nonlinear solid mechanics problems using 4-node tetrahedral
elements. Int J Numer Methods Eng 78:324353
30. Liu GR, Nguyen-Thoi T, Lam KY (2008) A novel Alpha finite
element method (FEM) for exact solution to mechanics problems
using triangular and tetrahedral elements. Comput Methods Appl
Mech Eng 197:38833897
31. Liu GR, Nguyen-Xuan H, Nguyen TT, Xu X (2009) A novel
weak form and a superconvergent alpha finite element method
for mechanics problems using triangular meshes. J Comput Phys
228(11):39114302
32. Liu GR, Zhang J, Lam KY (2008) A gradient smoothing method
(GSM) with directional correction for solid mechanics problems.
Comput Mech 41:457472
33. Liu GR, Xu XG (2008) A gradient smoothing method (GSM) for
fluid dynamics problems. Int J Numer Methods Fluids 58:1101
1133
34. Xu XG, Liu GR (2008) An adaptive gradient smoothing method
(GSM) for fluid dynamics problems. Int J Numer Methods Fluids
doi:10.1002/fld.2032
35. Xu XG, Liu GR, Lee KH (2009) Application of gradient smoothing method (GSM) for steady and unsteady incompressible flow
problems using irregular triangles. Int J Numer Methods Fluids
(submitted)
36. Lucy LB (1977) A numerical approach to the testing of the fission hypothesis. Astron J 82(12):10131024
37. Gingold RA, Monaghan JJ (1977) Smoothed particle
hydrodynamicstheory and application to non-spherical
stars. Mon Not R Astron Soc 181:375389
38. Liu MB, Liu GR, Lam KY (2003) Constructing smoothing functions in smoothed particle hydrodynamics with applications.
J Comput Appl Math 155(2):263284
39. Fulk DA, Quinn DW (1996) An analysis of 1-D smoothed particle hydrodynamics kernels. J Comput Phys 126(1):165180
40. Frenkel D, Smit B (2002) Understanding molecular simulation:
from algorithms to applications. Academic Press, New York
M.B. Liu, G.R. Liu

41. Hoogerbrugge PJ, Koelman J (1992) Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics. Europhys Lett 19:155
42. Groot RD (1997) Dissipative particle dynamics: Bridging the
gap between atomistic and mesoscopic simulation. J Chem Phys
107(11):4423
43. Gingold RA, Monaghan JJ (1982) Kernel estimates as a basis
for general particle method in hydrodynamics. J Comput Phys
46:429453
44. Hu XY, Adams NA (2006) Angular-momentum conservative
smoothed particle dynamics for incompressible viscous flows.
Phys Fluids 18:101702
45. Swegle JW, Hicks DL, Attaway SW (1995) Smoothed particle
hydrodynamics stability analysis. J Comput Phys 116(1):123
134
46. Morris JP (1996) Analysis of Smoothed Particle Hydrodynamics
with Applications. Monash University
47. Monaghan JJ (1982) Why particle methods work. SIAM J Sci
Stat Comput 3(4):422433
48. Monaghan JJ (1985) Particle methods for hydrodynamics. Comput Phys Rep 3:71124
49. Monaghan JJ (1992) Smooth particle hydrodynamics. Annu Rev
Astron Astrophys 30:543574
50. Johnson GR (1996) Artificial viscosity effects for SPH impact
computations. Int J Impact Eng 18(5):477488
51. Johnson GR, Beissel SR (1996) Normalized smoothing functions for SPH impact computations. Int J Numer Methods Eng
39(16):27252741
52. Randles PW, Libersky LD (1996) Smoothed particle hydrodynamics: some recent improvements and applications. Comput
Methods Appl Mech Eng 139(1):375408
53. Chen JK, Beraun JE, Carney TC (1999) A corrective smoothed
particle method for boundary value problems in heat conduction.
Int J Numer Methods Eng 46:231252
54. Chen JK, Beraun JE (2000) A generalized smoothed particle hydrodynamics method for nonlinear dynamic problems. Comput
Methods Appl Mech Eng 190:225239
55. Liu MB, Liu GR, Lam KY (2003) A one-dimensional meshfree
particle formulation for simulating shock waves. Shock Waves
13(3):201211
56. Liu MB, Liu GR (2006) Restoring particle consistency in
smoothed particle hydrodynamics. Appl Numer Math 56(1):19
36
57. Liu MB, Xie WP, Liu GR (2005) Modeling incompressible flows
using a finite particle method. Appl Math Model 29(12):1252
1270
58. Batra RC, Zhang GM (2004) Analysis of adiabatic shear bands
in elasto-thermo-viscoplastic materials by modified smoothedparticle hydrodynamics (MSPH) method. J Comput Phys
201(1):172190
59. Fang JN, Owens RG, Tacher L, Parriaux A (2006) A numerical study of the SPH method for simulating transient viscoelastic
free surface flows. J Non-Newton Fluid 139(12):6884
60. Fang JN, Parriaux A (2008) A regularized Lagrangian finite
point method for the simulation of incompressible viscous flows.
J Comput Phys 227(20):88948908
61. Fang JN, Parriaux A, Rentschler M, Ancey C (2009) Improved
SPH methods for simulating free surface flows of viscous fluids.
Appl Numer Math 59(2):251271
62. Dyka CT, Ingel RP (1995) An approach for tension instability in smoothed particle hydrodynamics (SPH). Comput Struct
57(4):573580
63. Dyka CT, Randles PW, Ingel RP (1997) Stress points for tension
instability in SPH. Int J Numer Methods Eng 40(13):23252341
64. Randles PW, Libersky LD, (2000) Normalized SPH with stress
points. Int J Numer Methods Eng 48(10):14451462

65. Vignjevic R, Campbell J, Libersky L (2000) A treatment of zeroenergy modes in the smoothed particle hydrodynamics method.
Comput Methods Appl Mech Eng 184(1):6785
66. Dilts GA (1999) Moving-Least-Squares-particle hydrodynamics.
I: Consistency and stability. Int J Numer Methods Eng 44:1115
1155
67. Dilts GA (2000) Moving least square particle hydrodynamics ii: conservation and boundaries. Int J Numer Methods Eng
48:15031524
68. Bonet J, Kulasegaram S (2000) Correction and stabilization
of smooth particle hydrodynamics methods with applications
in metal forming simulations. Int J Numer Methods Eng
47(6):11891214
69. Chen JS, Pan C, Wu CT, Liu WK (1996) Reproducing kernel particle methods for large deformation analysis of nonlinear structures. Comput Methods Appl Mech Eng 139:195227
70. Liu WK, Chen Y, Jun S, Chen JS, Belytschko T, Pan C, Uras RA,
Chang CT (1996) Overview and applications of the reproducing
kernel particle methods. Arch Comput Methods Eng 3(1):380
71. Rabczuk T, Belvtschko T, Xiao SP (2004) Stable particle methods based on Lagrangian kernels. Comput Methods Appl Mech
Eng 193(12):10351063
72. Belytschko T, Krongauz Y, Dolbow J, Gerlach C (1998) On
the completeness of the meshfree particle methods. Int J Numer
Methods Eng 43(5):785819
73. Zhang GM, Batra RC (2004) Modified smoothed particle hydrodynamics method and its application to transient problems. Comput Mech 34(2):137146
74. Belytschko T, Guo Y, Liu WK, Xiao SP (2000) A unified stability
analysis of meshless particle methods. Int J Numer Methods Eng
48:13591400
75. Fulk DA (1994) A numerical analysis of smoothed particle hydrodynamics. Air Force Institute of Technology
76. Monaghan JJ (2005) Smoothed particle hydrodynamics. Rep
Prog Phys 68(8):17031759
77. Hernquist L (1993) Some cautionary remarks about smoothed
particle hydrodynamics. Astrophys J 404(2):717722
78. Century Dynamics Incorporated (1997) AUTODYN release
notes version 3.1. AUTODYN interactive non-linear dynamic
analysis software
79. Benz W (1990) Smooth particle hydrodynamics: a review. The
numerical modelling of nonlinear stellar pulsations, problems
and prospects. Kluwer Academic, Boston
80. Morris JP (1996) A study of the stability properties of smooth
particle hydrodynamics. Publ Astron Soc Aust 13(1):97102
81. Omang M, Borve S, Trulsen J (2005) Alternative kernel functions
for smoothed particle hydrodynamics in cylindrical symmetry.
Shock Waves 14(4):293298
82. Jin HB, Ding X (2005) On criterions for smoothed particle hydrodynamics kernels in stable field. J Comput Phys 202(2):699
709
83. Capuzzo-Dolcetta R, Di Lisio R (2000) A criterion for the choice
of the interpolation kernel in smoothed particle hydrodynamics.
Appl Numer Math 34(4):363371
84. Cabezon RM, Garcia-Senz D, Relano A (2008) A one-parameter
family of interpolating kernels for smoothed particle hydrodynamics studies. J Comput Phys 227(19):85238540
85. Monaghan JJ, Lattanzio JC (1985) A refined particle method for
astrophysical problems. Astron Astrophys 149(1):135143
86. Johnson GR, Stryk RA, Beissel SR (1996) SPH for high velocity
impact computations. Comput Methods Appl Mech Eng 139(1
4):347373
87. Swegle JW, Attaway SW, Heinstein MW, Mello FJ Hicks DL
(1994) An analysis of smoothed particle hydrodynamics. Sandia
National Labs., Albuquerque
88. Liu GR, Quek SS (2003) The finite element method: a practical
course. Butterworth-Heinemann, Bristol
69
89. Chen JK, Beraun JE, Jih CJ (1999) An improvement for tensile
instability in smoothed particle hydrodynamics. Comput Mech
23(4):279287
90. Chen JK, Beraun JE, Jih CJ (1999) Completeness of corrective smoothed particle method for linear elastodynamics. Comput Mech 24(4):273285
91. Hernquist L, Katz N (1989) TREESPHa unification of SPH
with the hierarchical tree method. Astrophys J Suppl Ser
70(2):419446
92. Sod GA (1978) A survey of several finite difference methods
for systems of nonlinear hyperbolic conservation laws. J Comput Phys 27:131
93. Fletcher CAJ (1991) Computational techniques for fluid dynamics, 1: fundamental and general techniques. Springer, Berlin
94. Agertz O, Moore B, Stadel J, Potter D, Miniati F, Read J, Mayer
L, Gawryszczak A, Kravtosov A, Nordlund A, Pearce F, Quilis
V, Rudd D, Springel V, Stone J, Tasker E, Teyssier R, Wadsley J,
Walder R (2007) Fundamental differences between SPH and grid
methods. Mon Not R Astron Soc 380:963978
95. Schussler M, Schmitt D (1981) Comments on smoothed particle
hydrodynamics. Astron Astrophys 97(2):373379
96. Espanol P (1998) Fluid particle model. Phys Rev E 57(3):2930
2948
97. Wang L, Ge W, Li J (2006) A new wall boundary condition in
particle methods. Comput Phys Commun 174(5):386390
98. Revenga M, Zuniga I, Espanol P (1998) Boundary models in
DPD. Int J Mod Phys C 9(8):13191328
99. Liu MB, Meakin P, Huang H (2007) Dissipative particle dynamics simulation of pore-scale flow. Water Resour Res 43:W04411.
doi:10.1029/2006WR004856
100. Revenga M, Zuniga I, Espanol P (1998) Boundary models in
DPD. Int J Mod Phys C 9(8):13191328
101. Revenga M, Zuniga I, Espanol P (1999) Boundary conditions in dissipative particle dynamics. Comput Phys Commun
121(122):309311
102. Willemsen SM, Hoefsloot HCJ, Iedema PD (2000) No-slip
boundary condition in dissipative particle dynamics. Int J Mod
Phys C 11(5):881890
103. Duong-Hong D, Phan-Thien N, Fan X (2004) An implementation of no-slip boundary conditions in DPD. Comput Mech
35(1):2429
104. Crespo AJC, Gomez-Gesteira M, Dalrymple RA (2007) Boundary conditions generated by dynamic particles in SPH methods.
Comput Mater Continua 5(3):173184
105. Gong K, Liu H (2007). A new boundary treatment for smoothed
particle hydrodynamics. In: Asian and Pacific coasts 2007, Nanjing, China
106. Hieber SE, Koumoutsakos P (2008) An immersed boundary
method for smoothed particle hydrodynamics of self-propelled
swimmers. J Comput Phys 227(19):86368654
107. Klessen R (1997) GRAPESPH with fully periodic boundary conditions: fragmentation of molecular clouds. Mon Not R Astron
Soc 292(1):1118
108. Randles PW, Libersky LD (2005) Boundary conditions for a dual
particle method. Comput Struct 83(1718):14761486
109. Yildiz M, Rook RA, Suleman A (2009) SPH with the multiple boundary tangent method. Int J Numer Methods Eng
77(10):14161438
110. Takeda H, Miyama SM, Sekiya M (1994) Numerical simulation
of viscous flow by smoothed particle hydrodynamics. Prog Theor
Phys 92(5):939960
111. Morris JP, Fox PJ, Zhu Y (1997) Modeling low Reynolds number incompressible flows using SPH. J Comput Phys 136(1):214
226
112. Monaghan JJ (1994) Simulating free surface flows with SPH.
J Comput Phys 110(2):399406
70
113. Campbell PM (1989) Some new algorithms for boundary value
problems in smooth particle hydrodynamics. Technical Report,
Mission Research Corp, Albuquerque, NM
114. Libersky LD, Petschek AG, Carney TC, Hipp JR, Allahdadi
FA (1993) High strain Lagrangian hydrodynamics: a threedimensional SPH code for dynamic material response. J Comput
Phys 109(1):6775
115. Rapaport DC (2004) The art of molecular dynamics simulation.
Cambridge University Press, Cambridge
116. Liu MB, Meakin P, Huang H (2007) Dissipative particle dynamics simulation of fluid motion through an unsaturated fracture and
fracture junction. J Comput Phys 222(1):110130
117. Tartakovsky AM, Meakin P (2005) Simulation of unsaturated
flow in complex fractures using smoothed particle hydrodynamics. Vadose Zone J 4(3):848855
118. Liu MB, Meakin P, Huang H (2007) Dissipative particle dynamics simulation of multiphase fluid flow in microchannels and microchannel networks. Phys Fluids 19(3):033302
119. Colagrossi A, Landrini M (2003) Numerical simulation of interfacial flows by smoothed particle hydrodynamics. J Comput Phys
191(2):448475
120. Owen JM, Villumsen JV, Shapiro PR, Martel H (1998) Adaptive
smoothed particle hydrodynamics: methodology, II. Astrophys J
Suppl Ser 116(2):155209
121. Shapiro PR, Martel H, Villumsen JV, Owen JM (1996) Adaptive
smoothed particle hydrodynamics, with application to cosmology: methodology. Astrophys J Suppl Ser 103(2):269330
122. Fulbright MS, Benz W, Davies MB (1995) A method of
smoothed particle hydrodynamics using spheroidal kernels. Astrophys J 440(1):254262
123. Liu MB, Liu GR (2005) Meshfree particle simulation of micro
channel flows with surface tension. Comput Mech 35(5):332
341
124. Liu MB, Liu GR, Lam KY (2006) Adaptive smoothed particle hydrodynamics for high strain hydrodynamics with material
strength. Shock Waves 15(1):2129
125. Swegle JW, Attaway SW (1995) On the feasibility of using
smoothed particle hydrodynamics for underwater explosion calculations. Comput Mech 17(3):151168
126. Attaway SW, Hendrickson BA, Plimpton SJ, Gardner DR,
Vaughan CT, Brown KH, Heinstein MW (1998) A parallel contact detection algorithm for transient solid dynamics simulations
using PRONTO3D. Comput Mech 22(2):143159
127. Campbell J, Vignjevic R, Libersky L (2000) A contact algorithm
for smoothed particle hydrodynamics. Comput Methods Appl
Mech Eng 184(1):4965
128. Drumm C, Tiwari S, Kuhnert J, Bart HJ (2008) Finite pointset
method for simulation of the liquid-liquid flow field in an extractor. Comput Chem Eng 32(12):29462957
129. Huang H, Dyka CT, Saigal S (2004) Hybrid particle methods in
frictionless impact-contact problems. Int J Numer Methods Eng
61:22502272
130. Li Y, Liu GR, Luan MT, Ky Dai, Zhong ZH, Li GY, Han X (2007)
Contact analysis for solids based on linearly conforming radial
point interpolation method. Comput Mech 39(4):537554
131. Limido J, Espinosa C, Salauen M, Lacome JL (2007) SPH
method applied to high speed cutting modelling. Int J Mech Sci
49(7):898908
132. Mehra V, Chaturvedi S (2006) High velocity impact of metal
sphere on thin metallic plates: a comparative smooth particle hydrodynamics study. J Comput Phys 212(1):318337
133. Parshikov AN, Medin SA (2002) Smoothed particle hydrodynamics using interparticle contact algorithms. J Comput Phys
180(1):358382
134. Parshikov AN, Medin SA, Loukashenko II, Milekhin VA (2000)
Improvements in SPH method by means of interparticle contact
M.B. Liu, G.R. Liu
135.
136.
137.
138.
139.
140.
141.
142.
143.
144.
145.
146.
147.
148.
149.
150.
151.
152.
153.
154.
155.
algorithm and analysis of perforation tests at moderate projectile

velocities. Int J Impact Eng 24(8):779796
Plimpton S, Attaway S, Hendrickson B, Swegle J, Vaughan C,
Gardner D (1998) Parallel transient dynamics simulations: algorithms for contact detection and smoothed particle hydrodynamics. J Parallel Distrib Comput 50(12):104122
Prochazka PP, Kravtsov AN, Peskova S (2008) Blast impact on
structures of underground parking. Undergr Spaces Des Eng Environ Asp 102:1119
Sauer RA, Li SF (2008) An atomistically enriched continuum
model for nanoscale contact mechanics and Its application to
contact scaling. J Nanosci Nanotechnol 8(7):37573773
Schieback C, Burzle F, Franzrahe K, Neder J, Henseler P, Mutter D, Schwierz N, Nielaba P (2009) Computer simulations of
complex many-body systems. In: High performance computing
in science and engineering 08, Stuttgart, Germany
Seo S, Min O, Lee J (2008) Application of an improved contact algorithm for penetration analysis in SPH. Int J Impact Eng
35(6):578588
Vignjevic R, Campbell J (1999) A penalty approach for contact in smoothed particle hydrodynamics. Int J Impact Eng
23(1):945956
Vignjevic R, De Vuyst T, Campbell JC, Source C (2006) A frictionless contact algorithm for meshless methods. Comput Model
Eng 13(1):35
Belytschko T, Neal MO (1991) Contact-impact by the pinball
algorithm with penalty and Lagrangian methods. Int J Numer
Liu MB, Liu GR, Lam KY, Zong Z (2003) Smoothed particle hydrodynamics for numerical simulation of underwater explosion.
Comput Mech 30(2):106118
Liu MB, Liu GR, Lam KY, Zong Z (2001) A new technique
to treat material interfaces for smoothed particle hydrodynamics. In: Computational mechanicsnew frontiers for new millennium. Elsevier Science, Amsterdam
Monaghan JJ (2000) SPH without a tensile instability. J Comput
Phys 159(2):290311
Vignjevic R, De Vuyst T, Campbell J (2002) The use of a
homogeneous repulsive force for contact treatment in sph. In:
WCCM V, Vienna, Austria
Balsara DS (1995) von Neumann stability analysis of smoothed
particle hydrodynamics suggestions for optimal algorithms.
Belytschko T, Xiao S (2002) Stability analysis of particle methods with corrected derivatives. Comput Math Appl 43(35):329
350
Bonet J, Kulasegaram S (2001) Remarks on tension instability of
Eulerian and Lagrangian corrected smooth particle hydrodynamics (CSPH) methods. Int J Numer Methods Eng 52:12031220
Lanson N, Vila JP (2007) Renormalized meshfree schemes,
I: consistency, stability, and hybrid methods for conservation
laws. SIAM J Numer Anal 46(4):19121934
Randles PW, Petschek AG, Libersky LD, Dyka CT (2003) Stability of DPD and SPH. In: Meshfree methods for partial differential
equations. Springer, Berlin
Sigalotti LDG, Lopez H (2008) Adaptive kernel estimation and
SPH tensile instability. Comput Math Appl 55(1):2350
Gray JP, Monaghan JJ, Swift RP (2001) SPH elastic dynamics.
Comput Methods Appl Mech Eng 190(49):66416662
Beissel S, Belytschko T (1996) Nodal integration of the elementfree Galerkin method. Comput Methods Appl Mech Eng 139(1
4):4974
Melean Y, Sigalotti LDG, Hasmy A (2004) On the SPH tensile
instability in forming viscous liquid drops. Comput Phys Commun 157(3):191200

156. Benz W (1988) Applications of smooth particle hydrodynamics (SPH) to astrophysical problems. Comput Phys Commun
48(1):97105
157. Frederic AR, James CL (1999) Smoothed particle hydrodynamics calculations of stellar interactions. J Comput Appl Math
109:213230
158. Hultman J, Pharayn A (1999) Hierarchical, dissipative formation
of elliptical galaxies: is thermal instability the key mechanism?
Hydrodynamic simulations including supernova feedback multiphase gas and metal enrichment in cdm: structure and dynamics
of elliptical galaxies. Astron Astrophys 347:769798
159. Thacker RJ, Couchman HMP (2001) Star formation, supernova
feedback, and the angular momentum problem in numerical cold
dark matter cosmogony: halfway there. Astrophys J 555(1):L17
L20
160. Monaghan JJ, Lattanzio JC (1991) A simulation of the collapse and fragmentation of cooling molecular clouds. Astrophys
J 375(1):177189
161. Berczik P (2000) Modeling the star formation in galaxies
using the chemo-dynamical SPH code. Astrophys Space Sci
271(2):103126
162. Lee WH, Kluzniak W (1999) Newtonian hydrodynamics of the
coalescence of black holes with neutron stars, II: tidally locked
binaries with a soft equation of state. Mon Not R Astron Soc
308(3):780794
163. Lee WH (2000) Newtonian hydrodynamics of the coalescence
of black holes with neutron stars, III: irrotational binaries with a
stiff equation of state. Mon Not R Astron Soc 318(2):606624
164. Senz DG, Bravo E, Woosley SE (1999) Single and multiple detonations in white dwarfs. Astron Astrophys 349:177188
165. Monaghan JJ (1990) Modelling the universe. Astron Soc Aust
Proc 8(3):233237
166. Hong JM, Lee HY, Yoon JC, Kim CH (2008) Bubbles alive.
ACM Trans Graph 27(3)
167. Monaghan JJ, Kocharyan A (1995) SPH simulation of multiphase flow. Comput Phys Commun 87:225235
168. Hu XY, Adams NA (2006) A multi-phase SPH method for
macroscopic and mesoscopic flows. J Comput Phys 213(2):844
861
169. Garg R, Narayanan C, Lakehal D, Subramaniam S (2007) Accurate numerical estimation of interphase momentum transfer in
Lagrangian-Eulerian simulations of dispersed two-phase flows.
Int J Multiph Flow 33(12):13371364
170. Liu J, Koshizuka S, Oka Y (2005) A hybrid particle-mesh
method for viscous, incompressible, multiphase flows. J Comput
Phys 202(1):6593
171. Password F (2003) Cosmological smoothed particle hydrodynamics simulations: a hybrid multiphase model for star formation. Mon Not R Astron Soc 339(2):289311
172. Ritchie BW, Pa Thomas (2001) Multiphase smoothed-particle
hydrodynamics. Mon Not R Astron Soc 323(3):743756
173. Hu XY, Adams NA (2009) A constant-density approach for incompressible multi-phase SPH. J Comput Phys 228(6):2082
2091
174. Hu XY, Adams NA (2007) An incompressible multi-phase SPH
method. J Comput Phys 227(1):264278
175. Dalrymple RA, Rogers BD (2006) Numerical modeling of water
waves with the SPH method. Coast Eng 53(23):141147
176. Crespo AJC, Gomez-Gesteira M, Dalrymple RA (2007) 3D SPH
simulation of large waves mitigation with a dike. J Hydraul Res
45(5):631642
177. Shao SD, Gotoh H (2004) Simulating coupled motion of progressive wave and floating curtain wall by SPH-LES model. Coast
Eng J 46(2):171202
178. Gotoh H, Shao SD, Memita T (2004) SPH-LES model for numerical investigation of wave interaction with partially immersed
breakwater. Coast Eng 46(1):3963
71
179. Shao SD (2006) Incompressible SPH simulation of wave breaking and overtopping with turbulence modelling. Int J Numer
Methods Fluids 50(5):597621
180. Gomez-Gesteira M, Dalrymple RA (2004) Using a threedimensional smoothed particle hydrodynamics method for wave
impact on a tall structure. J Waterw Port C 130(2):6369
181. Crespo AJC, Gomez-Gesteira M, Carracedo P, Dalrymple RA
(2008) Hybridation of generation propagation models and SPH
model to study severe sea states in Galician Coast. J Mar Syst
72(14):135144
182. Crespo AJC, Gomez-Gesteira M, Dalrymple RA (2008) Modeling dam break behavior over a wet bed by a SPH technique.
J Waterw Port C 134(6):313320
183. Gotoh H, Sakai T (2006) Key issues in the particle method for
computation of wave breaking. Coast Eng 53(23):171179
184. Issa R, Violeau D (2008) Modelling a plunging breaking solitary
wave with eddy-viscosity turbulent SPH models. Comput Mater
Continua 8(3):151164
185. Khayyer A, Gotoh H, Shao SD (2008) Corrected incompressible SPH method for accurate water-surface tracking in breaking
waves. Coast Eng 55(3):236250
186. Monaghan JJ, Kos A, Issa N (2003) Fluid motion generated by
impact. J Waterw Port C 129(6):250259
187. Panizzo A (2005). SPH modelling of underwater landslide generated waves. In: Proceedings of the 29th international conference
on coastal engineering 2004, Lisbon, Portugal
188. Qiu LC (2008) Two-dimensional SPH simulations of landslidegenerated water waves. J Hydraul Eng ASCE 134(5):668671
189. Rogers BD, Dalrymple RA (2005) SPH modeling of breaking
waves. In: Proceedings of the 29th international conference on
coastal engineering 2004, Lisbon, Portugal
190. Shao SD (2009) Incompressible SPH simulation of water entry
of a free-falling object. Int J Numer Methods Fluids 59(1):91
115
191. Shao SD, Ji CM, Graham DI, Reeve DE, James PW, Chadwick
AJ (2006) Simulation of wave overtopping by an incompressible
SPH model. Coast Eng 53(9):723735
192. Violeau D, Buvat C, Abed-Meraim K, de Nanteuil E (2007) Numerical modelling of boom and oil spill with SPH. Coast Eng
54(12):895913
193. Yim SC, Yuk D, Panizzo A, Di Risio M, Liu PLF (2008) Numerical simulations of wave generation by a vertical plunger using
RANS and SPH models. J Waterw Port C 134(3):143159
194. Zou S, Dalrymple RA (2005) Sediment suspension modeling by
smoothed particle hydrodynamics. In: Proceedings of the 29th international conference on coastal engineering 2004, Lisbon, Portugal
195. Bulgarelli UP (2005) The application of numerical methods for
the solution of some problems in free-surface hydrodynamics.
J Ship Res 49(4):288301
196. Iglesias AS, Rojas LP, Rodriguez RZ (2004) Simulation of antiroll tanks and sloshing type problems with smoothed particle hydrodynamics. Ocean Eng 31(89):11691192
197. Kim Y (2007) Experimental and numerical analyses of sloshing
flows. J Eng Math 58(14):191210
198. Lohner R, Yang C, Onate E (2006) On the simulation of flows
with violent free surface motion. Comput Methods Appl Mech
Eng 195(4143):55975620
199. Rhee SH, Engineer L (2005) Unstructured grid based Reynoldsaveraged Navier-Stokes method for liquid tank sloshing. J Fluid
Eng 127:572
200. Souto-Iglesias A, Delorme L, Perez-Rojas L, Abril-Perez S
(2006) Liquid moment amplitude assessment in sloshing type
problems with smooth particle hydrodynamics. Ocean Eng
33(1112):14621484
72
201. Bursik M, Martinez-Hackert B, Delgado H, Gonzalez-Huesca A
(2003) A smoothed-particle hydrodynamic automaton of landform degradation by overland flow. Geomorphology 53(12):25
44
202. Zhu Y, Fox PJ, Morris JP (1999) A pore-scale numerical model
for flow through porous media. Int J Numer Anal Methods
23(9):881904
203. Zhu Y, Fox PJ (2002) Simulation of pore-scale dispersion in
periodic porous media using smoothed particle hydrodynamics.
204. Tartakovsky AM, Meakin P (2006) Pore scale modeling of immiscible and miscible fluid flows using smoothed particle hydrodynamics. Adv Water Resour 29(10):14641478
205. Cleary PW, Prakash M (2004) Discrete-element modelling and
smoothed particle hydrodynamics: potential in the environmental
sciences. Philos Trans R Soc A 362(1822):20032030
206. Ghazali JN, Kamsin A (2008) A real time simulation and modeling of flood hazard. In: 12th WSEAS international conference
on systems, Heraklion, Greece
207. Kipfer P, Westermann R (2006) Realistic and interactive simulation of rivers. In: Graphics interface 2006, Quebec, Canada
208. Bui HH, Fukagawa R, Sako K (2006) Smoothed particle hydrodynamics for soil mechanics. Taylor & Francis, London
209. Bui HH, Fukagawa R, Sako K, Ohno S (2008) Lagrangian meshfree particles method (SPH) for large deformation and failure
flows of geomaterial using elastic-plastic soil constitutive model.
Int J Numer Anal Methods 32(12):15371570
210. Bui HH, Sako K, Fukagawa R (2007) Numerical simulation of
soil-water interaction using smoothed particle hydrodynamics
(SPH) method. J Terramech 44(5):339346
211. Gallati M, Braschi G, Falappi S (2005) SPH simulations of the
waves produced by a falling mass into a reservoir. Nuovo Cimento C 28(2):129140
212. Herrera PA, Massabo M, Beckie RD (2009) A meshless method
to simulate solute transport in heterogeneous porous media. Adv
Water Resour 32(3):413429
213. Hui HH, Fukagawa R, Sako K (2006) Smoothed particle hydrodynamics for soil mechanics. Terramechanics 26:4953
214. Laigle D, Lachamp P, Naaim M (2007) SPH-based numerical investigation of mudflow and other complex fluid flow interactions
with structures. Comput Geosci 11(4):297306
215. Maeda K, Sakai H (2007) Seepage failure analysis with evolution
of air bubbles by SPH. In: New frontiers in Chinese and Japanese
geotechniques, proceedings of the 3rd Sino-Japan geotechnical
symposium, Chongqing, China
216. McDougall S, Hungr O (2004) A model for the analysis of rapid
landslide motion across three-dimensional terrain. Can Geotech
J 41(6):10841097
217. McDougall S, Hungr O (2005) Dynamic modelling of entrainment in rapid landslides. Can Geotech J 42(5):14371448
218. Moresi L, Muhlhous H, Dufour F (2001) An overview of numerical methods for Earth simulations
219. Morris JP, Zhu Y, Fox PJ (1999) Parallel simulations of porescale flow through porous media. Comput Geotech 25(4):227
246
220. Pastor M, Haddad B, Sorbino G, Cuomo S, Drempetic V (2009)
A depth-integrated, coupled SPH model for flow-like landslides
and related phenomena. Int J Numer Anal Methods 33(2)
221. Sakai H, Maeda K (2006) Seepage failure of granular ground
accounting for soil-water-gas interaction. In: Geomechanics and
geotechnics of particulate media, proceedings of the international
symposium on geomechanics and geotechnics of particulate media, Ube, Yamaguchi, Japan
222. Cleary PW (1998) Modelling confined multi-material heat and
mass flows using SPH. Appl Math Model 22(12):981993
M.B. Liu, G.R. Liu

223. Jeong JH, Jhon MS, Halow JS, Van Osdol J (2003) Smoothed
particle hydrodynamics: applications to heat conduction. Comput
Phys Commun 153(1):7184
224. Jiang F, Sousa ACM (2006) SPH numerical modeling for
ballistic-diffusive heat conduction. Numer Heat Transfer B, Fundam 50(6):499515
225. Jiang FM, Sousa ACM (2006) SPH numerical modeling for
ballistic-diffusive heat conduction. Numer Heat Transfer B, Fundam 50(6):499515
226. Rook R, Yildiz M, Dost S (2007) Modeling transient heat transfer
using SPH and implicit time integration. Numer Heat Transfer B,
Fundam 51(1):123
227. Sousa ACM, Jiang FM (2007) SPH as an inverse numerical tool
for the prediction of diffusive properties in porous media. Diffus
Solids Liq Heat Transfer, Microstruct Prop 553:171189
228. Gutfraind R, Savage SB (1998) Flow of fractured ice through
wedge-shaped channels: smoothed particle hydrodynamics and
discrete-element simulations. Mech Mater 29(1):117
229. Oger L, Savage SB (1999) Smoothed particle hydrodynamics for
cohesive grains. Comput Methods Appl Mech Eng 180(1):169
183
230. Ji SY, Li H, Shen HT, Wang RX, Yue QJ (2007) A hybrid Lagrangian-Eulerian numerical model for sea-ice dynamics.
Acta Oceanol Sin 26:1224
231. Ji SY, Shen HT, Wang ZL, Shen HH, Yue QJ (2005)
A viscoelastic-plastic constitutive model with Mohr-Coulomb
yielding criterion for sea ice dynamics. Acta Oceanol Sin
24(4):5465
232. Schafer C, Speith R, Kley W (2007) Collisions between equalsized ice grain agglomerates. Astron Astrophys 470(2):733739
233. Shen HT, Su JS, Liu LW (2000) SPH simulation of river ice dynamics. J Comput Phys 165(2):752770
234. Wang RX, Ji SY, Shen HT, Yue QJ (2005) Modified PIC method
for sea ice dynamics. China Ocean Eng 19(3):457468
235. Tartakovsky A, Meakin P (2005) Modeling of surface tension
and contact angles with smoothed particle hydrodynamics. Phys
Rev E 72(2):26301
236. Nugent S, Posch HA (2000) Liquid drops and surface tension with smoothed particle applied mechanics. Phys Rev E
62(4):49684975
237. Apfel RE, Tian Y, Jankovsky J, Shi T, Chen X, Holt RG, Trinh E,
Croonquist A, Thornton KC, Sacco JA (1997) Free oscillations
and surfactant studies of superdeformed drops in microgravity.
Phys Rev Lett 78(10):19121915
238. Fang HS, Bao K, Wei JA, Zhang H, Wu EH, Zheng LL (2009)
Simulations of droplet spreading and solidification using an improved SPH model. Numer Heat Transfer A, Appl 55(2):124
143
239. Lopez H, Sigalotti L, Di G (2006) Oscillation of viscous drops
with smoothed particle hydrodynamics. Phys Rev E 73(5):51201
240. Melean Y, Sigalotti LD (2005) Coalescence of colliding van der
Waals liquid drops. Int J Heat Mass Transfer 48:40414061
241. Melean Y, Sigalotti LDG (2005) Coalescence of colliding van der
Waals liquid drops. Int J Heat Mass Transfer 48(1920):4041
4061
242. Wang W, Huang Y, Grujicic M, Chrisey DB (2008) Study of
impact-induced mechanical effects in cell direct writing using
smooth particle hydrodynamic method. J Manuf Sci E, Trans
ASME 130(2)
243. Zhang MY, Zhang H, Zheng LL (2007) Application of smoothed
particle hydrodynamics method to free surface and solidification
problems. Numer Heat Transfer A, Appl 52(4):299314
244. Zhang MY, Zhang H, Zheng LL (2008) Simulation of droplet
spreading, splashing and solidification using smoothed particle hydrodynamics method. Int J Heat Mass Transfer 51(13
14):34103419

245. Zhang MY, Zhang H, Zheng LL (2009) Numerical investigation
of substrate melting and deformation during thermal spray coating by SPH method. Plasma Chem Plasma Process 29(1):5568
246. Zhou GZ, Ge W, Li JH (2008) A revised surface tension model
for macro-scale particle methods. Powder Technol 183(1):2126
247. Liu MB, Liu GR, Lam KY, Zong Z (2003) Meshfree particle simulation of the detonation process for high explosives in shaped
charge unlined cavity configurations. Shock Waves 12(6):509
520
248. Liu MB, Liu GR, Zong Z, Lam KY (2003) Computer simulation of high explosive explosion using smoothed particle hydrodynamics methodology. Comput Fluids 32(3):305322
249. Liu MB, Liu GR, Lam KY (2002) Investigations into water mitigation using a meshless particle method. Shock Waves
12(3):181195
250. Liu MB, Liu GR (2004) Smoothed particle hydrodynamics: some
recent developments in theory and applications. J Beijing Polytech Univ 30:6171
251. Liu MB, Liu GR, Lam KY, Zong Z (2003) Computer simulation of shaped charge detonation using meshless particle method.
Fragblast 7(3):181202
252. Alia A, Souli M (2006) High explosive simulation using multimaterial formulations. Appl Therm Eng 26(10):10321042
253. Bromm V, Yoshida N, Hernquist L (2003) The first supernova
explosions in the universe. Astrophys J 596(2):L135L138
254. Busegnies Y, Francois J, Paulus G (2007) Unidimensional SPH
simulations of reactive shock tubes in an astrophysical perspective. Shock Waves 16(45):359389
255. Liu MB, Liu GR, Lam KY (2003) Comparative study of the real
and artificial detonation models in underwater explosions. Electron Model 25(2):113124
256. Mair HU (1999) Review: hydrocodes for structural response to
underwater explosions. Shock Vib 6(2):8196
257. Swegle JW (1992) Report at Sandia National Laboratories
258. Libersky LD, Petschek AG (1991) Smooth particle hydrodynamics with strength of materials. In: Advances in the freeLagrange method including contributions on adaptive gridding
and the smooth particle hydrodynamics method, proceedings of
the next free-Lagrange conference, Jackson Lake Lodge, Moran,
WY, USA
259. Zhou CE, Liu GR, Ky Lou (2007) Three-dimensional penetration
simulation using smoothed particle hydrodynamics. Int J Comput
Methods 4(4):671691
260. Benz W, Asphaug E (1995) Simulations of brittle solids using smooth particle hydrodynamics. Comput Phys Commun
87(1):253265
261. Cleary P, Ha J, Alguine V, Nguyen T (2002) Flow modelling in
casting processes. Appl Math Model 26(2):171190
262. Cleary PW, Ha J, Ahuja V (2000) High pressure die casting simulation using smoothed particle hydrodynamics. Int J Cast Met
Res 12(6):335355
263. Chen JS, Pan C, Roque C, Wang HP (1998) A Lagrangian reproducing kernel particle method for metal forming analysis. Comput Mech 22(3):289307
264. Cleary PW, Ha J (2000) Three dimensional modelling of high
pressure die casting. Int J Cast Met Res 12(6):357365
265. Cleary PW, Prakash M, Ha J (2006) Novel applications of
smoothed particle hydrodynamics (SPH) in metal forming.
J Mater Process Technol 177(13):4148
266. Cleary PW, Prakash M, Ha J, Stokes N, Scott C (2007) Smooth
particle hydrodynamics: status and future potential. Prog Comput
Fluid Dyn 7(24):7090
267. Ha J, Cleary PW (2000) Comparison of SPH simulations of high
pressure die casting with the experiments and VOF simulations
of Schmid and Klein. Int J Cast Met Res 12(6):409418
73
268. Ha J, Cleary PW (2005) Simulation of high pressure die filling of

a moderately complex industrial object using smoothed particle
hydrodynamics. Int J Cast Met Res 18(2):8192
269. Hu W, Yao LG, Hua ZZ (2007) Parallel point interpolation
method for three-dimensional metal forming simulations. Eng
Anal Bound Elem 31(4):326342
270. Prakash M, Cleary PW, Grandfield J, Rohan P, Nguyen V (2007)
Optimisation of ingot casting wheel design using SPH simulations. Prog Comput Fluid Dyn 7(24):101110
271. Ala G, Francomano E, Tortofici A, Toscano E, Viola F (2007)
Corrective meshless particle formulations for time domain
Maxwells equations. J Comput Appl Math 210(12):3446
272. Ala G, Francomano E, Tortorici A, Toscano E, Viola F (2006)
Smoothed particle electromagnetics: a mesh-free solver for transients. J Comput Appl Math 191(2):194205
273. Ala G, Francomano E, Tortorici A, Toscano E, Viola F, di Ingegneria Elettrica D, e delle Telecomunicazioni E (2007) A meshfree particle method for transient full-wave simulation. IEEE
Trans Magn 43(4):13331336
274. Borve S, Omang M, Trulsen J (2001) Regularized smoothed particle hydrodynamics: a new approach to simulating magnetohydrodynamic shocks. Astrophys J 561(1):8293
275. Dolag K, Bartelmann M, Lesch H (1999) SPH simulations of
magnetic fields in galaxy clusters. Astron Astrophys 348(2):351
363
276. Francomano E, Tortorici A, Toscano E, Ala G, Viola F (2009) On
the use of a meshless solver for PDEs governing electromagnetic
transients. Appl Math Comput 209(1):4251
277. Jiang FM, Oliveira MSA, Sousa ACM (2006) SPH simulation
of transition to turbulence for planar shear flow subjected to a
streamwise magnetic field. J Comput Phys 217(2):485501
278. Meglicki Z (1994) Verification and accuracy of smoothed particle
magnetohydrodynamics. Comput Phys Commun 81(12):91
104
279. Pimenta LCA, Mendes ML, Mesquita RC, Pereira GAS (2007)
Fluids in electrostatic fields: an analogy for multirobot control.
IEEE Trans Magn 43(4):17651768
280. Price DJ, Monaghan JJ (2004) Smoothed particle magnetohydrodynamics, I: algorithm and tests in one dimension. Mon Not R
Astron Soc 348(1):123138
281. Price DJ, Monaghan JJ (2004) Smoothed particle magnetohydrodynamics, II: variational principles and variable smoothinglength terms. Mon Not R Astron Soc 348(1):139152
282. Price DJ, Monaghan JJ (2004) Smoothed particle magnetohydrodynamics: some shocking results. Astrophys Space Sci
292(1):279283
283. Antoci C, Gallati M, Sibilla S (2007) Numerical simulation
of fluid-structure interaction by SPH. Comput Struct 85(11
14):879890
284. Prakash M, Cleary PW, Ha J, Noui-Mehidi MN, Blackburn H,
Brooks G (2007) Simulation of suspension of solids in a liquid
in a mixing tank using SPH and comparison with physical modelling experiments. Prog Comput Fluid Dyn 7(2):91100
285. Anghileri M, Castelletti LML, Tirelli M (2005) Fluid structure interaction of water filled tanks during the impact with the
ground. Int J Impact Eng 31(3):235254
286. Chikazawa Y, Koshizuka S, Oka Y (2001) A particle method
for elastic and visco-plastic structures and fluid-structure interactions. Comput Mech 27(2):97106
287. Guilcher PM, Ducrozet G, Doring M, Alessandrini B, Ferrant P
(2006) Numerical simulation of wave-body interactions using a
modified SPH solver. In: Proceedings of the sixteenth international offshore and polar engineering conference, San Francisco,
CA
288. Hosseini SM, Amanifard N (2007) Presenting a modified SPH algorithm for numerical studies of fluid-structure interaction problems. IJE Trans B, Appl 20:167178
74
289. Muller M (2004) Interactive blood simulation for virtual surgery
based on smoothed particle hydrodynamics. Technol Health Care
12(1):2531
290. Hieber SE (2004) Remeshed smoothed particle hydrodynamics
simulation of the mechanical behavior of human organs. Technol
Health Care 12(4):305314
291. Tanaka N, Takano T (2005) Microscopic-scale simulation of
blood flow using SPH method. Int J Comput Methods 2(4):555
568
292. Tsubota K, Wada S, Yamaguchi T (2006) Simulation study on
effects of hematocrit on blood flow properties using particle
method. J Biomech Sci Eng 1(1):159170
293. Rosswog S, Wagner P (2002) Towards a macroscopic modeling
of the complexity in traffic flow. Phys Rev E 65(3):36106
294. Benson DJ (1992) Computational methods in Lagrangian and
Eulerian hydrocodes. Comput Methods Appl Mech Eng 99(2
3):235394
295. Walters WP, Zukas JA (1989) Fundamentals of shaped charges.
Wiley, New York
296. Anderson CE (1987) An overview of the theory of hydrocodes.
Int J Impact Eng 5(14):3359
297. Anghileri M, Castelletti LML, Invernizzi F, Mascheroni M
(2005) A survey of numerical models for hail impact analysis
using explicit finite element codes. Int J Impact Eng 31(8):929
944
298. Attaway SW, Heinstein MW, Swegle JW (1993) Coupling of
smooth particle hydrodynamics with the finite element method.
In: IMPACT-4: structural mechanics in reactor technology conference, Berlin, Germany
299. Batra RC, Zhang GM (2008) Modified Smoothed Particle Hydrodynamics (MSPH) basis functions for meshless methods, and
their application to axisymmetric Taylor impact test. J Comput
Phys 227(3):19621981
300. Brown K, Attaway S, Plimpton S, Hendrickson B (2000) Parallel strategies for crash and impact simulations. Comput Methods
Appl Mech Eng 184:375390
301. De Vuyst T, Vignjevic R, Campbell JC (2005) Coupling between meshless and finite element methods. Int J Impact Eng
31(8):10541064
302. Fahrenthold EP, Koo JC (1997) Energy based particle hydrodynamics for hypervelocity impact simulation. Int J Impact Eng
20(15):253264
303. Fawaz Z, Zheng W, Behdinan K (2004) Numerical simulation
of normal and oblique ballistic impact on ceramic composite armours. Compos Struct 63(34):387395
304. Fountzoulas CG, Gazonas GA, Cheeseman BA (2007) Computational modeling of tungsten carbide sphere impact and penetration into high-strength-low-alloy (HSLA)-100 steel targets.
J Mech Mater Struct 2(10):19651979
305. Groenenboom PHL (1997) Numerical simulation of 2D and
3D hypervelocity impact using the SPH option in PAMSHOCK(TM). Int J Impact Eng 20(15):309323
306. Hayhurst CJ, Clegg RA (1997) Cylindrically symmetric SPH
simulations of hypervelocity impacts on thin plates. Int J Impact
Eng 20(15):337348
307. Hayhurst CJ, Clegg RA Livingstone IA, Francis NJ (1996) The
application of SPH techniques in Autodyn-2D to ballistic impact problems. In: 16th international symposium on ballistics,
San Francisco, USA
308. Hiermaier S, Schafer F (1999) Hypervelocity impact fragment
clouds in high pressure gas numerical and experimental investigations. Int J Impact Eng 23(1):391400
309. Ioan A, Brizzolara S, Viviani M, Couty N, Donner R, Hermundstad O, Kukkanen T, Malenica S, Temarel P (2007) Comparison of experimental and numerical impact loads on ship-like sections. In: Advancements in marine structures. Taylor & Francis,
London
M.B. Liu, G.R. Liu

310. Ji XY, Li YC, Huang XC, Chen G (2007) Coupled FE-SPH simulation of truncated-conical projectiles penetrating steel sheets. In:
7th international conference on shock & impact loads on structures, Beijing, China
311. Johnson AF, Holzapfel M (2003) Modelling soft body impact on
composite structures. Compos Struct 61(12):103113
312. Johnson GR (1994) Linking of Lagrangian particle methods to
standard finite element methods for high velocity impact computations. Nucl Eng Des 150(23):265274
313. Johnson GR, Beissel SR, Stryk RA (2000) A generalized particle
algorithm for high velocity impact computations. Comput Mech
25(2):245256
314. Johnson GR, Peterson EH, Stryrk RA (1993) Incorporation of an
sph option into the epic code for a wide range of high velocity
impact computations. Int J Impact Eng 14:385394
315. Johnson GR, Stryk RA, Beissel SR, Holmquist TJ (2002) An
algorithm to automatically convert distorted finite elements into
meshless particles during dynamic deformation. Int J Impact Eng
27(10):9971013
316. Kikuchi M, Miyamoto B (2006) Numerical simulation of impact crush/buckling of circular tube using SPH method. Key Eng
Mater 306308:697702
317. Kitsionas S, Whitworth AP (2007) High-resolution simulations
of clump-clump collisions using SPH with particle splitting. Mon
Not R Astron Soc 378(2):507524
318. Lam KY, Shen YG, Gong SW (2001) A study of axial impact of
composite rods using SPH approach. Shock Vib 8(5):303312
319. Lee M, Yoo YH (2001) Analysis of ceramic/metal armour systems. Int J Impact Eng 25(9):819829
320. Libersky LD, Randles PW, Carney TC, Dickinson DL (1997) Recent improvements in SPH modeling of hypervelocity impact. Int
J Impact Eng 20(6):525532
321. Lukyanov AA, Reveles JR, Vignjevic R, Campbell J (2005) Simulation of hypervelocity debris impact and spacecraft shielding
performance. In: Proceedings of the 4th European conference on
space debris
322. Ma S, Zhang X, Qiu XM (2009) Comparison study of MPM and
SPH in modeling hypervelocity impact problems. Int J Impact
Eng 36(2):272282
323. Michel Y, Chevalier JM, Durin C, Espinosa C, Malaise F, Barrau JJ (2006) Hypervelocity impacts on thin brittle targets: experimental data and SPH simulations. Int J Impact Eng 33(1
12):441451
324. Nizampatnam LS, Horn WJ (2008) Investigation of multimaterial bird models for predicting impact loads. In: Proceedings
of the Asme turbo expo, Berlin, Germany
325. Park YK, Fahrenthold EP (2005) A kernel free particle-finite element method for hypervelocity impact simulation. Int J Numer
326. Rabczuk T, Eibl J (2006) Modelling dynamic failure of concrete
with meshfree methods. Int J Impact Eng 32(11):18781897
327. Randles PW, Carney TC, Libersky LD, Renick JD, Petschek AG
(1994) Calculation of oblique impact and fracture of tungsten
cubes using smoothed particle hydrodynamics. Int J Impact Eng
17(4):661672
328. Rosen P, Popescu V, Hoffmann C, Irfanoglu A (2008) A highquality high-fidelity visualization of the September 11 attack
on the World Trade Center. IEEE Trans Vis Comput Graph
14(4):937947
329. Shintate K, Sekine H (2004) Numerical simulation of hypervelocity impacts of a projectile on laminated composite plate targets by means of improved SPH method. Compos Part A, Appl
Sci 35(6):683692
330. Stellingwerf RF, Wingate CA (1993) Impact modeling with
smooth particle hydrodynamics. In: Conference: smooth particle
hydrodynamics in astrophysics OAT workshop, Trieste, Italy

331. Taylor EA (2001) Simulation of hollow shaped charge jet impacts onto aluminium whipple bumpers at 11 km/s. Int J Impact
Eng 26(110):773784
332. Vignjevic R, Lepage S, De Vuyst T (2005) Simulation of high
velocity impacts on thin metallic targets II (discrete elements).
WIT Trans Eng Sci 49:575582
333. Wang P, Shao JL, Qin CS (2009) Effect of loading-wave-front
width on micro-jet from aluminum surface. Acta Phys Sin
58(2):10641070
334. Zhang QM, Long RR, Liu ZF, Huang FL (2008) SPH simulations
of hypervelocity impact of al spheres on multi-plate structures.
Int J Mod Phys B 22(911):16041611
335. McCarthy MA, Xiao JR, Petrinic N, Kamoulakos A, Melito V
(2004) Modelling of bird strike on an aircraft wing leading edge
eade from fibre metal laminates, part 1: material modelling. Appl
Compos Mater 11(5):295315
336. Brown K, Attaway S, Plimpton S, Hendrickson B (2000) Parallel strategies for crash and impact simulations. Comput Methods
Appl Mech Eng 184(2):375390
337. Bicknell GV, Gingold RA (1983) On tidal detonation of stars by
massive black holes. Astrophys J 273:749760
338. Fulbright MS, Benz W, Davies MB (1995) A method of
smoothed particle hydrodynamics using spheroidal kernels. Astrophys J 440:254
339. Chin GL (2001) Smoothed particle hydrodynamics and adaptive
smoothed particle hydrodynamics with strength of materials. National University of Singapore, Singapore
340. Hiermaier S, Konke D, Stilp AJ, Thoma K (1997) Computational simulation of the hypervelocity impact of Al-spheres on
thin plates of different materials. Int J Impact Eng 20(15):363
374
341. Bangash MYH (1993) Impact and explosion. Blackwell Scientific, Oxford
342. Cole RH (1948) Underwater explosions. Princeton University
Press, Princeton
343. Fickett WW (1966) Flow calculations for pulsating onedimensional detonations. Phys Fluids 9:903916
344. Mader CL (1998) Numerical modeling of explosives and propellants. CRC Press, New York
345. Libersky LD, Randles PW (2006) Shocks and discontinuities in
particle methods. Shock Compress Condens Matter 845:1089
1092
346. Ma S, Zhang X, Lian YP, Zhou X (2009) Simulation of high explosive explosion using adaptive material point method. Comput
Model Eng 39(2):101123
347. Miyoshi H (2008) Numerical simulation of shaped charges using the SPH solver: jet formation and target penetration. Explos
Shock Wave Hypervelocity Phenom Mater II
348. Tanaka K (2008) Numerical studies of explosive welding by
SPH. Explos Shock Wave Hypervelocity Phenom Mater II
566:6164
349. Wataru K, Akiko M (2006) Explosion simulation by smoothed
particle hydrodynamics. Comput Methods 2:13971403
350. Xu JX, Liu XL (2008) Analysis of structural response under blast
loads using the coupled SPH-FEM approach. J Zhejiang Univ Sci
A 9(9):11841192
351. Liu MB, Liu GR, Lam KY (2001) Simulation of the explosive
detonation process by using SPH methodology. In: Computational fluids and solid mechanics. Elsevier Science, Amsterdam
352. Dobratz BM (1981) LLNL explosive handbook. Lawrence Livermore National Laboratory, Livermore
353. Chiesum JE, Shin YS (1997) Explosion gas bubbles near simple
boundaries. Shock Vib 4(11):1125
354. Karniadakis GE, Beskok A, Aluru A (2005) Microflows and
nanoflows: fundamentals and simulation. Springer, Berlin
355. Bird GA (1994) Molecular gas dynamics and the direct simulation of gas flows. Oxford University Press, Oxford
75
356. Beskok A, Karniadakis GE (2001) Microflows: fundamentals

and simulation. Springer, Berlin
357. Pan L, Ng TY, Xu D, Liu GR, Lam KY (2002) Determination of
temperature jump coefficient using the direct simulation Monte
Carlo method. J Micromech Microeng 12(1):4152
358. Pan LS, Liu GR, Khoo BC, Song B (2000) A modified direct
simulation Monte Carlo method for low-speed microflows. J Micromech Microeng 10:2127
359. Pan LS, Liu GR, Lam KY (1999) Determination of slip coefficient for rarefied gas flows using direct simulation Monte Carlo.
J Micromech Microeng 9(1):8996
360. Duong-Hong D, Wang JS, Liu GR, Chen YZ, Han J, Hadjiconstantinou NG (2007) Dissipative particle dynamics simulations of electroosmotic flow in nano-fluidic devices. Microfluid
Nanofluid 3:17
361. Hoover WG, Hoover CG, Kum O, Castillo VM, Posch HA, Hess
S (1996) Smooth particle applied mechanics. Comput Methods
Sci Eng 2:6572
362. Hoover WG, Pierce TG, Hoover CG, Shugart JO, Stein CM,
Edwards AL (1994) Molecular dynamics, smoothed-particle applied mechanics, and irreversibility. Comput Math Appl 28(10
12):155174
363. Hoover WG, Hoover CG (2003) Links between microscopic and
macroscopic fluid mechanics. Mol Phys 101(11):15591573
364. Espanol P (1997) Fluid particle dynamics: a synthesis of dissipative particle dynamics and smoothed particle dynamics. Europhys Lett 39(6):605610
365. Espanol P, Revenga M (2003) Smoothed dissipative particle dynamics. Phys Rev E 67(2):026705
366. Liu MB, Liu GR, Lam KY (2003) Computer simulation of
flip-chip underfill encapsulation process using meshfree particle
method. Int J Comput Eng Sci 4(2):405408
367. Li Q, Cai TM, He GQ, Hu CB (2006) Droplet collision and coalescence model. Appl Math Mech 27(1):6773
368. Tartakovsky AM, Meakin P (2005) A smoothed particle hydrodynamics model for miscible flow in three-dimensional fractures
and the two-dimensional Rayleigh Taylor instability. J Comput
Phys 207(2):610624
369. Tartakovsky AM, Meakin P, Scheibe TD, Eichler West RM
(2007) Simulations of reactive transport and precipitation with
smoothed particle hydrodynamics. J Comput Phys 222(2):654
672
370. Tartakovsky AM, Redden G, Lichtner PC, Scheibe TD, Meakin
P (2008) Mixing-induced precipitation: experimental study and
multiscale numerical analysis. Water Resour Res 44(6)
371. Tartakovsky AM, Meakin P, Ward AL (2009) Smoothed particle hydrodynamics model of non-Aqueous phase liquid flow and
dissolution. Transp Porous Media 76(1):1134
372. Brackbill JU, Kothe DB, Zemach C (1992) A continuum method
for modeling surface tension. J Comput Phys 100(2):335354
373. Morris JP (2000) Simulating surface tension with smoothed particle hydrodynamics. Int J Numer Methods Fluids 33(3):333
353
374. Lo YME, Shao S (2002) Simulation of near-shore solitary wave
mechanics by an incompressible SPH method. Appl Ocean Res
24(5):275286
375. Shao S (2006) Incompressible SPH simulation of wave breaking
and overtopping with turbulence modelling. Int J Numer Methods Fluids 50(5)
376. Shao S, Ji C (2006) SPH computation of plunging waves using
a 2-D sub-particle scale (SPS) turbulence model. Int J Numer
Methods Fluids 51(8)
377. Shao SD (2005) SPH simulation of solitary wave interaction with
a curtain-type breakwater. J Hydraul Res 43(4):366375
378. Frank J, Reich S (2003) Conservation properties of smoothed
particle hydrodynamics applied to the shallow water equation.
BIT 43(1):4155
76
379. Ata R, Soulaimani A (2005) A stabilized SPH method for inviscid shallow water flows. Int J Numer Methods Fluids 47(2):139
159
380. Rodriguez-Paz M, Bonet J (2005) A corrected smooth particle hydrodynamics formulation of the shallow-water equations.
Comput Struct 83(1718):13961410
381. Idelsohn SR, Onate E, Pin FD (2004) The particle finite element method: a powerful tool to solve incompressible flows with
free-surfaces and breaking waves. Int J Numer Methods Eng
61(7):964989
382. Oger G, Alessandrini B, Ferrant P (2005) Capture of air cushion
effects in a wedge water entry SPH simulation. In: Proceedings
of the fifteenth international offshore and polar engineering conference
383. Callati M, Braschi G, Falappi S (2005) SPH simulations of the
waves produced by a falling mass into a reservoir. Nuovo Cimento C 28(2):129140
384. Oger G, Doring M, Alessandrini B, Ferrant P (2006) Twodimensional SPH simulations of wedge water entries. J Comput
Phys 213(2):803822
385. Delorme L, Iglesias AS, Perez SA (2005) Sloshing loads simulation in LNG tankers with SPH. In: International conference on
computational methods in marine engineering, Barcelona
386. Greco M, Landrini M, Faltinsen OM (2004) Impact flows and
loads on ship-deck structures. J Fluids Struct 19(3):251275
387. Liu MB, Liu GR, Zong Z (2008) An overview on smoothed particle hydrodynamics. Int J Comput Methods 5(1):135188
388. Xu QX, Shen RY (2008) Fluid-structure interaction of hydrodynamic damper during the rush into the water channel. J Hydrodyn
20(5):583590
389. Ataie-Ashtiani B, Shobeyri G (2008) Numerical simulation of
landslide impulsive waves by incompressible smoothed particle
hydrodynamics. Int J Numer Methods Fluids 56(2):209232
M.B. Liu, G.R. Liu

390. Scanlon BR, Tyler SW, Wierenga PJ (1997) Hydrologic issues in
arid, unsaturated systems and implications for contaminant transport. Rev Geophys 35(4):461490
391. Nativ R, Adar E, Dahan O, Geyh M (1995) Water recharge and
solute transport through the vadose zone of fractured chalk under
desert conditions. Water Resour Res 31(2):253261
392. Schwartz LM, Martys N, Bentz DP, Garboczi EJ, Torquato S
(1993) Cross-property relations and permeability estimation in
model porous media. Phys Rev E 48(6):45844591
393. Huang H, Meakin P, Liu MB (2005) Computer simulation
of two-phase immiscible fluid motion in unsaturated complex
fractures using a volume of fluid method. Water Resour Res.
doi:10.1029/2005WR004204
394. Huang H, Meakin P, Liu MB, McCreery GE (2005) Modeling of multiphase fluid motion in fracture intersections and
fracture networks. Geophys Res Lett 32:L19402. doi:10.1029/
2005GL023899
395. Snyder LJ, Stewart WE (1966) Velocity and pressure profiles
for Newtonian creeping flow in regular packed beds of spheres.
AIChE J 12(1):167173
396. Sussman M, Smereka P, Osher S (1994) A level set approach for
computing solutions to incompressible two-phase flow. J Comput
Phys 114(1):146159
397. Hirt CW, Nichols BD (1981) Volume of fluid (VOF) method for
the dynamics of free boundaries. J Comput Phys 39(1):201225
398. Unverdi SO, Tryggvason G (1992) A front-tracking method
for viscous, incompressible, multi-fluid flows. J Comput Phys
100(1):2537
399. Zhu Y, Fox PJ (2001) Smoothed particle hydrodynamics model
for diffusion through porous media. Transp Porous Media
43(3):441471

Smoothed Particle Hydrodynamics (SPH) : An Overview and Recent Developments

Uploaded by

Copyright:

Available Formats

Smoothed Particle Hydrodynamics (SPH) : An Overview and Recent Developments

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Smoothed Particle Hydrodynamics (SPH) : An Overview and Recent Developments

Uploaded by

Copyright:

Available Formats

Arch Comput Methods Eng (2010) 17: 2576

Smoothed Particle Hydrodynamics (SPH): an Overview

Abstract Smoothed particle hydrodynamics (SPH) is a

G.R. Liu is SMA Fellow, Singapore-MIT Alliance.

problems in engineering and science. It plays a valuable role

applications in many types of complicated problems. The

M.B. Liu, G.R. Liu

Smoothed Particle Hydrodynamics (SPH): an Overview and Recent Developments

the governing equations of the classical Newtonian hydrodynamics.

5. SPH is suitable for problems where the object under

discontinuous SPH (DSPH) [55]. Liu et al. also proposed

M.B. Liu, G.R. Liu

2 SPH Approximation Techniques

where f is a function of the position vector x, and (x x )

Smoothed Particle Hydrodynamics (SPH): an Overview and Recent Developments

where h is the smoothing length defining the influence or

A smoothing function W is usually chosen to be an even

This condition is also termed as unity condition since the

The third condition is the compact condition

[f (x) + f (x)(x x) + r((x x)2 )]

f (x) = f (x) + r(h2 ).

It is clear that SPH kernel approximation of an arbitrary field

where the divergence in the integral is operated with respect

where is a constant related to the smoothing function for a

Using (5) and (9), (8) becomes

the following equation is obtained

is the unit vector normal to the surface S. Since the

It is clear that the differential operation on a function is

M.B. Liu, G.R. Liu

Fig. 1 Schematic illustration of the scenarios in which the support

Fig. 2 Schematic illustration of the scenarios in which the support

2.3 Particle Approximation

each associated with field properties. These particles can be

Smoothed Particle Hydrodynamics (SPH): an Overview and Recent Developments

one particle i is associated with a fixed lumped volume Vi

where N is the total number of particles within the influence

where the gradient W in the above equation is evaluated at

The particle approximation is, however, related to some

The above two identities may be substituted into the integral

One of the good features for the above two equations is

M.B. Liu, G.R. Liu

f1 f2  = f1 f2 .

Hence, an SPH approximation of the sum of functions

It is clear that the SPH approximation operator is a linear

For convenience, the SPH approximation operator  is

x and time t are the independent variables. The total stress

Substituting the SPH approximations for a function and

where v ij = v i v j . Equation (28) is a set of commonly

This is another approach to obtain density directly from the

Smoothed Particle Hydrodynamics (SPH): an Overview and Recent Developments

This normalization property ensures that the integral of

The dimension of the compact support is defined by the

from a global operation to a local operation. This will

This property makes sure that as the smoothing length

M.B. Liu, G.R. Liu

smoothing functions have been used. Lucy in the original

where d is 5/4h, 5/h2 and 105/16h3 in one-, two- and

where d is 1/ 1/2 h, 1/h2 and 1/ 3/2 h3 , respectively,

In one-, two- and three-dimensional space, d = 1/ h,

f (x) = f (x) + r(h2 ).

one particle i is associated with a fixed lumped volume Vi

f1 f2 = f1 f2 .

For convenience, the SPH approximation operator is

Wij, vj + fi,