An Improved Incompressible Smoothed Particle Hydrodynamics Method and Its Application in Free-Surface Simulations
An Improved Incompressible Smoothed Particle Hydrodynamics Method and Its Application in Free-Surface Simulations
An Improved Incompressible Smoothed Particle Hydrodynamics Method and Its Application in Free-Surface Simulations
December 2009
Rui Xu
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1
CONTENTS
2
CONTENTS
3
CONTENTS
4
CONTENTS
5
List of Figures
1.11 Incompressible flows simulated by the SPH method [47, 55, 117] 34
6
LIST OF FIGURES
3.3 Mirror particle and dummy particle method for a flat wall . . . 75
7
LIST OF FIGURES
4.8 The geometry and initial velocity field of vortex spin-down case. 92
4.15 Relative error for the horizontal velocity component, Re = 100 . 104
4.16 Relative error for the horizontal velocity component, Re = 1000 105
8
LIST OF FIGURES
5.11 Velocity and pressure field in the bluff body case, Re = 20 . . . . 124
5.16 Velocity and pressure field in the bluff body case, Re = 100 . . . 131
5.17 The velocity profiles along different cross sections, Re = 100 . . 132
5.18 The pressure profiles along different cross sections, Re = 100 . . 133
5.19 Drag Cd and lift Cl coefficients from code Saturne, Re = 100 . . . 134
6.1 Relative error profiles at the right end along y/D = 2.5 . . . . . 141
∂ϕ
6.2 Relative error profile for along x/D = 2.5 . . . . . . . . . . . 142
∂y
9
LIST OF FIGURES
6.12 The impulsive paddle simulations with different P emax limits . 154
10
LIST OF FIGURES
6.26 Wave propagation along the channel, wave height H = 0.05m . . 167
6.28 Wave propagation along the channel, wave height H = 0.1m . . 168
6.29 Free surface comparison between ISPH and SAWW, H = 0.1m . 169
11
List of Tables
12
Abstract
Several incompressible SPH (ISPH) methods have been proposed in the past
decade. Here only the widely-used projection method for the pressure and ve-
locity coupling is considered. Based on open source code SPHysics [121], the
stability and accuracy of three methods which enforce either a divergence-
free velocity field, density invariance, or their combination are tested here
through the standard Taylor-Green and spin-down vortex problems. It is
invariant ISPH, can maintain accuracy and stability but at a high compu-
tational cost. A new divergence-free ISPH approach is proposed here which
maintains accuracy and stability without increasing computational cost by
slightly shifting particles away from streamlines while correcting their hy-
13
Abstract
attempts to improve the accuracy have been conducted. The influence of the
Through the benchmark test cases of the lid-driven cavity and bluff body, the
algorithm and the code have been validated. For lid-driven cavity case, three
different Re values, Re = 100, 400 and 1000, have been simulated here. How-
ever, it has been shown that ISPH method provides a slower convergence
rate, compared to the finite volume method (FVM), which is consistent with
the first-order spatial convergence rate, while the pressure predictions from
ISPH match the converged finite-volume results very well. For the bluff body,
cases with two different Reynolds numbers are simulated here. The lift and
drag coefficients are calculated from both ISPH and the FVM code Saturne.
With Re = 20, both methods give almost identical predictions, but differences
appear for higher Reynolds condition, Re = 100.
The new ISPH approach has also been applied in free-surface flow simula-
tions. A viscous damping method is proposed to overcome the instability
on the free surface caused by truncated kernels. The range of maximum
Peclet number for the calculation of the damping viscosity is also suggested,
which is demonstrated in this thesis by the collapsing water column problem
and other test cases. The instabilities are effectively damped by increasing
kinematic viscosity for particles on and adjacent to the free surface. The
algorithm is validated against an analytical solution for the flow due to an
impulsively started plate and (almost) exact solutions for wet bed dam break
problems at zero and small times. The method is also validated for progres-
sive regular waves with paddle motion defined by linear theory. The waves
14
Abstract
propagate without decay. The highly accurate predictions show that increas-
ing surface viscosity locally by an order of magnitude does not reduce the
accuracy by damping the surface instabilities, generated by the inevitable
errors associated with truncated kernels. The quantitative validation is be-
lieved to be most thorough for this purpose. Good agreements are obtained
15
Declaration
No portion of the work referred to in the thesis has been submitted in sup-
port of an application for another degree or qualification of this or any other
16
Copyright
right”) and s/he has given The University of Manchester certain rights
to use such Copyright, including for administrative purposes.
Designs and Patents Act 1988 (as amended) and regulations issued un-
der it or, where appropriate, in accordance with licensing agreements
which the University has from time to time. This page must form part
of any such copies made.
17
Copyright
tion and commercialisation of this thesis, the Copyright and any Intel-
lectual Property and/or Reproductions described in it may take place is
available in the University IP Policy (see http://www.campus.manches-
ter.ac.uk/medialibrary/policies/intellectual-property.pdf), in any relevant
18
Acknowledgments
ester) gave me encouragement, guidance and support from the initial to the
final level enabled me to overcome many difficulties during three-year re-
search. Professor Dominique Laurence (The University of Manchester) in-
spired me with lots of creative ideas. Dr Charles Moulinec (CCLRC Dares-
bury Laboratory) did proof reading of this thesis and provided me enormous
help on work during PhD years. Dr Juan Uribe (The University of Manch-
ester) supported me on various computing problems. Thanks to the MACE
Scholarship from the School of Mechanical, Aerospace and Civil Engineer-
ing, University of Manchester. With this financial support, I can finish my
PhD research.
house with me during the first year in Manchester and cheered me up with
Italian and French jokes. I wish to thank my best friends the Yuan family,
the Moulinec family and Mr Nicolas Chini. They helped me get through the
difficult times, and provided all the emotional support, comaraderie, enter-
tainment, and caring. I cannot forget the football time with lots of friends,
Ali, Nicolas, Eric etc., and the drink sessions with Sarah, Maurice, Murat,
19
Acknowledgments
Armando and other friends. Also thanks to Pourya, Roberto and Penny who
Jiyun Du. They gave life to me, raised me, support me, teach me, and love
me. I am also grateful to my wife Haixia Xu. Your love, patience, under-
standing and encouragement have upheld me.
20
Chapter 1
1.1 Background
L uid mechanics has two philosophies for the description of fluid flows:
F Lagrangian and Eulerian. In Lagrangian concept, the observation
point moves with the fluid element, with the observer moving with a velocity
identical to the fluid element, while in the Eulerian concept, the observer
keeps a fixed position without moving, with all the flow quantities as func-
based approaches are the finite volume method (FVM), the finite difference
method (FDM) and the finite element method (FEM). Several commercial
21
Chapter 1. Introduction and Thesis Objectives
Several Lagrangian methods, such as the vortex method [16], the Smoothed
Particle Hydrodynamics (SPH) method [67, 33], the Finite Point or Finite
Pointset Method (FPM) [83, 106, 107, 109], the Radial Basis Function Method
(RBFM) [50, 51, 79], the Finite Volume Particle Method (FVPM) [39, 80],
have been introduced during the last decades. The incompressible SPH
(ISPH) method as one of the SPH topics has been heavily studied due to
its accurate and efficient numerical performance. The aims of this thesis
are to investigate the existing ISPH methods, improving the accuracy and
efficiency of these methods and applying the ISPH method for accurate free-
surface simulations. The vortex method can date back to 1930s, when Rosen-
head [97] calculated Kelvin-Helmholtz instabilities by hand. And the mod-
ern developments have been promoted by Chorin [13], Leonard [64] and Re-
hbach [92]. As an extensively and deeply studied method, the discussion and
investigation to the vortex method will not be deployed in this thesis. More
details about this method can be found in textbooks [16, 69]. Here, a brief in-
troduction to FVPM, FPM and SPH will be presented. The reasons studying
Incompressible SPH (ISPH) method is explained in this chapter.
22
Chapter 1. Introduction and Thesis Objectives
fore, the author consider it is necessary to give a brief introduction about this
method.
FVPM was first introduced by Hietel et al. [39]. Recent developments have
been performed by Keck and Hietel [52, 53], Nestor et. al [80] and Teleaga
[104]. In FVPMs, the fluid is represented by a set of particles, which in
turn are associated with normalised, overlapping, compactly supported ker-
nel functions. A test function is defined related to the particle volume. The
particles are viewed as discrete volumes to which the integral form of the
governing equations applies. The interaction between particles is calculated
from the flux rate between neighbouring particles. The conservation law is
ensured. In contrast to the standard SPH method, the FVPM is conservative
In 2000, Schick [98] applied anisotropic kernel functions and variable kernel
supports in the FVPM. Keck and Hietel [53] simulated an inviscid vortex
advection problem by applying the projection method [12] in FVPM. In 2008,
23
Chapter 1. Introduction and Thesis Objectives
Figure 1.1: Density and velocity graphs around an oscillating circle [105].
24
Chapter 1. Introduction and Thesis Objectives
(a) Re=40
(b) Re=100
Figure 1.2: Pressure contours and velocity vectors around a circular cylinder
in the fluid field [81].
In the Finite Pointset Method (FPM) [83, 106, 107, 109], all the points in
the computing domain are only interpolation bases, without any relation to
the physical properties of the fluid, such as density and volume. The first
FPM was only recently proposed in [83], but the promising potential has
25
Chapter 1. Introduction and Thesis Objectives
been shown through its high accuracy prediction with disordered point dis-
26
Chapter 1. Introduction and Thesis Objectives
tions can be enforced analytically by the Taylor expansion without the help
of extra particles outside of the simulation domain, and the accuracy of this
method can be theoretically increased with extra computing cost by consider-
ing more terms in the Taylor expansion for each particle. The stability of the
is obtained. All those computing procedures make FPM very time consum-
ing.
27
Chapter 1. Introduction and Thesis Objectives
Figure 1.5: Evolution of an initially circular fluid patch: FPM solution of the
pressure field at times = 0.2, 0.4, 0.6, 2.0 s [109].
28
Chapter 1. Introduction and Thesis Objectives
[71]. The use of SPH has since widely expanded in fluid dynamics and solid
mechanics [76] although its original applications were in astrophysics. The
basic concept of SPH is that continuous media are represented by discrete
particles with volume, density and mass. The particles have a kernel func-
tion to define their range of interaction, and the hydrodynamic variable fields
are approximated by integral interpolations. Meshes are not needed in the
simulation, which is a major advantage of SPH over Eulerian methods for
complex geometries.
The original purpose of the SPH method was to perform astrophysical sim-
ulations. In astrophysics, non-linear interaction between cosmic objects and
non-linear coupling between different interactions, such as electromagnetic,
gravitational, as well as the disordered distribution of cosmic objects make
SPH was also used for the star collisions and galaxy formulation predictions
[84].
29
Chapter 1. Introduction and Thesis Objectives
30
Chapter 1. Introduction and Thesis Objectives
(b)
(a)
Figure 1.8: Impact of a water column (SPH fluid) on a thin plate: FEM shell
(left), SPH shell (right). (a) elastic plate; (b) elastoplastic plate [87].
tic/elastoplastic plate was simulated, and the coupling between SPH fluid
particles and finite element points was achieved, shown in Fig. 1.8. This
work also demonstrated that by only using the SPH method, the fluid-body
interaction can be predicted accurately without the help of FEM for the
structure change. In 2008, Das and Cleary [21] used the SPH method to
31
Chapter 1. Introduction and Thesis Objectives
(a)
(b)
(c)
Figure 1.9: Fracture process for the plate during collision (left: colored by
von Mises stress and right: colored by damage). (a) time = 60µs; (b) time
= 65µs; (c) time = 100µs [21].
32
Chapter 1. Introduction and Thesis Objectives
Although the SPH method is still used for astrophysical simulations, its ap-
plications in fluid simulations are most notable. From compressible flows to
incompressible flows, the SPH method shows reliable performance. In 2007,
Monaghan [74] successfully simulated a shock tube with a SPH Riemann
solver, and the density and velocity jump of the gas in the shock tubes were
well predicted, shown in Fig. 1.10. Moreover, Welton [118] incorporated el-
ements of SPH to extract mean quantities from the particles, including the
mean pressure gradient, in simulations for compressible turbulent flows.
Figure 1.10: Velocity and density profiles in shock tubes. Dot points : SPH
simulations; lines: analytical results. [74]
Despite the good performance of the SPH method in compressible flow simu-
lations, the SPH method reveals some difficulties in incompressible flow sim-
ulations, and this is now a main research focus. After the first application
of the SPH method in incompressible flow simulations in [72], simulations
of incompressible flows with the SPH method turned to be a critical area in
CFD. Fig. 1.11 shows the simulations of the internal flow, free-surface flow
and turbulence phenomena with the SPH method [47, 55, 117]. Compared to
33
Chapter 1. Introduction and Thesis Objectives
Figure 1.11: Simulations of incompressible flows by the SPH method. (a) flow
through an orifice conduit [55]; (b) benchmark dam-break case with parallel
computations [47]; (c) simulation of fish pass with turbulence models [117].
mesh-based methods, the meshless SPH method is more suitable for violent
deformation, such as violent change of free surfaces. Also, similar to its vor-
34
Chapter 1. Introduction and Thesis Objectives
[12, 17, 44, 100, 63]. The recent alternative approach of imposing the kine-
matic constraint of a constant volume for each fluid particle through non-
thermodynamic pressure is also potentially competitive as shown through
comparisons with WCSPH in [25]. In this work, only the projection-based
ISPH methods are studied.
flows causes problems, such as sound wave reflection at the boundaries and
the high sound speed leading to a small time step resulting in high comput-
ing cost. Small fluctuations in fluid density cause big errors on the pressure
field and numerical instability due to the use of the state equation [63].
cess that converged when density changes were below a specified tolerance
was used. A similar approach was used in [58], where a pressure Poisson
35
Chapter 1. Introduction and Thesis Objectives
equation was solved instead of a penalty-like method with source terms pro-
field and a curl-free field respectively. Shao and Lo [100] used an incompress-
ible method, similar to that in [58], to describe the free surface in dam-break
flow. Colin et al. [14] proposed an improved Laplacian operator in a method
similar to [17]. In 2007 and 2008, Lee et al. [62, 63] pointed out that a truly
ISPH method could improve the accuracy of the SPH method. In 2007, differ-
tion. A new ISPH algorithm is proposed in this thesis. The comparison of the
new algorithm to the existing ones is conducted. The accuracy and stability
of ISPH methods based on the projection algorithm is also discussed.
VOF researchers Hirt and Nichols [40] said, it is the flaw of the Eulerian ap-
proach, fluid elements are convected through the grid interface with the av-
36
Chapter 1. Introduction and Thesis Objectives
eraged velocity instead of its own, which make it difficult to describe the free
surface. This is the reason why Fox and McDonald [30], in their introductory
fluid mechanics text book, state "Clearly the type of analysis depends on the
problem. Where it is easy to keep track of identifiable elements of mass, we
use a method of description that follows the particles".
to simulate the solitary wave propagating along the beach [75]. Meanwhile,
instead of simulating only one phase in two-phase flows, Colagrossi and Lan-
drini [15], with WCSPH method, treated two-dimensional interfacial flows
with different fluids separated by sharp interfaces. In 2004, Rogers et al.
[96] used a combined Large Eddy Simulation (LES) type scheme or subparti-
cle scale (SPS) scheme to capture coherent turbulent structures in breaking
waves. Dalrymple and Rogers [18] applied a LES model to wave propagation
and interaction with coastal defence. In 2007, Violeau and Issa [116] mod-
eled the free-surface flow with the consideration of complex turbulence mod-
els. To better capture the physical process of the violent free-surface change,
Ferrari et al. [28] performed simulations involving millions of particles, run-
ning on modern massively parallel supercomputers. Although WCSPH pre-
dicted some highly transient flows quite well, notably dam break flows, pres-
sures were extremely noisy and the method highly dissipative. An important
development was made by Vila [114] who introduced a Riemann formulation
between interacting particles, reducing pressure noise markedly. It is how-
ever not clear here whether ISPH or Riemann-based SPH will finally win the
The ISPH method with properties of noise-free pressure field and bigger time
step limits is also preferable in free-surface flow simulations. In 2002 [66]
37
Chapter 1. Introduction and Thesis Objectives
and 2003 [100], Lo and Shao applied a projection-based ISPH method in the
[15] was also adopted in the ISPH method by Hu and Adams [44]. In 2008,
Lee et al. [63] applied a so-called Truly Incompressible SPH (TISPH) method
to simulate a 2-D dam break. In [54], the ISPH method similar to the one
introduced in [63] is also compared with the counterpart in [66, 100], with
previous one showing obvious improvement over the later method. But the
free-surface predictions are still noisy in [54]. The ISPH method is also ap-
plied here for free-surface flow simulations. The problem in ISPH method
for free-surface prediction is first analyzed here, and one possible solution
is proposed. Rigorous validation is conducted against analytical or highly
In this thesis, the basic SPH methodologies and incompressible SPH meth-
ods, including the new ISPH approach, will be introduced in Chapter 2. The
structure of the incompressible SPH code will be presented in Chapter 3.
In Chapter 4, the comparison between existing ISPH methods and the new
ISPH method is conducted, and the accuracy of the new ISPH method will
be investigated. In Chapter 5, the validation of the method and the code is
38
Chapter 1. Introduction and Thesis Objectives
lar bluff body in a channel. And the application of the new ISPH approach in
free-surface flows is presented in Chapter 6. Finally, conclusions are drawn.
39
Chapter 2
form, defined by Eq. 2.1 and 2.2. Eq. 2.3 is also solved for the particle
positions. SPH deals with operators for first and second derivative terms,
and the pressure gradient as well as the viscous terms and Laplacian opera-
tors are presented in the following. Incompressibility here is enforced in the
projection method by a pressure Poisson equation [12].
∇·u=0 (2.1)
Du 1
= − ∇P + ν∇2 u + F (2.2)
Dt ρ
Dr
=u (2.3)
Dt
40
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
position.
For a variable A, in SPH formalism, the value of A at a point r, with r = (x, y),
is written as a convolution product of the variable A with the Dirac δ function
Z
A(r) = A(r0 )δ(| r − r0 |)dr0. (2.4)
Z
A(r) ≈ A(r0 )ωh (| r − r0 |)dr0, (2.5)
Ω
can be written as
X
A(ri ) ≈ Vj A(rj )ωh (rij ). (2.6)
j
where Vj is the volume of particle j, and subscript j stands for the neighbor-
ing particle; rij is the distance between particle i and j. A quintic kernel is
used for all the SPH interpolation in this thesis [77]. It is written as
41
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
(3 − q)5 − 6 (2 − q)5 + 15 (1 − q)5 ,
0 ≤ q ≤ 1;
(3 − q)5 − 6 (2 − q)5 ,
7
1 ≤ q ≤ 2;
ωij = (2.7)
478π
(3 − q)5 , 2 ≤ q ≤ 3;
0,
q ≥ 3,
the initial particle distance, c is a constant, equal to 1.3 for all simulations
in this thesis. Some other kernels can be found in Appendix A.
or
X
∇A(r) ≈ Vj A(rj )∇ωh (rij ), (2.10)
j
X
∇ · A(r) ≈ Vj A(rj ) · ∇ωh (rij ). (2.11)
j
42
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
To solve the Navier-Stokes equations, the gradient and the divergence op-
erators are needed. It has been shown in the previous part that the basic
gradient or divergence operator can be obtained directly from SPH interpo-
where α is a variant scalar. Substituting Eq. 2.12 with Eq. 2.10, the SPH
X h i
n n n−1
α ∇A ≈ Vj (α A)j − nα A i αj ∇ωij (2.13)
j
X
∇Ai ≈ Vj (Ai + Aj )∇ωij (2.14)
j
X
∇Ai ≈ − Vj (Ai − Aj )∇ωij (2.15)
j
In [4], Bonet and Lok presented a correction method to preserve the angular
momentum in gradient or divergence operators. In 2007, Oger et al. [82]
analyzed the accuracy of different operators, and pointed out that with Eq.
2.15 the accuracy can be improved up to second order with the normaliza-
43
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
tion of the kernel function. The same normalization of the kernel was also
proposed by Vila [114]. The expression for the kernel normalization can be
expressed as
∇Wij = L(r)∇ωij . (2.16)
where r = (x, y), ∇Wij is substituted into Eq. 2.15 for ∇ωij as a normalized
kernel first derivative, and
∂ωij ∂ωij −1
P P
j Vj (xj − x) ∂x j
Vj (xj − x)
∂y
L(r) =
P ∂ωij P ∂ωij .
(2.17)
Vj (yj − y) Vj (yj − y)
j ∂x j ∂y
with the normalization, the high accuracy obtained by the normalization still
places Eq. 2.15 in a superior position to Eq. 2.14.
In 1997, Morris et al. [77] deduced the laminar viscosity term by combining
SPH and FDM approaches. The viscous term is presented in Eq. 2.18. In
this approximated viscous term, only the first derivative of the kernel func-
44
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
tion is used. The calculation of first derivative can be used not only in the
divergence or gradient calculation, but also for the viscous term. This saves
the computing time. One of the viscous term expressions used in this work
is identical to the one suggested by Morris [77],
where m is the mass of the particle; ρ is the fluid density; µ is the dynamic
An approximate Laplacian operator with the same format is used for scalars.
Here is the example of the pressure.
45
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
X rij · ∇ωij
Γαβ = (rij )α (rij )β . (2.22)
j
rij2
And α and β are the coordinate directions. λ is the diffusion coefficient, which
is assumed to be spatially variant. The Laplacian operator can be obtained
Since the first time the SPH method was proposed, three widely-used wall
boundary conditions have been used, the dummy particle method [63, 100],
the mirror particle method [77, 72] and the boundary force method [18].
It is very difficult to build some boundary conditions, such as the homoge-
neous Neumann boundary for the pressure, with the boundary force method.
Therefore, the dummy particle and the mirror particle methods are normally
used in incompressible SPH simulations. Both of the boundary methods
46
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
have been applied in the incompressible SPH code. The comparison has been
In the dummy particle method, several extra layers of particles are placed
inside the wall, with the same pressure and velocity values as corresponding
normal direction solid wall particles. The dummy particle method is applied
Fig. 2.1 (a). For the inner corner, dummy particles for the two walls, which
are perpendicular to each other, overlap with the same treatment for the
flat wall. Special treatment is done by using the averaged value for dummy
particles at the diagonal position, Da in Fig. 2.1 (b). For the curved wall,
dummy particles are located along the opposite direction of the wall normal.
Dummy particles are easy to implement. They are set up at the first time
step, without update as the time marches, which is not the case in mirror
particle method. The discussion about wall boundaries with dummy particle
method, especially about the inner/outer corner and curved wall, will not be
extended here. More detail about the dummy particle method can be found
in [63, 100]. Dummy particles are only used for the flat wall boundary in this
47
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
Figure 2.1: The special treatment for dummy particles at the inner or outer
corner and the curved wall. Blue circles: fluid particles; red circles: wall
particles; yellow circles: dummy particles; purple circles: dummy particles
for the special treatment at corners.
48
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
Figure 2.2: Mirror particle method for the static wall boundary. uf and um
are the velocities for the fluid particle and its corresponding mirror particle
respectively; pf and pm are the pressures for the fluid particle and its corre-
sponding mirror particle respectively; w is the wall velocity, and here w = 0.
The mirror particle method is firstly proposed by Morris et al. in 1997 [77]
the same pressure as their corresponding fluid particles, and their velocity
49
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
where ui and uif are the velocities for fluid particle i and its corresponding
mirror particle respectively, w is the wall velocity. At a corner, the special
treatment, similar to the counterpart in the dummy particle method, is
performed, shown in Fig. 2.2.
In [5, 27, 44], it was shown that increasingly irregular particle distribu-
tions exhibit increasing numerical errors in results. Fang and Parriaux [27]
pointed out that an ill-conditioned matrix in the linear system could appear,
with increasing non-uniformity of particle distribution.
In [76], Monaghan pointed out that the tensile instability in SPH results
in the clustering of particles. The clustering is particularly noticeable in
materials which have an equation of state which can give rise to negative
pressures, but it can occur in the fluid where the pressure is always positive.
It is a particular problem in solid body computations where the instability
may corrupt physical fragmentation by numerical fragmentation which, in
some cases, is so severe that the dynamics of the system is completely wrong.
50
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
0.6
0.4
0.2
y/D
-0.2
-0.4
-0.6
-0.6 -0.4 -0.2 0 0.2 0.4 0.6
x/D
Here, the author takes this instability problem as the effect of the kernel
together at stagnation points in the fluid field, shown in Fig. 2.3. However,
the flaw of the smoothing kernel introduces interpolation error. Fig. 2.4 (a)
shows the contour graph of the 2-D quintic kernel value [77], and 2.4 (b)
its first derivative. From the cross-section profile of the first derivative, at
x = 0, Fig. 2.4 (c), it can be seen that when particles are getting close to
each other within a certain distance range, the interaction between them is
not increased, but reduced. This non-physical behavior of the kernel func-
51
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
0.2
0.3
0.2
∂f/∂y
0
f
0.1
0
-0.2
-2 -2
-1 -1
-2 -2
0 0
Y X -1 -1
1 1
0 0
Y X
2 2
1 1
2 2
(a) (b)
0.3
0.2
0.1
∂f/∂y
-0.1
-0.2
-0.3
-2 -1 0 1 2
y
(c)
Figure 2.4: (a) contour graph for 2-D quintic kernel; (b) contour graph for the
first derivative, ∂f
∂y
, of 2-D quintic kernel; (c) profile for the first derivative of
quintic kernel, ∂f
∂y
, at cross section x = 0.
Although it is reported that the quadratic kernel can reduce this non-physical
behavior [18], its low-order accuracy and second derivative property, which
strongly influences the stability of the SPH method [77], limits its applica-
52
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
by Chaniotis et al. [9]. In this thesis, a new method, based on the projec-
tion method with slight particle shifting, is introduced to overcome particle
clustering without relying on the background uniform mesh. The shifting
method is first introduced by Nestor et al. [80] in the context of Finite Vol-
ume Particle Method (FVPM), where the particles are shifted based on the
local particle distribution.
time marching scheme is applied, where both the density and mass of par-
ticles are constant. Particle positions, rni , are updated with velocity uni to
positions r∗i ,
r∗i = rni + 4tuni (2.24)
53
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
where, Fni is a body force. This intermediate velocity u∗i can be projected onto
4t n+1
u∗i = un+1
i + ∇pi . (2.26)
ρ
2.26. Therefore,
4t n+1
un+1
i = u∗i − ∇pi . (2.28)
ρ
un+1 + uni
i
rn+1
i = rni + 4t . (2.29)
2
In an alternative formulation, Lee et al. [63] derived the velocity from the
pressure correction with the particles at rni before advancing their positions
with the total velocity, termed as truly incompressible SPH. This formulation
is also applied here and produces results identical to [17] for the test cases
investigated. An interesting point is that Eq. 2.27 is obtained from the conti-
nuity equation, ∇ · u = 0, therefore, the divergence operator in the Laplacian
54
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
should be the same as the one used in the velocity divergence calculation,
and the gradient operator should be the same as the pressure gradient cor-
rection part in Eq. 2.28, which means the compatibility for the operators.
However, Cummins and Rudman [17] proved that the compatible Laplacian
operator causes the checker-board effect in the pressure field, similar to the
this case, the density is not kept constant in the course of the simulation. An
intermediate velocity, u∗i , is calculated without the pressure gradient term
as before,
u∗i = uni + ν∇2 uni + Fni 4t.
(2.30)
culating the pressure field pn+1 through a Poisson equation with a velocity
divergence on the right hand side (R.H.S.), the pressure field is obtained by
solving a Poisson equation with a relative density difference on R.H.S., as
55
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
where ρ0 and ρ∗ are respectively the initial and temporal fluid density of each
particle, calculated as
X
ρi = mj ωij . (2.33)
j
4t n+1
un+1
i = uni − ∇p . (2.34)
ρ∗
scheme is used.
un+1 + uni
i
rn+1
i = rni + 4t (2.35)
2
Eq. 2.32 and Eq. 2.27 are actually equivalent, if one assumes the continuity
equation for the flow to be valid at the intermediate position, r∗i , with
1 dρ
(∇ · u)r∗ =− . (2.36)
i ρ dt r∗i
Eq. 2.32 can be obtained by substituting Eq. 2.36 into Eq. 2.27. It will be
shown in §4.1 and 4.2 in Chapter 4 that the numerical performance of both
Poisson equations, Eq. 2.27 and Eq. 2.32, are quite different.
In [44], Hu and Adams point out that if only a divergence-free velocity field
56
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
Similar to Shao and Lo’s method [100], to keep the density constant, repre-
sented as ρni = ρ0i , the particle positions are adjusted. First, intermediate
velocities and particle positions are obtained by
∗,n+1/2 4t
ui = uni + ν∇2 uni + Fni (2.37)
2
∗,n+1/2
ri∗,n+1 = rni + ui 4t (2.38)
∗
4t2 ρni − ρi∗,n+1
1
∇· ∇p = (2.39)
2 ρ i ρni
∗
1
where ∇p is the intermediate pressure gradient to adjust the particle
ρ i
position, with ρi = ρ0i . If we define
X
σi = ωij (2.40)
j
then
X
ρi = mj ωij = mi σi . (2.41)
j
57
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
With the relation ρni = ρ0i = mi σi0 and ρi∗,n+1 = mi σi∗,n+1 , Eq. 2.39 can be
written as
∗
4t2 σi0 − σi∗,n+1
1
∇· ∇p = . (2.42)
2 ρ i σi0
After the pressure is calculated, the particle positions are adjusted through
∗
1
the pressure gradient ∇p following
ρ i
∗
4t2
1
rn+1
i = ri∗,n+1 − ∇p (2.43)
ρ i 2
∗,n+1/2 4t
∗,n+1/2 ∗,n+1/2
ui∗,n+1 = ui + ν∇2 ui + Fi . (2.45)
2
58
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
4t n+1
un+1
i = ui∗,n+1 − ∇pi (2.46)
2ρ
where, the pressure field is obtained from the pressure Poisson equation,
shown as
1 n+1 2
∇· ∇p = ∇ · ui∗,n+1 , (2.47)
ρ i 4t
with ρi = ρ0i . An interesting point for this method is that the pressure and
corresponding gradient fields are almost doubled due to the half-time-step
operation. In the second half step, the value of ui∗,n+1 is actually the same as
that in ISPH_DF. Comparing Eq. 2.27 and 2.47, the R.H.S., is doubled in Eq.
2.47. Therefore, the pressure and corresponding gradient fields are doubled
in ISPH_DFDI, which will be shown in §4.2 in Chapter 4. A noteworthy
59
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
as shown in Fig. 2.5 (a). It is observed in the simulations that this non-
uniform particle distribution hinders the convergence of linear solvers. Also,
on the one hand, if the streamlines direct the particles towards each other,
taking Taylor-Green vortices as an example, Fig. 2.3, the particles will clus-
ter; on the other hand, the non-physical behavior of the kernel, introduced in
Section 2.3, will weaken the interaction between particles when particle dis-
tances are within a close range. Then particles will continue clustering if the
inertial motion is big enough, which is quantified by the Reynolds number,
and
UD
Re = (2.48)
ν
this kernel flaw, also causes the instability in the projection-based ISPH.
The following method stabilizes the accurate ISPH_DF method. The parti-
cle distribution is effectively well maintained, shown in Fig. 2.5 (b). The
pressure field is calculated as in ISPH_DF. And the particles are advanced,
shifted slightly, and accordingly, the hydrodynamic variables corrected by
60
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
0
where, φ is a general variable; i and i are the particle’s old and new position
respectively; δrii0 is the distance vector between the particle’s new and old
position. The method with particle shifting is called ISPH_DFS. Here, only
the first two terms are applied in simulations, giving an order consistent
with the Laplacian operator. Higher order accuracy may be achieved with
additional terms. With the position shifting and interpolation, the particles
can jump from one streamline to another, and the particle clustering can be
avoided; the error caused by highly-distorted particle spacings is reduced.
0.6 0.6
0.4 0.4
0.2 0.2
y/D
y/D
0 0
-0.2 -0.2
-0.4 -0.4
-0.6 -0.6
-0.6 -0.4 -0.2 0 0.2 0.4 0.6 -0.6 -0.4 -0.2 0 0.2 0.4 0.6
x/D x/D
(a) (b)
The idea of shifting particle positions has been proposed previously in the
context of the Finite Volume Particle Method (FVPM) by Nestor et al. [80].
This allows the use of particle velocities with a correction, u0 , added to con-
serve fluid momentum. Different from FVPM, the adjustment of particle
61
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
velocity here introduces sink or source terms to the total momentum. The
position shift is applied to particles, and the use of a second-order Taylor ex-
pansion interpolates hydrodynamic variable values at the new position, as
given in Eq. 2.49. Similar to the velocity correction equation, u0 , in [80],
but modifying the particle shifting magnitude, α, in relation to the particle
convection distance and the particle size, the position shift reads
Mi
X r̄i2
Ri = nij (2.51)
r2
j=1 ij
considered; rij is the distance between particle i and particle j; r̄i is the aver-
age particle spacing in the neighborhood of i, and
Mi
1 X
r̄i = rij ; (2.52)
Mi j=1
nij is the unit distance vector between particles i and j. The number of
neighboring particles Mi is the same as that in the kernel interpolation for
62
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
particles inside of the fluid, while for particles i close to the free surface,
only particles which have smaller distance to the particle i than the distance
of particle i to the free surface are considered. The summation of nij actu-
ally represents the anisotropy of the particle spacings. r̄i2/r 2 is used here as
ij
The shifting distance should be large enough to prevent instability and small
enough not to cause inaccuracy due to the Taylor expansion correction. Val-
ues of C within the range 0.01 - 0.1 satisfy these criteria for these test cases
and the dependence on Reynolds number was not observed. A value of 0.04
is generally used here. In [80], an average particle spacing, as defined in Eq.
2.52, is used as the shifting magnitude α, considering the influence of parti-
cle size. However, it is observed in simulations that when shifting distances
are much larger than convection distances, large numerical error appears in
the fields of hydrodynamic variables, even with the Taylor expansion update.
To avoid this an upper limit on shifting distance α is simply set as Umax 4t,
updated at each time step. Note also that the shifting distance is always
much less than the smoothing length h.
63
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
• Calculate the pressure from the Poisson equation, as shown in Eq. 2.27;
This small change in ISPH_DF makes this method much more robust, and
the efficiency of simulations is improved in comparison with ISPH_DFDI.
It should be pointed out here that it is not a strictly conservative method,
similar to LFPM [27]. However it will be shown in the test cases in§4.3.3 that
smooth prediction for the hydrodynamic field [63]. By applying the projec-
tion algorithm in SPH, the pressure field is smoother [63]. However, the
64
Chapter 2. The SPH Methodology, Existing and New ISPH Methods
concluded here:
ever, the instability is set in the simulations when the particle spacing
is highly distorted. Another approach, called ISPH_DI here, was intro-
duced by Shao and Lo [100], and the incompressibility is achieved by
maintaining the particle volume constant. In 2007, Hu and Adams pro-
65
Chapter 3
pressure field is calculated from a state equation without solving the pres-
sure Poisson equation. This code is an open-source SPH code that has been
released in 1 Aug 2007, developed jointly by researchers at the Johns Hop-
kins University (U.S.A.), the University of Vigo (Spain), the University of
operators shown in §2.1.2 are built. Two linear solvers, namely Conjugated
Gradient (CG) and Bi-Conjugated Gradient (Bi-CGSTAB), are programmed.
The predictor-corrector time marching scheme is applied.
66
Chapter 3. The Incompressible SPH Code
the particle distance to all the other particles by which the status, inside or
outside of the support domain, is decided, the operations will be (N − 1) · N
for a N particles case. In the code SPHysics, square cells with a edge size
equal to the kernel support size are placed in the simulation domain [112].
For each particle, the linked list is built up with the consideration of particles
in neighbouring cells, which highly reduced the searching time. It has been
proved that this method can reduce the operation number from (N − 1) · N to
NlogN for N particles [70].
Figure 3.1: Neighbouring list in the incompressible SPH code. The neigh-
bouring particles are those located in the adjacent cells. In 2-D, only NW, N,
NE and E cells are used for neighbouring particle searching.
Another trick to speedup the building of the linked list is to search only
particles in northwest (NW), north (N), northeast (NE), east (E) and polar
(P) cells, shown in Fig.3.1. As shown in Chapter 2, all the SPH operators are
symmetric or anti-symmetric. The calculation for different operator is only
67
Chapter 3. The Incompressible SPH Code
needed once for each pair of particles. The influences of particles in west (W),
southwest (SW), south (S) and southeast (SE) cells to particles in P cell have
been considered when building linked list for particles in W, SW, S and SE
cells. By only sweeping through the NW, N, NE and E cells, the computing
cost is minimized.
A similar linked list building procedure can be easily extended to 3-D. Dif-
ferent linked list techniques can be found in [113].
After particles are moved to the intermediate position r∗ by Eq. 2.24 in Chap-
ter 2, an intermediate velocity field u∗ is calculated as shown in Eq. 2.25 in
Chapter 2. If the viscous operator in [77] is used, u∗ is written as
68
Chapter 3. The Incompressible SPH Code
Similar expressions can be obtained for the operator defined in [99]. The
pressure gradient is not considered here, while it can be included for u∗ cal-
culation, and the pressure calculation for time n + 1 P n+1 from the Poisson
solver will be the correction to P n [17].
2.20 and Eq. 2.15 respectively. The pressure Poisson equation can be written
as
X Vj Pin+1 − Pjn+1 rij · ∇ωij
1 X ∗ ∗
2 = V j uj − uj · ∇ωij . (3.2)
j
ρj (rij2 + η 2 ) 4t j
Vj rij · ∇ωij
A(i, j) = −2 ; (3.3)
ρj (rij2 + η 2 )
X
A(i, i) = − A(i, j); (3.4)
j
1 X
Vj u∗j − u∗j · ∇ωij .
B(i) = (3.6)
4t j
69
Chapter 3. The Incompressible SPH Code
4t X
un+1 = u∗i − Vj Pjn+1 − Pin+1 · ∇ωij .
i (3.7)
ρ j
And at the end the particle position at time n + 1 rn+1 is updated following
Eq. 2.29.
70
Chapter 3. The Incompressible SPH Code
1
Bi-CGSTAB with Jacobi preconditioner
CG with Jacobi preconditioner
Normalized residual
0.01
0.0001
found in the Appendix B. In this thesis, the Bi-CGSTAB linear solver with
the Jacobi Preconditioner is used in all the simulations.
71
Chapter 3. The Incompressible SPH Code
can be decided by the relative difference of results from two consecutive ex-
In Weakly Compressible SPH (WCSPH) method, the choice of the time step
is based on three limitations:
All these three time-step constraints were adopted by Lee et al. in the
Truly incompressible SPH method [63]. However only the first two condi-
tions are considered to determine the time step in the most of ISPH litera-
ture as the third condition is actually limiting the distance between particles
[77, 46] which can be avoided by the shifting of the particle position shown
in §2.2.4 or simply overcome by the algorithm, such as ISPH_DI [100] and
ISPH_DFDI [44]. Even with ISPH_DF [17] algorithm, the particle cluster-
ing does not depend on the time step, which has been tested in simulations of
Taylor-Green vortices. In [17], Cummins and Rudman stated that the time-
step constraint is different depending on the resolutions and the viscosity.
Also, it is presented in [100] that the viscous diffusion constraint is used as a
stability condition in the explicit finite difference method simulating viscous
flows, and for simulations with high resolution or large viscosity, the diffu-
sion constraint is more stringent than the CFL condition. In [44], in addition
72
Chapter 3. The Incompressible SPH Code
to the CFL and diffusion condition, Hu and Adams were also considering the
surface tension condition. In this work, both the CFL and diffusion condi-
tions are considered as time step criteria.
CFL condition
h
4tCF L ≤ σ (3.8)
uref
which is proportional to the initial particle spacing δr; uref can be the numer-
ical speed of sound c0 for WCSPH and the maximum fluid velocity umax for
ISPH.
h2
4tν ≤ ϑ min (3.9)
νE,i
where, ϑ is a constant, and here it is set as 0.2; νE,i is the effective viscosity
h2
for the particle i, and the minimum value for is taken into account. The
νE,i
effective viscosity νE,i is the kinematic viscosity ν in this work, as all the
flows are considered as laminar flows.
73
Chapter 3. The Incompressible SPH Code
It has been introduced in §2.1.3 that the wall boundaries can be achieved
through the mirror particle or dummy particle method. With the mirror
particle method, the Neumann boundary can be easily built up through al-
gebraic operations in the matrix. If particle jm is the corresponding mirror
particle of the fluid particle j as shown in Fig. 3.3, the discretized Poisson
... + P (j) · a(i, j) + ... + P (i) · a(i, i) + ... + P (jm ) · a(i, jm ) + ... = b(i). (3.10)
... + P (j) · [a(i, j) + a(i, jm )] + ... + P (i) · a(i, i) + ... = b(i). (3.12)
A similar linear system can be also obtained for dummy particle boundary,
... + P (j) · a(i, j) + ... + P (i) · a(i, i) + ... + P (k) · [a(i, k) + a(i, kd1 ) + a(i, kd2 )] = b(i).
(3.13)
Tests for linear solver and the wall boundary condition are conducted with
the Laplacian operator in [77]. A 1-D channel with periodic boundaries in
d2 P
horizontal direction is set up. The Poisson equation dy 2
= 1 is solved with
dP
a homogeneous Neumann wall boundary, dy
= 0 at y = 1, and a Dirich-
let boundary, P = 0 at y = 0. The analytical solution for this problem
74
Chapter 3. The Incompressible SPH Code
1.0 × 10−7 for the normalized residual. The Neumann wall boundary is set
up with both dummy particle and mirror particle methods. In Fig. 3.4 the
Poisson equation is solved with the same Laplacian operator Eq. 2.20 in
§2.1.2 in Chapter 2 with a vertical resolution of 41 particles. The mirror
particle method provides more accurate predictions for the wall boundaries
than the dummy particle method. The non-uniform distribution of particles
in the practical simulation may however complicate the accuracy behavior
of the mirror-particle boundary method. The convergence is tested with dif-
ferent vertical resolutions, 40, 50 and 100 particles, and compared with the
Figure 3.3: Mirror particle and dummy particle wall boundary. Black balls
are fluid particles, blue wall particles, red mirror particles or dummy parti-
cles.
75
Chapter 3. The Incompressible SPH Code
analytical P=0.5*y*y-y
-0.1 mirror particle boundary
dummy particle boundary
-0.2
-0.3
P
-0.4
-0.5
in the linked list for periodic boundaries. Therefore for particles close to
different boundary ends, the distances between them are calculated from
subtracting their coordinate difference from the domain length. For example,
the horizontal distance 4x for particles, i and j, at the two ends respectively
76
Chapter 3. The Incompressible SPH Code
analytical results
-0.1 vertical resolution 40
vertical resolution 50
vertical resolution 100
-0.2
-0.3
P
-0.4
-0.5
77
Chapter 3. The Incompressible SPH Code
The code structure is detailed in Fig. 3.7. The geometry of the case is defined
in a separate file, datclass.f. The geometric data is read into the solver. Af-
ter choosing the ISPH solvers, ISPH_DI, ISPH_DFDI and ISPH_DFS, and
the boundary methods, mirror particle and dummy particle boundary, the
solver is ready for simulations. The ISPH_DF solver can be simply obtained
by commenting out the particle shifting part in the ISPH_DFS solver. For
the ISPH_DI, ISPH_DFDI and ISPH_DFS solvers, the corresponding vis-
cous term and intermediate velocity/density will be calculated, and then the
discretized Laplacian operator and the right hand side (R.H.S.), dependent
on the solvers, are calculated, and the linear problem is solved. The veloc-
ity field is corrected with the pressure field obtained from the linear solver,
and particle positions are updated. Differently, particle position are shifted
and correction for hydrodynamic values are done in ISPH_DFS simulations.
For ISPH_DFDI solver, an internal iteration for the converged particle posi-
tion is conducted, a tolerance for the relative density change, defined as the
R.H.S. in Eq.2.39, is normally 1% or 0.5%, and the maximum internal iter-
ation normally is 20. The pressure and velocity field are corrected in a way
similar to ISPH_DF.
Four incompressible SPH algorithms are applied with the help of open-source
78
Chapter 3. The Incompressible SPH Code
79
Chapter 3. The Incompressible SPH Code
tion and the viscous diffusion condition. Both dummy-particle and mirror-
particle methods can be activated in the code to build up wall boundary
condition. It has been shown that mirror-particle boundary method pro-
vides more accurate predictions for the wall boundaries than dummy particle
method. The periodic boundaries are also available in the code.
80
Chapter 4
methods have been proposed in the past decade. This thesis is con-
cerned only with the widely-used projection method for the pressure and ve-
using the relative density difference as the right hand side (R.H.S.) in the
pressure Poisson equation. But this method is inaccurate and noisy. The
81
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
the full time step into two half time steps: 1) Within the first half time step,
the constant particle density is kept by an internal iteration, where the par-
ticles are shifted in a way similar to that in ISPH_DI method; 2) During the
second half time step, the divergence-free velocity field is obtained by solving
In this chapter, the accuracy and stability of the four projection-based ISPH
methods are tested through two academic cases, Taylor-Green vortices and
vortex spindown flows. Comparisons in terms of accuracy, stability and com-
puting time are carried out.
82
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
1
P = − e2κt [cos(4πx) + cos(4πy)] . (4.2)
4
U is the velocity scale, equal to 1.0 m/s here; kinematic viscosity ν are 0.1 m2 /s,
0.01 m2 /s and 0.001 m2/s in three runs with three different Reynolds numbers,
2
Re = 10, Re = 100, Re = 1000; κ = − 8π
Re
is the decay rate of the velocity field;
u and v are the horizontal and vertical velocity components respectively. Dif-
ferent resolutions are used to investigate the spatial accuracy of ISPH_DFS,
which will be shown in Chapter 4. The Reynolds number is calculated by
UD
Re = (4.3)
ν
where, D is the length of the unit square side. Because of the existence of
uations, the particle spacings are highly compressed in one direction, but
stretched in the other roughly normal direction, with ISPH_DF, shown in
83
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
0.6
0.4
0.2
y/D
-0.2
-0.4
-0.6
-0.6 -0.4 -0.2 0 0.2 0.4 0.6
x/D
particle spacings are compressed in one direction, and stretched in the other
direction roughly normal, as shown in Fig.4.1. When particles cluster to-
gether as happens at the four points (±0.25 m, ±0.25 m) in Fig.2.3, the error
caused by the kernel flaw, presented in §2.1.4, will increase with decrease
84
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
in the particle spacing. The clustering cannot be avoided for high Reynolds
0.5 0.5
0.4 0.4
0.3 0.3
0.2 0.2
0.1 0.1
u/U
v/U
0 0
-0.1 -0.1
-0.2 -0.2
-0.3 -0.3
-0.4 -0.4
-0.5 -0.5
-0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4
y/D x/D
(a) (b)
0.2
0.15
0.1
0.05
P/(ρU 2)
-0.05
-0.1
-0.15
-0.2
-0.25
-0.6 -0.4 -0.2 0 0.2 0.4 0.6
y/D
(c)
85
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
0.5
0.4 0.4
0.3
0.2 0.2
0.1
u/U
v/U
0 0
-0.1
-0.2 -0.2
-0.3
-0.4 -0.4
-0.5
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.2 0 0.2 0.4
y/D x/D
(a) (b)
0.4
0.2
P/(ρU 2)
-0.2
-0.4
-0.4 -0.2 0 0.2 0.4
y/D
(c)
86
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
(a) (b)
2
P/(ρU 2)
-2
(c)
87
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
reason, for cases with high Reynolds numbers Re = 100 and Re = 1000 results
from ISPH_DF method are not obtainable and not plotted. It can be observed
in Fig.4.2 that all four projection-based ISPH methods can provide very good
prediction for both velocity and pressure fields with Re = 10. When the
Reynolds number increases to 100, ISPH_DF could not stably simulate the
flow development, while the other three methods could continue simulations.
In Fig.4.3, it is shown that ISPH_DFDI and ISPH_DFS predict both velocity
and pressure field accurately; ISPH_DI underpredicts the velocities at x =
±0.25 m and y = ±0.25 m. Increasing Reynolds number to 1000, ISPH_DFS
provides smooth velocity and pressure profiles, but slightly underpredicts
Fig.4.5, 4.6 and 4.7 present the maximum velocity magnitude, umax , decaying
against time for Re = 10, Re = 100 and Re = 1000 respectively. For the low
Reynolds case, Re = 10, all the methods stably simulate the flow develop-
ment, and ISPH_DFS gives the most accurate prediction, shown in Fig.4.5.
Increasing Reynolds number to, or over 100, only ISPH_DI, ISPH_DFDI and
ISPH_DFS are presented, due to stability reasons. From Fig.4.6 and 4.7, it
can be seen that ISPH_DI does not provide accurate prediction. ISPH_DFDI
could simulate the flow with Re = 100. But with Re = 1000, and the same
particle resolution, 40 × 40, certain numerical noise appears on the profile at
the beginning of the flow development, shown in the enlarged part in Fig.4.7.
Results with Re = 1000 were also not reported in [44]. Because ISPH_DFS
88
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
0
10
-1
10
umax/U
10-2
10-3
-4
10
0 0.2 0.4 0.6 0.8 1
time/T
In Table 4.1 , the computing costs for the four methods are listed. ISPH_DFDI
is the most time-consuming method. The computing cost of ISPH_DFS method
All the four methods could manage to simulate the flow development with
a low Reynolds number, Re = 10. ISPH_DF cannot stably simulate higher
Reynolds cases. Although ISPH_DI could keep well distributed particle spac-
ings, the numerical noise generated by this method contaminates the results
89
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
100
-1
10
U max/U
10-2
-3
10
-4
10
0 2 4 6 8 10
time/T
90
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
Re = 10 Re = 100 Re = 1000
ISPH_DF 1,103 s unstable unstable
ISPH_DI 1,091 s 1,250 s 1,323 s
ISPH_DFDI 4,000 s 5,352 s 5,625 s
ISPH_DFS 1,213 s 1,281 s 1,387 s
Table 4.1: Computing costs comparison for Taylor-Green vortices with a res-
olution of 40 × 40. Physical time = 2.0 s. The same time step is used for all
the methods for a given Reynolds number.
with Re = 1000, small numerical noise appears on the profile, shown in Fig.
4.7. ISPH_DFS efficiently stabilizes all simulations, and supplies predictions
with less numerical noise for a given resolution.
u = U(y − 0.5)
(4.4)
v = U(0.5 − x)
91
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
ISPH methods. The Reynolds number calculation is defined in the same way
as for the case Taylor-Green vortices, Eq.4.3. Also, D, U and ν are the same
as those in Taylor-Green vortices.
The accuracy, stability and computing time are compared among all the four
ISPH methods. During simulations, it is found that particle spacings cannot
be well maintained if ISPH_DF is used. It is the same as the findings in [27]
that the irregular particle distribution hinders the convergence of the linear
0.8
0.6
Y/D
0.4
0.2
Figure 4.8: The geometry and initial velocity field of vortex spin-down case.
Black solid line: wall; black lines with arrow: streamlines; circles: fluid
particles in ISPH simulations.
92
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
with a resolution of 160 × 160 and are used for the comparison. With Re = 10,
all four methods manage to provide accurate predictions, shown in Fig.4.9.
However, when the Reynolds number increases to 100, the simulation be-
came unstable with ISPH_DF. As with the Taylor-Green vortex simulations,
particle spacings are strongly distorted, which makes simulations quite un-
stable.
Fig.4.9 and 4.12 present the simulation results with Re = 10 situation us-
ing all four ISPH methods. All the methods provide very good prediction for
the velocity decay. Note that for ISPH_DFDI pressure is used as an inter-
93
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
0.4 0.4
0.2 0.2
v/U
u/U
0 0
-0.2 -0.2
-0.4 -0.4
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
y/D x/D
(a) (b)
0.01
-0.01
-0.02
-0.03
P/(ρU )
2
-0.04
-0.05
-0.06
-0.07
-0.08
-0.09
-0.1
0 0.2 0.4 0.6 0.8 1
y/D
(c)
94
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
0.4 0.4
0.3 0.3
0.2 0.2
0.1 0.1
u/U
v/U
0 0
-0.1 -0.1
-0.2 -0.2
-0.3 -0.3
-0.4 -0.4
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
y/D x/D
(a) (b)
0.4
0.3
0.2
0.1
P/(ρU 2)
-0.1
-0.2
-0.3
-0.4
0 0.2 0.4 0.6 0.8 1
y/D
(c)
95
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
0.6 0.6
0.4 0.4
0.2 0.2
u/U
v/U
0 0
-0.2 -0.2
-0.4 -0.4
-0.6 -0.6
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
y/D x/D
(a) (b)
0.8
0.6
0.4
0.2
P/(ρU 2)
-0.2
-0.4
-0.6
-0.8
0 0.2 0.4 0.6 0.8 1
y/D
(c)
96
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
Re = 10 Re = 100 Re = 1000
ISPH_DF 913 s unstable unstable
ISPH_DI 934 s 997 s 1,233 s
ISPH_DFDI 1,262 s 5,339 s 7,679 s
ISPH_DFS 965 s 1,067 s 1,395s
Table 4.2: Computing costs comparison for vortex spin-down with a resolu-
tion of 40 × 40. Physical time t = 2.0s. The same time step is used for all the
methods under the same Reynolds situation.
In Table 4.2, the computing costs of all four methods are listed. Increas-
97
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
-1
10
U max/U
10-2
-3
10
-4
10
0 0.5 1 1.5
t/T
98
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
-1
10
U max/U
-2
10
-3
10
0 2 4 6 8 10
t/T
99
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
100
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
where, φex,max and φISP H,max are respectively the exact and simulated max-
imum values of a general variable, φ. To better investigate the error, an
estimated relative error εL2 , based on relative L2 norm, is also calculated,
given by
np 21
X
(φex,i − φISP H,i)2 vi
i=1
εL2 =
np
,
(4.6)
X
2
φex,i vi
i=1
101
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
error is
np 12
X
(φex,i − φISP H,i )2 vi
i=1
eL2 =
np
,
(4.7)
X
vi
i=1
The viscous and Laplacian operators are divided into two groups here, the
one presented in [77] called lower-order operators here and the one proposed
by Schwaiger [99], named higher-order operators. Fig. 4.15 and 4.16 show
time variations of estimated error with the viscous and Laplacian operators
in [77] for Re = 100 and Re = 1000 respectively. Based on the relative L2
norm, seeing Eq. 4.6, two maximum relative error values, 0.045% and 2.5%,
for the finest meshes are obtained respectively at Re = 100 and Re = 1000
102
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
is also estimated for both higher-order [99] and lower-order [77] operators,
shown in Fig. 4.17. An estimated absolute L2-norm value, eL2 , is calculated,
at t = 2.0s. The estimated spatial convergence speed is close to first order for
ISPH_DFS with both two groups of viscous and Laplacian operators, with
By using Eq. 2.21 in §2.1.2, the accuracy of viscous and Laplacian operators
is improved [99]. The convergence speed may also be possibly improved.
However, with only first-order Taylor expansion in the correction stage, Eq.
2.49, the convergence speed is limited. To give an order consistent with or
higher than the improved operators, a simple expansion for the second-order
gradient operator, Eq. 2.15 with the normalised kernel Eq. 2.16, is applied
here as
!
X X X
∇(∇φ) ' Vj Vk (φk − φj )∇Wjk − Vl (φl − φi )∇Wij ∇Wij , (4.8)
j k l
103
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
0.03
0.02
ε
0.01
0
0 1 2 3 4 5
time/T
(a)
0.03
0.02
εL2
0.01
0
0 1 2 3 4 5
time/T
(b)
Figure 4.15: Estimated relative error profiles for the horizontal velocity com-
ponent, u, in Taylor-Green vortices with different resolutions, Re = 100. (a)
estimated error based on Eq.4.5; (b) estimated error based on L2 norm. The
time is normalised by T , and T = U/D.
— = 40 × 40; - - - = 80 × 80; · · · = 120 × 120; -·- = 160 × 160; – – = 200 × 200.
104
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
0.15
0.1
ε
0.05
0 5 10 15
time/T
(a)
0.25
0.2
0.15
εL2
0.1
0.05
0
0 5 10 15
time/T
(b)
Figure 4.16: Estimated relative error profiles for the horizontal velocity com-
ponent, u, in Taylor-Green vortices with different resolutions, Re = 1000. (a)
estimated error based on Eq.4.5; (b) estimated error based on L2 norm. The
time is normalised by T , and T = U/D.
— = 40 × 40; - - - = 80 × 80; · · · = 120 × 120; -·- = 160 × 160; – – = 200 × 200; -··-
= 240 × 240.
105
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
0.01
eL2
0.001
0.01
h/D
106
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
where, ∇Wij and ∇Wjk are normalised kernels. By considering the whole
1
φi0 = φi + (∇φ)i · δrii0 + δrii0 ∇(∇φ)δrii0 + O(δrii3 0 ). (4.9)
2
In Fig.4.18, the L2-norm based errors are plotted for different correction sit-
uations, and both corrections, including the first derivative terms in Taylor
expansion and including the second derivative terms in Taylor expansion,
gave almost the same prediction, more accurate than the simulation with-
out correction stage. Therefore, more terms in the correction stage do not
increase the accuracy, producing more or less similar predictions. The error
in the simulation without correction is still bounded. The reason why the
bounded error is obtained even without the Taylor-expansion correction is
possibly that the shifting of particle positions is more like “a random walk”;
that is, an extra diffusion is introduced, which makes the algorithm still sta-
ble even without correction.
107
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
0.0001
L2 error
1e-08
0 5 10
time (s)
predictions for these flows. The numerical noise can be extremely high.
108
Chapter 4. ISPH Method Comparisons and Accuracy Study for ISPH_DFS
expansion is necessary.
109
Chapter 5
I
N
two benchmark test cases. The lid-driven cavity and bluff body, nor-
mally squares or circular cylinders, are often used as benchmark test cases
to assess new codes and algorithms. Here, the lid-driven cavity and the
flow around a circular cylinder, at various Reynolds numbers, are simulated.
In this Chapter, except where stated otherwise the ISPH method used is
110
Chapter 5. Validation by Benchmark Test Cases
The lid-driven cavity problem has long been used as a test or validation case
for new codes or new solution methods. A good set of data for comparison is
[32], where the data are listed out in a table for different Reynolds number
situations. The two-dimensional geometry is shown in Fig.5.1. Although the
lid-driven cavity case is often simulated as a steady case, in simulations,
the fluid is accelerated by the lid at the top of the cavity to a steady state.
Simulations with fully-developed flows are compared with numerical data in
[32] and STAR_CD results.
UD
Re = . (5.1)
ν
U and D are set to 1 m/s and 1 m for the convenience in the normalization
and the calculation of Reynolds number. Three simulations with different
Reynolds numbers, Re = 100, Re = 400 and Re = 1000, are simulated.
To obtain convergence towards the data set in [32], high resolutions are used
in high Reynolds number situations. Tab.5.1 shows different resolutions
111
Chapter 5. Validation by Benchmark Test Cases
Figure 5.1: Geometry of lid driven cavity case. The graph shows the stream-
line in the case with Re = 400, simulated by ISPH_DFS, with a resolution of
61 × 61.
112
Chapter 5. Validation by Benchmark Test Cases
In [32], there is no reference data for pressure field. Therefore, the lid-driven
cavity cases are simulated by STAR-CD as well to provide comparison ob-
jects. The mesh convergence tests are considered for finite volume runs. The
convergence criterion is set as 10−5 . The PISO algorithm is used for pres-
sure and velocity coupling. Second-order differencing schemes are used for
convection terms in the calculation.
As a widely-studied test case, the lid driven cavity is used to further vali-
In Fig.5.2, 5.3 and 5.4, the velocity profiles in the middle of the domain,
x/D = 0.5 m and y/D = 0.5 m, are presented for each Reynolds number,
Fig.5.5, 5.6 and 5.7 present pressure profiles along the diagonal line, from
(0.0 m, 0.0 m) to (1.0 m, 1.0 m). With a resolution of 160 × 160, STAR-CD
can achieve almost the same accuracy for velocity as [32], and correspond-
ing pressure results are used for validation. ISPH_DFS produces almost
113
Chapter 5. Validation by Benchmark Test Cases
0.4
1
0.8
0.2
0.6
u/U
v/U
0.4 0
0.2
-0.2
0
-0.2
-0.4
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
y/D x/D
(a) (b)
1.1
0.9
0.2
0.8
0.7
0.6 0
y/D
v/U
0.5
0.4
-0.2
0.3
0.2
0.1
-0.4
0
-0.1
-0.4 -0.2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
u/U x/D
(a) (b)
114
Chapter 5. Validation by Benchmark Test Cases
0.4
1
0.8 0.2
0.6 0
y/D
v/U
0.4
-0.2
0.2
-0.4
-0.6
0 0.5 1 0 0.2 0.4 0.6 0.8 1
u/U x/D
(a) (b)
identical results for Re = 100 for all resolutions. For the higher Reynolds
numbers, results can be seen to converge with increasing resolution, with a
close agreement with 240 × 240 resolution.
With the same resolution, ISPH_DFS predicts flows less accurately than the
finite volume method, which is probably caused by the first-order accuracy
limitation of ISPH_DFS, while the finite volume code STAR-CD can provide
115
Chapter 5. Validation by Benchmark Test Cases
2.5
1.5
P/(ρU 2)
0.5
0 0.5 1
X
Figure 5.5: Pressure profile along the square diagonal, from (0.0, 0.0) to (1.0,
1.0), in lid-driven cavity, Re = 100.
= STAR-CD with 160 × 160 resolution; – – – = ISPH_DFS with 41 × 41 reso-
lution; · · · = ISPH_DFS with 81×81 resolution; — = ISPH_DFS with 161×161
resolution.
116
Chapter 5. Validation by Benchmark Test Cases
1.2
0.8
0.6
P/(ρU )
2
0.4
0.2
-0.2
0 0.2 0.4 0.6 0.8 1
X
Figure 5.6: Pressure profile along the square diagonal, from (0.0, 0.0) to (1.0,
1.0), in lid-driven cavity, Re = 400.
= STAR-CD with 160 × 160 resolution; – – – = ISPH_DFS with 81 × 81 res-
olution; · · · = ISPH_DFS with 161 × 161 resolution; — = ISPH_DFS with
241 × 241 resolution.
117
Chapter 5. Validation by Benchmark Test Cases
0.8
0.6
P/(ρU 2)
0.4
0.2
-0.2
0 0.2 0.4 0.6 0.8 1
X
Figure 5.7: Pressure profile along the square diagonal, from (0.0, 0.0) to (1.0,
1.0), in lid-driven cavity, Re = 1000.
= STAR-CD with 160 × 160 resolution; – – – = ISPH_DFS with 81 × 81 res-
olution; · · · = ISPH_DFS with 161 × 161 resolution; — = ISPH_DFS with
241 × 241 resolution.
118
Chapter 5. Validation by Benchmark Test Cases
in a channel
4
Y/D
5 10 15
X/D
When the fluid passes over a body, separation will happen depending on the
body shape and the incident flow. Depending on the flow and fluid conditions,
separation may cause vortex shedding, which occurs as vortices are created
at the back of the body and become detached from each side periodically,
shown in Fig. 5.8. This phenomenon happens in many engineering problems,
such as buildings, bridges, offshore pipelines and even airfoils at high angles
of attack etc . If the frequency of the shedding vortex matches the resonance
frequency of the structure, the structure will resonate and the structure’s
movement can become self-sustaining, which may cause structure failure.
119
Chapter 5. Validation by Benchmark Test Cases
number is calculated as
UD
Re = , (5.2)
ν
Four different initial particle sizes, δr =0.025, 0.02, 0.01 and 0.005, were
used for cases with Re = 20, and three different initial particle sizes,δr =0.025,
0.02 and 0.015, were used for cases with Re = 100. Particle numbers in each
cases are listed in Table 5.2.
The geometry of the case is shown in Fig. 5.9. To allow the flow to become
fully developed, the channel is much longer for the case with Re = 100. The
geometry parameters in different Reynolds situations are listed in Table 5.3.
120
Chapter 5. Validation by Benchmark Test Cases
The prediction of the velocity and pressure field from ISPH method are com-
pared with code Saturne, which is a finite volume code developed in R & D
in Electricite de France (EDF) [1]. The same geometry size and boundary
conditions are used in F.V. simulations. Three different meshes with differ-
ent resolutions, 16858 cells, 67353 cells and 255711 cells, are used for the case
with Re = 20, while another three different meshes with different resolu-
tions, 15914 cells, 62094 cells and 251293 cells, are used for the other case
with Re = 100. One of the meshes with a resolution of 67353 in the case with
Re = 20 is shown in Fig. 5.10. In simulations by code Saturne, the SIMPLE
121
Chapter 5. Validation by Benchmark Test Cases
Re=20
Fig. 5.11 shows the velocity and pressure contour in the case with Re =
20 at the time when the flow is fully accelerated. The new ISPH approach
provides very smooth prediction for both velocity and pressure field. Fig.
5.12 and 5.13 present the velocity and pressure profiles, at Re = 20, along the
cross sections, x/D = 2.5, 5.0, 7.5, 10.0, 12.5 and 30.0. Both the results from
ISPH with an initial particle distance δr = 0.01 and code Saturne with a grid
number of 255711 are plotted. The vertical coordinates are normalised by
122
Chapter 5. Validation by Benchmark Test Cases
the cylinder diameter D, the velocity by the bulk/average velocity U and the
1
pressure by ρU 2 . To better illustrate the profiles change along the channel,
2
excluding the profile at x/D = 2.5, normalised variable values are shifted
with the cross-section coordinates x, for example, the normalised horizontal
velocity u/U is changed by u/U + x, for x = 1.0. It can be observed in Fig.
5.12 (a) that due to the blockage effect of the cylinder, the horizontal velocity
is greatly reduced in the middle of the channel, and at the upstream of the
cylinder the vertical-velocity profiles show both positive and negative peaks
at upper and bottom halves of the channel. The positions of the positive and
negative peaks will be swapped downstream of the cylinder due to the flow
recirculation. The blockage effect is reduced as the fluid moves along the
channel. Both code Saturne and ISPH give almost identical predictions for
the velocity and pressure field.
Other noteworthy variables to predict in this bluff body case are the lift and
drag forces. They are calculated as
Z Z
Fd = pnx − (τxx · nx + τxy · ny ) (5.3)
Γ Γ
Z Z
Fl = pny − (τyx · nx + τyy · ny ) (5.4)
Γ Γ
where nx and ny are two components of the normal direction at the cylinder
∂ux ∂uy
surface; τxy is the viscous stress, τxy = µ + . The forces are usually
∂y ∂x
normalised by the bulk/average velocity U as
Fd
Cd = (5.5)
0.5ρU 2 D
Fl
Cl = (5.6)
0.5ρU 2 D
123
Chapter 5. Validation by Benchmark Test Cases
U
0.18
0.167143
4 0.154286
0.141429
0.128571
0.115714
Y/d
0.102857
0.09
2 0.0771429
0.0642857
0.0514286
0.0385714
0.0257143
0 0.0128571
0
0 5 10
X/d
(a)
V
0.05
0.0428571
4 0.0357143
0.0285714
0.0214286
0.0142857
Y/d
0.00714286
0
2 -0.00714286
-0.0142857
-0.0214286
-0.0285714
-0.0357143
0 -0.0428571
-0.05
0 5 10
X/d
(b)
P
0.03
0.028
0.026
4 0.024
0.022
0.02
Y/d
0.018
0.016
2 0.014
0.012
0.01
0.008
0.006
0 0.004
0.002
0 5 10
X/d
(c)
Figure 5.11: Velocity and pressure field in the bluff body case, Re = 20, initial
particle size δr = 0.01, at the time when the fluid fully developed. (a) Contour
of the horizontal velocity U; (b) Contour of the vertical velocity V ; (c) Contour
of the Pressure P .
124
Chapter 5. Validation by Benchmark Test Cases
Saturne
5 ISPH
3
y/d
-1
0 2 4 6 8
u/U
(a)
Saturne
5 ISPH
3
y/d
-1
0 1 2 3 4 5 6 7
v/U
(b)
Figure 5.12: The velocity profiles along different cross sections, x/D = 2.5,
5.0, 7.5, 10.0, 12.5 and 30.0, in the channel, Re = 20. The vertical coordinates
y and velocity components u and v are normalised by the cylinder diameter
D and the bulk velocity U. To better present the profiles along the channel,
excluding the profile at x/D = 2.5, normalised velocity values u/U and v/U
are shifted with the cross-section coordinates x.
(a) Profile of the horizontal velocity u; (b) Profile of the vertical velocity v.
125
Chapter 5. Validation by Benchmark Test Cases
Saturne
5 ISPH
3
y/d
-1
0 1 2 3 4 5 6 7
2
p/(0.5 ρ U )
Figure 5.13: The pressure profiles along different cross sections, x/D = 2.5,
5.0, 7.5, 10.0, 12.5 and 30.0, in the channel, Re = 20. The vertical coordinates
1
y and pressure P are normalised by the cylinder diameter D and ρU 2 re-
2
spectively. To better present the profiles along the channel, excluding the
1
profile at x/D = 2.5, normalised pressure values p/ ρU 2 is shifted with the
2
cross-section coordinates x.
126
Chapter 5. Validation by Benchmark Test Cases
6
dx=0.01
dx=0.025
4
Cd/Cl
variance in the interpolation, the results of drag and lift force are not as
smooth as counterparts in the finite-volume simulations, shown in Fig. 5.14.
After the fluid is fully accelerated, the time-averaged values are calculated
for the coefficients. In the case with Re = 20, the values for Cd with different
resolutions and extrapolated values, based on Richardson extrapolation, are
plotted in Fig. 5.15. The extrapolated values, 4.640 for ISPH and 4.642 for
Saturne, are calculated.
Re=100
In the case with Re = 100, three different initial particle spacings, δr = 0.025,
0.02 and 0.015, are tested here. As different shedding frequencies predicted
127
Chapter 5. Validation by Benchmark Test Cases
4.7
4.6
4.5
Cd
ISPH
4.4 Richardson extrapolation
4.3
4.2
50 100 150 200
1/dx
(a) Drag coefficients change against increasing resolution 1/dx
4.66
4.63
Cd
4.62
4.61
4.6
4.59
4.58
10000 1e+05
n
(b) Drag coefficients change against increasing mesh number n
Figure 5.15: Drag coefficients Cd from both ISPH and code Saturne as a
function of resolution, Re = 20. dx in (a) is the initial particle size. An
extrapolated value has been calculated.
128
Chapter 5. Validation by Benchmark Test Cases
and different time steps used in both ISPH and FVM, it is difficult to com-
pare contour graphs at a certain time. In Fig. 5.16, the contours of velocity
and pressure fields from ISPH are plotted, and it is shown in the graphs
that the new ISPH approach provides smooth predictions for the velocity
and pressure fields. Fig. 5.17 and 5.18 show velocity and pressure profile
comparisons between FVM and ISPH. Due to the anisotropic and chang-
ing particle distribution in ISPH, calculating the time-averaged value is not
achieved here. Therefore, all the velocity and pressure values within several
shedding periods along the cross sections are plotted in the same graph as
time-averaged velocity profiles from code Saturne. It can be seen from Fig.
5.17 and 5.18 that ISPH method gives similar velocity and pressure predic-
tions to code Saturne. The maximum and minimum values for Cd and Cl
from code Saturne are also presented in Fig. 5.19. The extrapolated values
for maximum and minimum drag coefficients are 2.829 and 2.673, while the
lift coefficients are ±1.215. Although the averaged values are quite close to
the predictions by code Saturne it is difficult for ISPH to provide smooth pro-
files for the drag and lift coefficients as particles moving around the cylinder
causes inaccuracy and noise in the force calculations, as shown in Fig. 5.20.
A clear vortex shedding period, T ' 1.54s, is obtained from the finest reso-
lution, as presented in the zoomed part in Fig. 5.20. Comparing with the
counterpart from code Saturne, T = 1.68s, there is 8% relative difference.
fD
The Strouhal numbers, St = U
, are calculated for both methods based on
the bulk velocity, U, and cylinder diameter, D. They are listed in Tab. 5.4.
Because the maximum and minimum drag coefficients are difficult to obtain
from ISPH simulations, and it is difficult to run another simulation with a
higher resolution on a single CPU, the averaged drag coefficient, Cd = 2.34, is
129
Chapter 5. Validation by Benchmark Test Cases
Table 5.4: Strouhal numbers, drag and lift coefficients from ISPH and
FVM, Re = 100. FVM simulations are conducted with open source code,
Code_Saturne [1].
calculated for ISPH simulation; the maximum and minimum lift coefficients
from the ISPH simulation with δr = 0.015 are approximately ±0.87. The rel-
ative differences between results from ISPH and code Saturne are 15% for
the drag coefficient and 28.4% for lift coefficient.
It should be noted here that accurate predictions are obtained in [81] by us-
ing local refinement around the bluff body and fixing all particles or particles
around the bluff body with the help of an Arbitrary Lagrangian-Eulerian
(ALE) [41] type algorithm. The numerical error introduced by the irreg-
ular particle distributions, which appear in SPH simulations with moving
particles, is avoided in [81]. Also the local refinement around the cylinder
improves the accuracy of force predictions. However, without an ALE-type
ISPH algorithm, if particle positions are fixed and resolutions are refined
ever, wall boundary methods with smaller numerical errors are needed in
the future.
130
Chapter 5. Validation by Benchmark Test Cases
U
0.9
0.828571
0.757143
4 0.685714
0.614286
0.542857
0.471429
Y/d
0.4
2 0.328571
0.257143
0.185714
0.114286
0.0428571
-0.0285714
0
-0.1
0 5 10
X/d
(a)
V
0.4
0.342857
0.285714
4 0.228571
0.171429
0.114286
0.0571429
Y/d
0
2 -0.0571429
-0.114286
-0.171429
-0.228571
-0.285714
-0.342857
0
-0.4
0 5 10
X/d
(b)
P
0.65
0.607143
0.564286
4 0.521429
0.478571
0.435714
0.392857
Y/d
0.35
2 0.307143
0.264286
0.221429
0.178571
0.135714
0.0928571
0 0.05
0 5 10
X/d
(c)
Figure 5.16: Velocity and pressure field in the bluff body case, Re = 100,
initial particle size δr = 0.015, at time = 196.65 s. (a) Contour of the horizontal
velocity U; (b) Contour of the vertical velocity V ; (c) Contour of the pressure
P.
131
Chapter 5. Validation by Benchmark Test Cases
3
y/d
0
0 1 2 3 4 5 6 7
u/U
(a)
3
y/d
0
0 1 2 3 4 5 6
v/U
(b)
Figure 5.17: The time-averaged velocity profiles along different cross sec-
tions, x/D = 2.5, 5.0, 7.5, 10.0 and 12.5, in the channel, Re = 100. The vertical
coordinates y and velocity components u and v are normalised by the cylinder
diameter D and the bulk velocity U. To better present the profiles along the
channel, normalised velocity values u/U and v/U are shifted with the cross-
section coordinates x/0.5. Red dots: F.V. method; black dots: ISPH method.
(a) Profiles of the horizontal velocity u; (b) Profiles of the vertical velocity v.
132
Chapter 5. Validation by Benchmark Test Cases
0.8
0.6
y/D
0.4
0.2
0
0 5 10 15 20 25
2
P/(0.5ρU )
Figure 5.18: The pressure profiles, from left to right, along different cross
sections, x/D = 2.5, 5.0, 7.5, 10.0 and 12.5, in the channel, Re = 100. The ver-
tical coordinates y and pressure P are normalised by the cylinder diameter
1
D and ρU 2 respectively. To better present the profiles along the channel
2
profiles are shifted. Red dots: F.V. method; black dots: ISPH method.
133
Chapter 5. Validation by Benchmark Test Cases
2.85
2.8
Cd max
Cd max / Cd min
Cd min
2.75 Extraploted Cd max
Extraploted Cd min
2.7
2.65
2.6
0 50000 1e+05 1.5e+05 2e+05 2.5e+05 3e+05
Grid Resolution
(a) Maximum and minimum drag coefficients change against increas-
ing grid number
Cl max
Cl max / Cl min
Cl min
0 Extraploted Cl max
Extraploted Cl min
-1
-2
0 50000 1e+05 1.5e+05 2e+05 2.5e+05 3e+05
Grid Resolution
(b) Maximum and minimum lift coefficients change against increas-
ing mesh number n
Figure 5.19: Drag Cd and lift Cl coefficients from code Saturne with different
resolutions, Re = 100. An extrapolated value has been calculated.
134
Chapter 5. Validation by Benchmark Test Cases
1
Cl / Cd
0
-1
-2
Cd
-3 Cl
-4
0 25 50 75 100 125 150
time (s)
(a)
4
1
Cl / Cd
-1
-2
Cd
-3 Cl
-4
0 25 50 75 100 125 150
time (s)
(b)
4
1
Cl / Cd
-1
-2
Cd
-3 Cl
-4
0 25 50 75 100 125 150
time (s)
(c)
Figure 5.20: Drag Cd and lift Cl coefficients from ISPH with different resolu-
tions, Re = 100. (a) dx = 0.025; (b) dx = 0.02; (c) dx = 0.015.
135
Chapter 5. Validation by Benchmark Test Cases
Through the benchmark test cases, lid-driven cavity and bluff body, the al-
gorithm and the code have been validated. Both cases are simulated by both
ISPH and FVM, with the commercial software STAR-CD or code Saturne [1].
Three different Re situations, Re = 100, 400 and 1000, have been simulated
for the lid driven cavity. However, it has been shown that the ISPH method
provides a slower convergence rate, compared to second-order of FVM, as
shown by the velocity profile comparison in Fig. 5.4. The pressure predic-
tions from ISPH match the converged finite-volume results very well.
For the bluff body, cases with two different Reynolds numbers are simulated
here. The vortex shedding happened at Re = 100. The lift and drag coef-
ficients are calculated on both ISPH and FVM with code Saturne [1]. At
Re = 20, both methods give almost identical predictions. However, differ-
ences appear at higher Reynolds number, Re = 100. The ISPH prediction of
force coefficients is quite noisy with low resolutions. The vortex-shedding pe-
riod is only obtained from the highest resolution in ISPH. The time-averaged
drag coefficient is calculated from ISPH simulation with the highest resolu-
tion. Relative differences between ISPH and code Saturne, 8%, 15% and
28.4%, are obtained for vortex-shedding period, drag and lift coefficients re-
spectively. Wall boundary methods with less numerical errors are needed in
the future.
136
Chapter 6
ISPH_DFS Method in
Free-Surface Flow Simulations
as undertaken for internal flows. The following phenomena can occur in one
simulation or in different simulations separately:
137
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
3. wave propagation.
We use the analytical solution for impulsive plate motion with zero grav-
ity [86] which shows asymptotic-like behavior at the fluid-solid intersection.
The dam-break flow has received considerable attention but has not been
fully exploited. For the wet-bed case we use the highly accurate solution
for flow acceleration at zero time, with a singularity in the free surface at
the lower gradient discontinuity [102]. With an analytic initial condition
a highly accurate high-order boundary integral method for potential flow
[24] is available for small times showing jet-like mushroom behavior [102].
a piston-type paddle for several periods and surface profiles, velocities and
pressures are compared against accurate stream-function theory [94]. This
is undertaken for small and moderate waves for which linear wavemaker
theory is realistic. To the authors’ knowledge, this is the first time ISPH
simulations of free-surface flows have been validated for both pressure and
velocity fields.
In this Chapter, without special declaration, the ISPH method used is ISPH_DFS
method introduced in Chapter 2.
138
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
may have larger values. This misjudgment for the free-surface boundary will
introduce certain error in the results, although this appears to be insignifi-
cant in this work. To keep the free-surface prediction uncontaminated by the
artificial movement of particles, the shifting in §2.2.4 is not applied for the
free-surface particles.
also reported in [26]. Even with the improved Laplacian operator, Eq. 2.21,
Schwaiger has shown a large relative error to the analytical value at the free
139
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
surface, almost 70% [99]. This error strongly distorts the free surface, and
The same numerical test as that in [99] has been repeated here, where the
singularity in the denominator of the relative error and to keep the same
geometry as in [99], a uniform 50 × 50 particle distribution is used within a
unit square with a dimension 2.0 ≤ x/D ≤ 3.0 and 2.0 ≤ y/D ≤ 3.0. The
relative error, defined as
| φSP H − φexact |
ε= , (6.1)
φexact
where φSP H is the value obtained from SPH calculation, and φexact is the
analytical value, has been plotted in Fig. 6.1. Fig. 6.1 shows the same
results as in [99]: the relative errors for the LS operator on the free surface
are around 70%, while the relative errors for the LM operator on the free
surface reaches 4000%; at the row just below the free surface, the relative
error for LS operator rapidly reduces to around 4%, while LM operator still
gives a high relative error value, around 500%. In this work, the LS operator
is used for viscous and Laplacian terms. A similar test is also conducted here
for gradient operator introduced in §2.1.2 on the same particle distribution.
The test function is ϕ = y m , with m =1, 2 and 3. The relative error for
∂ϕ
along the line x/D = 2.5 is plotted in Fig. 6.2. The maximum relative
∂y
error, located at the free surface, is around 1% with the gradient operator for
140
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
2
10
101
0
10
Relative Error
-1
10
10-2
-3
10
10-4
-5
10
2.9 2.91 2.92 2.93 2.94 2.95 2.96 2.97 2.98 2.99 3
X/D
Figure 6.1: Relative error profiles at the right end along y/D = 2.5 for test
function ϕ = xm + y m with the Laplacian operator in [99] and that in [77],
where ϕ is the general function; m are equal to 2, 3 and 4; D is the domain
length, equal to 1. Profiles for the Laplacian operator in [99] are black, while
profiles for the Laplacian operator in [77] are red. The same symbols are
used for the same test functions ϕ. Circles : test function with m = 2; dia-
monds : test function with m = 3; squares : test function with m = 4.
the test function, and for function ϕ = y, the accuracy has reduced to the
machine precision.
Because the pressure value on the free surface is constant, conveniently zero,
the Poisson equation is not solved for the free-surface particles, and the er-
ror introduced by the Laplacian operator on the free surface thus thus does
141
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
-1
10
10-3
-5
10
Relative Error
-7
10
10-9
-11
10
10-13
-15
10
∂ϕ
Figure 6.2: Relative error profile for along x/D = 2.5 for test function ϕ =
∂y
y m with the gradient operator introduced in §2.1.2, where ϕ is the general
function; m are equal to 1, 2 and 3; D is the domain length, equal to 1. Circles
: test function with m = 1; diamonds : test function with m = 2; squares :
test function with m = 3.
142
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
numbers considered, the relative errors are no larger than 7 × 10−5 , and ex-
actly 0 on the free surface. However, the relative errors could be greater with
the non-uniform particle distribution. In further tests, particles above 0.8
are shifted in both x and y coordinates with small random numbers between
±0.5dx and between ±0.1dx , where dx is the mean particle spacing. Fig. 6.4
(a) and 6.4 (b) show the relative error values along the channel for the solu-
tion of the Poisson equation 4ϕ = 1.0. It can be observed that the maximum
error around the free surface reaches almost 17% in Fig. 6.4 (a) and 1.5%
in Fig. 6.4 (b). It is also reported in [2], [90] and [99] that the disordered
particle distribution will complicate the convergence behavior of operators.
This lack of convergence caused by truncated-kernel errors on the free sur-
face will result in noisy free-surface evolution which will be demonstrated
later. Therefore, suppressing this free-surface instability is desirable but of
143
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
∆ϕ = 1
6e-05 ∆ϕ = y
2
∆ϕ = y
3
5e-05 ∆ϕ = y
Relative Error
4e-05
3e-05
2e-05
1e-05
0
0 0.2 0.4 0.6 0.8 1
y/D
Figure 6.3: Relative errors in the solution for 4ϕ = y m in a 1-D periodic open
dϕ
channel with boundary conditions ϕ = 1.0 at y = 1.0 and = 10.0 at y = 0.0
dy
with a uniform particle distribution, where m = 0, 1, 2, 3 in the tests.
It has been observed in simulations that the noisy free surface can be smoothed
by artificially increasing viscosity, which is shown in Fig.6.5 as below, where
the smooth free surface is predicted in the mud-flow simulation with higher
viscosity, 10−3 m2 /s. Here we are trying to increase the viscosity of the free-
surface particles and particles adjacent to the free surface, within a distance
of twice of the initial particle spacing, in order to obtain the same effect
as the high viscosity used in mud-flow simulations without strongly influ-
encing the free-surface predictions. The global maximum Peclet number
umax dx
P emax = , with umax as the maximum global velocity and dx as the
νd
initial particle spacing, is introduced here to calculated the higher viscos-
umax dx
ity νd around the free surface; that is, νd = . the The effect on the
P emax
144
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
0.2
0.15
Relative Error
0.1
0.05
0
0 0.5 1
y/D
(a)
0.02
0.015
Relative Error
0.01
0.005
0
0 0.5 1
y/D
(b)
145
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
0.8
y/a
0.6
0.4
0.2
Figure 6.5: High viscosity reduces noise in the collapse of a mud column with
a dimension of a × 2a, where a is the dimension unit.
tank with a dimension of 4a × 4a. The pressure field and free surface are
compared for cases with and without free-surface damping. The geometry of
the case is shown in Fig. 6.6. Fig. 6.7 shows the free surfaces at the same
time with and without free-surface damping. It can be observed that with
damping the predicted free surface is smooth compared with the simulation
without viscosity damping. It should be recognized that higher viscosity not
only limits the unbounded numerical error to a small value, but also damps
small-scale physical motion of the fluid. However, by choosing the damping
146
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
water column shown above. The maximum fluid velocity at the free sur-
√
face, umax , is estimated as 4ga [72, 63]. The maximum Peclet number,
umax dx
P emax = , is calculated through the free surface evolution for each
ν
case. If the surface viscosity is too high, the artificial diffusion will influence
the free-surface prediction, while values which are too low cannot provide
enough artificial damping for the truncated-kernel error. It has been shown
through numerical experiments that the viscosity value should be limited
to satisfy the maximum local Peclet number approximately in the range,
7 ≤ P emax ≤ 150 shown in Fig. 6.8, and a value around 30.0 is used in this
work. Through the P emax limit the estimated maximum free-surface velocity
and the particle size, the damping viscosity can be calculated for each case.
147
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
2.5
1.5
Y/a
0.5
0
0 1 2 3 4
X/a
(a)
2.5
1.5
Y/a
0.5
0
0 1 2 3 4
X/a
(b)
2.5
1.5
Y/a
0.5
0
0 1 2 3 4
X/a
(c)
148
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
To reduce the influence from free-surface damping, the viscosity value for
particles adjacent to the free surface has been reduced to half or a quarter of
free-surface viscosity.
It may be noted that a variable damping viscosity can also be calculated from
the local velocity, particle size and a specified local Peclet number. The same
dam break case is also simulated by specifying the local Peclet number of 10.
The same effect as constant damping viscosity has been observed, shown in
Fig. 6.9. In the results presented below a constant value has been used for
simplicity.
In the numerical experiment, the water column height at the left side, H,
and the water front edge, X, are recorded, comparing with experimental
results from [57] and WCSPH results from [116] with the k − turbulence
model. In Fig. 6.10, it is shown that both ISPH and WCSPH obtained quite
similar results, both overpredict the water front although the column height
predictions match the experiment very well. The free-surface damping did
not strongly influence the predictions of the water column height and the
water front edge. Lots of reasons may cause the overprediction of the water
front. It should be pointed out that both ISPH and WCSPH suffer from
the truncated kernel error, which may be an explanation for the difference
between simulations and experiments.
The similar numerical scheme was also applied in [89, 6], where the arti-
ficial dissipation term was used to diffuse discontinuities with the smooth-
ing scale similar to that resolved by the numerical method. In the previous
ISPH free-surface simulations in [63], instabilities were not apparent but the
k − ε turbulence model was used which effectively increases viscosity above
149
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
P P
1000 1000
916.667 916.667
3 833.333 3 833.333
750 750
666.667 666.667
583.333 583.333
500 500
2 2
Y/a
Y/a
416.667 416.667
333.333 333.333
250 250
166.667 166.667
1 83.3333
1 83.3333
0 0
0 0
0 1 2 3 4 0 1 2 3 4
X/a X/a
(a) (b)
P P
1000 1000
916.667 916.667
3 833.333 3 833.333
750 750
666.667 666.667
583.333 583.333
500 500
2 2
Y/a
Y/a
416.667 416.667
333.333 333.333
250 250
166.667 166.667
1 83.3333
1 83.3333
0 0
0 0
0 1 2 3 4 0 1 2 3 4
X/a X/a
(c) (d)
P P
1000 1000
916.667 916.667
3 833.333 3 833.333
750 750
666.667 666.667
583.333 583.333
500 500
2 2
Y/a
Y/a
416.667 416.667
333.333 333.333
250 250
166.667 166.667
1 83.3333
1 83.3333
0 0
0 0
0 1 2 3 4 0 1 2 3 4
X/a X/a
(e) (f)
Figure 6.8: Graphs of the free-surface profiles and pressure contours in the
collapse of water column, t ' 0.6s. (a) P emax = 7.0, dx/a = 0.04; (b) P emax =
150.0, dx/a = 0.04; (c) P emax = 7.0, dx/a = 0.025; (d) P emax = 150.0, dx/a =
0.025; (e) P emax = 7.0, dx/a = 0.0125; (f) P emax = 150.0, dx/a = 0.0125.
150
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
P
1000
916.667
3 833.333
750
666.667
583.333
500
2
Y/a
416.667
333.333
250
166.667
1 83.3333
0
0
0 1 2 3 4
X/a
Figure 6.9: Graphs of the free surface and pressure contour in the collapse
of water column, t ' 0.6s, dx/a = 0.0125. The damping viscosity is calculated
from the local maximum Peclet number limit, which is 10 for free-surface
particles, and 20 for particles close to the free surface.
1 4
Expt. Expt.
ISPH ISPH
0.8 WCSPH WCSPH
3
0.6
H/2a
X/a
0.4
2
0.2
0 1
0 1 2 3 4 5 0 1 2 3
1/2 1/2
t(2g/a) t(2g/a)
151
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
its molecular value. For almost inviscid fluids the effect on the accuracy of
An analytical solution for the free surface was derived by Peregrine [86] to
first order in time, assuming zero gravity. The free-surface elevation, η, is
written as
2Ut πx
η=− ln tanh (6.2)
π 4h
This transient free-surface flow has been further generalized by Roberts [95]
The free surfaces around the paddle, with two different paddle velocities,
U = 0.2m/s and U = 1.0m/s, are simulated for a water depth of 0.5m. Two
152
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
movement of the paddle, the number of solid particles on the bed is actually
reduced according to the change of the bed length.
Fig. 6.12 shows the effect of P emax and Fig. 6.13 and 6.14 show the free-
surface development with two different constant paddle velocities, U = 0.2m/s
and U = 1.0m/s for different resolutions and P emax = 30.0. With smaller
P emax limit, the free-surface predictions are less smooth than those with
larger P emax limit, shown in Fig. 6.12. Generally, the ISPH predicts the free
surface well with the limit P emax = 30.0. At the paddle, a jet is generated,
but the free surface becomes scattered at the crest. At the crest, only very
few particles are present in the SPH interpolation, therefore the truncated-
kernel error will be higher here than elsewhere on the free surface. What
is more, the value of ∇ · r cannot accurately identify the free-surface parti-
cles, as shown in Fig. 6.13 and Fig. 6.14. Some particles inside of the free
surface may also be taken as free-surface ones. The error caused by this
153
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
Pemax = 150
1.8 Pemax = 7
Analytical results
1.6
y/d
1.4
1.2
0.8
0 0.5 1 1.5 2
x/d
Figure 6.12: The free-surface simulations due to the impulsive paddle mo-
tion with different P emax limits, with a resolution of dx = 0.0125, a constant
velocity U = 1.0m/s, at time 0.04s, 0.08s, 0.12s, and 0.16s, from left to right
respectively.
mis-identification is also another reason for the scattered free surface at the
crest of the jet.
In this simulation case, a higher resolution is also used in both cases shown
in Figs. 6.13 and 6.14. More accurate predictions are obtained through the
higher resolution simulations.
The dam break is a widely used test case for impulsively started, rapid evolv-
ing free-surface flows. Experiments have been undertaken but results are
affected by the withdrawal of the plate which has a similar time scale to
the formation of the initial surface structures [102, 49]. For the dam break
154
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
1.8
1.6
y/a
1.4
1.2
0.8
0 0.5 1 1.5 2
x/a
Figure 6.13: The free-surface simulations due to the impulsive paddle mo-
tion, with a constant velocity U = 0.2m/s, P emax = 30.0, at time 0.2s, 0.4s, 0.6s
and 0.8s, from left to right respectively. — = the analytical results [86]; 4 =
the ISPH simulation with a resolution of dx = 0.0125; ◦= the ISPH simulation
with a resolution of dx = 0.00625.
1.8
1.6
y/a
1.4
1.2
0.8
0 0.5 1 1.5 2
x/a
Figure 6.14: The free-surface simulations around the impulsive paddle, with
a constant velocity U = 1.0m/s, P emax = 30.0, at time 0.04s, 0.08s, 0.12s, and
0.16s, from left to right respectively. — = the analytical results [86]; 4 = the
ISPH simulation with a resolution of dx = 0.0125; ◦= the ISPH simulation
with a resolution of dx = 0.00625.
155
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
Figure 6.15: Geometry of dam break with wet bed. The board separating two
water columns is removed or dissolved instantaneously.
with a wet bed, Stansby et al. [102] showed that the pressure distribution at
t = 0.0 s and the resulting flow acceleration may be determined accurately by
the solution of Laplace’s equation for pressure. They also computed time evo-
lution of the highly distorted free surface for small time assuming potential
flow and solving by a highly accurate boundary integral method [24]. These
solutions for t = 0.0 s and small times provide ideal tests for the ISPH algo-
rithm developed here. The initial geometry is shown in Fig. 6.15. The gate
separating the two columns is dissolved or removed immediately at t = 0.0 s.
The fluid density ρ and kinematic viscosity ν are 1000kg/m3 and 1.0×10−6 m2 /s.
In the investigation of the initial stage of dam break with the wet bed, the
geometry parameters, h1 , h2 , x1 and x2 , are set as 1.0m, 0.1m, 2m and 2m
respectively. For t = 0.0 s, the Laplacian equation for pressure is solved in
finite difference form on a uniform mesh with different mesh sizes of 0.025m,
0.02m, 0.01m, 0.005m and 0.0025m, with the 2nd-order differencing. The lin-
156
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
same resolutions are used for ISPH simulations. The horizontal accelera-
tion component ax and the angle, θ, between the acceleration vector and the
horizontal, from both methods are compared. The accelerations along the
free-surface are plotted with results from both methods.
For the free-surface evolution within a small time, the free-surface profiles
at different time are digitised from [102], and used as validation of ISPH
predictions.
For the initial condition (t = 0 s), Fig. 6.16 presents the pressure contours
from both the Laplacian solution with the second-order finite difference method
[102] and ISPH with a resolution of 0.01m. The results are almost identical to
each other. To compare further the results with each other, two cross-section
pressure profiles, at y = 0.0 and y = 0.5, are plotted, as shown in Fig. 6.17.
With a regular particle distribution, similar to a uniform structured mesh in
finite difference or finite volume methods, the ISPH and Laplacian solvers
provide an almost identical prediction.
Fig. 6.18 presents the acceleration profiles from both the ISPH and Lapla-
cian solvers along the free surface, with Fig. 6.18 (a) showing horizontal
acceleration values ax at the free surface along x = 2.0m, and Fig. 6.18 (b)
showing vertical acceleration values ay at the free surface along y = 0.1m.
The ISPH and Laplacian solvers provide almost identical predictions. ISPH
157
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
1.5
1
500
1000
10 00
1500
50
0
2000
15
2500 00
3000
3000
20
00
Y/h1
350
35000
1 50
0.5
0
40
25
00
00
50 0
0
6000 5 50
500
055 50
1000 10 200000
600 00 00
40
6500 0
00 4540
3 00
7000
7500 65
3000
00
500
75 0 7 0 00 00
0
25
0
-0.5
0 0.5 1 1.5 2 2.5 3 3.5 4
X/h1
10000
4200
3000
6000
P (pa)
2400
P (pa)
1200
2000
600
0 0
0 1 2 3 4 0 0.4 0.8 1.2 1.6 2
x (m) x (m)
(a) (b)
Figure 6.17: Pressure profiles at two different cross sections, (a) y = 0.0 m
and (b) y = 0.5 m, at t = 0.0s.
158
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
35 22
20
30
18
16
25
14
12
ax (m/s2)
ay (m/s2)
20
10
8
15
6
10 4
5 0
-2
0
0 0.2 0.4 0.6 0.8 1 2 2.2 2.4 2.6 2.8 3
y (m) x (m)
(a) (b)
Figure 6.18: The fluid acceleration comparison along the free surface at x =
2.0m and y = 0.1m with the same initial mesh/particle size, dx = dy = 0.025m.
Solid line with square symbols : ISPH; Dashed line : finite-difference Lapla-
cian solver. (a) Profile of horizontal acceleration values ax at the free surface
x = 2.0m; (b) Profile of vertical acceleration values ay at the free surface
y = 0.1m.
70 42
65
40
60
55 38
ax (m/s2)
50
θ( )
o
36
45
40 34
35
32
30
25 30
-100 0 100 200 300 400 500 -100 0 100 200 300 400 500
-1 -1
1/dx (m ) 1/dx (m )
(a) (b)
159
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
6.19 (a), the horizontal accelerations from both methods are increasing with
the increase of the resolution, and in Fig. 6.19 (b), the angle θ, between ac-
celeration vector and the horizontal direction are converging towards about
45o . The accelerations from ISPH are slighter lower but identical predictions
at a singularity by different schemes are not expected, although they are
For small time evolution, shown in Fig. 6.20 and 6.21, the predictions of the
free surface from ISPH with P emax = 7.0, P emax = 150.0 are compared with
digitised data from [102]. It is shown in Fig. 6.20 that ISPH with P emax = 7.0
damps some feature of the small-scale motion on the free surface, which
is indicated by the loss of accurate mushroom-jet prediction, while in Fig.
6.21 ISPH with P emax = 150.0 provides small noise distorting the capture
of the small-scale motion. By increasing the P emax limit to 30.0, the ISPH
prediction generally matches the results based on the nonlinear potential-
flow theory very well. However, at t = 0.08s in Fig. 6.22, the cap of the
mushroom-shape jet slightly deviates from the prediction in [102], which
may be caused by the truncated-kernel error along the free surface. It should
be noted that in [102] a tanh profile was used to define the initial condition
to avoid sharp corners as required for the boundary integral method. This
only changed the SPH distribution by one or two particles and had negligible
160
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
0.8
0.6
Y/h1
0.4
0.2
Figure 6.20: Free surface comparison at the initial stage of dam break with
wet bed. From the left to the right, the solid lines and groups of symbols
show the free surface at t = 0.024s, t = 0.04s, t = 0.066s and t = 0.08s. —
: digitised free-surface elevations from [102]; ◦ : free-surface particles from
the ISPH simulation with P emax = 7.
0.8
0.6
Y/h1
0.4
0.2
Figure 6.21: Free surface comparison at the initial stage of dam break with
wet bed. From the left to the right, the solid lines and groups of symbols
show the free surface at t = 0.024s, t = 0.04s, t = 0.066s and t = 0.08s. —
: digitised free-surface elevations from [102]; ◦ : free-surface particles from
the ISPH simulation with P emax = 150.
161
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
0.8
0.6
Y/h1
0.4
0.2
Figure 6.22: Free surface comparison at the initial stage of dam break with
wet bed. From the left to the right, the four solid lines and groups of symbols
show the free surface at t = 0.024s, t = 0.04s, t = 0.066s and t = 0.08s. — :
digitalised free-surface elevations from [102]; ◦ : free-surface particles from
the ISPH simulation with P emax = 30.
the channel
In this work, regular waves with wave heights, H = 0.05m and H = 0.1m,
in water of depth d = 0.5m are simulated with the ISPH method. Waves
are generated by the piston-type paddle in a rectangular wave tank with a
length of 18m. Based on linearized wavemaker theory [23], the paddle motion
follows
S
u= ωcos(ωt). (6.3)
2
162
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
H 2 (cosh 2kd − 1)
= , (6.4)
S sinh 2kd + 2kd
k = 2π/l (6.5)
The angular wave frequency ω in Eq. 6.3 has the relation with the wavenum-
ber as
Even with small wave heights waves show slightly non-sinusoidal form and
comparisons are made with highly accurate stream function theory of Rie-
163
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
The wave conditions are listed in Table 6.1. The water depth d is 0.5 m, the
wave tank length L is 18m. The fluid density ρ is 1000kg/m3, the kinematic
viscosity ν 10−6m2 /s, the initial particle spacing δr is 0.005 m.
To absorb the wave reflection from the wall at the end of the open channel, an
U = U0 f (x) (6.7)
α is a coefficient, equal to 2.0; δx0 is the damping zone length, equal to 3.0m;
x0 is the damping zone starting point, equal to 15.0m, x0 = L − δx0 , and L is
164
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
0.8
0.6
f(x-x0)
0.4
0.2
0
0 0.2 0.4 0.6 0.8 1
dx/δx0
Figure 6.24: Profile of exponential damping function used for wave absorp-
tion. α = 2.0; δx0 = 3.0; x0 = 0.0; dx = x − x0 .
the channel length. A similar exponential scheme has been used by Larsen
and Dancy [61] for Boussinesq wave modelling. It is shown in Fig. 6.24 that
when the particles are getting close to the right end of the wave tank, the
Fig. 6.25 shows the wave propagation case with different damping viscosi-
ties. The wave with a wave height of 0.05m and a wave length of 1.5m is
165
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
as smooth as with P emax = 7.0. However, since it has been observed in earlier
simulations that using too small P emax number limit can introduce too much
numerical diffusion, P emax = 30.0 has again been chosen for the following
wave propagation cases.
a) 1.5
1
Y/(2d)
0.5
-0.5
0 1 2 3 4 5
X/(2d)
b) 1.5
1
Y/(2d)
0.5
-0.5
0 1 2 3 4 5
X/(2d)
Figure 6.25: Damping investigation for wave propagation. (a) P emax = 7; (b)
P emax = 150.
Fig. 6.26 and 6.28 show the wave propagation along the open channel for two
different wave heights, H = 0.05m and H = 0.1m, with snapshots of the free
surface at different times. It is observed from the graphs that the first wave
decayed soon after it is generated by the paddle. The second wave decayed
gradually during the propagation process. But the other waves propagate
along the channel without decaying. It can be observed that the damping
zone effectively reduces the reflection from the right end of the channel under
166
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
1.5 +
++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ ++
1 +++++++++ ++++ +++ + +++ + + + + + + + + + + ++ + + + + + + + + + + + + + + + + + + + + + + + + ++ +
+
-0.5 ++ + +++ ++++++ + + +++++ +++++++ + +++++ +++++ ++++++ +++++++++++++ ++++++++++++++
-1 +++ + ++++ ++ ++++++++++ ++ +++ ++++ ++++++ +++ ++ +++++ + +++ ++++ ++ +++++
+ ++ + ++ + + + ++ +
0 2 4 6 8 10 12 14 16
X/(2d)
Figure 6.26: Wave propagation along the channel, from the top to the bottom,
t = 1.70625s, t = 5.3625s, t = 9.01875s, t = 12.675s, t = 16.33125s, t = 19.9875s,
the wave height H = 0.05m, the wave length L = 1.5m, the water depth
d = 0.5m, P emax = 30.0. + : the free-surface particles in the simulation.
these situations, with very small reflection apparent in the last profiles in
Fig. 6.27 and 6.29 show comparisons of free-surface predictions from ISPH
and stream function theory (code SAWW) at t = 19.99s and t = 9.75s for
the two different wave heights, H = 0.05m and H = 0.1m, respectively. It
can be observed that for the larger wave height, waves are noticeably non-
sinusoidal but a generally good agreement between ISPH and SAWW is ob-
Fig. 6.30 and 6.31 present the total pressure, non-hydrostatic pressure and
velocity profiles predicted by both ISPH and SAWW, below the crest and
167
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
0.7
0.6
y/(2d)
+ + + + + + + + + +
+ + + +
0.5 + +
+ + +
+ + + + + + + + +
+
+ +
0.4
0 1 2 3 x/(2d) 4 5 6
+++
1 ++++++++ +++++ + + + + ++++++ +++++ +++++++++++++++++++++++++++++++++
+++ +++
0.5 ++++++++ +++++++++ ++++++ + + ++ ++++++++ +++++++++++++++++++++++++
y/(2d)
0 5 10 15
x/(2d)
Figure 6.28: Wave propagation along the channel, from the top to the bottom,
t = 0.78s, t = 3.12s, t = 4.68s, t = 6.24s, t = 7.8s, t = 9.36s, t = 12.48s, the
wave height H = 0.1m, the wave length L = 3.0m, the water depth d = 0.5m,
P emax = 30.0. + : the free-surface particles in the simulation.
168
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
0.8
0.6
++ +++ ++ ++ +++
y/(2d)
+ ++ + +
+ + + ++++ + +
+ + + + +
++ + + + + + + +
+ +
+++ +++
+ ++ + + ++ ++
+ + +++ ++
0.4
0.2
0 5 x/(2d) 10 15
trough, with positions given on the figure caption. The non-hydrostatic pres-
sure is given by
p = P − ρg(η − y), (6.9)
169
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
0.5 0.5
0.45 0.45
0.4 0.4
0.35 0.35
0.3 0.3
y (m)
y (m)
0.25 0.25
0.2
0.2
0.15
0.15
0.1
0.1
0.05
0.05
0
0
0 1000 2000 3000 4000 -0.1 0 0.1
P (Pa) u (m/s)
(a) (b)
0.55
0.5
0.45
0.4
0.35
0.3
y (m)
0.25
0.2
0.15
0.1
0.05
-0.05
-300 -200 -100 0 100 200 300
P nh (Pa)
(c)
Figure 6.30: Comparing velocity and pressure predictions from both ISPH
with P emax = 30.0 and SAWW [8].Wave height H = 0.05m; Wave length l =
1.5m; Water depth d = 0.5m; t = 19.9875s.
(a) Pressure profiles below wave crest and trough; (b) Horizontal velocity
profiles below wave crest and trough; (c) Non-hydrostatic pressure profiles
below wave crest and trough.
— : prediction from SAWW below the wave crest ; − − −: prediction from
SAWW below the wave trough ; ◦ : prediction from ISPH along below the
wave crest, x = 1.61m ; : prediction from ISPH line below the wave trough,
x = 2.36m.
170
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
0.7 0.7
0.6 0.6
0.5 0.5
0.4 0.4
y (m)
y (m)
0.3 0.3
0.2 0.2
0.1 0.1
0 0
-0.1 -0.1
-2000 0 2000 4000 6000 8000 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4
P (Pa) u (m/s)
(a) (b)
0.7
0.6
0.5
0.4
y (m)
0.3
0.2
0.1
-0.1
-400 -300 -200 -100 0 100 200
P nh (Pa)
(c)
Figure 6.31: Comparing velocity and pressure predictions from both ISPH
with P emax = 30.0 and SAWW [8].Wave height H = 0.1m; Wave length L =
3.0m; Water depth d = 0.5m; t = 9.75s
(a) Pressure profiles below wave crest and trough; (b) Horizontal velocity
profiles below wave crest and trough; (c) Non-hydrostatic pressure profiles
below wave crest and trough.
— : prediction from SAWW below the wave crest ; − − −: prediction from
SAWW below the wave trough ; ◦ : prediction from ISPH below the wave
crest, x = 5.943m; : prediction from ISPH below the wave trough, x =
7.443m.
171
Chapter 6. ISPH_DFS Method in Free-Surface Flow Simulations
ity at the surface and, with sufficient particle resolution and surface P emax
numbers below about 30, accurate flow predictions can be made. Surface
P emax numbers may also be reduced by increasing particle resolution without
changing viscosity but this would increase computational times considerably.
The algorithm with artificially increased surface viscosity is tested for the
impulsive flows of the impulsively started plate and the wet bed dam break
for which there are analytical and highly accurate solutions. It is also shown
that regular wave propagation can be accurately simulated without decay.
Surface profiles, total pressures and velocities are simulated accurately, al-
172
Chapter 7
7.1 Conclusions
I
N
methods, ISPH_DF [17, 63], ISPH_DI [100], ISPH_DFDI [44] have been
investigated. Success has been achieved in several different applications.
However, it has been discussed that due to the kernel flaw and the stream-
line movement of particles in the SPH method , ISPH_DF [17, 63] is not a
stable method; the simulations with ISPH_DI [100] are quite stable, but the
inaccurate prediction has dropped the advantage to use incompressible SPH;
ISPH_DFDI [44] can provide accurate predictions but with much higher com-
puting expense and doubled pressure and pressure gradient fields.
Based on the ISPH_DF [17, 63], the ISPH_DFS is a stable algorithm without
sacrificing efficiency. Particles are shifted slightly across streamlines, avoid-
ing particle spacing distortion and the error resulting from particle cluster-
ing, and their hydrodynamic properties are adjusted by Taylor series inter-
173
Chapter 7. Conclusions and Future Work
polation. This small change to ISPH_DF method has stabilized the ISPH
method with very little increase of computing expense. Through the test
cases, Taylor-Green vortices, vortex spindown, it has been shown that that
ISPH_DFS method is much more efficient than ISPH_DFDI method [44];
the accuracy is improved over ISPH_DI method [100]; the stability is highly
order [82], the influence of the viscosity and Laplacian operators is investi-
gated. It is shown in §5.1 that although the operator proposed by Schwaiger
[99] has improved the accuracy performance, the spatial convergence speed
is not enhanced. Similar convergence speeds have been obtained for both
types of the viscous and Laplacian operators [77, 99]. This is most possibly
caused by the particle convection method – particles move after the velocity
at the new time step is updated, but will have the same velocity at the new
place as at the old position. This delayed velocity updating may be one of
main reasons causing slow spatial convergence.
Two benchmark cases, lid-driven cavity and bluff body in the channel, were
simulated. The results are compared with Ghia’s data [32] and the finite vol-
ume simulations from the commercial software STAR-CD and code Saturne
[1]. It has been shown that for low Reynolds number, ISPH method gave
quite good prediction as finite volume method did, while the low convergence
174
Chapter 7. Conclusions and Future Work
and lift forces were calculated in the bluff-body cases. It has been shown
that the movement of the particles around the circular cylinder has resulted
in the noisy force values. This noise may disable or place difficulty on the
prediction of the vortex-shedding frequency under the situation of low reso-
lutions. Although the increase of the resolution can provide higher accuracy
without considering the computing expense, a better approach for the simu-
lation of fluid on the wall is needed, especially for the accurate simulation of
fluid-body interaction.
distributed particles, while this small error will be amplified as the parti-
cle distribution is distorted, and the free surface could be quite noisy with
ISPH method. One treatment has been introduced here by slightly increas-
ing the viscosity on the free surface. This has been proved to be an effective
way to reduce the noise on the free surface without strongly influencing the
free-surface prediction.
Three test cases, initial-stage investigation to the dam break, impulsive pad-
dle and wave propagation, have been studied. Good agreements are obtained
for the studied cases. The thorough quantitative validation has been con-
ducted for pressure and velocity predictions in the context of the SPH method
through the wave propagation test case. It is shown that for the free-surface
prediction, or for the velocity and pressure field, this ISPH method can pro-
vide accurate predictions.
175
Chapter 7. Conclusions and Future Work
In this work, all the simulated cases are with simple geometries. In this
way, the complexity caused by the geometries is avoided, and the boundary
complicated geometries.
At the end, the truncated-kernel error around free surface has been in-
vestigated in this work. A numerical approach to overcome the instability
caused by the truncated-kernel error on the free surface has been proposed.
However, it should admitted here that this approach is mainly a numerical
method. A better way is needed.
176
Appendix A
Here, only a brief instruction about code SPHysics will be presented to facili-
tate the further illustration about the incompressible SPH code, which is de-
veloped based on this WCSPH code SPHsysics. More details about SPHysics
can be obtained from the user guide [121].
Smoothing kernel
the distance between particle i and its neighbouring particle j increase. The
smoothing length h is proportional to the initial particle distance dr, nor-
mally h = c · dr, and decides the support size of the particle i, Li . For differ-
ent kernels in a manner of linear algebraic equations, Li = αh with different
177
Chapter A. The SPHysics Code
ω(q) = γe(−q )
2
(A.1)
where, γ is a coefficient, equal to 1/ (πh2 ) in 2D, and 1/ π 3/2 h3 in 3D.
3 2 3 3
q − q+ 0≤q≤2
ω (q) = γ 16 4 4 (A.2)
0
q≥2
3 2 3 3
1 − q + q 0≤q≤1
2 4
1
ω (q) = γ (2 − q)3 1≤q≤2 (A.3)
4
0
q≥2
1 − q (2q + 1) 0 ≤ q ≤ 2
4
ω (q) = γ 2 (A.4)
0
q≥2
178
Chapter A. The SPHysics Code
Momentum equation
In the momentum equation Eq. 2.2, one strong effect is viscosity. The viscous
effect can be simulated by the artificial viscosity [71] or by the approximated
laminar viscosity [77]. Excluding the artificial viscosity and the laminar
viscosity, the Sub-Particle Scale (SPS) approach, which was first introduced
by Gotoh et al. [34] in the context of Moving Particle Semi-implicit method
(MPS), to model turbulence was also developed in the code SPHysics.
The artificial viscosity was first proposed by Monaghan [71]. Due to its sim-
plicity, it is widely used in the SPH simulations. If we define the viscous
X
ψi = mj Πij ∇ωij (A.5)
j
−ζcij µij
uij · rij < 0
ρij
Πij = (A.6)
0
uij · rij > 0
and
huij · rij
µij = (A.7)
rij2 + η 2
where, rij = ri −rj ; uij = ui −uj ; cij is the averaged speed of sound for particle
ci +cj ρi +ρj
i and j, cij = 2
; ρij is the averaged density for particle i and j, ρij = 2
;
179
Chapter A. The SPHysics Code
SPS model, similar to Large Eddy Simulation, is used to consider the effect
of fluid scale smaller than the particle size and simulate the turbulence effect
in the fluid, seeing Gotoh et al. [34]. In 2004, Gotoh et al. [35] extended the
SPS model into incompressible SPH. To include the sub-particle effect in the
fluid, the momentum equation is changed into
Du 1 1
= − ∇P + ν∇2 u + F + ∇τ (A.8)
Dt ρ ρ
where F is the body force; τ is the SPS stress tensor. With the Boussinesq’s
hypothesis, SPS stress tensor is modeled as
2 2 2 2
τij = ρ 2νt Sij − kδij − CI ∆ δij | Sij | (A.9)
3 3
2
where, νt is the eddy viscosity, νt = [min (Cs 4l2 )] | S |, CI = 0.0066, Sij is the
stress tensor element, Cs is the Smagorinsky constant, with a range of 0.1 -
0.2, and a value of 0.12 is used in SPHysics code. 4l is the particle-particle
spacing, and it is used as the implicit filter here, equivalent to the grid size
in mesh-based method; | S |= (2Sij Sij )1/2 .
locity, is used for the calculation of coefficient. To satisfy the CFL condition,
the time step is limited to a very small value, which slows down the simula-
tion efficiency. Compressibility causes problems with sound wave reflection
180
Chapter A. The SPHysics Code
at the boundaries. However, due to its simplicity, the state equation is still
where, γ = 7, B = c20 ρ0 /γ, with ρ0 being the reference density and c0 being
r
∂P
the speed of sound at the reference density and c0 = c (ρ0 ) = .
∂ρ
With substituting the pressure gradient term and SPS stress gradient term
with Eq. 2.14, the viscous term with Eq. 2.19, and considering the SPS
model, the momentum equation, Eq. A.8 can be discretized as
du X Pi Pj X 4mj µrij · ∇ωij
' − mj 2
+ 2 ∇ωij + 2
u
2 ) ij
dt i j
ρi ρj j
(ρi + ρj ) (r ij + η
X τi τj
+ F+ mj 2
+ 2 ∇ωij (A.11)
j
ρi ρj
Continuity equation
dρ
+ ρ∇ · u = 0. (A.12)
dt
Therefore, the particle density can be updated with Eq. 2.15 as the formula-
tion
dρ X
= −ρ∇ · u = mj uij · ∇ωij . (A.13)
dt i j
The right hand side of Eq. A.13 is also briefly written as D in the later part.
181
Chapter A. The SPHysics Code
Also the particle density can be calculated from Eq. 2.6, which is
X
ρ= mj ωij . (A.14)
j
For the accuracy consideration, both of the expression suffer from the trun-
cated kernel error at the free surface. For the efficiency consideration, Eq.
A.13 needs the first derivative of the kernel function, but Eq. 2.41 the kernel
value, which means extra memory and computing cost will be put in because
Time marching
(also known as the modified Euler method) [72], Verlet and Beeman [3], are
used. The Predictor-Corrector and Verlet scheme are introduced here. Oth-
ers can be found in [31].
Predictor-Corrector
Predicting stage:
n+1/2 ∆t n
ui = uni + F (A.15)
2
n+1/2 ∆t n
ρi = ρni + D (A.16)
2
n+1/2 ∆t n
ri = rni + u (A.17)
2
182
Chapter A. The SPHysics Code
Correcting stage:
n+1/2 ∆t n+1/2
ui = uni + F (A.18)
2
n+1/2 ∆t n+1/2
ρi = ρni + D (A.19)
2
n+1/2 ∆t n+1/2
ri = rni + u (A.20)
2
n+1/2
un+1
i = 2ui − uni (A.21)
n+1/2 n+1/2
ρi = 2ρi − ρni (A.22)
n+1/2 n+1/2
ri = 2ri − rni (A.23)
Verlet Scheme
un+1
i = uin−1 + 2∆tFn ; (A.24)
ρn+1
i = ρin−1 + 2∆tD n ; (A.25)
∆t2 n
rn+1
i = rni + ∆tuni + F . (A.26)
2
un+1
i = uni + ∆tFn ; (A.27)
ρn+1
i = ρni + ∆tD n ; (A.28)
183
Chapter A. The SPHysics Code
∆t2 n
rn+1
i = rni + ∆tuni + F . (A.29)
2
184
Appendix B
Linear Solvers
CG is the most popular iterative solver. It is widely used to solve the Pressure
Poisson equation in incompressible simulations. The CG method is an algo-
r0 = B − AX0 (B.1)
185
Chapter B. Linear Solvers
P 0 = r0 (B.2)
Iterate the following conjugate-gradient root search process until the conver-
gence criterion is satisfied,
rTk rk
αk = (B.3)
PTk APk
Xk+1 = Xk + αk Pk (B.4)
rTk+1rk+1
βk = (B.7)
rTk rk
k =k+1 (B.9)
the code, is simply listed below. The same as that in CG solver, the linear
system is simply written as AX = B with an initial guessed value X0 and a
preconditioner K, and the residual is
r0 = B − AX0 (B.10)
186
Chapter B. Linear Solvers
r̄0 is an arbitrary vector, and the dot product (r0 , r̄0 ) satisfies
ρ0 = α = ω0 = 1; (B.12)
v0 = p0 = 0. (B.13)
Then the following conjugated-gradient root searching process with the resid-
ual smoothing can be conducted until the convergence of the linear system is
achieved.
ρi = (r̄0 , ri−1 ) ; β = (ρi /ρi−1 ) (α/ωi−1) ; (B.14)
vi = Ay; (B.17)
α = ρi / (r̄0 , vi ) ; (B.18)
t = Az; (B.21)
Xi = Xi−1 + αy + ωi z. (B.23)
187
Chapter B. Linear Solvers
ri = s − ωi t, (B.24)
and the root searching will restart again from Eq. B.14.
Although this linear solver demands more calculation during one iteration
than CG solver, its robustness and efficiency on the solution of linear systems
generated in ISPH algorithms, shown by the numerical experiment in §3.3,
have given itself more advantages than the CG solver.
Convergence criterion
The convergence of the linear solver is achieved when the iteration number
reaches the maximum iteration number, or
k r i k2
≤ε (B.25)
k r 0 k2
where k · k2 is the l2 -norm; subscript i and 0 are for the ith iteration and
the initial value respectively; ε is the linear solver tolerance, ε = 10−5. The
discussion about the convergence of linear system can be gained from [29].
188
Appendix C
Publication lists
During three years PhD, the author did not only focus on the meshless SPH
method, but also studied some other topics, which provides great help for the
research. The publications during the PhD time are listed as follow.
Journal Publications:
189
Chapter C. Publication lists
Conference Publications:
P.K. Stansby, A.C. Hunt, R. Xu, P.H. Taylor, A.G.L. Borthwick, T. Feng, D.R.
Laurence, Wave overtopping from focused wave groups: experiments and
190
Bibliography
191
BIBLIOGRAPHY
192
BIBLIOGRAPHY
1993.
[22] S. F. Davis. Flux difference splittings and limiters for the resolution of
193
BIBLIOGRAPHY
for steep unsteady water waves. In Methods for Fluid Dynamics II (ed.
K. W. Morton & M. J. Baines), pp. 671 - 679. Oxford University Press,
1986.
194
BIBLIOGRAPHY
[35] H. Gotoh, S. Shao, and T. Memita. SPH-LES model for numerical in-
vestigation of wave interaction with partially immersed breakwater.
Coastal Engineering Journal, 46 (1): 39 - 63, 2004.
6662, 2001.
195
BIBLIOGRAPHY
[40] C. Hirt and B. Nichols. Volume of fluid (VOF) method for the dynamics
[45] Y. Imaeda, T. Tsuribe and S.-I. Inutsuka. Gas accretion from the el-
liptic gas disk to the binary system. ERCOFTAC SIG SPHERIC IIIrd
International workshop, Lausanne, Switzland, 2008.
196
BIBLIOGRAPHY
126, 2008.
[48] H. Jasak and A. Gosman, Automatic resolution control for the finite-
volume method, part 1 : A-posteriori error estimates, Numerical Heat
ator Models for Free-Surface Fluid Flows. Ph.D. thesis, Kyoto Univer-
sity, 2008.
197
BIBLIOGRAPHY
[59] J. Kuhnert and S. Tiwari. Finite pointset method based on the projec-
tion method for simulations of the incompressible Navier-Stokes equa-
198
BIBLIOGRAPHY
1985.
199
BIBLIOGRAPHY
[79] D. Stevens, H. Power, M. Lees and H. Morvan. The use of PDE cen-
tres in the local RBF Hermitian method for 3D convective-diffusion
problems. Journal of Computational Physics, 228: 4606 - 4624, 2009.
200
BIBLIOGRAPHY
201
BIBLIOGRAPHY
1576–1581, 1971.
sis applied to viscous free surface fluid flow. International Journal for
Numerical Methods in Fluids, 7: 953 - 984, 1987.
1981.
202
BIBLIOGRAPHY
[100] S. Shao and E. Lo. Incompressible SPH method for simulating Newto-
nian and non-Newtonian flows with a free surface. Advances in Water
Resources, 26: 787 - 800, 2003.
[102] P. K. Stansby, A. Chegini and T.C. D. Barnes. The initial stages of dam-
break flow. Journal of Fluid Mechanics, 374: 407 - 424, 1998.
203
BIBLIOGRAPHY
[106] S. Tiwari and J. Kuhnert. Finite Pointset Method Based on the Pro-
[111] H. A. Van Der Vorst. Bi-CGSTAB: A fast and smoothly converging vari-
ant of Bi-CG for the solution of non-symmetric linear systems. S. J. Sci.
Stat. Comput., 13: 631 - 644, 1992.
204
BIBLIOGRAPHY
[115] D. Violeau and R. Issa. Modelling turbulent free surface flows with
pressible SPH (ISPH) Based on the Projection Method and A New Ap-
proach. Journal of Computational Physics, 228 (18): 6703-6725, 2009.
205
BIBLIOGRAPHY
http://wiki.manchester.ac.uk/sphysics/index.php/Main_Page. Lastest
modification: 28 January 2010.
206