GaoM PDF

NUMERICAL SIMULATION OF LIQUID SLOSHING IN
RECTANGULAR TANKS USING CONSISTENT PARTICLE

METHOD AND EXPERIMENTAL VERIFICATION

GAO MIMI
(B.ENG., SHANGHAI JIAO TONG UNIVERSITY, CHINA)

A THESIS SUBMITTED
FOR THE DEGREE OF DOCTOR OF PHILOSOPHY
DEPARTMENT OF CIVIL ENGINEERING
NATIONAL UNIVERSITY OF SINGAPORE
2011

i
Acknowledgement
The author would like to express her sincerest gratitude and appreciation to the
following people for their invaluable guidance, advice and encouragement,
Professor Koh Chan Ghee, for his professional guidance, wisdom, advice and
continual support throughout the duration of my PhD. There was a time when my
world was filled with darkness and I even thought of giving up. It was Professor
Kohs great support and encouragement that helped me go through the dark days and
accomplish this thesis. His guidance and help are greatly appreciated.
Professor T Balendra, for his invaluable advice and patient help in the study. I am
very grateful for his understanding and support in the whole PhD study.
Associate Professor Ang Kok Keng and Professor Choo Yoo Sang, for their
suggestions during the qualifying examination which helped me greatly to better
define the research focus.
Research Fellow, Dr. Luo Chao for many useful discussions and for his help in the
experiments. The thought provoking discussions with him have contributed to the
success of the numerical model.
All the staff in the structural engineering laboratory for their kind assistance in
providing technical and logistic support for the experimental work.
Dr. Zhang Zhen, Dr. Teng Mingqing and Dr. Duan Wenhui for their insightful
discussions and advice.
Finally, to my parents, sisters and brother for their unconditional encouragement and
support, without which this thesis would not have been completed successfully.

ii

iii
Table of Contents
Acknowledgement .......................................................................................................... i
Table of Contents ..........................................................................................................iii
Summary ...................................................................................................................... vii
List of Figures ............................................................................................................... ix
List of Tables ............................................................................................................. xvii
Nomenclature .............................................................................................................. xix
Chapter 1 Introduction ................................................................................................... 1
1.1 Overview .................................................................................................... 2
1.2 Sloshing in membrane LNG tank .............................................................. 2
1.3 Study of liquid sloshing ............................................................................. 4
1.4 Research scope and objectives ................................................................... 6
1.5 Organization of the thesis .......................................................................... 7
Chapter 2 Literature Review ........................................................................................ 11
2.1 Research works involving mainly experimental study ............................ 12
2.2 Analytical study of liquid sloshing .......................................................... 14
2.3 Numerical study of liquid sloshing .......................................................... 16
2.3.1 Mesh-based methods ............................................................................ 16
2.3.2 Meshless methods ................................................................................ 24
2.4 LNG and LNG sloshing ........................................................................... 28
2.4.1 LNG and its carrier system .................................................................. 28
2.4.2 Sloshing phenomena ............................................................................ 30
Chapter 3 Formulation of Consistent Particle Method ................................................ 37
3.1 Introduction .............................................................................................. 37
3.2 Moving particle semi-implicit method..................................................... 38

iv
3.2.1 Governing equations ............................................................................ 39
3.2.2 MPS formulation .................................................................................. 42
3.2.3 Modeling of incompressibility ............................................................. 44
3.2.4 Boundary conditions ............................................................................ 44
3.2.5 Drawbacks of MPS .............................................................................. 46
3.3 CPM based on Taylor series .................................................................... 47
3.3.1 Introduction .......................................................................................... 47
3.3.2 Approximation of gradient and Laplacian by Taylor series ................ 50
3.3.3 Main features of CPM .......................................................................... 54
3.3.4 Performance test of the Laplacian based on Taylor series ................... 63
3.4 Concluding remarks ................................................................................. 67
Chapter 4 Numerical Simulation of Incompressible Free Surface Flows by CPM ..... 83
4.1 Introduction .............................................................................................. 83
4.2 Benchmark examples ............................................................................... 84
4.2.1 Hydrostatic pressure in a static tank .................................................... 84
4.2.2 Dam break with 2 / =
w
L d .................................................................. 86
4.3 Parametric study of CPM ......................................................................... 87
4.3.1 Influence of weighting functions in weighted least-square solution ... 88
4.3.2 Influence of influence radius ............................................................... 90
4.3.3 Influence of particle sizes .................................................................... 92
4.3.4 Influence of time step........................................................................... 93
4.3.5 Computational cost .............................................................................. 94
4.4 Numerical simulation of free oscillation of liquid ................................... 95
4.5 Numerical simulation of violent fluid flows with breaking ..................... 97
4.5.1 Free oscillation of liquid in a container with large amplitude ............. 98

v
4.5.2 Dam break with 5 . 0 / =
w
L d ............................................................... 99
4.5.3 Dam break with obstacle .................................................................... 103
4.6 Concluding remarks ............................................................................... 105
Chapter 5 Liquid Sloshing in Rectangular Tanks: Experimental Study and CPM
Simulation .................................................................................................................. 131
5.1 Introduction ............................................................................................ 131
5.2 Experimental setup................................................................................. 132
5.2.1 Experimental facilities ....................................................................... 132
5.2.2 Water Tank......................................................................................... 133
5.2.3 Shake table ......................................................................................... 133
5.2.4 Wave probes....................................................................................... 133
5.2.5 Pressure sensor ................................................................................... 134
5.2.6 Displacement transducer .................................................................... 135
5.2.7 High speed camera ............................................................................. 135
5.2.8 Other considerations .......................................................................... 135
5.3 Sloshing experiments and comparison with CPM solutions.................. 136
5.3.1 Experiments of sloshing waves in high-filling tank .......................... 138
5.3.2 Experiments of sloshing waves in low-filling tank............................ 150
5.3.3 Experiments with sloshing wave impact on the tank ceiling ............. 153
5.4 Concluding remarks ............................................................................... 154
Chapter 6 Conclusions and Future Research ............................................................. 189
6.1 Conclusions ............................................................................................ 189
6.2 Future work ............................................................................................ 191
References .................................................................................................................. 193
Appendix A: CD for animation files and explanation note ....................................... 207

vi

vii
Summary
The use of numerical simulation has made an enormous impact on the study of free
surface motion of incompressible liquid such as liquid sloshing. Simulating this
complex problem has many important applications, ranging from coastal protection
and offshore structure design to LNG/oil sloshing on vessels. Furthermore, animated
wave motion has great potential in modern movies and computer games where violent
liquid motion is featured.
In this context, conventional mesh-based numerical methods have met difficulties in
simulating waves involving discontinuity of liquid motion (e.g. wave breaking). Even
with some free-surface capturing techniques incorporated, such as marker-and-cell
and volume of fluid, mesh-based methods suffer from the problem of numerical
diffusion. This is mainly due to the discretization of advection terms in the Navier-
Stokes equation in Eulerian formulation. In addition, tracking of free surface requires
complex and time consuming algorithm to update the time varying nonlinear
boundary.
In recent years, a new generation of computational methods known as meshless
(mesh-free) methods has been shown to outperform conventional mesh-based method
in dealing with discontinuous fluid motion. Lagrangian meshless methods called
particle methods have shown very good potential in dealing with large-amplitude free
surface flows, moving interfaces and deformable boundaries. The problem of
numerical diffusion does not arise in particle methods. Nevertheless, in many of the
existing particle methods such as Smoothed Particle Hydrodynamics (SPH) method
and Moving Particle Semi-Implicit (MPS) method, the approximation of partial
differential operators requires a pre-defined kernel function. Accuracy is not

viii
necessarily satisfactory when the particle distribution is irregular. In particular, these
particle methods tend to give severe and spurious pressure fluctuation.
In this thesis, a new particle method addressing the above-mentioned problems is
proposed for 2D large amplitude free-surface motion. Called the Consistent Particle
Method (CPM), it eliminates the use of kernel function which is somewhat arbitrarily
defined. The required partial differential operators are approximated in a way
consistent with Taylor series expansion. A boundary particle recognition method is
applied to help define the changing liquid domain. The incompressibility condition of
free surface particles is enforced by an adjustment scheme. With these improvements,
the CPM is shown to be robust and accurate in long time simulation of free surface
flow particularly for the smooth pressure solution without spurious fluctuation.
The CPM is applied to study different 2D free surface flows, i.e. free oscillation of
water in static tank, dam break in tank with different water depth-to-height ratios, dam
break with obstacle. In the simulation of both gentle and violent free surface motion,
the CPM outperforms the original MPS method in both particle distribution and
pressure solution.
An important free surface problem, 2D liquid sloshing in rectangular tanks is then
studied experimentally and numerically by CPM. A series of sloshing experiments are
carried out making use of a hydraulic-actuated shake table. Standing waves in high
filling tanks, traveling waves in low filling tanks and breaking waves in a closed tank
are well simulated by CPM in terms of free surface profiles and pressure fields. The
CPM solution of pressure history shows tremendous improvement compared with
MPS results. In all cases considered, the CPM solutions of free surface elevation and
pressure are in very good agreement with the experimental results.

ix
List of Figures
Figure 1-1. An example of prismatic LNG tank (Photo: Business Wire website*) ...... 9
Figure 2-1. Effect of liquid density and viscosity in sloshing simulation scanned from
Lee et al. (2007b) ......................................................................................................... 35
Figure 2-2. The number of LNG ships as at year 2006 scanned from Foss (2007) ..... 35
Figure 2-3. The LNG ship orders by containment system scanned from Foss (2007) 35
Figure 2-4. Typical permissible filling levels scanned from Lloyds Register (2008) 36
Figure 3-1. Schematic of a typical reference particle with its neighbor particles ........ 73
Figure 3-2. Algorithm of the MPS method .................................................................. 73
Figure 3-3. (a) The arc boundary recognition method; (b) Example of angle list for
particle A and B ........................................................................................................... 74
Figure 3-4. Incompressibility adjustment of free surface particles.............................. 74
Figure 3-5. Candidate list and neighbor list generation ............................................... 75
Figure 3-6. Flowchart of CPM ..................................................................................... 76
Figure 3-7. A center point surrounded by 24 irregularly spaced points ............... 77
Figure 3-8. Comparison of Laplacian using CPM (Eq. (3-31) ), MPS (Eq. (3-16) ) and
SPH (Eq. (3-39) ) (L
0
=0.1, Test function
4 4
) , ( y x y x + = ) .......................................... 77
Figure 3-9. Convergence test of Laplacian for regular points (Test function
4 4
) , ( y x y x + = ) ............................................................................................................ 77
Figure 3-10. Convergence test of Laplacian in MPS for irregular points (Test function
4 4
) , ( y x y x + = ) ............................................................................................................ 78
Figure 3-11. Convergence test of Laplacian in SPH for irregular points (Test function
4 4
) , ( y x y x + = ) ............................................................................................................ 78
Figure 3-12. Convergence test of Laplacian in CPM for irregular points (Test function
4 4
) , ( y x y x + = ) ............................................................................................................ 78

x
Figure 3-13. An example of irregular nodes (1156 in total) in x-y domain [0.8, 4.2] x
[0.8, 4.2] ....................................................................................................................... 79
Figure 3-14. Analytical and numerical result of Laplacian (Test function
) sin( ) , ( xy y x = ) ........................................................................................................... 80
Figure 3-15. Difference of Laplacian values with analytical result (Test function
) sin( ) , ( xy y x = ) ........................................................................................................... 81
Figure 3-16. Error of Laplacian of different numerical algorithm (Test function
) sin( ) , ( xy y x = ) ........................................................................................................... 82
Figure 4-1. Schematic view of the initial particle distribution for static water tank . 109
Figure 4-2. Time history of hydrostatic pressure at point A by MPS ........................ 109
Figure 4-3. Time history of hydrostatic pressure at point A with d=0.2 m scanned
from Khayyer and Gotoh (2008)................................................................................ 109
Figure 4-4. Time history of hydrostatic pressure at point A with d=0.2 m scanned
from Khayyer and Gotoh (2009)................................................................................ 110
Figure 4-5. Comparison of time histories of hydrostatic pressure at point A ............ 110
Figure 4-6. Particle distribution at t=5 s. ................................................................... 111
Figure 4-7. Hydrostatic pressure field of the whole tank of water at t=5 s................ 111
Figure 4-8. Geometry and initial particle distribution of the dam break example ..... 111
Figure 4-9. Comparison of dam break simulation using MPS with experimental results
.................................................................................................................................... 112
Figure 4-10. Pressure field of the dam break example by MPS method ................... 113
Figure 4-11. Pressure field of the dam break example by CPM ................................ 113
Figure 4-12. Comparison between different weighting function ............................... 114
Figure 4-13. Comparison of CPM solutions using different weighting functions with
experiments ................................................................................................................ 114

xi
Figure 4-14. Effect of weighting functions used in CPM .......................................... 115
Figure 4-15. Effect of influence radius used in CPM ................................................ 115
Figure 4-16. Comparison of the leading edge location in CPM solution with published
results---experimental: Hirt and Nichols (1981); Lagrangian FEM: Ramaswamy and
Kawahara (1987) ........................................................................................................ 116
Figure 4-17. Effect of particle size used in CPM ....................................................... 116
Figure 4-18. Dam break profiles using different particle sizes in CPM .................... 117
Figure 4-19. Dam break profiles using different time step t ................................... 117
Figure 4-20. CPU time for different operations (a) CPU time of CPM; (b) Fraction
over the total time of CPM; (c) Fraction over the total time of MPS ........................ 118
Figure 4-21. (a) A schematic view of the tank; (b) Initial particle distribution ........ 118
Figure 4-22. Comparison of time histories of surface elevation amplitude with
published results......................................................................................................... 119
Figure 4-23. Comparison of time histories of surface elevation amplitude results by
CPM and MPS ........................................................................................................... 119
Figure 4-24. Pressure fields at different time instants for MPS and CPM simulation
.................................................................................................................................... 120
Figure 4-25. Comparison of the free oscillation of liquid for large amplitude .......... 121
Figure 4-26. Pressure contours of the free oscillation of liquid by CPM .................. 121
Figure 4-27. Pressure contours of the free oscillation of liquid for larger amplitude by
CPM ........................................................................................................................... 122
Figure 4-28. Configuration of the tank in dam break example and the positions of the
water depth probes and pressure sensor (Fekken, 1998) ........................................... 122
Figure 4-29. Wave profiles of the dam break example by CPM ............................... 123

xii
Figure 4-30. Comparison of pressure at Point P1 with published results (Fekken, 1998)
.................................................................................................................................... 124
Figure 4-31. Pressure contour of dam break by CPM at (a) t=0.7 s, (b) t=1.475 s and
(c) t=2.95 s ................................................................................................................. 124
Figure 4-32. Comparison of water heights at the four points with published results
(Fekken, 1998) ........................................................................................................... 125
Figure 4-33. Pressure contours of the dam break example by CPM ......................... 126
Figure 4-34. Velocity field of the fluid particles of the dam break example by CPM
.................................................................................................................................... 127
Figure 4-35. Geometry and definition of the dam break with obstacle ..................... 128
Figure 4-36. Initial particle distribution of the dam break with obstacle................... 128
Figure 4-37. Graphical comparisons of the dam break behavior ............................... 129
Figure 5-1. The experimental setup ........................................................................... 159
Figure 5-2. Tank with pressure sensor mounted ........................................................ 159
Figure 5-3. Experimental apparatus and working principle ....................................... 159
Figure 5-4. Definition of parameters for liquid sloshing in a rectangular tank ......... 160
Figure 5-5. Fixing tools of rectangular tank on the shake table ................................. 160
Figure 5-6. Schematic Diagram and picture of wave probe ...................................... 160
Figure 5-7. (a) Wave probe mounted to an adjustable stand (b) Calibration of the
wave probes ............................................................................................................... 161
Figure 5-8. Calibration results of wave probe 1 ........................................................ 161
Figure 5-9. Calibration results of wave probe 2 ........................................................ 161
Figure 5-10. Calibration results of wave probe 3 ...................................................... 162
Figure 5-11. Experimental installation of pressure sensor......................................... 162

xiii
Figure 5-12. (a) Devices used for the calibration of pressure sensor (b) calibration of
the pressure sensor ..................................................................................................... 162
Figure 5-13. Calibration result of the pressure sensor ............................................... 163
Figure 5-14. Displacement transducer and its experimental installation ................... 163
Figure 5-15. High speed camera set-up ..................................................................... 163
Figure 5-16. Securing technique ................................................................................ 164
Figure 5-17. Displacement signal of shake table (5mm amplitude) .......................... 164
Figure 5-18. Comparison of free surface elevation at Point P
1
( 0 . 1 /
0
= ) ............ 164
Figure 5-19. Particle distribution simulated by MPS method and comparison with
experimental results ................................................................................................... 165
Figure 5-20. Pressure contours simulated by MPS method ....................................... 166
Figure 5-21. Particle distribution of MPS simulation with arc method and IA of free
surface particles ......................................................................................................... 167
Figure 5-22. Pressure contours of MPS solution with arc method and IA of free
surface particles ......................................................................................................... 167
Figure 5-23. Comparison of pressure history at Point P
2
( 0 . 1 /
0
= ) ..................... 168
Figure 5-24. Particle distribution of CPM simulation ............................................... 169
Figure 5-25. Pressure contours of CPM solution ....................................................... 169
2
( 0 . 1 /
0
= ) ..................... 170
2
( 0 . 1 /
0
= ) ..................... 170
1
( 0 . 1 /
0
= ) ............ 170
Figure 5-29. Free-surface elevation vs. time ( 0 . 1 /
0
= ) by CPM .......................... 171
1
( 1 . 1 /
0
= ) ............. 172
2
( 1 . 1 /
0
= ) ...................... 172

xiv
1
( 1 . 1 /
0
= ) ............. 173
2
( 1 . 1 /
0
= ) ...................... 173
0
= ) by CPM .......................... 174
1
( 9 . 0 /
0
= ) ............ 175
2
( 9 . 0 /
0
= ) ..................... 175
1
( 9 . 0 /
0
= ) ............ 175
2
( 9 . 0 /
0
= ) ..................... 176
1
( 583 . 0 /
0
= ) ........ 176
2
( 583 . 0 /
0
= ) ................. 176
1
( 583 . 0 /
0
= ) ........ 177
2
( 583 . 0 /
0
= ) ................. 177
Figure 5-43. Maximum and minimum free surface elevations vs. excitation frequency
.................................................................................................................................... 178
Figure 5-44. Maximum and minimum hydrodynamic pressure vs. excitation
frequency.................................................................................................................... 178
Figure 5-45. Maximum and minimum free surface elevations under different
excitation amplitudes ................................................................................................. 179
Figure 5-46. Maximum and minimum hydrodynamic pressure under different
excitation amplitudes ................................................................................................. 179
Figure 5-47. Maximum and minimum free surface elevations vs. filling depths ...... 180
Figure 5-48. Maximum and minimum hydrodynamic pressure vs. filling depths..... 180
Figure 5-49. Standing Wave at the initial stage, t =2.75 s ........................................ 181
Figure 5-50. Traveling wave starts to form, t =4.80 s................................................ 181
Figure 5-51. Multi-crested waves traveling, t =6.10 s .............................................. 181

xv
Figure 5-52. Hydraulic run-up (before formation of bore), t =6.75 s ....................... 181
Figure 5-53. Formation of bore, t =6.90 s ................................................................. 182
Figure 5-54. Bore splits into multi-crested traveling waves, t =7.10 s ..................... 182
Figure 5-55. Multi-crested traveling waves, t =7.50 s .............................................. 182
Figure 5-56. Hydraulic run-up, t =7.85 s ................................................................... 182
Figure 5-57. Comparison of the free surface elevation at Point P
1
............................ 183
Figure 5-58. Comparison of the pressure history at Point P
2
..................................... 183
Figure 5-59. Pressure contours for different time instants by CPM .......................... 184
Figure 5-60. Velocity field of the fluid particles by CPM ......................................... 185
Figure 5-61. Free surface profiles at different time instants. Left: experiments; Right:
CPM simulation. ........................................................................................................ 186
Figure 5-62. Time history of free surface elevations at the right wall ....................... 187
Figure 5-63. Pressure contours for different time instants by CPM .......................... 187
Figure 5-64. Computed pressure history at Point P
2
.................................................. 188
1
.................................................. 188

xvi

xvii
List of Tables
Table 2-1. Material properties of LNG in comparison with water .............................. 33
Table 3-1. Summary and comparison of particle methods .......................................... 71
Table 3-2. Parameters for performance test of Laplacian ............................................ 71
Table 3-3. Parameters for performance test of Laplacian ............................................ 71
Table 4-1. Numerical examples studied in Chapter 4 ................................................ 107
Table 4-2. Parameters used in the study of the effect of least-square weighting
function ...................................................................................................................... 107
Table 4-3. CPU time (s) and fraction per time step in CPM...................................... 107
Table 5-1. Parameters of the tank and liquid for high-filling sloshing ...................... 157
Table 5-2. Parameters of the tank and liquid for low-filling sloshing ....................... 157
Table 5-3. Parameters of the tank and liquid for high-filling sloshing with breaking
.................................................................................................................................... 157

xviii

xix
Nomenclature
A amplitude of external excitation
A matrix used in Eq. (3-21)
5 1
a a coefficient used in Eq. (3-28)
B width of the tank
5 1
b b coefficients used in Eq. (3-28)
5 1
c c coefficients used in Eq. (3-28)
D number of space dimensions
d water depth
5 1
d d coefficients used in Eq. (3-28)
5 1
e e coefficient sused in Eq. (3-28)
E residual error vector used in Eq. (3-21)
E 2-norm of E
f function value vector of fluid particles
) , ( y x f differential function of two variables x and y
g body force vector
H height of the tank
h coordinate difference in x direction
k coordinate difference in y direction
k (subscript) the k-th time step
L length of tank
0
L initial particle distance
m mode number
N number of total particles

xx
n particle number density
*
n temporary particle number density
0
n initial constant particle number density
p pressure
) , ( y x P a particle P with coordinate ) , ( y x
R radius of circle around a particle
r particle distance
r particle coordinate vector
*
r temporary particle coordinate vector
e
r cut-off radius of influence area
t , t time, time step
v particle velocity vector
*
v temporary particle velocity vector
max
V maximum velocity
w kernel function
parameter used in Eq. 3-19
coefficient used in Eq. (3-37)
free surface elevation
coefficient used in Eq. (3-5)
kinematic viscosity of fluid
density of fluid
) , ( y x differential function of two variables x and y
excitation frequency
0
first natural frequency

Chapter 1 Introduction
1
Global growing needs for energy are constantly driving the demand for energy
sources such as natural gas. One critical part of the natural gas supply is transportation.
When a natural gas source is near the market, it can be transported by pipeline.
However, when long distance supply is required, the gas needs to be converted to
liquid state for transportation flexibility, which is liquefied natural gas (LNG). LNG is
made by cooling natural gas to a temperature of approximately minus 163 degrees
Centigrade. At this low temperature, natural gas becomes a liquid and its volume is
reduced by more than 600 times. LNG is easier to store and transport than its gaseous
form since it takes up much less space. These advantages call for the need for more
LNG carriers designed to meet harsher operational requirements.
LNG carriers have usually been operated in fully loaded condition or with a
minimum filling of liquid cargo during ballast voyage. Recently however, there has
been growing demand for membrane type LNG tanks that can operate with cargo
loaded to any filling level. The sloshing induced loads in the tanks at these partial
filling levels is the main concern for vessels operated in this manner. Thus, a better
physical understanding and numerical modeling of sloshing waves in the partially
filled tanks is crucial for the designing of the tank structures and developing
mitigation measures and devices to reduce the undesirable effects of sloshing in the
LNG carriers and storage tanks. The research findings will then greatly enhance the
operational flexibility in LNG transport and delivery as well as safety and cost
effectiveness of LNG vessel design.
2
1.1 Overview
From the viewpoint of competitive energy sources, electricity generation is
increasingly dependent on gas which acts as a flexible contributor. The trend comes
from the worldwide desire to reduce dependence on nuclear power, coal and oil
energy for economic and strategic reasons. Currently the overwhelming majority of
gas is supplied by pipelines in the world. Due to geographical constraints, there is a
lack of connectivity between pipelines or between countries. As an attractive
alternative for flexibility and strategic reasons, gas transportation through vessels
becomes a highly desired approach to trade gas all over the world. The conversion of
natural gas into a liquid state makes the worldwide gas trade more economical and
convenient. LNG is made by cooling natural gas to a temperature of approximately
minus 163 degrees Centigrade. LNG makes the long-distance delivery possible,
especially for some regions where pipeline transport is not accessible. It is a natural
result for energy industry to design and develop more LNG carriers which works
under harsher operational requirements. LNG is expected to play an increasing role in
the natural gas industry and global energy markets in the next several decades.
1.2 Sloshing in membrane LNG tank
To meet the growing demand, many LNG ships and terminals have been
proposed and built. For example, Singapore is building LNG terminal costing about
S$1.5 billion (Energy Market Authority, 2006). There is a trend towards the use of
membrane tanks instead of the self supporting storage systems, mainly due to the fact
that membrane tanks utilize the hull shape more efficiently. Space utilization is an
important consideration on sea vessels. Generally, a membrane-type tank is of
prismatic shape (see Figure 1-1 for example) which is more compatible with the shape
3
of ship hull than other shapes such as sphere. The width and height of a LNG tank are
typically around 30m.
LNG carriers are usually operated under the fully loaded condition or with a
minimum filling level for the tank cooling-down purpose during the ballast voyage.
The typical filling level of the LNG tank is greater than 95 percent or less than 10
percent of the tank height. To overcome this constraint, it would be desired to have
large membrane type LNG tanks that can be operated under any filling level. This
demand has stemmed from the emergence of a spot energy market for LNG
containers. The sloshing load at filling levels other than the fully loaded or ballast
condition is the main concern for vessels operated in this manner.
When strong sloshing occurs, liquid moves against the sides of the container
with gradually increased surface elevation. The large liquid movement creates highly
localized impact pressure on tank walls. If the excitation frequency is near or equal to
the natural sloshing frequency, the high dynamic pressures due to resonance may
damage the tank walls.
In this context, an academically challenging and practically important aspect is
the sloshing-induced loads in the membrane-type tanks. Particularly for partially filled
tanks, sloshing can cause both high loads and fatigue upon the containment system
and the hull structure of the transportation vessel. Besides, severe sloshing poses a
potential threat to the stability of ship motion, thereby restricting the operational
flexibility in LNG transport and delivery. Hence research into liquid sloshing is of
great importance to the energy industry.
4
1.3 Study of liquid sloshing
Liquid sloshing has been a research subject attracting much attention over the
last few decades. There has been a considerable amount of work on the study of liquid
sloshing. Most of the early computational studies on liquid sloshing problems were
based on linear wave theory, assuming that the free surface elevation is small
(Housner, 1957; Abramson, 1996). However, the linear theory will result in big
errors in the time-history response when the external excitation is large or near the
natural frequency of liquid sloshing. Thus in the last few decades, researchers began
to use fully nonlinear wave theory to numerically study and simulate the liquid
sloshing in containers. The numerical study of nonlinear liquid sloshing has been
actively performed since 1970s. Different numerical methods based on mesh such as
finite difference, finite element and finite volume method were applied in the studies
(Wu et al., 1998; Koh et al., 1998; Chen and Nokes, 2005).
Mesh-based methods, however, encounter problems when the sloshing amplitude
becomes large especially when there are possibilities of free surface breaking, since
discontinuity of the domain in such cases can not be modeled in a mesh-based method
without the help of other free surface capturing approaches. Although some free
surface capturing methods such as Volume of Fluid (Hirt and Nicholls, 1981) and
Level Set (Osher and Sethian, 1988) can be used to improve mesh-based methods,
they require complex computer programming to solve extra boundary equations in
order to capture the time varying free surface and update the computational mesh.
Furthermore, problems of numerical diffusion arise owing to the discretization of the
advection terms in the Navier-Stokes equation in mesh-based method using Eulerian
grids.
5
Recently there is a growing interest in developing the next generation
computational methods, namely meshless methods as alternatives to conventional
mesh-based methods. Meshless methods in a Lagrangian description are also named
as particle methods. Particle methods have some outstanding advantages which are
not possessed by mesh-based methods and are expected to out-perform the
conventional mesh-based method in some aspects. For example, particles have a
natural ability to represent the coalescence and fragmentation behavior of breaking
waves, especially when they are used to simulate free-surface sloshing. The numerical
diffusion problem in the conventional mesh-based methods using the Eulerian
formulation does not arise in particle methods in which the Lagrangian formulation is
adopted. Hence, as a potential algorithm to simulate breaking wave phenomenon,
particle methods deserve more research.
Nevertheless, in most of the existing particle methods such as Smoothed Particle
Hydrodynamics (SPH) and Moving Particle Semi-Implicit method (MPS), the
approximation of partial differential operators is determined by a pre-defined
weighting function. The approximation error can be large when the particle
distribution is irregular. As a result, existing particle methods such as MPS and SPH
suffer from pressure fluctuations especially in fluid problems with long time
simulation.
In this study, a new particle method named Consistent Particle Method (CPM) is
developed. The required partial differential operators are computed in a way that is
consistent with Taylor series expansion. A predictor-corrector algorithm is used to
solve the coupled equations efficiently. The Poisson equation of pressure is solved in
the context of incompressible flow. A boundary particle recognition method is applied
6
to help define the changing liquid domain. The proposed method shows better
performance both in the accuracy and stability of the scheme compared with the
original MPS.
1.4 Research scope and objectives
A better understanding and numerical modeling of sloshing waves is vital to the
design of LNG carriers and other similar engineering applications where free surface
motion is the main concern. The proposed research mainly addresses a major
challenge in such problems, i.e. accurate simulation of nonlinear behavior of sloshing
in tank including possible wave overturning and breaking. The key question is how to
predict the maximum sloshing motion and maximum hydrodynamic pressure for a
given set of external excitations.
The first objective of this thesis is to develop a numerical model and solution
strategy suitable for simulating liquid sloshing motion in a moving tank, with specific
attention on the potential application to LNG carriers. To achieve this objective, a
new particle method in the simulation of free surface flow problem will be proposed
and verified by various numerical examples. The numerical simulations by the new
particle method will be carried out to investigate the differences in sloshing induced
loads on the tank at various filling conditions.
The second objective is to conduct experimental study for partial verification of
the numerical model, making use of a shake table facility available in the Structural
Engineering Laboratory of National University of Singapore. The experimental results
will be compared with the numerical simulation results.
7
1.5 Organization of the thesis
This thesis contains seven chapters as follows.
Chapter 1 introduces the background of LNG and motivation for studying LNG
sloshing. An overview of liquid sloshing problems and studies is presented. Based on
that, the scope and objectives of this research are defined.
Chapter 2 of the dissertation covers a detailed literature review in the field of
liquid sloshing in containers. A summary of the state-of-the-art accomplishments to
date is given, including applications and limitations of different numerical methods.
The work done on liquid sloshing motion by conventional mesh-based numerical
method is mostly confined to a sloshing wave without breaking. Particle methods
without mesh are found to be more robust in dealing with large amplitude sloshing
with possibility of wave breaking. The advantages and disadvantages of different
meshless methods are discussed.
Chapter 3 gives a detailed description of a new method called Consistent Particle
Method (CPM). An existing particle method called MPS is introduced first. The
limitations of the MPS method are discussed and demonstrated through numerical
examples. There are mainly three improvements in the proposed new method
compared with other conventional particle method such as MPS and SPH. Firstly, the
discretization of derivatives in the Poisson equation of pressure and gradient model in
CPM is based on Taylor series expansion. Secondly, a more effective free surface
boundary recognition method is introduced, which can greatly improve the stability of
the pressure field. Lastly, the problem of imposing an incompressibility condition on
8
the free surface particles is addressed. Numerical examples are presented in this
chapter to demonstrate the capability of the CPM.
In Chapter 4, six numerical examples are tested to show the performance of
CPM compared with the existing particle method MPS. The numerical solutions of
different free surface flow problems, such as flow due to dam collapse without and
with obstacle, are presented based on the proposed CPM and are compared with
reference solutions or experimental data. Numerical results using the proposed CPM
shows the significant improvement in the pressure field compared with MPS solution.
Chapter 5 gives an experimental verification of the numerical modeling with
focus on the liquid sloshing in rectangular tank. Water sloshing under different
excitations is studied experimentally and numerically. The proposed CPM is again
found to be capable of simulating free surface flows problems. The sloshing wave
patterns in rectangular tanks under different filling depths are studied in this chapter.
The effects of external excitation frequencies and amplitudes are also investigated.
Finally, liquid sloshing at high-filling level with impact on the tank ceiling is studied
and simulated.
Lastly, Chapter 6 presents conclusions and suggestions for future work.
A CD is attached containing computer animation files for some selected
numerical examples. A brief explanation note is given in Appendix A.
9

Figure 1-1. An example of prismatic LNG tank (Photo: Business Wire website*)
*: http://www.businesswire.com/multimedia/home/20081217005080/en/1735488/ExxonMobil-
Technology-Yields-World%E2%80%99s-Largest-LNG-Carrier

10

Chapter 2 Literature Review
11
Free surface flow problems, such as wave breaking near shores and moving ships,
green water on ship decks, liquid sloshing in containers (e.g. LNG tankers) and
interaction of waves with floating structures, have received considerable attention of
researchers over the past few decades (Huijsmans et al., 2004; Greco et al., 2005;
Lohner et al., 2006). Violent free surface flow has a profound impact on offshore and
marine structures (Armenio, 1997; Soulaimani and Saad, 1998; Apsley and Hu, 2003;
Idelsohn et al., 2004). In this thesis, the sloshing phenomenon can be defined as the
highly nonlinear motion of the free surface in a moving partially filled tank. Liquid
sloshing generates dynamic loads on the structure of the tank and thus is an issue of
great concern in the design of membrane-type LNG vessels (Tveitnes et al., 2004).
Sloshing loads in liquid transportation tanks affect not only the structure of ships but
also their movement and stability on sea waves (Kim et al., 2003). This liquid sloshing
may cause loss of human lives, economic and environmental resources owing to the
unexpected failure of the vessels.
Liquid sloshing in storage tanks due to wind and earthquake is also a concern in
design. Various finite element schemes have been developed to study the seismic
response of liquid storage tanks by Brown (1982), Veletsos and Tang (1986) and
Rammerstofer et al. (1990). Balendra et al. (1982a, b) and Yi and Natsiavas (1990)
studied the mode shapes and natural frequencies of cylindrical storage tanks using the
finite element method for the liquid and tank wall.
The sloshing effect coupled with ship motion was studied by Kim et al. (2007)._
Godderidge et al. (2009) investigated the effect of compressibility of fluid and found
12
that the inclusion of fluid compressibility has a significant effect on the pressure
evolution of a sloshing flow. Kim et al. (2010) studied the fatigue strength of the
insulation system of MARK-III type LNG carriers. In contrast to the destructive effect
of liquid sloshing in transportation and storages tanks, liquid sloshing in a container, as
a vibration absorber, has been found to have a positive effect in suppressing structural
vibrations due to external loads. Called tuned liquid dampers (TLDs), they have been
used for tall buildings, long span bridges and offshore structures subjected to wind,
waves and earthquakes (Modi and Seto, 1997; Modi et al., 2003). Different tank shapes
are used in TLDs such as rectangular, cylindrical and U-shaped (Ibrahim, 2005). Koh
et al. (1994, 1995) and Shankar and Balendra (2002) proposed multiple TLDs tuned to
several natural frequencies of structures. Balendra et al. (1995, 1999) investigated the
vibration control effect of tuned liquid column dampers (TLCD) in various buildings.
An active control system involving a TLCD is developed by Balendra et al. (2001) for
the vibration control of a single-degree-of-freedom tower subjected to wind excitation.
2.1 Research works involving mainly experimental study
Due to the complexity of sloshing, there has been a considerable amount of work
carried out to understand the complex sloshing behavior and to design appropriate
devices to suppress it. One of the earlier experimental investigations of nonlinear, free-
surface standing waves was reported by Taylor (1953) who focused on the wave crest
in the center of a rectangular tank. Using different scaled model tanks, Strandberg
(1978) conducted experiments to investigate overturning of moving vehicles due to
liquid sloshing under different working conditions. It was found that the overturning
limit of a half filled tank vehicle could be reduced to half of that of a fully filled tank
vehicle where no liquid sloshing occurs.
13
Wang et al. (1996) experimentally studied the waves in a water-filled circular
tank excited by two shakers at opposite sides of the tank. The excitation frequency
used was near one of the natural frequencies of sloshing waves. Standing waves as
well as breaking waves were observed and investigated in their experiments. Pawell
(1997) conducted experimental studies on liquid sloshing interaction with cylindrical
tanks subjected to pitching excitation.
Tveitnes et al. (2004) carried out a series of experiments and proposed a
simplified load formula for estimation of sloshing load in preliminary design. La
Rocca et al. (2005) investigated the problem of sloshing waves of a two-liquid system
experimentally and theoretically. The Lagrangian variational approach was used in the
mathematical model by applying to the potential formulation of the fluid motion. The
theoretical solution of the mathematical model gave good agreement with their
experimental results. Rognebakke et al. (2005) conducted a series of tests on a scaled
tank. They investigated the high filling impacts by using pressure sensors and
analyzing the pressure data statistically. Akyildiz and Unal (2005) investigated the
pressure distributions at different locations and three-dimensional (3D) effects on
liquid sloshing experimentally. Sloshing in a rectangular tank at a scaled model was
studied for various filling levels. Romero et al. (2006) studied lateral sloshing forces
within scaled vehicle tanks under high filling levels. Recently Yan et al. (2009)
conducted experiments to investigate liquid sloshing in a relatively large size test tank
with Reuleaux triangle type cross-section. The experiments were performed on the
tank with and without laterally placed baffles under three fill conditions.
14
2.2 Analytical study of liquid sloshing
Although experimental study is most direct in obtaining maximum impact
pressure due to violent sloshing, some technical issues in the application of
experimental data to actual tank designs have not been completely resolved (Kim et al.,
2004). For example, the scaling law from model test data to real ships is not clear
(Lloyd's Register, 2005). Furthermore, experiments are time consuming and expensive,
particularly for large size tanks at real working condition. As an alternative method to
experimental approach, analytical methods have been developed for simulation of
liquid sloshing.
An analytical model allows understanding of sloshing mechanics and extensive
parametric study. Analytical formulation of liquid equations is well documented by
many researchers for tanks with various regular geometries. The general equation of
liquid motion in closed containers is often simplified by assuming the container rigid
and impermeable. Furthermore, the liquid is assumed ideal which is inviscid,
incompressible, and irrotational. Capillary or surface tension effects are ignored in a
gravitational field.
Stolbetsov (1967) studied the nonlinear sloshing in a rectangular tank due to
horizontal excitation theoretically. He used a perturbation technique and presented two
types of steady-state solutions. Ockendon and Ockendon (1973) proposed an
analytical scheme for resonant sloshing due to external horizontal and vertical
excitations. A third order asymptotic solution was derived mathematically. Kim et al.
(1996) developed an analytical solution of a partially filled rectangular tank under
horizontal and vertical ground excitation. Interaction of the liquid and flexible wall of
the tank was taken in to account. Their solutions of two-dimensional (2D) analysis
15
agreed well with those by other numerical methods. Gavrilova (2004) developed an
analytical solution of a cylindrical tank with rigid tank wall using the Bubnov-Galerkin
method. The coupling vibration frequency of the system was obtained.
Faltinsen et al. (2000, 2003) studied the liquid sloshing in a square-base tank in
frequency domain. The tank was forced under 3D arbitrary motions with frequency
close to the lowest natural frequency. In their work, steady-state waves were solved
using a Bubnov-Galerkin scheme combined with an asymptotic technique. A
quantitative comparison of the free surface elevation using their proposed theoretical
method with the experimental results was presented. A sensitivity study of the initial
condition was performed. In the work of Faltinsen et al. (2006) the same method was
applied in simulation 3D sloshing in a square base tank with emphasis on the swirling
waves. Analytical solutions were validated by experiments performed by them with a
cubic tank of dimensions 0.6m.
Analytical methods can predict the sloshing waves accurately when the motion is
not violent. But wave overturning and breaking cannot be studied analytically. Ideal
fluid has to be assumed, i.e. inviscid, irrotational and incompressible. In addition,
analytical study of liquid sloshing is often performed for simple boundary conditions.
For complex geometry and extensive fluid-structure interaction, analytical solutions
are difficult, if not impossible, to obtain.
Hence, for general problems, numerical methods are necessary in the study of
fluid problems involving free surfaces (i.e. liquid sloshing) with more complex liquid
properties and boundary conditions.
16
2.3 Numerical study of liquid sloshing
Recent advances in computational methods and computer power make it possible
for numerical methods to be applied to study large motion problems of free surface
flow. Traditional numerical methods are mesh based, such as Finite Difference
Method (FDM), Boundary Element Method (BEM), Finite Volume Method (FVM)
and Finite element Method (FEM). The main common feature of these methods is that
computation is based on a pre-defined mesh.
2.3.1 Mesh-based methods
There has been a considerable amount of work using mesh-based methods in the
simulation of liquid sloshing. Most of the early studies on liquid sloshing problems
were based on linear wave theory. Free surface elevation was assumed to be
sufficiently small so that the nonlinear effects could be neglected. Abramson (1996)
used a linear theory to simulate small amplitude sloshing in a container. Solaas and
Faltinsen (1997) adopted a perturbation theory to investigate sloshing in 2D tanks of
general shape. Linear theory is not accurate in the time-history response when external
excitation is large or near the natural frequency of liquid sloshing. Hence, in recent
years, researchers have begun to use fully nonlinear wave theory to numerically study
liquid sloshing in containers.
2.3.3.1. Boundary Element Method (BEM)
Among the numerical methods, the BEM is often used to analyze nonlinear free
surface problems. Faltinsen (1978) and Nakayama and Washizu (1981) simulated large
amplitude sloshing in 2D rectangular tanks using BEM. Grilli and Svendsen (1990)
examined the corner problems (corners of fluid domain were modeled by double-nodes)
17
and investigated accuracy in the BEM simulation of nonlinear wave flows. Koh et al.
(1998) proposed a coupled BEM-FEM scheme to study the fluid-structure interaction
during liquid sloshing in a 3D rectangular tank. The tank structure was modeled by
FEM while the fluid domain by BEM. Experiments were conducted to validate the
proposed numerical scheme. Good agreement in terms of the free surface elevations
and hydrodynamic pressures was obtained. Dutta and Laha (2000) analyzed the small
amplitude liquid sloshing using a low-order boundary element method. Gedikli and
Erguven (2003) adopted variational BEM to investigate the effect of a rigid baffle on
the natural frequencies of the liquid in a cylindrical tank. Zhang et al. (2004)
developed a fully nonlinear 3D numerical wave tank base on a higher order BEM in
the time domain. Numerical examples were presented to show the good performance
of their mixed-Eulerian-Lagrangian scheme. Huang et al. (2010) developed a time-
domain Green function based BEM to simulate liquid sloshing in tanks. Experiments
were conducted to validate the numerical simulation results.
An important feature of the BEM is that only the boundary has to be discretized in
order to carry out the integrations. Because the interior of a solution domain is not
discretized, there is much less approximation involved in representing the solution
variables, making data generation much easier (Fenner, 1983). However, the
coefficient matrix in BEM is generally fully populated with non-zero terms, and is not
symmetric. It has been pointed out by Bettess (1981) that for simple elements, the
FEM is more efficient than BEM in which the mesh is the FEM mesh with internal
nodes and elements removed.
2.3.3.2. Finite Difference Method (FDM)
18
The FDM is also widely used in the study of liquid sloshing problems. Chen et al.
(1996) developed an FDM to simulate large amplitude liquid sloshing in 2D container
due to seismic load. Chen and Chiang (2000) used time-independent FDM to study
sea-wave induced sloshing in a floating tank. The fluid was assumed to be inviscid,
incompressible and irrotational. The coupled interaction effect of sloshing fluid and
tank motion was investigated by the FDM. Kim et al. (2004) applied the FDM in
simulating violent sloshing flows in 2D and 3D prismatic tanks. The impact pressure
on tank ceiling was studied. Numerical solutions were compared with existing
experimental data for which favorable agreement was achieved. Frandsen and
Borthwick (2003) and Frandsen (2004) developed fully nonlinear FDM solutions
based on inviscid flow assumption. The sloshing motions were studied in 2D tanks
under both horizontal and vertical external excitations. Gu et al. (2005) coupled the
Level Set technique with the finite difference solver to study two-phase flow in a 3D
square tank. They found the Level Set method robust in tracking the free surfaces.
Valentine and Frandsen (2005) studied 2D sloshing waves in a rectangular tank under
horizontal excitation using FDM. The evolutions of sloshing of two-layer and three-
layer fluid systems were investigated. Chen and Nokes (2005) developed a novel time
dependent FDM for simulation of 2D sloshing motion in a tank. A fully nonlinear
model was developed where fluid viscosity was included. The numerical solution of
2D waves was compared with other published results and good agreement was
obtained. Lee et al. (2007c) studied the coupling effect of liquid sloshing in LNG tank
with ship motion. The FDM with SURF scheme was applied to simulate liquid
sloshing. Liu and Lin (2008) studied 3D liquid sloshing in rectangular tanks using
FDM. VOF was used to capture the free surface. Wu and Chen (2009) developed a 3D
time-independent FDM to study sloshing waves in a square-base tank under coupled
19
surge-sway motions. Five types of waves under various excitation angles and a wide
range of excitation frequencies were presented.
The most attractive feature of FDM is that it is relatively easy to implement. Its
basic form is, however, restricted to handle regular shapes and simple boundaries. To
handle complicated geometries, FEM is more straightforward to apply (Bathe, 1996).
2.3.3.3. Finite Element Method (FEM) and Finite Volume Method (FVM)
The FEM has been extensively used in the study and simulation of the liquid
sloshing problems. Ikegawa (1974) analyzed sloshing liquid under single component
of horizontal excitation using FEM. Nakayama and Washizu (1980) used FEM to
analyze nonlinear sloshing of liquid in a 2D rectangular tank subjected to pitching
excitations. Wu and Eatock Taylor (1994) developed 2D nonlinear solutions of the free
surface flows using FEM. Numerical examples were presented to validate the proposed
scheme. Nomura (1994) used the FEM to investigate the interaction between sloshing
viscous liquid and rigid body. Wu et al. (1998) gave a broad account of both 2D and
3D sloshing problems based on FEM. In their paper, the potential flow assumption,
where the viscosity of the fluid was neglected, was made.
Pal et al. (1999) developed a numerical model to study the dynamics of liquid
inside a thin-walled composite cylindrical tank based on FEM. The effect of tank
flexibility on the sloshing waves was investigated. Recently, Biswal et al. (2003)
developed a finite element formulation to investigate the vibration modes of liquid in a
liquid-filled cylindrical tank. But their results were limited to cylindrical tanks. Wave
breaking was not included in their works. Pal et al. (2003) studied the effect of the free
surface nonlinearity on the coupled response of cylindrical liquid-filled composite
container. FEM was used to model both the fluid and structure domain. A mixed-
20
Eulerian-Larangian approach was developed. Ideal fluid was assumed. The results
were provided for a cylindrical container using a 2D finite element approach.
Bermudez et al. (2003) used FEM to compute sloshing modes in a container with
an elastic baffle. Linear velocity potential formulation in the frequency domain was
adopted in their work. Turnbull et al. (2003) developed a -transformed FEM for the
simulation of 2D free surface wave motions based on potential theory. Forced sloshing
in a base-excited rectangular tank was studied with the proposed numerical algorithm.
Cho and Lee (2004) simulated a large amplitude liquid sloshing in 2D tanks using fully
nonlinear FEM. Sudharsan et al. (2004) studied the large amplitude sloshing in a 2D
rectangular tank using FEM. Fluid structural interaction was investigated by the
proposed scheme. Ideal fluid was assumed in these works.
Mitra and Sinhamahapatra (2005) studied the coupled slosh dynamics of liquid in
containers using pressure based FEM. The analysis was, however, restricted to linear
problems where small amplitude wave was assumed. Wang and Khoo (2005) studied
nonlinear sloshing in rectangular container under random excitations. FEM solutions
were obtained using the fully nonlinear potential wave theory. The spectra of random
waves and forces were investigated. The nonlinear effects of the random waves were
studied and typical nonlinear features of the waves were captured. Wang et al. (2006)
proposed a damping estimation scheme of liquid sloshing with small amplitude based
on FEM. The potential wave theory was used in their work. The proposed scheme is
not applicable to problems with large sloshing amplitude or large liquid viscosity.
Guillot (2006) developed a discontinuous Galerkin FEM and applied it to liquid
sloshing problems with shallow water depths. Eatock Taylor et al. (2008) proposed a
coupled FE and BE model to study nonlinear transient waves in numerical wave tanks.
21
A mixed Eulerian-Lagrangian formulation was implemented in quadratic iso-
parametric elements. Wave overturning was not captured in their numerical model.
Some of the above-mentioned works such as Nakayama and Washizu (1980),
Abramson (1996), Solaas and Faltinsen (1997), Dutta and Laha (2000), Gedikli and
Erguven (2003) and Wang et al. (2006) concentrated on small amplitude sloshing,
while others like Nakayama and Washizu (1981), Cho and Lee (2004), Wang and
Khoo (2005), Chen et al. (1996), Frandsen and Borthwick (2003) and Sudharsan et al.
(2004) were for relatively large amplitude motions without breaking.
The FVM is widely used in computational fluid dynamics where the solution
favors simpler and lower order approximation within each cell (Bucchignani et al.,
2004; Ahmadi et al., 2007; Greaves, 2007). Bucchignani (2004) studied 2D sloshing in
rectangular tank using FVM based on potential flow. Zhang et al. (2005) developed a
FVM code for the numerical simulation of free surface flow in a container. Recently,
Ming and Duan (2010) investigated liquid sloshing in rectangular tank using FVM
based on unstructured grid. A high order VOF method was adopted in their study to
capture the free surface.
Finite element and finite volume analysis can account for complex geometry with
relative ease. On the other hand, finite element solution can be time consuming and
expensive. In addition, if local nonlinear effects such as overturning and breaking
waves are considered, mesh-based methods like FEM, FVM and FDM meet
difficulties in simulating waves involving discontinuity of liquid motion. To solve this
problem, different interface capturing methods have been proposed by many
investigators. Updating of the free surface is a key factor in the identification of the
flow domain as well as in the application of free surface boundary conditions. The
22
error accumulated in the free surface tracking as time progresses may cause numerical
instability in the sloshing response (Fletcher, 1991; Chen et al., 1996).
2.3.3.4. Free surface capturing method used in mesh-based methods
The most well known approaches to capture the free surfaces are PIC (particle-in-
cell), MAC (marker-and-cell) and VOF (volume of fluid) methods. The PIC scheme
uses particles on the free surfaces and FDM to solve the governing equations (Harlow,
1963). Another similar approach, MAC method, is based on Lagrangian concepts and
can treat overturning waves and reentry inception with simple logic. Marker particles
which move with the fluid are used in MAC method to track the movement of free
surfaces (Harlow and Welch, 1965). MAC has been widely used to solve complex
computational fluid dynamic problems (J ohnson, 1996). Mikelis and J ournee (1984)
simulated 2D liquid sloshing using the MAC method. A series of experiments were
conducted on scaled tanks in their work. The measured pressures were compared with
their numerical transient solution and reasonable agreement was achieved. Armenio
and La Rocca (1996) studied the sloshing of water in rectangular containers with the
filling depths of liquid in shallow water hypotheses numerically with modified form of
MAC method. Experiments were carried out to verify the performance of the
numerical solutions.
Though PIC and MAC are flexible and robust, they are quite complicated in
programming and need additional storage required for locating the marker particles,
and the additional programming complexity to locate the cells containing the free
surface. It significantly increases the computational effort especially for 3D cases
(Griebel et al. 1998). More recently Hirt and Nicholls (1981) developed the VOF
method in capturing the moving boundaries of fluids. VOF method solves an
23
additional partial differential equation for the volume fraction at each time step besides
the conservation equations. This method can define sharp interfaces and is robust.
Nevertheless, tracking and reconstruction of free surfaces remains complicated and
difficult, especially in three dimensions (Qian et al., 2006). Osher and Sethian (1988)
proposed a Level-Set method to deal with moving boundaries. This method also needs
to solve an additional level-set function except for the conservation equations.
All of the above free surface capturing methods can properly compute the
instantaneous free surface displacement. However, they all require complex computer
programming in order to treat the time varying free surface boundary and update the
computational mesh. Furthermore, the problem of numerical diffusion arises owing to
the discretization of the advection terms in the Navier-Stokes equations in the mesh-
based methods using Eulerian grids.
In summary, mesh-based methods are efficient when the sloshing amplitude is
small, where several assumptions can be made to help solve the problem. Nevertheless,
mesh-based methods have met difficulties in simulating waves involving discontinuity
of liquid motion (e.g. wave breaking). These methods may suffer from mesh distortion
in problems with extremely large fluid motion if no additional effort of free surface
capturing scheme, such as VOF method, is introduced. Even with some free-surface
tracking techniques incorporated, mesh-based methods using Eulerian formulations
suffer from the problem of numerical diffusion. In addition, tracking of free surfaces
require a complex and time consuming algorithm to update the rapidly changing
nonlinear boundary. Furthermore, many aforementioned works about liquid sloshing
problems are limited to cylindrical tanks because most of the storage tanks are of
24
cylindrical shape. Little research has been reported on liquid sloshing in rectangular
tanks in the context of large motions.
2.3.2 Meshless methods
In recent years, a new generation of numerical methods has been developed, i.e.
meshless (mesh-free) methods that outperform conventional mesh-based methods in
dealing with discontinuous motion. Meshless methods have been found to have
advantages in dealing with problems for large-amplitude free surface flows, moving
interfaces, large motion and complex and deformable boundaries. Liquid sloshing with
the possibility of wave breaking belongs to the class of large motion problems with a
free surface. In particular, meshless methods have good potential in liquid sloshing
problems when the sloshing amplitude is large, especially with possible breaking of
free surface. The main idea of meshless methods is to obtain numerical solutions of
partial differential equations through a set of particles instead of meshes. Since no
numerical meshes are needed, multi-scale resolution problems can be modeled with
relative ease in meshless methods through adjusting the particle distribution densities
over the solution domain. Particles have a natural ability to represent the coalescence
and fragmentation behavior of breaking wave, especially when they are used to
simulate free-surface sloshing. In addition, the problem of numerical diffusion in
mesh-based methods does not arise in Lagrangian formulated particle methods. Hence,
with great potential to simulate breaking wave phenomena, meshless methods have
recently attracted increasing attention.
Onate et al. (1996) proposed a meshless method called the Finite Point Method
using a Lagrangian formulation to solve fluid mechanics problems. Idelsohn et al.
(2001, 2003) generalized this meshless idea with finite-element type approximations.
25
They later proposed a new method called the Particle Finite Element Method (PFEM)
using the same meshing technique (Idelsohn et al., 2004). The PFEM method was
robust to solve fluid-structure interaction problems with moving free surfaces.
However, extra effort was needed to generate a new mesh at each time step, and this
made the method inefficient in computation.
Another group of meshless methods using particles to describe the motion of the
fluid based on the Lagrangian formulation was developed. The first idea on this
approach was proposed by Monaghan (1988, 1994) for the treatment of astrophysical
hydrodynamic problems. The method was called Smooth Particle Hydrodynamics
(SPH). Many other researchers later generalized this method to solve fluid mechanic
problems. For example, Shao and Lo (2003) developed an improved SPH to simulate
Newtonian flows with a free surface. They tested the method by a typical 2D dam
break example. Souto-Iglesias et al. (2004) simulated sloshing in the anti-roll tanks
using SPH. Souto-Iglesias et al. (2006) investigated the sloshing moment amplitudes in
a rectangular tank under different rolling frequencies using SPH. Experiments were
conducted for validation of the numerical solutions. Very recently Fang et al. (2009)
proposed an improved SPH method for simulating free surface flows with viscous
fluids. Two enhanced variants were presented. Better accuracy and stability was
achieved compared with the original SPH in the work. Delorme et al. (2009)
investigated the impact pressure in shallow water sloshing under forced rolling motion
experimentally and numerically using SPH. Traveling waves and breaking waves
under resonant sloshing condition were observed. A modified SPH scheme was also
proposed to study the density effect to the impact pressure. Good agreement in terms
of free surface shape and global dynamic of the flow was shown between experimental
and numerical results.
26
In the SPH method, the incompressible flow condition is simulated as the limit of
the compressible NavierStokes equations with some stiff equation of state (Monaghan,
1988). The compressibility is adjusted by changing the numerical speed of sound in the
equation of state, which helps to control density fluctuation to some pre-defined ratio.
The incompressible limit is obtained by choosing a very large speed of sound in the
equation of state in SPH. However, the large value of the speed of sound restricts the
time step to be very small due to the CourantFriedrichsLewy condition (Courant et
al., 1967). Furthermore, the compressibility adjustment is case-dependent and it has to
be calibrated for every specific fluid dynamic problem.
Perrone and Kao (1975) and Liszka and Orkisz (1980) proposed a so called
Generalized Finite Difference (GFD) method for structural problems. Chew et al.
(2006) proposed an Arbitrary Lagrangian-Eulerian (ALE) formulation of GFD method
and applied it in vortex flow and flows past moving bodies. Tiwari and Kuhnert (2007)
presented a meshless method for two-phase incompressible flows named finite pointset
method. Nevertheless, neither of these methods has been applied in the study of
incompressible free surface flow problems.
A finite point Euler algorithm for aerodynamics was proposed by Batina (1992,
1993). Recently Mendez and Velazquez (2004) applied this method to simulate 2D
laminar unsteady flows. Recently, Liu et al. (2005) developed a Lagrangian particle
method named finite particle method (FPM) and applied it to modeling incompressible
free surface flow problems. Nevertheless, the same approach of imposing the
incompressibility condition as that in SPH was adopted in these studies.
Hence, a particle method with more straightforward modeling of incompressible
fluid was proposed by Koshizuka and Oka (1996) and Koshizuka et al. (1998). Called
the Moving Particle Semi-implicit (MPS) method, it adopts particle number density in
27
the computation of density parameter. The incompressibility condition is imposed by
solving the Poisson equation of pressure.
The MPS method has been applied in many fluid mechanic problems. In
comparison with the SPH method, the MPS method showed better stability (Ataie-
Ashtiani and Farhadi, 2006). Koshizuka et al. (1998) used MPS to simulate breaking
waves. Yoon et al. (1999) proposed a hybrid particle-gridless method for
incompressible flows based on MPS and verified their scheme by a liquid sloshing
problem in a 2D rigid tank. Zhang et al. (2006) improved the original MPS method and
applied it to convective heat transfer problems. Ataie-Ashtiani and Farhadi (2006)
improved the stability of MPS and used it to simulate free-surface flows. Shibata and
Koshizuka (2007) analyzed ship water impact on a deck. Recently, Lee et al. (2007a)
developed a coupled MPS and FEM method to analyze sloshing in tanks with elastic
thin walls that undergo large deflection. The results showed that their method could
solve current benchmark problems. Sueyoshi et al. (2008) used MPS to simulate wave-
induced nonlinear motion of a 2D floating body.
Actually the MPS method evolves from the well-known Marker and Cell (MAC)
method (Harlow and Welch, 1965). The predictor-corrector algorithm of MAC is
inherited by MPS for solving Poisson equation of pressure. The main difference is that
MAC evaluates the Laplace (or Laplacian) operator and gradient operators based on a
pre-defined regular mesh while MPS establishes an approximate model of Laplace and
gradient operators based on neighbor particles. Based on the transient diffusion model
and weight function, Koshizuka and Oka (1996) proposed discrete Laplace and
gradient operators to form the basis of MPS. Although subsequent researchers
28
(Khayyer and Gotoh, 2008, 2009; Zhang et al., 2006) made some modifications, the
basic idea of the approximation of these two operators are the same.
Nevertheless, similar to the SPH method, the MPS method employs a smoothing
function to approximate a function value and its derivatives. The approximation of
partial differential operators is determined with a pre-defined kernel function. The
accuracy of the approximation is affected by the particle distribution since the kernel
function relies on the inter-particle distance only. As a result, MPS suffers from
pressure fluctuations especially in fluid problems with long time simulation.
To address the above-mentioned problems, a new particle method named
Consistent Particle Method (CPM) is developed in which the discretizaiton of the
equations is based on derivatives consistent with a Taylor series expansion. A two-step
algorithm is used to solve the Navier-Stokes equations. The Poisson equation of
pressure is solved in the context of incompressible flow. A boundary particle
recognition method is applied to help determine the solution domain. Hence, the
pressure field is obtained with no spurious pressure fluctuation in the simulation of
both gentle and violent free surface flows. This is a great advantage of CPM
considering that the pressure solution has not been well simulated in conventional
particle methods. The detailed formulation of CPM will be explained in Chapter 3.
2.4 LNG and LNG sloshing
2.4.1 LNG and its carrier system
Natural gas is colorless, odorless, non-corrosive and non-toxic. Liquefied natural
gas (LNG) is natural gas that has been cooled to a liquid at -260F (-162.2C) at
atmospheric pressure. LNG reduces its volume by more 600 times, making it more
29
economical to store and transport (Foss, 2007). LNG is a liquefied hydrocarbon
mixture consisting predominantly of methane (CH
4
) with a small admixture of higher
hydrocarbons (primarily ethane). For many reasons it is sufficient to consider it to be
liquid methane (Webber et al., 2008). The material properties of LNG are listed in
Table 2-1 and compared with those of water.
Lee et al. (2007b) performed a parametric sensitivity study on the parameters of
sloshing loads in LNG tanks by using a computational fluid dynamics program. They
investigated the effect of parameters such as liquid viscosity and density,
compressibility and ullage pressure of gas and laminar or turbulence fluid model in the
study of sloshing load. They explored the sloshing impact pressure with the liquid
property of water and LNG. Figure 2-1 shows the effect of density and viscosity of
liquid in simulation of impact pressure of sloshing waves in the work of Lee et al.
(2007b). They concluded that the effects of viscosity and density ratio are
insignificant, as was also pointed out by Bass et al. (1985). In this context,
experimental sloshing tests using water in the laboratory can thus be considered a good
approximation of LNG sloshing.
Specially designed tankers (double-hulled ships) are used to transport LNG
between terminals. LNG is stored in a containment system within the inner hull at
atmospheric pressure and cryogenic temperature (-260F). There are mainly three
types of cargo containment systems used in modern standards, i.e. the spherical (Moss)
design, the membrane design and the structural prismatic design. The Moss and
membrane design systems are most commonly used historically (Foss, 2007). The
Moss design system has better resistance to sloshing forces compared with prismatic
membrane type. The tanks in the Moss design, however, do not fit the contour of the
30
ships hull, which results in a diseconomy. There is thus a trend towards the use of
membrane system in recent years (Foss, 2007). Figure 2-2 shows that 51 percent of
LNG ships were membrane design in year 2006. The rising trend towards membrane
design is evident in LNG ship orders shown in Figure 2-3, which shows that 85 percent
of LNG ships in orders from 2005 to 2010 are membrane design.
2.4.2 Sloshing phenomena
Earlier LNG vessels were relatively small and employed spherical tanks, hence
sloshing influences can be neglected. As pointed by Kim et al. (2004) and Chiu (2006),
the spherical tanks of Moss design system are more robust and resistant to the sloshing
effects of LNG, while the membrane systems (which are more efficient in hull space
usage) are particularly sensitive to sloshing loads. With the increase of the LNG
market, there is need to increase the size of tanks for more storage space. This will
essentially cause more severe sloshing impacts for partially filled tanks. There are two
reasons for violent fluid motion in a tank (Abramson, 2003). One is that the frequency
of regular ship motion caused by a sea state is close to the natural frequency of the free
surface. The other reason is that the fluid viscosity of LNG for suppressing sloshing
motions is very small.
Membrane tanks are potentially susceptible to damage due to strong sloshing
since no internal device is used to mitigate fluid motion in a relatively huge space. The
characteristics of the free surface in a membrane tank change with the fill level, as
described below (Lloyds Register, 2005).
(1) At the highest fill levels above about 90%H (H =maximum height of LNG
tank), the vapor space seeks the instantaneous highest point in the tank and this has a
major influence on the fluid motion near the surface.
31
(2) At high fill levels between 70%H and 90%H, the free surface takes the form
of a standing wave. The oscillatory motion of the surface produces generally vertical
velocity components and a shallow contact angle. A consequence is that air and vapor
may become trapped between the fluid and tank wall causing a cushioning effect of the
sloshing impact.
(3) At intermediate levels (less than 70%H), the free surface develops a traveling
wave or bore, which is characterized by a high velocity and high impact angle but with
low vapor entrapment. Such a wave can transfer a significantly higher impact force to
the walls of the tank than a standing wave. Furthermore, cushioning effects due to
vapor entrapment are considered to be less likely.
In this study, sloshing in membrane type system is considered. For simplicity a
rectangular-shape tank is used to represent the tank in membrane type system. The
flexibility of the tank wall, which is commonly assumed rigid in many research works
(Ibrahim, 2005; Lloyds Register, 2008; Chen et al., 2008), is neglected in present
study.
It is pointed out by Lloyds Register (2008) that, in the most recent LNG ships of
conventional size (up to 155,000m3 capacity) filling is permitted for 70%H and above
for higher filling levels (H is the internal tank height) or 10%L
T
and below for lower
filling levels (L
T
is the internal tank length), as illustrated in Figure 2-4. At mid-fill
height, ship and tank periods may be close to each other which generally lead to the
highest sloshing impacts. Mid-fill heights are therefore not permitted. Recently,
however, there has been growing demand for membrane-type LNG tanks that can
operate with cargo loaded to any filling level which allow for maximum operational
flexibility. The sloshing loads at filling levels other than the fully loaded or ballast
32
condition is the main concern for vessels operating in this manner. Hence, liquid
sloshing around mid-filling is mainly addressed in present study.

33

Table 2-1. Material properties of LNG in comparison with water
LNG*
Water ^
Boiling point 111.7 K 373.1 K
Heat of vaporization 510 kj/kg. K 2257 kj/kg. K
Dynamic viscosity 1x 10
-4
Pa.s 1x 10
-3
Pa.s
Critical temperature 190.4 K 647 K
Critical pressure 4.6 MPa 22.06 MPa
Density 425 kg / m
3
1000 kg / m
3

Phase at room temperature Gas
Liquid

Source:
*: Webber et al., (2008) and NIST Chemistry Webbook.
^: http://en.wikipedia.org/wiki/Properties_of_water; and Haynes, W. M. ed., (2010).

34

35

Figure 2-1. Effect of liquid density and viscosity in sloshing simulation scanned from Lee et al. (2007b)

Figure 2-2. The number of LNG ships as at year 2006 scanned from Foss (2007)

Figure 2-3. The LNG ship orders by containment system scanned from Foss (2007)
36

Figure 2-4. Typical permissible filling levels scanned from Lloyds Register (2008)

Chapter 3 Formulation of Consistent Particle Method
37
3.1 Introduction
Numerical modeling of highly nonlinear free surface flow is difficult because of
the changing boundary condition of free surface, i. e. not only the free surface forms
part of boundaries of the computational domain, but also the positions and velocities of
the free surface vary with time and are not known a priori. In addition, free surface
flow may involve unsteady fragmentation and coalescence processes when wave
breaking occurs.
Particle methods in which Lagrangian particles are used to represent the fluid
domain becomes appealing in the application of free surface flow problems due to their
natural ability in characterization of particle fragmentation and coalescence. Compared
with the conventional numerical methods such as FDM and FVM, particle methods
show some outstanding advantages which are not possessed by conventional mesh-
based methods.
According to the literature review in Chapter 2, we first focus on a promising
particle method, i.e. the moving particle semi-implicit (MPS) method, developed by
Koshizuka et al. (1995). The method is based on solving a Poisson equation of
pressure which is derived from Navier-Stokes equation. By detailed evaluation of the
discrete algorithms used in the MPS method, we have found a problem with MPS and
another particle method SPH. It is that they rely on a pre-defined kernel function in
approximation of the gradient and Laplace operators. Kernel function can cause
instability which is often encountered in particle methods (Ataie-Ashtiani and Farhadi,
2006). Besides, even with properly selected kernel function, the pressure solution of
38
MPS method presents severe fluctuation, as will be shown later in this chapter and
subsequent chapters.
A new particle method named Consistent Particle Method (CPM) is therefore
proposed to overcome this drawback in the MPS method. Instead of using pre-defined
kernel function, the approximations of derivatives in the governing equations are based
on Taylor series expansion. In this chapter, the MPS method is first introduced. The
new method, CPM, is then proposed with three key improvements over the original
MPS. The program codes of MPS and CPM are written in Matlab (v.7.0) and run on
computer equipped with Pentium 4
TM
CPU (duo 3.6 GHz).
3.2 Moving particle semi-implicit method
The moving particle semi-implicit (MPS) method was first developed by
Koshizuka et al. (1995). It is a Lagrangian method developed for simulating
incompressible fluid phenomenon based on Navier-Stokes equation. The fluid is
modeled as an assembly of interacting particles. Computational meshes are
unnecessary. The governing equations are discretized based on the interaction models
of gradient and Laplacian. The pressure field is obtained by solving Navier-Stokes
equation in the whole fluid domain. In addition, as no advection term is used in the
Lagrangian description in the method, there is no numerical diffusion problem which is
present in conventional mesh-based methods.
The MPS method has been applied to a wide range of problems such as
hydrodynamics (Koshizuka and Oka, 1996; Ataie-Ashtiani and Farhadi, 2006; Shibata
et al., 2009; Gotoh et al., 2005; Sueyoshi et al., 2008), wave breaking (Koshizuka et al.,
1998; Khayyer and Gotoh, 2008; Alam et al., 2007), wave impact (Khayyer and Gotoh,
39
2009), fluid-shell interaction (Lee et al., 2007a) and heat transfer problems (Zhang et
al., 2006). Various kernel functions were studied by Ataie-Ashtiani and Farhadi (2006)
to improve the stability of the MPS method. Different treatments of the boundary
conditions and methods of solving the Poisson equation of pressure were considered to
improve the accuracy and stability of the MPS method (Zhang et al., 2006, Khayyer
and Gotoh, 2008; Khayyer and Gotoh, 2009). Nevertheless, even with these reported
improvements, the problem of pressure fluctuation in long time simulation is still not
completely resolved.
The detailed analysis and improvements will be stated in the following section.
3.2.1 Governing equations
The governing equations of viscous Newtonian fluid are expressed by the
conservation laws of mass and momentum as follows (Lee et al., 2007a):
0
1
= + v
Dt
D
(3-1)
g v
v
+ + =
2
1
p
Dt
D
(3-2)
where is the density of the fluid, v the particle velocity vector, p the pressure, the
kinematic viscosity and g the body force (normally the gravity) .
The computational procedure of the MPS method comprises two steps. The first
step is the predictor step in which the velocity field is computed without enforcing
incompressibility. In the second step, namely the corrector step, incompressibility is
enforced through the calculation of Poisson equation of pressure.
Predictor-Corrector algorithm
40
(1) Predictor step. Determine the temporary particle velocities and positions by
considering only the body force and viscosity terms in Eq. (3-2).
t
k k k
+ + = ) (
2 *
g v v v (3-3)
t
k
+ =
* *
v r r (3-4)
where
k k
t t t =
+1
, ) , (
k k k
t x v v = , ) , (
k k k
t x = .
k
v and
k
r are the particle velocity
and position at time
k
t ;
*
v and
*
r are the temporary particle velocity and position,
respectively.
Incompressibility is not satisfied in this step and the fluid density that is
calculated based on the temporary particle position deviates from the constant density.
(2) Corrector step. Pressure computed from Poisson equation of pressure is used
to enforce incompressibility.
The momentum equation can be split into
g v
v v v v v v v v v
+ +
+
=
+
=
=
+ + 2
* * *
1 1
'
p
t t t Dt
D
k k k k
(3-5)
where
*
1
' v v v =
+ k
and
k
v v v =
* *
.
Substituting g v
v
+ =
2
*
t
from Eq. (3-3) into Eq. (3-5) we can get

1
1
'
+
+
k
k
p
t
v
(3-6)
The mass conservation equation can be split into

) ' ( ) (
)
'
(
1 1 1 1
* * *
1
* * *
1 1
v v v v v + = + =

+
=
+
=
=
+
+ +
k
k k k k
t t t Dt
D
(3-7)
41
The mass conservation equation at the temporary position can be expressed as
*
*
1
v =

Substituting the above equation into Eq. (3-7) gives
) ' (
' 1
v =
(3-8)
By combining Eqs. (3-6) and (3-8) the Poisson equation of pressure is obtained
(Koshizuka et al., 1998):

1
1
*
2
1
2
) (
+
+
+

=
k
k
k
t
p

(3-9)
By imposing the new density
1 + k
=
0
on the right side of Eq. (3-9),
incompressibility is satisfied at the new time step. In the MPS method, particle
number densities
*
n (corresponding to
*
) and
0
n (corresponding to
0
) are used in Eq.
(3-9) (Koshizuka and Oka, 1996; Koshizuka et al., 1998; Gotoh et al., 2005).
After applying appropriate discretization algorithm for the Laplace operator in Eq.
(3-9), a system of linear equations is obtained and can be solved using the existing
solver such as incomplete Cholesky decomposition conjugate gradient method
(Kershaw, 1978) and Gauss-Seidel method (Ataie-Ashtiani and Farhadi, 2006). In
present study the Gauss-Seidel method is adopted.
Once the pressure term is obtained, the new particle velocities can be computed
by Eq. (3-6) . Finally the new position of the particle is updated by
t
k k k
+ =
+ + 1 1
v r r (3-10)
42
The time step is controlled in the computation to satisfy the CourantFriedrichs
Lewy (CFL) condition (Courant et al., 1967; Ataie-Ashtiani and Farhadi, 2006) as
follows,

max
0
2 . 0
V
L
t (3-11)
where
0
L is the initial particle distance and
max
V the maximum particle velocity in the
computation.
3.2.2 MPS formulation
In MPS, the fluid is represented by a set of particles. Each particle carries a mass
m and moves at velocity v. A particle interacts with its neighbor particles covered by
an influence radius through a kernel function ) (r w , where r is the distance between
two particles (Figure 3-1).
The most commonly used kernel function applied in the MPS method is the one
proposed by Koshizuka and Oka (1996):

<
=
e
e
e
r r
r r
r
r
r w
0
0 1
) ( (3-12)
where
e
r is the influence radius of the influence area of each particle.
0
) 4 ~ 2 ( L r
e
=
is commonly used in the simulations (Ataie-Ashtiani and Farhadi, 2006; Gotoh et al.,
2005) where
0
L is the initial distance between two adjacent particles.
The particle number density for particle i is then defined as summation of the
weights of its neighbor particles (Koshizuka and Oka, 1996; Koshizuka et al., 1998),
|) (|
i j
i j
i
w n r r =
(3-13)
43
The gradient operator is modeled as a weighted average of the gradient vectors
between the center particle i and its neighbor particles. The gradient of the pressure is
expressed as (Koshizuka and Oka., 1996),
( )

(
(
=
i j
i j i j
i j
i j
i
w
p p
n
D
p |) (| ) (
| |
2
0
r r r r
r r
(3-14)
where D is the number of space dimensions and
0
n the initial constant particle number
density. As discussed by Koshizuka et al. (1998), the pressure value of any neighbor
particle can be used for
i
p if the distribution of neighbor particles is isotropic. In the
above equation, '
i
p is used instead of
i
p . The value of '
i
p is computed by (Koshizuka
et al., 1998)
} 0 |) (| | { ) min( ' =
i j j i
w j for p p r r (3-15)
By using the minimum value '
i
p among the neighbor particles within the
influence radius
e
r , the inter-particle forces will always be repulsive since '
i j
p p is
positive. This is good for the numerical stability as pointed out by Koshizuka et al.
(1998).
The Laplacian model is proposed by Koshizuka et al. (1998) based on the concept
of diffusion, which is in the form:
( ) | |
=
i j
i j i j i
w p p
n
D
p |) (| ) (
2
0
2
r r
(3-16)
where

=
V V
dV r w dV r w r ) ( / ) (
2
(3-17)
In Eq. (3-16), the physical property (such as pressure) of the center particle is
distributed to its neighbor particles using a weight functionw. Since diffusion is a
44
linear process, the distributions from one particle to the others can be superimposed.
General concepts of these models can be found in the work of Koshizuka et al. (1998)
and Gotoh et al. (2005). The discretization of viscosity components in the viscosity
term is the same as that of the Laplacian term shown in Eq. (3-16) (Gotoh et al., 2005).

3.2.3 Modeling of incompressibility
In the MPS method, density is represented by particle number density defined in
Eq. (3-13). Incompressibility is satisfied by keeping the particle number density
constant. By applying particle number density in the right hand side of Eq. (3-9) and
imposing
0
1
n n
k
=
+
(Koshizuka et al., 1998), the Poisson equation of pressure becomes,

0
0
*
2
1
2
) ( n
n n
t
p
k
=
+

(3-18)
The computational procedure of MPS is summarized in the flowchart shown in
Figure 3-2.
3.2.4 Boundary conditions
3.2.4.1 Wall boundaries
In the MPS method, wall boundaries are also modeled by particles. They balance
the pressure of the inner fluid particles and prevent them from penetrating the wall
(Ataie-Ashtiani et al., 2008). A solid wall boundary is represented by one line of fixed
particles. The Poisson equation of pressure is solved on these particles. Furthermore,
in order to ensure the particle number density is not truncated near the solid boundaries,
several lines of dummy particles are placed outside the wall boundaries. Velocities are
45
always zero at the solid wall particles and dummy particles so as to retain the non-slip
boundary condition.
There are different ways to place the dummy particles. The first one is to use
fixed particles which are spaced according to the initial particle configuration (Gotoh
and Sakai, 1999). Lo and Shao (2002) used another method involving the so called
ghost particle, where mirror particles of inner fluid particles adjacent to the wall
were used. The fixed particle approach is simpler than the ghost particle approach
since the ghost particles keep changing their positions every time step. Additional
computational cost is required in the ghost particle method to generate dummy
boundaries at every time step. Hence in the present study we adopt the fixed particle
approach.
In the fixed particle approach, most researchers applied the same pressure for the
dummy particles as that of the solid wall particles (Lo and Shao, 2002; Ataie-Ashtiani
et al., 2008). As well in ghost particle method, they applied the same pressure for the
ghost particle as that of the respective fluid particle (Colagrossi and Landrini, 2003;
Molteni and Colagrossi, 2009). However, as pointed out by Delorme et al. (2009), a
pressure difference between dummy particles and solid wall particles (or ghost
particles and the respective fluid particles) should be introduced to take into account
the effect of the gravity field and other body forces. Here in this study, we use the
fixed particle approach to place dummy particles and the method of applying pressure
for the dummy particles that is discussed by Delorme et al. (2009).
When moving wall boundaries is encountered, the motion of walls is
implemented by updating the positions of wall particles at every time step. For liquid
sloshing problems in Chapter 5, the tank motion is imposed by this approach.
46
3.2.4.2 Free surface
Since no particle exists outside the free surface, the particle number density of the
free surface particles decrease on this boundary. Thus if a particle satisfies Eq. (3-19),
it is considered to be free surface particle (Koshizuka et al., 1998).
( )
0
*
n n
i
< (3-19)
where is a numerical parameter which is usually between 0.8 and 0.99. It is pointed
out by Koshizuka and Oka (1996) that the simulation results are not sensitive to the
value of in this range. Fragmentation and coalescence of the fluid can be simulated
easily with this simple boundary recognition method. The Dirichlet boundary condition
of pressure is applied to the free surface particles.
3.2.5 Drawbacks of MPS
Although the MPS method can simulate and represent well the fluid general
motion, it produces large pressure fluctuation in the solutions (Khayyer and Gotoh,
2009). Khayyer and Gotoh (2008) pointed out that the spurious pressure fluctuations
seem to be unavoidable in the MPS method for incompressible fluid problems. They
proposed a corrected MPS algorithm to ensure the momentum conservation of the
particles. Their algorithm is further revised by introducing a new source term for the
Poisson pressure equation and allowing slight compressibility in their recent work
(Khayyer and Gotoh, 2009). They demonstrated the improved performance through
several wave breaking cases. However, even with all these improvements, the pressure
fluctuation does not vanish especially when large amplitude particle motion appears.
One main reason of the pressure fluctuation in the MPS method is that free
surface particles are not accurately recognized. Since free surface particles are detected
47
by the number density using Eq. (3-19), some in-domain fluid particles with zero or
low pressure may be falsely considered as free surface particles when the fluid motion
are violent. In this case pressure fluctuation occurs due to incorrect pressure solutions
of the whole domain. The pressure fluctuation in turn produces zero or low pressure
particles resulting in a vicious cycle.
Another reason of pressure fluctuation lies in the approximation of the Laplacian
and gradient model in MPS method. This is a common problem for particle methods
such as SPH and MPS in which the partial differential operators are approximated by a
pre-defined kernel function. The approximation of these operators is highly dependent
on the kernel selected. The accuracy of approximation is not guaranteed when the
particle distribution is irregular, as will be demonstrated in Sect. 3.3.4.
This thesis mainly focuses on overcoming the above-mentioned two problems of
the pressure fluctuation in fluid flows. A new particle method named Consistent
Particle Method (CPM) is proposed in this study. Three main improvements are
included, the most important one dealing with the partial differential operators
approximated based on Taylor series expansion.
3.3 CPM based on Taylor series
3.3.1 Introduction
As discussed in the previous section, although the MPS method has been shown
to simulate the profile of free surface flow, it has some serious drawbacks for pressure
solution. The pressure history shows spurious high fluctuation and sometimes even
negative pressure values appear.
48
In the two popular particles methods MPS and SPH, the discretization of partial
derivatives such as Laplace and gradient operators rely on a kernel function defined
only according to inter-particle distance. The main problem of using kernel function is
that the Laplacian model may become unstable once the distribution becomes irregular.
The consistency of the Laplacian and gradient model in conventional particle method
such as SPH and MPS is not guaranteed. To overcome the shortcoming, a new method
based on Taylor series expansion is proposed.
The method of approximating partial derivatives based on Taylor series expansion
is the core algorithm in the generalized finite difference method (GFD). GFD method
is a meshless method which is proposed in the early 1970s for structural problems. One
of the early investigations was by Perrone and Kao (1975) who studied large deflection
response of a flat membrane. Liszka and Orkisz (1980) presented a FIDAM (Finite
Difference at Arbitrary Mesh) code based on GFD using moving least-square
interpolation, and applied it in structural mechanics. Nevertheless this method has not
received much attention possibly because the method does not have the favorable
properties such as symmetric positive and well conditioned matrix (Luo and Haussler-
Combe, 2002). Luo and Haussler-Combe (2002) developed a GFD based on
minimizing global residual and applied in solving solid mechanics problems. Gavete et
al. (2003) improved the GFD by using variable radius of influence. The results
obtained were compared with other meshless methods. Nevertheless, no application of
the method to real problems was presented. Recently Chew et al. (2006) proposed an
Arbitrary Lagrangian-Eulerian (ALE) formulation of GFD method. They applied the
improved scheme to vortex flow and flows past moving bodies. To our knowledge,
however, GFD method has not been used in the study of incompressible free surface
flows. The main reason may be that the GFD suffers from ill-conditioning of matrix if
49
there is poor spatial arrangement of the support particles, which is very likely to occur
in free surface flow problems with large particle motions and severe irregularity of
particle distribution. Hence the application of GFD in large motion fluid problems with
free surface is limited.
By incorporating the incompressibility model of MPS using number density in
GFD, we may simultaneously overcome the problem of ill-conditioning in GFD
method and improve the approximation accuracy not achieved by MPS method. In this
context, we propose a new particle method which combines the merits of GFD in
accurate approximation of the gradient and Laplace operators and MPS in obtaining
even particle distribution through the modeling of incompressibility. Since in the
proposed method the discretiaztion of partial differential operators are consistent with
Taylor series expansion, we name it Consistent Particle Method (CPM).
Our algorithm in approximation of differential operators is similar to GFD
developed by Perrone and Kao (1975) and Liszka and Orkisz (1980) for structural
mechanics. And it resembles the finite point methods of Onate et al. (1996, 2000) and
the finite point solver of Mendez and Velazquez (2004). Mendez and Velazquez (2004)
pointed out that the evolution of these particle methods (GFD, finite point solver, finite
point method and our proposed CPM) leads to a closer resemblance between them.
Nevertheless, the main difference between our work and the previous research is that
we impose strict incompressibility of free surface flows by computing particle number
density and solving Poisson equation of pressure. A free surface recognition method
and free surface incompressibility adjustment approach are applied to help define the
boundary and impose accurate boundary condition. In this sense, the drawbacks of
50
low accuracy in some particle methods (such as SPH and MPS) and ill-condition
matrix in GFD due to poor spatial particle distribution are resolved.
A brief summary and comparison of particle methods is listed in Table 3-1. As
shown in the table, our proposed CPM does not require mesh and pre-defined kernel
function. Incompressible fluid is modeled by solving Poisson pressure equation. By
incorporating the arc free surface recognition method it can be applied to simulate free
surface flow problems.
In this section we present the development of the discretization equations based
on Taylor series similar to GFD. The method is first validated through numerical tests
of Laplace operator evaluation. The performance of the proposed CPM method in
modeling free surface flow will be demonstrated in the next Chapter.
3.3.2 Approximation of gradient and Laplacian by Taylor series
For any differentiable function ) , ( y x f , the Taylor series expansion around a
point ) , (
0 0
y x P can be expressed in the form (Gavete et al., 2003)
) ( ,
2
1
, ,
2
1
, , ) , (
3
0
2
0 0
2
0 0 0
r O f k hkf f h kf hf f y x f
yy xy xx y x
+ + + + + + = (3-20)
where
0
x x h = ,
0
y y k = , ) , (
0 0 0
y x f f = ,
x
y x f
f
x
=
) , (
,
0 0
0
,
y
y x f
f
y
=
) , (
,
0 0
0
,
2
0 0
2
0
) , (
,
x
y x f
f
xx
=
y x
y x f
f
xy

=
) , (
,
0 0
2
0
,
2
0 0
2
0
) , (
,
y
y x f
f
yy
= ,
2 2
k h r + = .
For the reference particle ) , (
0 0
y x P with its neighbor particles j (j=1,2N), if the
values of the function ) , (
j j
y x f are all known at the neighbor particles, we may
truncate the Taylor series to the second order and approximate the derivatives of the
51
function at the reference particle ) , (
0 0
y x P by solving a system of linear equations
(Liszka and Orkisz, 1980; Chew et al., 2006).
Writing Eq. (3-20) for each of the neighbor particles, we derive a set of linear
equations
| |{ } { } 0 = f f A D (3-21)
with
| |
(
(
(
(
(
(
(
(
=
2 2
2
2 2 2
2
2 2 2
2
1 1 1
2
1 1 1
2
1
2
1
2
1
2
1
2
1
2
1
N N N N N N
k k h h k h
k k h h k h
k k h h k h

A (3-22)
{ } { }
=
0
0 2
0 1
0
0
0
0
0
,
,
,
,
,
,
f f
f f
f f
f
f
f
f
f
D
N yy
xy
xx
y
x
f f (3-23)
where
0
x x h
j j
= ,
0
y y k
j j
= and N is the total number of particles in the influence
area of particleP.
The minimum number of neighbor particles required in Eq. (3-21) to determine
{ } f D vector is five since we are seeking solutions of five unknowns. With exactly five
neighbor particles, however, the matrix in Eq. (3-22) is highly susceptible to ill-
conditioning due to irregular spatial arrangement of these particles. To overcome this
problem, more than five particles can be used in order to improve the accuracy of the
approximation (Liszka and Orkisz, 1980). In this case, the over-determined linear
equations is solved by minimizing the residual error vector
52
| |{ } { } f f A E = D (3-24)
We consider the following norm of E with weight factors incorporated (Gavete et
al., 2003)
| | | |
=
+ + + + + =
N
j
j yy j xy j j xx j y j x j j
w f k f k h f h f k f h f f
1
2
0
2
0 0
2
0 0 0
, , , , , E (3-25)
where
j
w is the weighting function used in the weighted least-square solution of the
equation. Incorporating the weight function
j
w in Eq. (3-25) is equivalent to solving
the over-determined Eq. (3-21) using the weighted least-square approach. Different
least-square weighting functions may be selected (Liszka and Orkisz, 1980; Gavete et
al., 2003; Chew et al., 2006).
By minimizing the norm

{ }
0 =
f
E
D
(3-26)
a set of five equations with five unknowns is obtained as follows
(
(
(
(
(
(
(
(
(
(
(
(

0
0
0
0
0
4
2
3
2
2 2
2
3
2
2
2
3
2 2 2 2
3
2 2 2 2 2
2 2
2
3
2
4
2
2
2
3
2
3
2 2 2
2
2 2 2 2
2
2 2 2
3
2 2 2 2
,
,
,
,
,
4 2 4 2 2
2 2
4 2 4 2 2
2 2
2 2
yy
xy
xx
y
x
j
j
j j
j
j j
j
j
j
j j
j
j j
j j j j
j j
j j j j j j j
j j
j
j j
j
j
j
j j
j
j
j
j
j j j j
j j
j j j j j j
j j
j j j j
j
j j j j j j
f
f
f
f
f
k
w
k h
w
k h
w
k
w
k h
w
k h
w k h w
k h
w k h w k h w
k h
w
k h
w
h
w
k h
w
h
w
k
w k h w
k h
w k w k h w
k h
w k h w
h
w k h w h w

53

|
|
|
|
|
|
|
|
|
.
|
\
|
=

2 2
2 2
2
2
0
2
2
2
0
2
2
2
0
2
2
2
0
2
2
0
2
j
j
j
j j
j j j j j j j
j
j
j
j j
j j j j j
j j j j j
k
w f
k
w f
k h w f k h w f
h
w f
h
w f
k w f k w f
h w f h w f
(3-27)
All the elements in the first matrix in Eq. (3-27) depend only on the relative
coordinates between the reference particle and its neighbor particles. The main task in
this method is to solve the derivatives in Eq. (3-27) for all the particles in the domain.
As pointed out by Chew et al. (2006), for stationary boundary problems where all
nodes are fixed, one only needs to invert the matrix in Eq. (3-27) once to obtain the
derivatives. However, for moving boundary problems, for example, fluid flows with
free surface, the matrix needs to be solved for every time step when the particles
change positions with time.
The solution of Eq. (3-27) can be expressed as,

|
|
|
|
|
|
|
|
|
.
|
\
|
(
(
(
(
(
(

2 2
2 2
,
,
,
,
,
2
2
0
2
2
2
0
2
2
2
0
2
2
2
0
2
2
0
2
5 4 3 2 1
5 4 3 2 1
5 4 3 2 1
5 4 3 2 1
5 4 3 2 1
0
0
0
0
0
j
j
j
j j
j j j j j j j
j
j
j
j j
j j j j j
j j j j j
yy
xy
xx
y
x
k
w f
k
w f
k h w f k h w f
h
w f
h
w f
k w f k w f
h w f h w f
e e e e e
d d d d d
c c c c c
b b b b b
a a a a a
f
f
f
f
f
(3-28)
where the parameters from
5 1
a a to
5 1
e e are solved from Eq. (3-27).
Hence from Eq. (3-28) the gradient operator is expressed as
( ) ) ( )
2 2
(
0
1
2
5 4
2
3 2 1
2
0
f f
k
a k h a
h
a k a h a w f
j
N
j
j
j j
j
j j j x

(
(
+ + + + =

=
(3-29)
54
( ) ) ( )
2 2
(
0
1
2
5 4
2
3 2 1
2
0
f f
k
b k h b
h
b k b h b w f
j
N
j
j
j j
j
j j j y

(
(
+ + + + =

=
(3-30)
The Laplace operator is obtained as

( )
) (
2
) ( ) (
2
) ( ) ( ) (
0
2
5 5 4 4
1
2
3 3 2 2 1 1
2
0
2
f f
k
e c k h e c
h
e c k e c h e c w f
j
j
j j
N
j
j
j j j
(
(
(
+ + + +
+ + + + + =

=
(3-31)
3.3.3 Main features of CPM
As mentioned before, the method of computing Laplace and gradient operators
based on Taylor series is the core algorithm in the GFD method. The main drawback
of the GFD method is the possibility of obtaining ill-conditioned matrix due to poor
particle distribution. Since the element in the matrices of Eq. (3-27) is only dependent
on the geometrical positions of the reference particle and its neighbors, the distribution
of the particles in the domain greatly affects the stability of the algorithm. Ill-
conditioned matrix can arise from a variety of causes, such as the total number of
neighbor particles less than the number of unknown derivatives in Eq. (3-21) is
obtained, extremely close/separation between some particles, and large number of
particles aligning in a nearly collinear manner. To overcome this problem, many
different approaches have been explored (Luo and Haussler-Combe, 2002; Gavete et
al., 2003).
The selection of particles in the neighborhood of the reference particle is a crucial
factor affecting the accuracy of approximation. The most commonly used criterion in
particle methods is to select the neighbor particles according to their distance from the
reference particle (Monaghan, 1994; Shao and Lo, 2003; Koshizuka et al., 1998). It is
55
simple but fails very often in the GFD method due to the irregular density of particles
(Liszka and Orkisz, 1980). Perrone and Kao (1975) presented an eight segments
method, where the domain around the reference particle is divided equally into eight
quadrants. The nearest particle to the reference particle in each quadrant is selected.
Liszka and Orkisz (1980) found it too complicated and suggested a simple criterion of
four quadrants. Gavete et al. (2003) also adopted the same approach as Liszka and
Orkisz (1980).
Although the quadrant criterion can improve the performance of the
approximation, it cannot completely avoid ill-conditioned matrices. Gavete et al. (2003)
suggested variable influence radius for each particle. The residual at each point is
obtained first. If ill-conditioned matrix (i.e. large residual value) for any particle is
detected, the radius of that particle is increased to capture more neighbor particles. Luo
and Haussler-Combe (2002) proposed an improved method based on minimizing
global residual with respect to particle parameters. However, they found their results in
their method are sensitive to the number of neighbor particles or the size of domain of
influence.
Nevertheless, no matter what criteria of particle selection are used, the ill-
conditioned matrix cannot be fully avoided due to the irregular density of particles. A
possible solution is to generate regular density of particles as used in MPS method
where incompressibility is employed by keeping the particle number density constant.
Hence, the proposed method is motivated by combining the merits of GFD and MPS in
consistent representation of differential operators and producing even particle
distribution through the idea of number density, respectively.
56
3.3.3.1 Laplacian and gradient model
According to the derivations in Sect. 3.3.2, the Poisson equation of pressure (see,
Eq. (3-18)) can be written as

( )
0
0 *
2
2
5 5 4 4
1
2
3 3 2 2 1 1
2 2
) (
) (
2
) ( ) (
2
) ( ) ( ) (
(
(
(
+ + + +
+ + + + + =

=
t
p p
k
e c k h e c
h
e c k e c h e c w p
i j
j
j j
N
j
j
j j j i
(3-32)
Incompressibility of the particles inside the fluid domain is satisfied through
keeping the density constant. In the MPS method particle number density is used
(Koshizuka et al., 1998). In the proposed CPM, the same density computation method
is adopted. The Gauss-Seidel method is used to solve the discretized linear equation
system of Eq. (3-32). The viscosity term in Eq. (3-3) is also approximated based on the
same algorithm in Eq. (3-31).
The gradient model of pressure based on the solution of the first derivatives in Eq.
(3-27) is
( ) ) ( )
2 2
(
1
2
5 4
2
3 2 1
2
i j
N
j
j
j j
j
j j j x i
p p
k
a k h a
h
a k a h a w p
(
(
+ + + + =

=
(3-33)
( ) ) ( )
2 2
(
1
2
5 4
2
3 2 1
2
i j
N
j
j
j j
j
j j j iy
p p
k
b k h b
h
b k b h b w p
(
(
+ + + + =

=
(3-34)
In the MPS method, the gradient operator defined in Eq. (3-14) with
i
p replaced
by '
i
p is shown to be good for the stability of the numerical algorithm (Koshizuka and
Oka 1996; Koshizuka et al., 1998). Similarly in our proposed CPM algorithm, a
modified gradient model is used as follows
57
( ) ) ' ( )
2 2
(
1
2
5 4
2
3 2 1
2
i j
N
j
j
j j
j
j j j x i
p p
k
a k h a
h
a k a h a w p
(
(
+ + + + =

=
(3-35)
( ) ) ' ( )
2 2
(
1
2
5 4
2
3 2 1
2
i j
N
j
j
j j
j
j j j iy
p p
k
b k h b
h
b k b h b w p
(
(
+ + + + =

=
(3-36)
where } 0 |) (| | { ) min( ' =
i j j i
w j for p p r r .
The use of
i
p' improves the stability of the algorithm by imposing the repulsive
forces between any pair of particles (Koshizuka et al., 1998). This feature is very
important in the algorithm since it reduces the possibility of local accumulation of fluid
particles at any time step and hence greatly improves the stability of the algorithm. The
stability improvement of the gradient model proposed in Eqs. (3-35) and (3-36) is very
significant in moving boundary problems, such as fluid flows with free surface.
The parameter
j
w in the above Laplacian and gradient model is the weighting
function used in Eq. (3-25) for the weighted least-square solution of the derivatives. As
discussed by researchers using GFD method for the same purpose, three weighting
functions are considered: quartic spline and cubic spline functions used by Gavete et al.
(2003) and inverse distance weighting coefficients proposed by Liszka and Orkisz
(1980). The explanation and comparison of these weighting functions will be presented
in Chapter 4.
3.3.3.2 Predictor-corrector steps
The same two-step algorithm introduced in Sect. 3.2.1 is used in the CPM method.
The algorithm is commonly used in the solution of the Navier-Stokes equation using
different numerical method, such as PFEM (Idelsohn et al., 2003) and incompressible
SPH (Shao and Lo, 2003).
58
3.3.3.3 Boundary surface recognition
In free surface flow problems, a key consideration is how to recognize free
surface boundaries particles. In mesh-based methods like FEM and FDM, the
boundary node labels are known as a priori. However, in particle methods, a new
particle distribution is obtained at each time step. The boundary particles should be
recognized before the boundary conditions can be implemented. Hence, in free surface
flow problems, the identification of free surface particles is of great importance in
obtaining the pressure field accurately. The MPS method uses density threshold in Eq.
(3-19) to differentiate the free surface particles. However, as discussed in Sect. 3.2.5,
this simple threshold method is not effective in recognizing free surface particles
accurately. In addition, the pressure solution is greatly affected by the parameter as
well as the influence radius
e
r defined in the algorithm. Since the pressure values of
free surface particles are theoretically zero and hence known, the total number of
unknowns in the Poisson equation of pressure excludes the numbers of free surface
particles. The density threshold method in recognition of free surface particle will
cause great pressure errors when some inner fluid particles are falsely recognized as
free surface particles.
Considering this problem, several boundary recognition algorithms may be used
to define the boundaries from a collection of particles. Edelsbrunner et al. (1983)
proposed a 2D alpha-shape method and later Edelsbrunner and Mucke (1994)
developed a 3D alpha-shape method. Idelsohn et al. (2003, 2004) and Idelsohn and
Onate (2006) used the alpha-shape method in their proposed PFEM. Dilts (2000)
introduced another boundary recognition method named arc method. It is tested to be
fast and can yield the exact solution. Detailed investigation of the efficiency of the arc
59
method can be found in a reference (Dilts, 2000). The arc method is used in this study
in the CPM simulations since it is easier to implement compared with the alpha-shape
method.
The basic idea of arc method can be illustrated in Figure 3-3(a). If any arc of the
circle (with radius R) around a center particle is not covered by any circle of its
neighbors, the center particle is treated as a free surface boundary particle. As shown
in the figure, the intersection arc of the center particle covered by the circle of
neighbor particle is recorded according to the position of neighbor particle in the form
of angle interval of the two interaction points. After that, all the angle intervals covered
by the neighbor particles are listed and sorted in ascending order by the left endpoints
(either clockwise or anti-clockwise). The angle intervals for the center particle A and B
are shown in Figure 3-3(b) for example. The sorted list is then scanned in sequence,
looking for gaps between two neighbor angle intervals. If any gap is found, the center
particle is considered as a free surface boundary particle. As can be seen in the figure,
particle A is recognized as a free surface particle and particle B in-domain fluid particle.
A more detailed discussion of the arc method is documented by Dilts (2000).
3.3.3.4 Incompressibility adjustment (IA) of the free surface particles
The incompressibility condition of fluid particles is imposed by Eq. (3-18), which
corrects the temporary particle number densities to the constant
0
n . The number
density of free surface particles is less than that of the in-domain fluid particles due to
the absence of particles above the free surface. Hence the correction of number density
to
0
n is not applicable for the free surface particles since the neighbor particle
distributions are truncated. The truncation depends on the radius of influence. Khayyer
and Gotoh (2009) noted this problem and proposed a solution method where a
60
modified source term is used in the right hand side of Eq. (3-18). They compute
Dt Dn / from the derivatives of kernel function w since number density n is a
function of w. But the source term in this formulation varies depending on the kernel
function used.
Here we proposed a scheme to impose incompressibility condition for the free
surface particles explicitly. Since the incompressibility condition is satisfied for all the
fluid particles, if the two free surface particles tend to get very close at some time
instant, their velocities are modified (as explained below) so that they will not collide.
This will greatly improve the stability of the numerical algorithm.
As shown in Figure 3-4, for the two free surface particles i and j , their velocities
are decomposed into the components normal and tangential to their coordinate
vector
i j
r r . Then we check whether their relative tangential velocity is greater than a
criterion shown below
( )
( )
0
L t
i j
i j t
j
t
i
>

r r
r r
v v (3-37)
where
t
i
v and
t
j
v are the tangential velocity components of the two particles and
0
L the
initial particle distance. The parameter is chosen around 0.1-0.3 with the selection of
time step t according to the condition in Eq. (3-11).
If the condition in Eq. (3-37) is satisfied for any pair of free surface particles, their
tangential velocity components will be adjusted as follows

2
' '
t
j
t
i
t
j
t
i
v v
v v
+
= = (3-38)
where '
t
i
v and '
t
j
v are the adjusted tangential velocity components of the two particles.
61
In this way, incompressibility of free surface particles is enforced. The
momentum sum of the two particles is still maintained. The proposed algorithm is
verified to be able to improve the distribution of the free surface particles so that
collision will not occur. It hence improves the stability of the algorithm.
3.3.3.5 Selection of neighbor particles
In particle methods, a particle interacts with its neighbor particles in a restricted
area within an influence radius
e
r (Figure 3-1). The influence radius is usually much
smaller than the solution domain size. For instance, Oger et al. (2006) have pointed out
that an influence radius covers 20 neighbor particles provide acceptable results in SPH
simulation of a 2D problem. The number of neighbor particles is a small fraction of the
total particle number in the solution domain. Identification of the neighbor particles for
a center particle requires the computation of distances between all pairs of particles. If
N is the total number of particles in a solution domain, the computational operations
for searching neighbor particles are of the order
2
N . The computational cost in
selection of neighbor particles is then dominant in the total computation time.
Therefore, an efficient neighborhood search algorithm is of great importance in fluid
problems with a large number of particles.
Viccione et al. (2008) compared two strategies of neighbor list generation, named
Verlet list method and cell-lined list method, and presented a sensitivity study on the
efficiency of the two methods. Koshizuka et al. (1998) proposed a candidate list
method, which is in principle the same as the Verlet list method in the work of
Viccione et al. (2008), to improve the efficiency of MPS method. The candidate list
method is shown to be efficient in decreasing the computational time in MPS code.
62
Here in this thesis we adopted the candidate list method of Koshizuka et al. (1998)
because of the simplicity of the approach.
The basic principle of the candidate list method is shown in Figure 3-5. For a
given problem with N particles, two lists are generated. One is for the neighbor
particles within the influence radius
e
r and the other is for the candidate particles
within the distance of h r
e
+ , where h is an additional distance. The candidate list is
updated everyk time steps by searching all N particles. The neighbor list is updated
every time step by searching all the candidate particles in candidate list. The
parameter k and h are selected properly depending on the motion of the problem so
that all neighbor particles within the influence radius should be captured in the
candidate list. Koshizuka et al. (1998) have shown that by proper selection of k , the
computational operations of list generation are reduced from the order of
2
N to the
order of
5 . 0 5 . 1
M N where M (<<N) is the number of neighbor particles.
3.3.3.6 Distinctive features of CPM
In summary, there are three distinctive features of the proposed CPM which
distinguish CPM from other particle methods.
Firstly, the method establishes discrete equations from derivatives obtained based
on Taylor series. The derivatives are solved from the function values of a collection of
particles in a weighted least-square approach. The Laplacian and gradient model are
established based on the derivatives. The improvement in accuracy of Laplacian model
in CPM will be presented in the following section (performance test).
Secondly, a more efficient method (arc method) for free surface particle
recognition is adopted in CPM. The arc method can produce more accurate boundary
63
particle recognition results compared with the density threshold method in MPS. The
performance of the arc method will be shown in the sloshing examples in Chapter 5.
Finally, the incompressibility condition of free surface particles is adjusted by an
IA scheme. A velocity modification scheme for the free surface particles is adopted to
impose incompressibility condition for them explicitly. The accuracy and stability
improvements of the arc method and IA scheme will be demonstrated in the following
two chapters through numerical examples.
The flowchart of CPM is shown in Figure 3-6, with the three distinctive features
highlighted.
3.3.4 Performance test of the Laplacian based on Taylor series
(1). Laplacian at a fixed point
Based on the algorithm described in Sect. 3.3.2, a simple performance test is
presented here to test the accuracy and convergence of the Laplacian model by Taylor
series solution. As shown in Figure 3-7, irregularly distributed points are generated by
adding random noise (uniform distribution) to the coordinates of 24 regularly spaced
points. The random positions are within a radius of 0.1L
0
around the original position
of regular grid. The operation of randomization is repeated for 50 times and thus 50
groups of irregular distribution are obtained.
For a given function ) , ( y x , the Laplacian of this function at the center point (dot
shown in Figure 3-7) can be evaluated by Eq. (3-16) in MPS method and Eq. (3-31) in
Taylor series method. Different discretization equations of Laplacian are proposed by
many researchers in SPH method (Shao and Lo, 2003; Lee et al., 2008). Here for
64
comparison a typical one shown in the following equation (Shao and Lo, 2003) is
investigated.

( )
+

+
=
|
|
.
|
\
|

i j
ij
ij i ij j i
j i
j
i
r
r w
m
2 2 2
) ( ) (
8 1
r
(3-39)
where
j
m is the mass of particlej ,
i
and
j
are the densities of particlei and j and
is usually equals to
e
r 1 . 0 where
e
r is the influence radius of particle i .
According to Eqs. (3-16) and (3-39), Laplacian models of MPS and SPH rely on
a kernel function ) (
ij
r w . In order to compare the performance of these two
discretization models, the same kernel function (Shao and Lo, 2003) is used as follows
( )
>
<
|
|
.
|
\
|

|
|
.
|
\
|
|
|
.
|
\
|
+
|
|
.
|
\
|
=
e
e e
e e
e
e e e
ij
r r
r r r
r
r
r
r r
r
r
r
r
r
r w
0
5 . 0 2 2
7
10
5 . 0 0 6 6 1
7
40
3
2
3 2
2
(3-40)
Therefore 50 corresponding values of Laplacian by each of the three methods can
be evaluated based on the parameters listed in Table 3-2. A test function
4 4
) , ( y x y x + = is selected in the study.
As shown in Table 3-2, since the tested function is
4 4
) , ( y x y x + = , the exact
value of Laplacian at the point with coordinate (1, 1) is 24. The computed Laplacian
values using the three methods are shown in Figure 3-8. The plot of Laplacian using
CPM is very close to the exact value. On the contrary, the results of MPS and SPH
method have prominent fluctuation. The solution using Taylor series in CPM is more
accurate than those of the MPS and SPH methods.
65
Besides accuracy, convergence is another performance which needs to be tested.
Generally speaking, for any numerical method, a finer discretization can help improve
the accuracy of simulation. For the same test function
4 4
) , ( y x y x + = , a convergence
study of Laplacian using the three methods is performed with different initial point
distance
0
L . The convergence test of the three methods in regularly distributed points
(without imposing random noise for the coordinates) is first studied with
0
L
decreasing from 0.7 to 0.001, as shown in Figure 3-9. One can see in the figure that all
the Laplacian value obtained from approximation schemes in the three methods
converge to 24 when
0
L becomes smaller. The convergence speed of MPS and SPH
seems better than CPM when regular points are used.
The convergence performance in irregular point distribution is then investigated
with 2 . 0
0
= L , 1 . 0
0
= L , 01 . 0
0
= L , and 001 . 0
0
= L . Small random noises are imposed
to the neighbor point coordinates. The convergence test results using MPS and SPH
are shown in Figure 3-10 and Figure 3-11 respectively. One may observe in the two
figures that the Laplacian approximation algorithms in the MPS and SPH methods
produce non-convergent solutions since a smaller particle size 001 . 0
0
= L results in
larger fluctuation of the Laplacian values. In contrast, Figure 3-12 shows the results of
convergence performance using the CPM method. The algorithm based on Taylor
series presents a good convergence trend: when the particle size is smaller, the plot
curve of Laplacian converges to the exact value. Therefore, we can conclude that the
approximation algorithms of Laplacian in MPS and SPH method produce convergent
solutions only when the particle distribution is regular. The convergence performance
is poor when the particles are irregularly distributed. In contrast, the approximation
66
algorithm in CPM presents good and convergent solutions even for irregular particle
distribution.
(2). Laplacian of a 2D surface
Instead of calculating Laplacian at a fixed point with coordinates (1, 1), a more
general case is studied here. The parameters used in this case are shown in Table 3-3.
A more complicated test function ) sin( ) , ( xy y x = is used in this case. The analytical
(exact) Lapalcian of the test function is ) sin( ) ( ) , (
2 2 2
xy y x y x + = , which is a
smoothed surface in the x-y domain.
In order to obtain numerical results of the Laplacian by the three methods, a set of
regular-grid points (34 x 34) in the domain is first generated. The irregular grid points
are then generated by imposing small random noise (random position within a radius
of 0.06
0
L around the regular point), as can be seen in Figure 3-13. Based on the point
positions and corresponding function values at each point, Eqs. (3-16), (3-31) and
(3-39) are used respectively in the three methods to evaluate the Laplacian at the points.
The comparison of three numerical solutions with the exact Laplacian solution in
the x-y domain ] 2 . 4 , 8 . 0 [ ] 2 . 4 , 8 . 0 [ is shown in Figure 3-14. It is shown that the
surfaces of Laplacian using MPS and SPH in general agree with the exact solution.
However, the surfaces are not smooth and some fluctuations of the Laplacian values
appear in certain points. While the surface using Eq. (3-31) in CPM is very smooth and
agree very well with the analytical solution.
In order to clearly show the differences of the numerical solutions with the exact
solution, the error (i.e. exact numerical ) surfaces of Laplacian are plotted in Figure
3-15. It can be seen that the error of CPM solution is much smaller than the errors of
67
SPH and MPS. The global error of the numerical results in Figure 3-15 is estimated by
the following equation and the values for MPS, SPH and CPM are 0.1323, 0.1041 and
0.01815, respectively.

2
2
2
2 2
) , (
) , ( ) , (
analytical i i
analytical i i numerical i i
y x
y x y x
err

= (3-41)
A series of particle distributions with random noise is then generated and the
global errors are plotted in Figure 3-16. One can see in the figure that the global error
of MPS and SPH approximation of Laplacian values are about 13% and 10% while
that of CPM is only about 1.8%. It can be concluded that CPM significantly out-
performs SPH and MPS in Laplacian approximation for irregular particle distribution.
3.4 Concluding remarks
In this chapter, a promising particle method based on Lagrangian formulation
named MPS is first studied. The method is based on solving a Poisson equation of
pressure which is derived from Navier-Stokes equation (Koshizuka and Oka, 1996). A
two-step algorithm is used to solve the coupled governing equations. At the predictor
step, temporary particle velocity and position are computed without considering the
pressure gradient. At the corrector step, pressure is computed by solving the Poisson
equation of pressure. As a driving source, the pressure gradient is then computed and
used to update the particles velocities and positions.
Nevertheless, in MPS (and also some other particle methods such as SPH), a pre-
defined kernel function is used to compute gradient and Laplace operators. Since the
kernel function is constructed artificially, an accurate solution cannot be guaranteed,
especially when the distribution of particles is irregular. Another reason which
deteriorates the performance of MPS method is that the free surface particles are not
68
accurately recognized at every time steps. Since the free surface particles are detected
by the number density, some in-domain fluid particles may be considered as free
surface particles, particularly when the fluid motion is violent.
In view of the shortcomings of the conventional particle method such as MPS, a
more robust particle method CPM is proposed. There are three key features in the
proposed method. Firstly, the method establishes discrete equations from derivatives
obtained consistent from Taylor series. The performance of the Laplacian model is
verified by a numerical test, and comparison with SPH and MPS method is made. It is
shown in the test that the approximation of Laplacian in CPM is significantly more
accurate than MPS and SPH. Furthermore, the algorithm in CPM shows better
convergence performance as the particle distance become smaller.
Secondly, a boundary recognition algorithm named arc method is adopted to
overcome the shortcomings of the density threshold method used in MPS. The density
threshold method cause great pressure errors when some inner fluid particles are
falsely recognized as free surface particles (Khayyer and Gotoh, 2009). With the
proposed new boundary recognition method, the problem is addressed.
And finally, the incompressibility conditions of the free surface particles are
adjusted by a simple scheme. The relative velocity components of any pair of neighbor
free surface particles are checked and adjusted if necessary to make sure they will not
collide. In this scheme, the momentum conservation of the two particles is still
maintained.
The accuracy and stability performance of the proposed CPM scheme will be
demonstrated in the following chapters through different free surface flow problems. A
69
parametric study of some numerical parameters in CPM such as influence radius and
particle sizes will also be carried out in the next chapter.

70

71
Table 3-1. Summary and comparison of particle methods
Particle
Methods
Require mesh
(back-
ground)?
Require
pre-defined
kernel?
Applied to
Free surface
flow?
Incompressi-
bility
Free surface
recognition
PFEM - Solve PPE*
Alpha-shape
method
SPH -
Slightly
compressible
State equation
( p )
MPS - Solve PPE
Number density
threshold
GFD - - - Solve PPE -
Finite
particle
- -
Slightly
compressible
State equation
( p )
Finite point
solver
- - -
Pseudo-
compressible
-
Finite point
method
- - - Solve PPE -
Finite
pointset
- - - Solve PPE -
CPM - - Solve PPE Arc method
PPE*: Poisson pressure equation (or Poisson equation of pressure)
Table 3-2. Parameters for performance test of Laplacian
Interaction radius
e
r
Initial distance
0
L
Test function
( ) , x y
Kernel function
( ) ,
e
w r r
0
1 . 2 L 0.1
4 4
y x + Eq. (3-40)

Table 3-3. Parameters for performance test of Laplacian
Interaction radius
e
r
Initial distance
0
L
Test function
( ) , x y
Kernel function
( ) ,
e
w r r
0
1 . 2 L 0.1 ) sin(xy Eq. (3-40)

72

73

Figure 3-1. Schematic of a typical reference particle with its neighbor particles

Figure 3-2. Algorithm of the MPS method

Initial configuration
End
Compute intermediate velocities
& particle positions
Compute particle number density
Solve Poisson equation of pressure
Update particle velocity and position
Continue?
No
Yes
74

(a)

Center Particle A Center Particle B
Neighbor
particle
Range of
angle covered
(clockwise)
Overlapped
with previous
arc?
Neighbor
particle
Range of
angle covered
(clockwise)
Overlapped
with previous
arc?
1 (0, 0.54) N 1 (0, 0.52) Y
2 (0.42, 0.97) Y 2 (0.27, 0.84) Y
3 (0.81, 1.37) Y 3 (0.64, 1.20) Y
4 (1.03, 1.58) Y

5 (1.35, 1.95) Y
6 (1.59, 2.15) Y
Found gap? Y (1.37~2) Found gap? N
Particle A: Free surface particle Particle B: In-domain particle
(b)
Figure 3-3. (a) The arc boundary recognition method; (b) Example of angle list for particle A and B

Figure 3-4. Incompressibility adjustment of free surface particles
75

Figure 3-5. Candidate list and neighbor list generation
76

Figure 3-6. Flowchart of CPM
Initial configuration
0 0
,v r
Compute temporal values
* *
,v r
Selection of neighbor particles
Compute particle number density
*
n
Free surface recognition by the
arc method (Sect. 3.3.3.3)

Solve Poisson equation of pressure
1 + k
p
Compute pressure gradient
1 +
k
p
Solve local matrix to get coefficients for Laplacian
and gradient (Sect. 3.3.2 and Sect. 3.3.3.1)
Yes
Yes
No
No
Check of
termination
Check of
output
Update velocities of next step
1 + k
v

Incompressibility adjustment of free
surface particles (Sect. 3.3.3.4)

Update positions of next step
1 + k
r

Output
End
77

Figure 3-7. A center point surrounded by 24 irregularly spaced points

Figure 3-8. Comparison of Laplacian using CPM (Eq. (3-31)), MPS (Eq. (3-16)) and SPH (Eq. (3-39))
(L
0
=0.1, Test function
4 4
) , ( y x y x + = )

Figure 3-9. Convergence test of Laplacian for regular points (Test function
4 4
) , ( y x y x + = )
78

Figure 3-10. Convergence test of Laplacian in MPS for irregular points (Test function
4 4
) , ( y x y x + = )

Figure 3-11. Convergence test of Laplacian in SPH for irregular points (Test function
4 4
) , ( y x y x + = )

Figure 3-12. Convergence test of Laplacian in CPM for irregular points (Test function
4 4
) , ( y x y x + = )
79

Figure 3-13. An example of irregular nodes (1156 in total) in x-y domain [0.8, 4.2] x [0.8, 4.2]

80

Figure 3-14. Analytical and numerical result of Laplacian (Test function ) sin( ) , ( xy y x = )
81

Figure 3-15. Difference of Laplacian values with analytical result (Test function ) sin( ) , ( xy y x = )
82

Figure 3-16. Error of Laplacian of different numerical algorithm (Test function ) sin( ) , ( xy y x = )

Chapter 4 Numerical Simulation of Incompressible Free Surface Flows by CPM
83
Chapter 4 Numerical Simulation of Incompressible
Free Surface Flows by CPM
4.1 Introduction
The CPM formulation presented in Chapter 3 is applied in numerical simulation
of free surface flow problems in this chapter. Two benchmark examples, i.e. pressure
computation in static tank and dam breaks in tank, are studied using CPM.
Comparison of CPM simulation results with MPS solutions is made in terms of free
surface motion and pressure contours. A parametric study is performed to investigate
the effects of different parameters used in CPM, such as the influence radius, least-
square weighting function and initial particle distance. In CPM, the discretization of
the governing equation is realized by the derivatives which are solved from the
equation based on Taylor series. The weighting function used in the weighted least-
square solution of the derivatives plays an important role in the stability and accuracy
of the proposed method. Three types of weighting functions are investigated. An
investigation of computational cost is performed following the parametric study.
Thereafter, different types of free surface flows are investigated and simulated using
the proposed CPM. More numerical examples are presented to show the capability of
the numerical algorithm and compare the results with published numerical and/or
experimental results. Both the MPS and CPM results are generated by the author
except some results with given references. All the numerical examples studied in this
chapter are summarized in Table 4-1.
84
4.2 Benchmark examples
4.2.1 Hydrostatic pressure in a static tank
Water in a static rectangular tank is the simplest free surface problem that can be
used to validate the performance of numerical algorithm since the velocities and
hydrostatic pressures of all the water particles are known. Starting from assumed
initial particle distribution, the two-step dynamic solution algorithm proposed in
Chapter 3 is used to simulate the variation of the particles in a static tank without
excitation. The velocities of the particles should approach to zero and the pressure of
particles is constant hydrostatic pressure when static equilibrium is achieved. In order
to check the accuracy of the pressure variation in the MPS method and CPM, the
pressure value at a fixed point is computed. The case was first studied by Khayyer and
Gotoh (2008) to show the ability of their proposed corrected MPS algorithm. Their
results are presented here for comparison.
Figure 4-1 shows a schematic view of the problem and the particle distribution in
the initial position. The dark color solid particles represent the fixed particles (solid
wall and dummy particles) where velocity is zero. The light color hollow particles are
the fluid particles. The water depth considered is 0.3 m. The pressure at a particle
(point A in the figure) located at the bottom of the tank is recorded. The initial particle
distance is 0.015 m and totally 41x 20 fluid particles are used in the simulation.
The time history of hydrostatic pressure at point A obtained by the MPS method
is shown in Figure 4-2. It can be seen that there is a significant pressure fluctuation
even in a static water tank. The fluctuation of pressure is so significant that the real
pressure history cannot be observed. The same observation and findings has been
85
reported by Khayyer and Gotoh (2008) for the original MPS method with water depth
0.2 m. The time history of hydrostatic pressure of their results is shown in Figure 4-3.
They proposed a corrected MPS algorithm (CMPS) to ensure the momentum
conservation of particles and showed that the CMPS reduced the pressure fluctuation,
as shown in Figure 4-3. In their more recent work (Khayyer and Gotoh, 2009), the
algorithm is further revised by introducing a new source term (MPS-HS) for the
Poisson pressure equation and allowing slight compressibility (MPS-WC). Figure 4-4
shows the time history of hydrostatic pressure at point A with these improvements in
the work of Khayyer and Gotoh (2009). Although the fluctuation of pressure history is
reduced with the combination of these improvements (Figure 4-4), the pressure
fluctuation does not vanish especially when large amplitude particle motion appears.
The hydrostatic pressure in the static tank is studied by CPM and the time history
of hydrostatic pressure at point A is shown in Figure 4-5. The dashed line represents
the pressure solution of CPM. It can be seen in the figure that the solutions using
CPM agree very well with the theoretical solution. In contrast, large pressure
fluctuation is observed in the original MPS method. The pressure fluctuation problem
is eliminated in CPM without resorting to relaxing incompressibility condition or
introducing any artificial source term.
Figure 4-6 shows the particle positions of the fluid domain at t=5 s using MPS
method and CPM. The fluid particles simulated by the MPS method show uneven
distribution especially near the free surface. This is due to its inability to satisfy the
incompressibility condition for the free surface particles.
Figure 4-7 shows the comparison of the hydrostatic pressure field at the same
time instant. The MPS method presents incorrect pressure field where some particles
86
with zero pressure (in dashed circle) appear inside the fluid domain. In the CPM
simulation, the pressure filed is generally smooth and presents a good solution of
hydrostatic pressure which increases linearly with respect to the particle depth.
4.2.2 Dam break with 2 / =
w
L d
A schematic view of the problem and initial particle distribution is shown in
Figure 4-8. It is a 2D simplification of water dam. The water located on the left is kept
still initially by the dam as a barrier. Once the barrier is suddenly removed at t=0 s,
the water collapses and starts to flow along the tank bottom. The example was solved
by Koshizuka and Oka (1996) both experimentally and numerically by the MPS
method. Later it became a benchmark problem to test the validity of Lagrangian
formulation in fluid flows. In this section, the results obtained using the MPS method
proposed in 2D domain are firstly presented. Viscosity and surface tension effects are
neglected.
The geometry ratios of the initial water depth over width and tank length are
2 / =
w
L d and 5 . 0 / =
T
L d , where d is the water depth and
w
L ,
T
L are the water width
and tank length. In this example, the water column is represented by 18 x 36 particles
with initial particle distance 0.008 m and the time step is 0.001 s.
Figure 4-9 shows MPS solution of the particle positions at different time instants.
Water runs onto the bottom surface until near 0.3 s when it impinges on the right
vertical wall. Breaking waves appear at about 0.6 s. As can be seen in the figure, the
MPS solution is in very good agreement with the experimental results in the change of
free surface profile with time.
87
Although the flow profile of dam break process is well simulated by MPS
method, the pressure field is not. Figure 4-10 shows that the pressure field is not
smooth at some time instances. Particles with zero pressure are detected not only on
the free surface but also (incorrectly) inside the fluid domain. For example, as we can
see in Figure 4-10, some particles inside the fluid domain has low pressure at time
t=0.2 s, 0.5 s and 0.9 s (dark particles in dashed circle) and are recognized falsely as
free surface particles for which pressure is set to be zero according to the criterion in
Eq. (3-19). Hence, the pressure contours at these time instants are not accurate any
more. Furthermore, by checking the pressure values of all the fluid particles, we found
that negative pressure value appeared at some time instant. This causes large pressure
fluctuation inside the fluid domain.
In contrast, the particle distribution and pressure solutions using CPM are shown
in Figure 4-11 at various time instants. The agreement with the experimental results in
Figure 4-9 both in the shape of the free surface and time development is excellent.
The pressure contours using the CPM method are shown to be much smoother than
the MPS solution in Figure 4-10. No zero pressure solution appears inside the fluid
domain. This mitigates pressure fluctuation inside the fluid domain.
4.3 Parametric study of CPM
There are several parameters in CPM which affect the stability and accuracy of
the simulation results. The first one is the weighting function used in the weighted
least-square solution of the derivatives based on Taylor series. Different least-square
weighting functions have been proposed (Gavete et al., 2003) and applied in the
simulation of problems such as solid mechanics and steady flows. The influence of
the weighting functions in free surface flows used in CPM is studied in this section.
88
Three different weighting functions are analyzed in the study of a simple dam
collapse example. Other parameters such as particle size (or initial particle distance)
and influence radius are also important in the proposed scheme. The effects of these
parameters will also be investigated in this section through the dam break example in
Sect. 4.2.2.
4.3.1 Influence of weighting functions in weighted least-square solution
In the solution of the derivatives based on Taylor series, a weighted least-square
approach is used. The weighting function used in Eq. (3-25) mainly helps to reduce
the distance effect of the particles since Taylor series expansion assumes the particles
involved are close. When the distance of the particles used in the Taylor series
expansion is too large, the accuracy of the second-order approximation decreases.
Hence, it is reasonable to give less weight for particles further away in the least-
square solution of Eq. (3-26).
As discussed by researchers who use weighting function in GFD method for the
same purpose, the stability of the algorithm is affected by the weighting coefficients
in the weighted least-square solution of the derivatives. We will investigate three
different weighting functions in CPM: the quartic spline and cubic spline functions
used by Gavete et al. (2003), and inverse distance weighting coefficients proposed by
Liszka and Orkisz (1980). They are listed as below.
Quartic Spline (Gavete et al., 2003)

>
+
=
e
e
e e e
r r
r r
r
r
r
r
r
r
r w
0
) ( 3 ) ( 8 ) ( 6 1
) (
4 3 2
(4-1)
Cubic Spline (Gavete et al. 2003)
89

>
< +
< +
=
e
e e
e e e
e
e e
r r
r r r
r
r
r
r
r
r
r r
r
r
r
r
r w
0
5 . 0 ) (
3
4
) ( 4 ) ( 4
3
4
5 . 0 0 ) ( 4 ) ( 4
3
2
) (
3 2
3 2
(4-2)
Inverse distance (Liszka and Orkisz, 1980)

>
=
e
e
r r
r r
r
r w
0
1
) (
3
(4-3)
All the three weighting functions are summarized in Table 4-2 and compared in
Figure 4-12. The effects of the three different weighting functions in CPM are
investigated using the same dam break example as in Sect. 4.2.2. A total of 18 x 36
fluid particles are used in the model. The initial particle distance is
0
L =0.008 mand
the radius of influence is selected as
0
5 . 2 L r
e
= .
The CPM solutions using the three least-square weighting functions are shown in
Figure 4-13. Agreeable wave profiles and smooth pressure contours can be obtained
using all these weighting functions. It is shown that the water profiles using cubic
spline weight function progress slightly faster than the experimental results, while the
solutions using quartic spline and inverse distance weighting functions agree well
with the experiments.
A sensitivity analysis is performed on the weighting functions, where the
proportion of number of free surface particles to the total fluid particles (N
f
) is plotted
in Figure 4-14. It can be seen that both the cubic spline and the inverse distance
weighting functions can produce stable results. The proportion of free surface
particles increases from less than 10% in the beginning of dam break to about 35%
90
when breaking wave occurs and gradually becomes stable at about 10%. When the
quartic spline weighting function is used, the number of free surface particles at some
time instants from 0.6 s to 1 s is abnormally high (up to more than 60% of total fluid
particles). This means that high fluctuation occurs so that most of the fluid particles
are recognized with zero pressure value. However, the fluctuation is suppressed after
2 s and the method still can run until the liquid calms down to be a tank of static
liquid.
By comparing the results of three weighting functions in Figure 4-13 and Figure
4-14, we found that the inverse distance weighting function is preferred. The time
history of the total number of free surface particle is smooth even when water hits the
tank wall and breaking occurs. The reason may be that the inverse distance weighting
function preserves a sharper increase when the particle distance decreases, while the
weights computed with quartic and cubic spline weighting functions increase much
more slowly when particle distance decreases (see Figure 4-12). Hence, very close
particle distance or even particle overlap is less likely to happen when the inverse
distance weighting function is used. In the subsequent numerical examples, the
inverse distance weighting function is therefore used.
4.3.2 Influence of influence radius
Influence radius
e
r is another important parameter used in CPM. It cannot be too
small, otherwise not enough particles are included to obtain the solutions of
derivatives. Since in Eq. (3-27) fives unknowns are to be determined for the center
particle, at least 5 neighbor particles are need. However, the distribution of the five
particles may not be satisfactory to avoid ill-conditioned matrix. Therefore more than
five particles are included for a stable and convergent solution. The influence radius
91
cannot be too big either. One reason is that it will increase the computational cost in
solving the derivatives from Eq. (3-27). Another important reason is that if
e
r is too
large, particles far away are unnecessarily included in the truncated Taylor series
expansion which assumes small distance from the center particle.
A sensitivity analysis is performed on the influence radius and the results of the
number of free surface particles are shown in Figure 4-15. The minimum influence
radius for regularly distributed particles should be
0
2L because it guarantees
minimum particle numbers needed in the solution of the five unknowns. Hence in this
study we investigate the different influence radii from
0
5 . 1 L to
0
1 . 3 L . It can be seen
in the figure that, numerical instability occurs when the influence radius is
0
5 . 1 L . In
this case, the number of free surface particles increases suddenly to a very high
number at about 0.2 s when the fluid has not hit the wall yet. In contrast, the number
of free surface particles is stable for influence radius from
0
6 . 1 L to
0
1 . 3 L .
Figure 4-16 indicates the location of leading edge of the collapsed water with
time using different influence radius. Experimental results presented by Hirt and
Nichols (1981) and published numerical solution using Lagrangian FEM by
Ramaswamy and Kawahara (1987) are also plotted in the same figure for comparison.
It can be seen that, the CPM solution for leading edge of water using influence radius
0
6 . 1 L is slightly slower than the experimental result, while the solution using
0
1 . 3 L
is slightly faster than the experimental result. The CPM solutions of leading edge
using influence radii from
0
8 . 1 L to
0
5 . 2 L are in good agreement with the experimental
results by Hirt and Nichols (1981). Hence, the influence radius needs not be too large
92
due to the reasons mentioned before. In the later numerical examples in this thesis,
either
0
1 . 2 L or
0
5 . 2 L will be used.
4.3.3 Influence of particle sizes
In particle methods, the solution domain is represented by particles. The initial
particle distance or the so-called particle size, is another parameter that needs to be
studied. The particle size is, to some extent, the same as the mesh size in mesh-based
method. The particle size should not be too large since too few fluid particles will
result in an inaccurate representation of liquid motion in the simulation. In theory, the
more particles used in representing the domain, the better accuracy of simulation
results may be obtained. However, the computational cost increases correspondingly.
Therefore, a reasonable particle size should be selected depending on the severity of
liquid motion.
In this section, a sensitivity study on the particle size is performed in CPM in
terms of the number of free surface particles. As shown in Figure 4-17, five different
particle sizes are used in the simulation, which are
0
L =0.004 m, 0.006 m, 0.008 m,
0.012 m and 0.016 m respectively. It is expected that decreasing the particle size will
result in an increase of the total number of the fluid particles. Hence the number of
free surface particles is normalized with respect to the total fluid particles. It is shown
in Figure 4-17 that, all the four particle sizes used in the simulation generate stable
results. The variations of proportions of free surface particles with time using particle
size from 0.004 m to 0.008 m show similar results. When particle size 0.012 m and
0.016 m are used, the proportions of free surface particles are higher than other cases.
This means that large particle size presents inaccurate representation of the liquid
motion, with large portion of particles recognized as free surface.
93
The profiles of dam break process using four different particle sizes are shown in
Figure 4-18 for different time instants. In the left-most column, solution based on the
particle size 0.016 m is too large to capture the detailed wave profile. It is shown that
the liquid motion is well represented using the particle size 004 . 0
0
= L m. The profile
of liquid motion is very similar to the solution based on a smaller particle size L
0
=
0.003 m.
We can conclude that for the dam break case in this study, the particle
size 004 . 0
0
= L m can be selected as a reasonable value. For other cases with different
geometries and liquid motion types, a similar convergence study performed here is
needed to help make a choice of suitable particle size. A relatively large particle size
should be used first to examine the overall liquid motion. A refined particle size will
then be used to model the details of changing free surface profile.
4.3.4 Influence of time step
Similar to the MPS method, a two-step predictor-corrector algorithm is used to
solve the governing equation. The time step is controlled by the CFL condition shown
in Eq. 3-11. In this section, the choice of time step is studied and discussed by using
the same dam break example.
The particle size 004 . 0
0
= L m is selected here according to previous section.
The maximum fluid velocity V
max
can be estimated as (Lee et al., 2008)
gH V 2
max
=
where H is the initial water height. From Eq. 3-11 we can obtain the maximum time
step is about t=0.00025 s.
94
Four different time steps are used in this study, which are 0.0002 s, 0.00025 s,
0.0005 s and 0.001 s. The dam break profiles using these time steps are shown in
Figure 4-19 for the three time instants. It can be seen that when larger time step are
used (0.0005 s and 0.001 s), the particle distributions at t=0.3 s are not well simulated.
The solutions diverge due to large particle motion at one time step. When the time
step is selected according to the CFL condition in Eq. 3-11, the simulation results are
very good. It is shown that the simulation results of water profiles using a smaller
time step t=0.00015 s are practically the same as those for t=0.00025 s. But the
computational cost using smaller time step is higher. Therefore the time step
satisfying the CFL condition in Eq. 3-11 is selected in this research.
4.3.5 Computational cost
Compared with MPS, the additional computational cost in CPM lies in the
solution of derivatives based on local particle coordinates (See Sect. 3.3.2 and Sect.
3.3.3.1). The computational cost incurred in this operation needs to be investigated.
Table 4-3 shows the CPU time per time step for the dam break example in Sect. 4.2.2
with different particle numbers (N). The computer used is equipped with Pentium 4
TM

CPU (duo 3.6 GHz). The CPU times of three main operations are listed: generation of
neighbor and candidate list (t
list
, see Sect. 3.3.3.5), weighted least-square solving of
derivatives (t
WLS
, see Sect. 3.3.2), solving Poisson equation of pressure (t
PPE
, see Sect.
3.3.3.1). The total time cost includes all operations in a time step (t
total
). The time
fraction of the MPS method for four particle numbers is also listed for comparison.
The CPU time and fractions in Table 4-3 are plotted in Figure 4-20 for better
presentation visualization of the results. It is shown that when the particle number is
small, i.e. N=1095, the most time consuming part (35%) of CPM is the solving of
95
derivatives based on Taylor series. The list generation operation only takes up 13%.
In MPS method, however, the list generation occupied 41% of the total time (Gotoh et
al., 2005). With the increase of particle numbers, the time used in least-square
solution of derivatives in CPM reduces from 35% to less than 10% while that of list
generation increases from 13% to about 30% of the total time. It can be concluded
that computational cost in the addition algorithm of solving derivatives based on local
particle coordinates in CPM is not significant in problems requiring large particle
number. The time used in solving Poisson equation of pressure in CPM increases with
increasing particle number and becomes the dominant part of the total computational
time when the particle number is large. This was also pointed out by Koshizuka et al.
(1998) in the case of MPS method.
4.4 Numerical simulation of free oscillation of liquid
After the parametric study performed in the previous section, the CPM can be
applied in simulation of various free surface flows. In general the motion of the free
surface fluid flows can be classified as fluid motion without wave breaking and with
breaking.
Many particle methods have been developed with emphasis on breaking waves,
such as SPH, FPM and MPS. However these methods frequently suffer from
difficulties in simulation of steady-state fluid motion even without breaking. The
pressure fluctuation problem in these particle methods makes the solution unstable if
long time simulation is performed, as demonstrated in Sect. 4.2.1 for MPS.
In this section, we will demonstrate the capability of CPM in simulation of fluid
motion without breaking, i.e. free oscillation of liquid in a container with small
96
amplitude. The schematic view of the container used for the analysis is illustrated in
Figure 4-21(a). A square liquid column of width b with initially a sinusoidal free
surface profile in a rectangular 2D tank is fixed on ground. The initial velocities of the
liquid particles are all zero. The analytical solution is available when the amplitude is
small (Radovitzky and Ortiz, 1998). The case was studied by many researchers
numerically (Ramaswamy et al., 1986; Ramaswamy, 1990; Radovitzky and Ortiz,
1998; Idelsohn et al., 2004). The initial free surface profile of the liquid column is as
follows
x
b
a x

sin ) 0 , ( = (4-4)
The initial particle distribution is shown in Figure 4-21(b). In this example we
use the same parameters as those adopted by Idelsohn et al. (2004). The initial
amplitude of free surface is a=0.1 unit and the square liquid column (in steady-state)
is of width b=1 unit. 50 x 53 fluid particles are used to model the water inside the tank.
Viscosity is assumed to be negligible.
The time variation of the free surface elevation at 2 / b x = is computed with
the MPS method and CPM. The comparison of the free surface elevation with
analytical results is plotted in Figure 4-22. The free surface elevation obtained using
CPM agrees well with the numerical solution of Idelsohn et al. (2004) using PFEM.
Compared with the analytical results, it is seen that both the two numerical solutions
exhibit some numerical viscosity. This has been also discussed by Idelsohn et al.
(2004). Furthermore, the solution of Idelsohn et al. (2004) using PFEM presents some
phase difference compared with the theoretical results, as shown in the figure. In
contrast, this phase difference in the solution using CPM method is smaller.
97
The solution of free surface elevation in the MPS method and CPM are
compared in Figure 4-23. The MPS solution exhibits fluctuation after t=0.5 s. After
t=2 s the fluctuation becomes so severe that the sinusoidal sloshing profile of free
surface cannot be clearly identified. This shows that MPS method fails to produce
long time simulation of the free oscillation of liquid.
Figure 4-24 compares the pressure fields at several time instants using MPS,
CPM and analytical solutions. Compared with MPS, the CPM solution is much better
in terms of free surface and agrees very well with analytical results. Before t=1.26 s
the MPS solution for free surface profile agrees with the CPM solution. However,
after that the free surface profiles in MPS solution become volatile, which makes
the recognition of free surface elevation difficult. The MPS solution of pressure
contours shows poor performance. There are particles with inaccurate pressure inside
the fluid domain in the MPS solution, whereas CPM gives smooth pressure solutions
correctly. The drawback of the MPS method in the pressure solution makes it difficult
to study the fluid motion for a long time simulation due to the pressure fluctuation
problem.
More cases of fluid motion without breaking will be studied in Chapter 5, which
focus on sloshing in moving rectangular tank.
4.5 Numerical simulation of violent fluid flows with breaking
Study of violent free surface flow with wave breaking is a challenging problem
in fluid dynamics. The fast changing free surface means that the geometry of the
solution domain changes with time and the free surface after some time becomes
discontinuous. The free surface may break in many parts, and some parts may join
98
together in one. The fragmentation and coalescence of fluid particles makes numerical
simulation difficult especially for conventional mesh-based methods. For this reason
many particle methods have been developed and applied in the study of such
problems. The most well known ones are SPH and MPS. As discussed before, in
conventional particle methods such as MPS, the solution of the pressure field is not
satisfactory.
In this section, we will demonstrate the ability of the proposed CPM in modeling
this kind of problem, especially with regards to improvement in pressure solutions.
4.5.1 Free oscillation of liquid in a container with large amplitude
The first case we considered is the same free oscillation of non-viscous water
column but with larger initial amplitude. Upon oscillation, the waves may break.
Linear analytical solution is not applicable for large wave amplitude.
The amplitude of the initial free surface is one third of the width of square water
column in the tank. The same problem is studied by Idelsohn et al. (2004) using
PFEM method. Their simulation result is presented here for comparison. In our
simulation 39 x 39 fluid particles are used.
Figure 4-25 compares the CPM solution of wave profiles at several time instants
with the PFEM solutions obtained by Idelsohn et al. (2004). The general free surface
profiles are in good agreement. Breaking particles are captured in both solutions.
Nevertheless, there is a large density variation in the PFEM solution of Idelsohn et al.
(2004) in Figure 4-25 (a). This large density variation makes differentiation between
breaking particles and non-breaking particles difficult and also affects the boundary
particle recognition. The error in identifying the free surface particles would, in turn,
99
affect the pressure solution in the Poisson equation of pressure, thereby further
deteriorating the accuracy of PFEM solution.
The CPM simulation solution is shown in Figure 4-25 (b). Breaking particles
are observed in the third picture. Compared with the PFEM solution, the CPM
solution does not present large density variation inside the fluid domain. As discussed
before, this helps improve the simulation accuracy.
The pressure contours at the respective time instants for CPM in Figure 4-25 are
shown in Figure 4-26. Again it is seen that the pressure distributions inside the fluid
are smooth with no abnormal pressure values in the tank. The pressure field is well
represented even when breaking particles appear.
In order to obtain more violent breaking phenomena, an even larger amplitude
free oscillation example is investigated. The initial amplitude of free surface is a=0.6
unit. A total of 27 x 27 fluid particles are used. The CPM solution of pressure fields at
different time instants is shown in Figure 4-27. It is shown that the larger initial
amplitude produce more violent breaking particles from t=0.42 s to 1.02 s. The
pressure contours are smooth even when the breaking wave occurs. Furthermore, the
free surface particles (with zero pressure) are accurately recognized in CPM. No
spurious pressure value in the fluid domain is observed.
4.5.2 Dam break with 5 . 0 / =
w
L d
To further demonstrate the performance of CPM in the pressure field solution for
free surface flow with breaking, a larger dam break example with 5 . 0 / =
w
L d and
186 . 0 22 . 3 / 6 . 0 / =
T
L d is considered. This dam break example is based on an
experiment performed by the Maritime Research Institute Netherlands (MARIN). A
100
schematic view of the example is shown in Figure 4-28. A large tank of 3.22 m in
length is used with an open roof. A volume of water with 0.6 m in height and 1.2 m in
length was released to flow at the initial time step.
Measurements were performed in the experimental study. As shown in Figure 4-
28, four vertical wave probes were used to measure the water elevations at positions
H1 to H4. One of the pressure sensors was located at position P1 to capture the
pressure value at that point. Details of the experimental study can be found in the
references (Fekken, 1998; Kleefsman et al., 2002).
In our numerical model, the initial particle distance is
0
L =0.01 m. The kernel
size
0
5 . 2 L r
e
= and time step t =0.001 s are adopted. Viscosity of the liquid is
neglected. The animations of CPM simulation of the dam breaking wave profiles,
pressure and vector fields are recorded in the attached CD with their respective file
names Dam_marin_w.avi, Dam_marin_p.avi and Dam_marin_v.avi.
The water profiles obtained in CPM simulation at different time instants are
shown in Figure 4-29. The water flows once the dam breaks at the initial time step
and hits the right wall of the tank in less than one second. The water runs up along the
right wall and then overturns and falls back (t=1.4 s) before hitting the water body in
the tank and generating a big splash (t=1.7 s to t=2.1 s). The splash falls into the tank
and the water gradually calms down after that.
Figure 4-30 shows the time history of the computed pressure at point P1 using
CPM superimposed on a scanned figure that presents the experimental result of
MARIN and numerical solution by Fekken (1998). In the figure, the dashed line
represents the experimental results. The pressure solution by CPM agrees reasonably
101
well with the experimental result before t=2 s. In terms of agreement with the
experimental results, the CPM solution is better than the numerical solution of Fekken
(1998) who used the software ComFlow developed by MARIN. From t=2 s to t=4 s
the pressure in CPM solution is higher than the experimental result.
There are two sharp pressure peaks in the pressure history, as shown in Figure 4-
30. The first one occurs at about t0.7 s when water from the collapsed dam hits the
right wall. The second and higher peak appears at t1.5 s when the overturned
breaking wave falls and hits the remaining water body which then exerts additional
pressure on the tank wall. Furthermore, there is a gentle peak at t3 s when the
overturned wave fully falls down to the tank bottom and induce a relatively large
undulation of water. The undulation hence causes a pressure increase at point P1 due
to the moderate water run-up around the wall. It is clear in the figure that, both these
three peaks are successfully captured in CPM. In contrast, the numerical result by
Fekken (1998) fails to capture the two peaks at t0.7 s and t3 s. The pressure
contours at the occurrence of these three pressure peaks are shown in Figure 4-31,
where high pressures around point P1 are observed.
Figure 4-30 also shows that there is some minor fluctuation of pressure value
after water hits the right wall and falls back to water body. It is caused by the impact
of breaking particles on wall or main water body. Some of the particles fly outside the
solution domain and later fall back to the domain with big velocities. These particles
will affect the local pressure distribution of the fluid particles when they fall down
and re-join the fluid. This may be mitigated with some special treatment of the flying
particles once they separate from the main solution domain. For example, a particle
slip method is adopted by Losasso et al. (2008) in which ballistic particle threshold is
102
maintained to differentiate solitary particles that are far from the main fluid body and
follow ballistic trajectories. Nevertheless in this study we do not focus on the breaking
particles and thus no special action is taken.
Figure 4-32 shows the water heights computed by CPM at the four measure
positions, and the results of MARIN and Fekken (1998) are again scanned for
comparison. It is shown that, the computed water heights at the four positions in
general agree well with the experimental results. The agreement is especially
excellent before and long after the water hits the tank wall. During the impact of water
on the wall or main water body, i.e. from t=1.5 s to t=2.2 s, some differences between
the CPM solution and the experimental results are observed at points H1 and H2. It
can be seen that the difference in the maximum height are less than 16% in CPM
solution, while in the result of Fekken (1998) the difference is up to 42% at point H2.
The agreement of our simulation results with experiments is thus much better than the
numerical solution of Fekken (1998).
The water profiles with pressure contours are presented in Figure 4-33 for some
representative time instants. It can be seen in the figure that the water hits the wall in
less than one second after release from the dam. Large pressure occurs during the
hitting of the water on the tank wall. Some water particles break away upon hitting
the wall. Some particles roll up along the wall and then overturn and fall back to the
main water body. The overturning wave generates high pressure in the fluid when
hitting the fluid particles, as shown in the figure for t=1.5 s. The pressure solution
approaches hydrostatic pressure field after t=2.5 s since then the water becomes gentle
and gradually calms down.
103
The pressure contours shown in the figures are smooth without any zero pressure
values inside the domain. The CPM is shown to be robust in simulation of the violent
fluid motion even when the fluid motions are complicated including possible wave
breaking. The capability of CPM shown in the numerical examples reveals the great
potential of the proposed particle method in the study and simulation of problems in
fluid dynamics.
The velocity fields after the dam break at different times in Figure 4-34 shows
the fast changing flow motion. At t=0.5 s the velocities of the particles at the left
corner are relatively small while those of the particles at the front of the water are
large. At t1.0 s when the water hits the right wall, the vertical velocity components
of the particles near the wall are large. From t=1.3 s to t=1.5 s overturning wave is
observed, for which the velocity vectors of the particles exhibit the direction of the
fluid motion. The overturning wave hits the main water body at t1.5 s and water
splash appears from t=1.7 s to t=2.1 s. Some local violent vortexes of the water are
observed at t=3.1 s and t=3.5 s from the velocity fields. The water gradually calms
down after t=3.5 s. The velocity results obtained prove the good performance of CPM
in the simulation of violent free surface problem with possible water breaking.
4.5.3 Dam break with obstacle
To further demonstrate the capability of the proposed CPM in simulation of the
violent breaking waves, another dam break example is presented in this section. A
more violent and complicated fluid motion is generated in this example, owing to the
presence of an obstacle in the tank. The dam break examples with and without
obstacle can be considered as a simple model of green water flow on the deck of a
ship (Fekken, 1998; Lee et al., 2002; Greco et al., 2004; Lohner et al., 2006).
104
The example was studied experimentally by Koshizuka et al. (1995). A
schematic view of the problem is shown in Figure 4-35. The tank geometry is the
same as the case studied in Sect. 4.2.2. In addition, an obstacle block with geometry
h h 2 , where 24 = h mm is used in the experiment, is located at the middle of the
tank. The same problem was also studied by Larese et al. (2008) using the PFEM
method. Their simulation results are presented here for comparison.
The initial particle distribution of water is shown in Figure 4-36. The initial
particle distance
0
L =0.0073 m and totally 36 x 18 fluid particles are used. The
influence radius is
0
5 . 2 L r
e
= and time step is 0005 . 0 = t . Viscosity of the liquid is
neglected. The obstacle block is modeled as part of the rigid tank wall.
Figure 4-37 shows the water profiles at several time instants after the dam breaks.
The left column shows the experimental results of the water profiles captured by
camera (Koshizuka et al., 1995). The middle column shows the numerical results by
Larese et al. (2008), who used the particle FEM method in their simulation. The right
column of Figure 4-37 presents the simulation results by the proposed CPM.
The experimental results show that the water collapses and, at about 0.1 s, water
hits the obstacle block, upon which the wave profile changes tremendously. The wave
runs up to a high elevation upon impact with the obstacle block. The big water splash
then hits the right wall and falls back to the right part of the tank. After that the water
flows to the left and gradually calms down.
One can see in the figure that, the simulation results using CPM agree well with
the experimental results and PFEM solutions until t=0.4 s. After that, the numerical
simulation results of the water profiles differ somewhat from the experiments. It is
105
shown in Figure 4-37 that, for the experimental results at t=0.5 s, the water particles
on the right side of the obstacle are higher than the obstacle. The numerical simulation
results of both CPM and PFEM show that the particles fall down towards the tank
bottom. The main reason is that there is air entrapped in the overturned waves in the
experiment. Since only one fluid is modeled in the numerical simulation, air
entrapment is not considered in either of the two numerical methods. However, in the
real experiments, air entrapped inside the overturned wave slow the fall of water
particles. Therefore a two-phase model is needed if numerical simulation is to be
improved.
Nevertheless, the CPM results agree with the experimental results well after the
main water body totally fall down to the tank. It can be seen that for t=1.0 s the CPM
solution agree with the experimental results much better than the PFEM solutions.
The animation of CPM solution for this case is recorded in the attached CD with file
name Dambreak_obstacle.avi.
In this chapter, six numerical examples are presented to examine the
performance of CPM in simulation of different free surface flow problems.
In the benchmark examples, the MPS method present volatile free surface
profiles and the solution of pressure field presents inaccurate result with zero pressure
inside the fluid domain. In contrast, the CPM solution presents much smoother and
accurate pressure contours without volatile free surface.
Different problems of free surface flows are investigated with the proposed CPM.
The pressure field is successfully presented in the CPM. No severe pressure
106
fluctuation appears in the simulation of both gentle fluid motions without breaking
and violent free surface flows when the fluid waves are complicated and with the
possibility of breaking. The CPM solutions agree very well with published numerical
and/or experimental results.
A parametric study is performed to investigate the effects of different
parameters in CPM using the dam break example. It is found that the inverse distance
weighting function is more suitable for the proposed CPM for better numerical
accuracy. The influence radius
0
1 . 2 L r
e
= or
0
5 . 2 L can be used in the simulation of
violent fluid motion. The initial particle distance should be chosen according to the
geometry and fluid motion types of the specified problem to be studied. The
computational cost of CPM is also investigated. It is found that the time used in the
addition algorithm of solving derivatives based on local particle coordinates in CPM
is not a significant portion when large particle numbers are used.

107
Table 4-1. Numerical examples studied in Chapter 4
Examples
Source of results and methods used for
comparison
Source of result Method
Hydrostatic pressure in a static tank
Khayyer and Gotoh
(2008, 2009)
Improved MPS
Dam break with 2 / =
w
L d
Koshizuka and Oka
(1996)
Experiment
Free oscillation of liquid with small
amplitude
Idelsohn et al. (2004) Analytical and PFEM
Free oscillation of liquid with large
amplitude
Idelsohn et al. (2004) PFEM
Dam break with 5 . 0 / =
w
L d Fekken (1998)
Finite volume method
+improved VOF
Dam break with obstacle
Koshizuka et al. (1995)
Larese et al. (2008)
Experiment
PFEM

Table 4-2. Parameters used in the study of the effect of least-square weighting function
Cases 1 2 3
Weighting function Eq. (4-1) Eq. (4-2) Eq. (4-3)
e
r
2.5 L
0
2.5 L
0
2.5 L
0

0
L
0.008 m 0.008 m 0.008 m

Table 4-3. CPU time (s) and fraction per time step in CPM
N t
List
(fraction) t
WLS
(fraction) t
PPE
(fraction) t
Total
(fraction)
CPM MPS CPM MPS CPM MPS CPM MPS
1095
0.052
(0.13)
(0.41)
0.135
(0.35)
/
0.098
(0.25)
(0.07)
0.389
(1)
(1)
1743
0.310
(0.28)
(0.31)
0.236
(0.22)
/
0.338
(0.31)
(0.16)
1.094
(1)
(1)
3471
0.987
(0.30)
(0.30)
0.593
(0.18)
/
1.203
(0.37)
(0.33)
3.275
(1)
(1)
5775
2.960
(0.29)
(0.28)
0.970
(0.09)
/
5.350
(0.52)
(0.49)
10.266
(1)
(1)
9000
9.500
(0.28)

3.050
(0.09)
/
20.01
(0.58)

34.322
(1)

10730
13.013
(0.26)

4.508
(0.09)
/
29.595
(0.60)

49.530
(1)

108

109

Figure 4-1. Schematic view of the initial particle distribution for static water tank
0
2000
4000
6000
8000
10000
0.0 1.0 2.0 3.0 4.0 5.0
Time (s)
P
r
e
s
s
u
r
e

(
P
a
)Original MPS
Theoretical solution

Figure 4-2. Time history of hydrostatic pressure at point A by MPS

Figure 4-3. Time history of hydrostatic pressure at point A with d=0.2 m scanned from Khayyer and
Gotoh (2008)
0.3 m
Point A
110

Figure 4-4. Time history of hydrostatic pressure at point A with d=0.2 m scanned from Khayyer and
Gotoh (2009)

2700
2900
3100
3300
3500
0.0 1.0 2.0 3.0 4.0 5.0
Time (s)
P
r
e
s
s
u
r
e

(
P
a
)Theoretical solution
CPM

Figure 4-5. Comparison of time histories of hydrostatic pressure at point A

111

Figure 4-6. Particle distribution at t=5 s.

Figure 4-7. Hydrostatic pressure field of the whole tank of water at t=5 s

Figure 4-8. Geometry and initial particle distribution of the dam break example
Original MPS CPM
Original MPS CPM
Zero
pressure
112

Figure 4-9. Comparison of dam break simulation using MPS with experimental results
Left two Columns: numerical results by MPS; Right Column: experiments (Koshizuka and Oka, 1996)

t=0.1 s t=0.2 s
t=0.2 s
t=0.3 s t=0.4 s
t=0.4 s
t=0.5 s t=0.6 s
t=0.6 s
t=0.7 s t=0.8 s
t=0.8 s
113

Figure 4-10. Pressure field of the dam break example by MPS method

Figure 4-11. Pressure field of the dam break example by CPM
t=0.1 s t=0.2 s t=0.3 s
t=0.4 s t=0.5 s t=0.6 s
t=0.7 s t=0.8 s t=0.9 s
t=0.1 s t=0.2 s t=0.3 s
t=0.4 s t=0.5 s t=0.6 s
t=0.7 s t=0.8 s t=0.9 s
Zero
pressure
Zero
pressure
Zero
pressure
114
r/r
0
1
2
3
4
5
0 0.2 0.4 0.6 0.8 1 1.2 1.4
e
w
(
r
)
Quartic Spline
Cubic Spline
Inverse distance

Figure 4-12. Comparison between different weighting function

Figure 4-13. Comparison of CPM solutions using different weighting functions with experiments
Quartic Spline
t=0.6 s
Cubic Spline
t=0.6 s
Inverse distance
t=0.6 s
Quartic Spline
t=0.8 s
Cubic Spline
t=0.8 s
Inverse distance
t=0.8 s
115
f
0
0.2
0.4
0.6
0.8
0.0 1.0 2.0 3.0 4.0
Ti me (s)
N
o
.

o
f

f
r
e
e

s
u
r
f
a
c
e

p
a
r
t
i
c
l
e
s
/
N
Quartic Spline
Cubic Spline
Inverse distance

Figure 4-14. Effect of weighting functions used in CPM
e 0
e 0
e 0
e 0
e 0
e 0
0
50
100
150
200
250
300
350
400
450
500
550
600
650
0.00 0.25 0.50 0.75 1.00
Ti me (s)
N
o
.

o
f

f
r
e
e

s
u
r
f
a
c
e

p
a
r
t
i
c
l
e
s
r =1.5 L
r =1.6 L
r =1.8 L
r =2.1 L
r =2.5 L
r =3.1 L

Figure 4-15. Effect of influence radius used in CPM
116

Figure 4-16. Comparison of the leading edge location in CPM solution with published results---
experimental: Hirt and Nichols (1981); Lagrangian FEM: Ramaswamy and Kawahara (1987)
f
0
0
0
0
0
0
0.2
0.4
0.6
0.8
0.0 0.5 1.0 1.5 2.0
Ti me (s)
N
o
.

o
f

f
r
e
e

s
u
r
f
a
c
e

p
a
r
t
i
c
l
e
s
/
N
L =4 mm
L =6 mm
L =8 mm
L =12 mm
L =16 mm

Figure 4-17. Effect of particle size used in CPM

117
(a) L
0
=0.016 m (b) L
0
=0.008 m (c) L
0
=0.006 m (d) L
0
=0.004 m (e) L
0
=0.003 m

Figure 4-18. Dam break profiles using different particle sizes in CPM

t=0.0002 s t=0.00025 s t=0.0005 s t=0.001 s
t=0.1 s
t=0.2 s
t=0.3 s
Figure 4-19. Dam break profiles using different time step t

t=0.4 s
t=0.7 s
t=0.8 s
t=1.0 s
118
List WLS PPE
0
10
20
30
40
50
1095 1743 3471 5775 9000 10730
Total parti cl e numbers (N)
C
P
U

t
i
m
e

(
s
)
t t t

(a)

List WLS PPE
0%
20%
40%
60%
80%
100%
1095 1743 3471 5775 9000 10730
T
i
m
e

f
r
a
c
t
i
o
n
t t t
List WLS PPE
0%
20%
40%
60%
80%
100%
1095 1743 3471 5775
T
i
m
e

f
r
a
c
t
i
o
n
t t t

(b) (c)
Figure 4-20. CPU time for different operations (a) CPU time of CPM; (b) Fraction over the total time
of CPM; (c) Fraction over the total time of MPS

(a) (b)
Figure 4-21. (a) A schematic view of the tank; (b) Initial particle distribution

x
y
b
b
CPM
CPM
MPS
119

Figure 4-22. Comparison of time histories of surface elevation amplitude with published results

Figure 4-23. Comparison of time histories of surface elevation amplitude results by CPM and MPS

120
MPS CPM Analytical
t=0.18 s

t=1.26 s

t=4.32 s

t=4.86 s

Figure 4-24. Pressure fields at different time instants for MPS and CPM simulation
121
(a)
(b)
Figure 4-25. Comparison of the free oscillation of liquid for large amplitude
(a) PFEM results by Idelsohn et al. (2004); (b) results by CPM

Figure 4-26. Pressure contours of the free oscillation of liquid by CPM

t=0 s t=0.6 s t=1.32 s t=1.88 s
t=0 s t=0.6 s t=1.32 s t=1.88 s
t=0 s t=0.6 s t=1.32 s t=1.88 s
Breaking
particles
122

Figure 4-27. Pressure contours of the free oscillation of liquid for larger amplitude by CPM

Figure 4-28. Configuration of the tank in dam break example and the positions of the water depth
probes and pressure sensor (Fekken, 1998)

t=0 s t=0.21 s t=0.42 s t=0.63 s
t=0.84 s t=1.02 s t=1.50 s t=1.86 s
123

Figure 4-29. Wave profiles of the dam break example by CPM
t=0.5 s t=1.0 s
t=1.3 s t=1.4 s
t=1.5 s t=1.7 s
t=2.0 s t=2.1 s
t=2.5 s t=3.1 s
t=3.5 s t=4.0 s
124

Figure 4-30. Comparison of pressure at Point P1 with published results (Fekken, 1998)

Figure 4-31. Pressure contour of dam break by CPM at (a) t=0.7 s, (b) t=1.475 s and (c) t=2.95 s
t=0.7 s
t=1.475 s
t=2.95 s
(a)
(b)
(c)
125

Figure 4-32. Comparison of water heights at the four points with published results (Fekken, 1998)

126

Figure 4-33. Pressure contours of the dam break example by CPM
t=0.5 s t=1.0 s
t=1.3 s t=1.4 s
t=1.5 s t=1.7 s
t=2.0 s t=2.1 s
t=2.5 s t=3.1 s
t=3.5 s t=4.0 s
127

Figure 4-34. Velocity field of the fluid particles of the dam break example by CPM
t=0.5 s t=1.0 s
t=1.3 s t=1.4 s
t=1.5 s t=1.7 s
t=2.0 s t=2.1 s
t=2.5 s t=3.1 s
t=3.5 s t=4.0 s
128

Figure 4-35. Geometry and definition of the dam break with obstacle

Figure 4-36. Initial particle distribution of the dam break with obstacle

129

Figure 4-37. Graphical comparisons of the dam break behavior
Left: experimental results by Koshizuka et al. (1995); Middle: PFEM results by Larese et al. (2008);
Right: CPM results
t=0 s
t=0.1 s
t=0.2 s
t=0.3 s
t=0.4 s
t=0.5 s
t=1.0 s
130

Chapter 5 Liquid Sloshing in Rectangular Tanks: Experimental Study and CPM Simulation
131
Chapter 5 Liquid Sloshing in Rectangular Tanks:
Experimental Study and CPM Simulation
5.1 Introduction
Sloshing is the motion of a liquid with free surface in a partially filled liquid
container. The phenomenon is often associated with many engineering problems, and
one particular concern is the transportation of LNG in membrane tanks on ships. The
basic analysis of liquid sloshing involves the estimation of hydrodynamic pressure
distribution, free surface elevations and the natural frequencies of the system. These
parameters have a significant effect on the dynamic stability and performance of the
moving containers. The high pressures and impact forces created can result in serious
implications such as structural deformation within the walls, fatigue of the membrane
layer or even structural failure in severe cases. Another repercussion could be the
capsizing of the ship, which is caused by the high pressure and as a result, a large
overturning moment to cause the instability.
Sloshing is related to the continuously changing free surface. It is difficult to
simulate the free surface, especially when nonlinear phenomena such as breaking
wave and overturned wave are considered. In this regard, sloshing experiment is
essential for understanding and verification of numerical simulating results.
Various experimental studies have been carried out to study dynamic sloshing in
rough sea conditions. The most significant experimental study is the one carried out
by Lloyds Register (2008). Damages to containment system in late 1970s and recent
NO96 vessels have motivated Lloyds Register to carry out experimental investigation
132
to find out the causes. It is concluded that LNG filling depths play a significant role in
affecting the structural integrity of cargo containment systems. Another experimental
study was carried out by Daewoo Shipbuilding & Marine Engineering Co Ltd (DSME)
which was referred to by Chen et al. (2008). DSME has found that high impact
pressure is generated at low filling depth of approximately 30%, which is in line with
Lloyds Registers findings.
In this chapter, water sloshing in rectangular tank under lateral excitation is first
investigated experimentally. Rectangular tanks with two different sizes are designed
to investigate the different aspects of sloshing behavior. Parameters such as water
depth and excitation frequencies are varied to study their effects on sloshing.
Numerical simulations of typical sloshing cases are then performed using the
proposed CPM and the results are compared with the experimental results.
5.2 Experimental setup
5.2.1 Experimental facilities
Figure 5-1 shows the experimental setup. A rectangular tank was mounted on a
shake table which was subjected to lateral displacement excitation. The tank made of
Plexiglas was partially filled with water. The free surface elevation of water was
measured using wave gauges. A pressure sensor was used to measure water pressure
on one side wall (see Figure 5-2 for details). A hydraulic actuator generates horizontal
motion of the shake table. Its displacement was measured by a displacement
transducer mounted on the shake table. All the output signals from the shake-table,
wave gauges and pressure sensor were sent to a digital oscilloscope where data can be
acquired and stored. The dynamic motion of the fluid was recorded by a high speed
133
camera located in front of the tank. The overall instrumentation scheme is illustrated
in Figure 5-3.
5.2.2 Water Tank
The water tank used in the sloshing experiments is schematically shown in
Figure 5-4. The tank with length L and height H was filled with water up to the filling
depth d. A wave probe was placed 20 mm from the right wall at Point P
1
to measure
the free surface elevation () near that corner. The pressure sensor was located at the
left corner of the tank bottom, which is P
2
shown in the figure.
5.2.3 Shake table
A 1.5 m x 1 m shake table driven by a hydraulic actuator (load capacity 25 kN
and stroke capacity 5 cm) is used in the sloshing experiments. The shake table can
only generate unidirectional lateral excitation (translation in one degree of freedom).
In order to ensure that the excitation from the shake table can be fully transmitted
to the tank, the rectangular tank has to be tightly fastened on the shake table. This is
done by anchoring the tank against two square steel tubes that are bolted tightly
against the shake table (Figure 5-5).
5.2.4 Wave probes
Three capacitance type wave probes (KENEK CHT4-60) were installed to
measure the free surface elevation of the sloshing waves. The wave probes worked
based on the principle of linear variation of capacitance with water surface change.
Capacitor is formed from insulated wire held taut by a supporting rod, with water
serving as a ground (see Figure 5-6). The capacitance of the wire changes according
134
to the water surface displacement. Hence, the wave height is obtained by detecting the
change in capacitance of wire. Care was to be taken to ensure that the wave probe is
installed vertically for accuracy.
Calibration is necessary to obtain conversion coefficient between voltage and
water level. Calibration is performed in a deep water bucket, by mounting wave probe
on an adjustable stand which allows it to be adjusted up and down as shown in Figure
5-7. The readings are taken at interval of 10 mm from 10 cm to 27 cm of the water
height. The calibration results are plotted as shown in Figure 5-8, Figure 5-9 and
Figure 5-10 for the three wave probes. As seen in the figures, the calibrated
coefficients of the wave probes 1 to 3 are 49.55 mm/v, 47.47 mm/v and 49.85 mm/v,
respectively.
5.2.5 Pressure sensor
A WIKA Model S-10 pressure sensor was installed on one side wall at about 20
mm from the bottom of the tank (Figure 5-11). It was used to measure the
hydrodynamic pressure response during sloshing. Care was taken to ensure water-
tightness at the connection point.
Calibration was carried out to obtain the different voltage readings corresponding
to the different water pressures. As shown in Figure 5-12, the pressure sensor was
mounted at the bottom of a cylindrical container for calibration. A bubble level was
used to make sure the cylindrical container was standing vertically. Pressure readings
were collected by varying the water level. The calibration result is plotted in Figure
5-13 which shows a calibrated coefficient of 103.87 mm/V. Considering that water
density is =1000 kg/m
3
and gravity g=9.78 m/s
2
, the calibrated coefficient for
pressure sensor is 1015.8 Pa/V.
135
5.2.6 Displacement transducer
The lateral displacement generated by the shake table was measured by using a
conductive plastic resistive element type potentiometer 13FLP25A shown in Figure
5-14. The calibration coefficient for the potentiometer is 10.29 mm/V.
5.2.7 High speed camera
As seen in Figure 5-15, a high speed camera (Nikon/Basler A504K) was set up at
about 2 m in front of the experimental tank for image capture and connected to a data
processor. The camera could record the wave profiles at a maximum frame rate of 500
fps and provided a graphical view of water sloshing. DC lighting was used to provide
the necessary brightness. The images taken were selectively digitized to derive the
free surface profiles.
5.2.8 Other considerations
Precautions were taken to ensure good accuracy of results during the
experiments. Firstly, the water surface was ensured to be stationary before the shake
table was started. Secondly, the tank was placed and fastened to align with the motion
of the shake table such that 2D water sloshing was obtained (Figure 5-16). The inner
surface of the tank was also made smooth to reduce boundary layer effect and
minimize 3D effects of sloshing waves. Thirdly, connecting cables were arranged so
as not to cross one another and sufficient lag lengths were reserved to avoid cable
being stretched when the shake table moved.
136
5.3 Sloshing experiments and comparison with CPM solutions
Liquid sloshing in a tank is a highly nonlinear phenomenon in which the free
surface profiles can be considered as a composition of several wave modes. The
superposition of the modes depends on the liquid depth, tank geometry and external
excitations.
For a rectangular tank, the natural frequency of the fluid can be estimated from
linear theories (Su et al., 1982) by the following equation,
) tanh( d
L
m
L
m
g
m

= (5-1)
where d is the depth of water in the tank, L the length of the tank, m the mode
number and g the acceleration due to gravity taken as 9.78 m/s
2
(gravitation in
Singapore, http://en.wikipedia.org/wiki/Gravity_of_Earth).
As seen from the above equation, there exist an infinite number of natural
frequencies. However, only the fundamental frequency (m=1) is significant for marine
engineering application (Su et al., 1982). In our experimental study only the first
natural frequency of the first sloshing mode is considered.
It should be noted that resonant sloshing is an ideal scenario difficult to achieve
in experiments due to the accuracy in the natural frequency of sloshing. Furthermore,
Hill (2003) and Chen and Nokes (2005) pointed out that the first natural frequency is
shifted due to the nonlinearity of sloshing and the effect of the viscosity of the fluid.
The accuracy of the calculation of the real natural frequency is further decreased by
the precision of the parameters used in Eq. (5-1) which were measured during
experiments. Therefore, in the following resonant sloshing case, the excitation
frequency is close to but is not necessarily exact the natural frequency of liquid
137
sloshing. Based on the estimation of the first natural frequency, a series of
experiments near the resonance frequency were done to investigate violent sloshing
phenomenon.
In our experiment study, the shake table displacement is sinusoidal and governed
by ) cos 1 ( t A x = , where A is the amplitude of excitation and the excitation
frequency. Initially we used the excitation displacement t A x sin = for t0, as many
investigators adopted it such as Faltinsen (1978) and Wu et al. (1998). The
corresponding velocity is t A u cos = in this case. For the shake table given an
initial condition with a non-zero velocity, there will be a big initial impulse occurring
in a very short time. Using ) cos 1 ( t A x = , the initial velocity t A u sin = is zero.
This is found to help reduce the initial impulse of the shake table. Ideally, the initial
conditions of the shake table should be 0
0
=
= t
x , 0
0
=
= t
x and 0
0
=
= t
x . A more
effective way is to introduce ramping function which reduces the influence of the
initial velocity and acceleration, as also pointed out by Faltinsen et al. (2000).
Nevertheless, in this study we adopted the displacement of shake table
) cos 1 ( t A x = for simplicity. One example of the recorded displacement of the
shake table is plotted and shown in Figure 5-17. The excitation amplitude is A =5 mm
and the excitation frequency =6.85 rad/s (1.09 Hz).
Tap water is used in the experiments and the property of water including
viscosity is considered in CPM simulation. The animations of water sloshing by CPM
as well as the experiments recorded by high speed camera of most of the cases
presented in this chapter (except cases for parametric study) are recorded in the
attached CD.
138
5.3.1 Experiments of sloshing waves in high-filling tank
Depending on the liquid filling height in the tank, the sloshing waves can be in
the form of a standing wave, a traveling wave or a bore (Wu et al., 1998). Armenio
and La Rocca (1996) have found that these three types of wave scenario may all
appear depending on the ratio of filling depth over tank length, i.e. L d / . In this
section, experiments are designed to observe the free-surface elevation under relative
high filling depth where 5 . 0 / L d since Chen et al. (2008) and Lloyds Register
(2005) have pointed out that the sloshing is most severe when the filling depth (d) is
from nearly 50% to 70%H, where H is the tank height.
In the study of sloshing in high-filling water tank, a tank with dimension 0.6 m
(L) x 0.6 m (H) x 0.3 m (B) is used. The filling depth of the tank is 0.3 m with the
ratio 5 . 0 / = L d . The parameters of the tank and liquid are shown in Table 5-1.The
first natural frequency is calculated by Eq. (5-1) as 6.85 rad/s.
With the parameters of the tank and liquid listed in Table 5-1 four different cases
are considered with different excitation frequencies. The frequency ratios of the
excitation frequency () to the natural frequency (
0
) of the sloshing system
examined in this section are
0
/ =1.0, 1.1, 0.9, 0.583 respectively. The frequency
ratios were also used by Wu et al. (1998) in their sloshing study by FEM. In all the
four cases, the excitation amplitudes are the same as A=5 mm.
5.3.1.1 Resonant water sloshing with
0
=1
When the frequency of the external excitation is close to the natural frequency of
the sloshing system, resonant sloshing occurs. The time history of the free surface
139
elevation at P
1
is shown in Figure 5-18. One can see in the figure that the wave
amplitude increases with time until steady state is reached.
The MPS method is first used to simulate the case and the time history of the free
surface elevation is plotted in Figure 5-18 for comparison. As can be seen in the
figure, the MPS method can reasonably simulate the overall free surface elevations
although there is some phase difference as compared with the experimental results.
This may be due to the small difference of the excitation frequency assumed in the
MPS method from the actual excitation frequency in the experiment. The overall
sloshing trend and wave profile can reasonably be modeled in the MPS method.
Typical wave profiles at two different time instants (t=2.10 s and t=4.35 s) are shown
in Figure 5-19 (b) in comparison with the experimental results presented in Figure
5-19 (a). The general sloshing phenomenon of liquid in the tank is represented fairly
well by MPS and agrees with the experimental results. It is seen in the figure that
there are some hot spots where particles are tightly spaced and some cold spots
where particles are loosely spaced (see zoom-in figures in Figure 5-19 (c)). These hot
spots indicate too high number density and cold spot too low number density,
implying that the incompressibility condition is not well enforced.
Nevertheless, as mentioned in Chapter 4, the pressure fluctuation in the MPS
method is severe. The pressure contours at the same time instants are shown in Figure
5-20 with local zoom in at the bottom corners. We can see that the pressure
distribution in the liquid is not smooth. A large number of particles with zero pressure
value occur inside the fluid domain below the free surface (dark points in the zoom-in
figures). This greatly affects the pressure field and thus inaccurate pressure
distribution is obtained.
140
As proposed in Chapter 3, there are several ways to improve the pressure
fluctuation in the original MPS method. The arc method for free surface particle
recognition is the first improvement (Sect. 3.3.3). Another improvement relates to
enforcing the incompressibility condition of the free surface particles (Sect. 3.3.3).
The same case of resonant sloshing is studied to demonstrate the effects of the
two improvements mentioned above. With the two improvements incorporated in the
original MPS method, the sloshing waves are better simulated. Figure 5-21 shows the
particle distributions at the same time instants are more regular compared with the
original MPS simulation results in Figure 5-19. The hot and cold spots practically
disappear in the figure. The pressure contours at these two time instants are plotted in
Figure 5-22. Compared with the results shown in Figure 5-20 using the original MPS
method, the pressure contours are considerably smoother. The problem of spurious
zero pressure inside the fluid domain also vanishes. It can reasonably be inferred that
the numerical solution of pressure contour represents the actual one generally well.
Nevertheless, even with the two improvements in the MPS method, the high
frequency pressure fluctuation still appears in the time signal of pressure. Figure 5-23
shows the time history of pressure at P
2
. The comparison between the original MPS
solution and experimental results is made in Figure 5-23 (a). The pressure fluctuation
obtained by the original MPS method is very large with amplitude of fluctuation up to
2-3 times of the actual pressure value. The fluctuation is so severe that it masks the
actual trend of the pressure variation due to sloshing. With the arc method for free
surface recognition used in the MPS method, the pressure fluctuation is greatly
reduced (Figure 5-23 (b)). The amplitude of fluctuation is reduced to less than the
actual pressure value. When the incompressibility adjustment (IA) of the free surface
141
particles is also introduced, the fluctuation of pressure is further reduced (Figure 5-23
(c)). The amplitude of the hydrodynamic pressure can be better estimated in this case.
Now the CPM as proposed in Chapter 3 is used to simulate the same case of
resonant sloshing. The particle distribution and pressure contours at the same time
instants are shown in Figure 5-24 and Figure 5-25, respectively. It can be seen that the
particle distribution in CPM simulation is quite regular even near the free surface. The
pressure contours are much smoother compared with those obtained by MPS in Figure
5-22.
The pressure history at P
2
is shown in Figure 5-26. The comparison of the
pressure solution with experimental result shows excellent agreement. Due to
enforcement of the incompressibility condition for all the fluid particles, there are
only some minor fluctuations of the pressure values when the sloshing amplitude
becomes large. It can be easily improved by imposing some artificial compressibility
of the fluid, as introduced by Hu and Kashiwagi (2004) and Khayyer and Gotoh
(2009). In our study we do not implement this since the fluctuation is so minor that it
can be neglected. A little phase difference is observed after 6 s. This difference is
likely due to the frequency ratio (
0
/ ) used in the experiment is not exactly the
same as the one assumed in the numerical model.
To show the great improvement achieved by CPM over the original MPS, the
pressure solutions are plotted together in Figure 5-27 in comparison with the
experimental result. The tremendous improvement of the pressure history in the
proposed CPM shows the capability of the new particle method in simulation water
sloshing.
142
The comparison of the free surface elevation at P
1
using the proposed CPM with
the experimental result is shown in Figure 5-28. Excellent agreement can be seen in
this figure. Again some phase difference after 6 s is also observed, similar to the
pressure history in Figure 5-26.
The free surface profiles at different time instants simulated using CPM are
shown in Figure 5-29. There is a very good agreement between numerical solutions
and experimental results. Standing waves are observed in the tank. The free surface
becomes higher and higher with time increase during the resonant sloshing until
steady state. Nonlinear sloshing wave with sharp crest and flat trough is observed. It
can be seen in Figure 5-29 that the pressure contours exhibits a reasonable trend and
without irregular pressures inside the fluid domain. The animations of CPM
simulation and experiments of this case are recorded in the attached CD with file
name Slosh_1_0.avi.
5.3.1.2 Water sloshing with
0
=1.1
Besides resonant sloshing, it is interesting to study waves under excitation
frequency near natural frequency. As pointed out by Wu et al. (1998), when the
external excitation frequency is slightly away from the natural frequency, modulated
waves can be obtained. This is also called beating phenomenon by some other
investigators (Pal et al., 1999; Faltinsen et al., 2000; Cho and Lee, 2004). In order to
show this, an excitation frequency ratio 1 . 1 /
0
= is considered in this section.
Similar to the resonant sloshing case, an experimental test was performed with
1 . 1 /
0
= and excitation amplitude A=5 mm. The MPS method is used first to
simulate the sloshing waves. The comparison of the time history of the free surface
143
elevation using the MPS method and experimental results is shown in Figure 5-30.
The beating phenomenon is captured by the wave probe in the experiment. The wave
amplitude is modulated with a clear envelop. The period of the envelope is
approximately 9.5 times of the period of the sloshing waves. The solutions of the free
surface using MPS method in general agree with the experimental results before t=4 s.
However, after that the beating waves are not successfully represented. The envelope
of the modulated waves cannot be observed. The comparison of the hydrodynamic
pressure measured at point P
2
is shown in Figure 5-31. It is again shown that the
pressure fluctuation using the original MPS method is very severe.
Next, the CPM is used. The solutions of the free surface and pressure are
presented in Figure 5-32 and Figure 5-33, respectively. It is clearly shown that the
numerical simulations give excellent agreement with the experimental data, in terms
of both the free surface elevation and pressure. In particular, Figure 5-32 shows that
the beating phenomenon is successfully captured in the CPM simulation. The period
of the envelope of the beating waves agree well with the experiments. In Figure 5-33,
the pressure history presents a reasonable trend without obvious spurious fluctuation
and matches with the experimental results very well. Some spikes in the experimental
results are shown in the figure. This is caused by the shake table for which the
hydraulic actuator motion may not be smoothly produced, i.e. it generates sudden
acceleration during each start of back and forth motion when the direction of motion
is changed.
The free surface profiles at different time instants simulated using CPM are
shown in Figure 5-34. It can be seen that standing waves are present in the rectangular
tank. Good agreement between numerical solutions and experimental results is found.
144
With the excitation ratio 1 . 1 /
0
= in this case, the sloshing wave is modulated at
around 7.95 s, when the free surface elevation becomes nearly flat. After that, the
wave again becomes higher until it comes to the end of the next beating period. It is
also shown in Figure 5-34 that the pressure contours are smooth at every time instant.
The animations of CPM simulation and experiments of this case are recorded in the
attached CD with file name Slosh_1_1.avi.
0
=0.9
Wu et al. (1998) have proved that for the small amplitude sloshing, the frequency
of the envelop of the beating waves is governed by
0
= and its period is
/ 2 . That means the excitation frequency ratios 1 . 1 /
0
= and 9 . 0 /
0
= would
present the same envelope of the beating waves. For large amplitude sloshing
however, where nonlinearity affects the natural frequency of the sloshing waves,
exactly the same envelopes of the two ratios cannot be obtained. But similar beating
waves and close beating frequency are expected.
Hence in this section, we studied the water sloshing waves under external
excitation frequency 9 . 0 /
0
= experimentally and numerically. The numerical
solutions using MPS method are shown in Figure 5-35 and Figure 5-36 for free
surface elevations and pressure history, respectively. Again the envelope of the
beating waves is not represented in the MPS method. The pressure fluctuation shown
in Figure 5-36 is large, similar to the results of excitation frequency ratio 1 . 1 /
0
= .
The simulation results using CPM are shown in Figure 5-37 and Figure 5-38.
Both the free surface elevation and the pressure solutions agree well with the
experimental results. The envelope of the beating waves is well simulated in the
145
proposed CPM. Similar to the case with 1 . 1 /
0
= , pressure spikes in the
experimental result are also shown in Figure 5-36 and Figure 5-38. The free surface
profiles at different time instants are very similar to the profiles shown in Figure 5-34
for the excitation ratio 1 . 1 /
0
= , and hence these as well as the animation are not
presented.
0
=0.583
The previous three cases study the sloshing waves in 50% filling tank with
excitation frequency close to the natural frequency. When the excitation frequency is
far away from the natural frequency, the sloshing effect is not significant. An
excitation frequency ratio 583 . 0 /
0
= used by Wu et al. (1998) is used in the
experiment.
The free surface elevation and pressure history based on the original MPS
method are shown in Figure 5-39 and Figure 5-40, respectively. The simulation
results obtained by the proposed CPM are shown in Figure 5-41 and Figure 5-42.
Pressure spikes in experimental results are also observed for 583 . 0 /
0
= . Both MPS
and CPM can give good agreement with the experiments in term of the free surface
elevations. But the CPM simulation gives smoother pressure history compared with
MPS method.
It can be seen from Figure 5-41 and Figure 5-42 that the sloshing waves under
excitation frequency 583 . 0 /
0
= are not significant. The hydrodynamic pressure
due to sloshing waves can be neglected compared with the hydrostatic pressure.
Hence sloshing under excitation frequency far away from the natural frequency is not
a concern in the study of liquid sloshing.
146
5.3.1.5 Parametric study of excitation frequencies at the same filling depth
To further demonstrate the effects of the external excitation frequencies, the
maximum free surface elevations of the sloshing waves for different excitation
frequencies ( 417 . 1 , 1 . 1 , 0 . 1 , 9 . 0 , 583 . 0 /
0
= ) are plotted in Figure 5-43. Numerical
solutions by CPM are also presented for comparison in the figure with dashed lines.
For all the cases, the filling depths of the tank are 50%, the excitation amplitudes are
A=5 mm.
In order to compare the extreme free surface elevations, the minimum free
surface elevations are plotted with their absolute values. In Figure 5-43 it can be seen
that the CPM solutions (dashed lines) of the extreme wave elevation agree with the
experimental results (solid lines). As expected, the sloshing amplitude of free surface
is large when the excitation frequency is close to the natural frequency.
It is interesting to note in Figure 5-43 that the absolute value of the minimum
elevation of the sloshing waves is much lower than the maximum elevation. This is a
result of wave nonlinearity, where sharp wave crest and flat trough appear in
nonlinear waves (Wu et al., 1998; Wang and Khoo, 2005). It is also seen in the figure
that the nonlinearity becomes more significant when the sloshing is close to resonance,
where the difference between the crest and trough is the largest.
Figure 5-44 plots the maximum and minimum hydrodynamic pressure with
different excitation frequencies. Again it can be seen that the CPM solutions agree
well with experimental results. As the figure shows, the absolute value of the
minimum hydrodynamic pressure is not significantly different from the maximum
147
hydrodynamic pressure. Nonlinearity of the sloshing waves is less significant in
pressure than in wave elevation.
5.3.1.6 Parametric study of excitation amplitudes at the same filling depth
In real situations, LNG carriers often encounter random waves of different
frequencies and amplitudes. Therefore, besides the study of the excitation frequencies,
different excitation amplitudes should also be considered in sloshing experiments.
Resonant sloshing with different excitation amplitudes is studied experimentally
and numerically. The filling depth of the tank is the same as before, i.e. 50% filling of
the tank height. Four different excitation amplitudes, namely 5 mm, 12.5 mm, 15 mm,
20 mm are used in the study. The extreme free surface and pressure are shown in
Figure 5-45 and Figure 5-46, respectively. The CPM solutions (dashed lines) are in
very good agreement with experimental results (solid lines).
From the maximum and minimum free surface elevations shown in Figure 5-45,
one can see that with increasing excitation amplitude, the maximum and minimum
free surface elevation increased almost linearly, with the increasing rate of the
minimum free surface elevation smaller than that of the maximum. This indicates that
as the excitation amplitude increases, nonlinearity of the sloshing wave increases and
gives flatter trough.
Figure 5-46 presents the maximum and minimum hydrodynamic pressure with
different excitation amplitudes. With larger excitation amplitude, the hydrodynamic
pressure increases in an approximately linear manner. The CPM solutions of extreme
pressure agree very well with the experimental results.
148
5.3.1.7 Parametric study of filling depths with the same excitation amplitude
As discussed before, the sloshing waves in high-filling level tanks exhibit the
standing wave types. The comparison of the sloshing waves under different filling
depths is worthy of investigation.
Resonant sloshing with excitation amplitude of 5 mm is studied in this section
for five typical filling depths 25%, 50%, 60%, 70%, and 75%. Figure 5-47 shows that
the maximum surface elevation occurs between 45% - 70% filling depths. This is in
line with the guideline of Lloyds Register (2008).
Figure 5-48 shows the normalized extreme hydrodynamic pressure versus the
ratio of filling depth. The computed solutions by CPM are agreeable with the
experimental data. It is interesting to see that the normalized extreme hydrodynamic
pressure decreases almost linearly as the filling depth increases. It implies that the
proportion of the hydrodynamic pressure decreases in the total pressure when the
hydrostatic pressure is large with high filling depth.
5.3.1.8 Discussion
In this section, water sloshing in tank with filling depths 50%H and above are
studied experimentally and numerically. The MPS method and the proposed CPM are
applied to simulate water sloshing waves. The numerical solutions are compared with
experimental results.
The solution of free surface elevation using the MPS method in general agrees
with the experimental results. Nevertheless, the MPS solution of pressure history at a
fixed position exhibits severe spurious fluctuation. In addition, the particle
distribution in MPS simulation shows hot spots with particles too densely spaced
149
and cold spots with particles too loosely spaced (see Figure 5-19 and Figure 5-20).
This indicates that the incompressibility condition is not satisfactorily enforced. The
proposed ways to improve fluctuation in the MPS method, i.e. the arc method for free
surface particle recognition and incompressibility adjustment of free surface particles,
are proved to some extent to be effective to improve particle distribution (without
obvious hot and cold spots) and help mitigate pressure fluctuation amplitude (see
Figure 5-21).
The proposed CPM is applied to simulate water sloshing with different excitation
conditions. The particle distribution in CPM simulation is regularly spaced with no
hot or cold spot observed even near the free surface (Figure 5-24). The pressure
contour of CPM solution is quite smooth. Compared with MPS results, the CPM
solutions agree better with experimental results in terms of time histories of free
surface elevation and especially pressure.
Different parameters which affect the sloshing waves are analyzed. It is
concluded that a filling depth of approximately 45%-70% of the tank height produces
the highest sloshing pressure and free surface elevation. The larger excitation
amplitude produces higher free surface elevation and hydrodynamic pressure.
Additionally, for various excitation frequencies, resonant frequency produces the
largest free surface elevation and pressure. The conclusions from the research findings
are consistent with other published results (Wu et al., 1998; Cho and Lee 2004;
Ibrahim, 2005). In addition, the experimental data acquired form a valuable database
which could be used to validate numerical models.
150
5.3.2 Experiments of sloshing waves in low-filling tank
Experiments and numerical simulations at partial filling levels in an LNG tank
were investigated by many researchers, indicating that the sloshing motion in an LNG
tank at the low-filling level is quite different from that experienced at the high-filling
level (Ibrahim, 2005). In the low-filling case, a more distinct wave motion is observed
in the tank compared with standing wave in high filling level. At the low-filling level,
where the filling height is less than 20 percent of tank width (for lateral motion) or
length (for longitudinal motion), a phenomenon known as a bore can be observed in
the tank (Huang and Hsiung, 1996). The phenomenon can have some practical
implications on the carriers. When tank motion is large, the front of the hydraulic run-
up becomes steeper, developing a breaking wave. If the bore hits the bulkhead before
breaking, large impact pressure can occur and damage the tank. The uniform velocity
of the traveling wave also results in a large drag force on the lower part of the pump
tower and its supporting system, which is an additional load on the system that must
be considered when designing storage tanks on the carriers. Armenio and La Rocca
(1996) have reported that, depending on the ratio of fluid depth to breadth of tank, a
combination of standing wave, traveling wave and a bore may exist. The same
observation was made by Wu et al. (1998) using a FEM simulation. Therefore it is of
interest to see the CPM performs in simulating sloshing waves in relatively shallow
filling depth.
In this study, a tank with dimension 0.6 m (L) x 0.6 m (H) x 0.3 m (B) is used.
The filling depth of the tank is 0.03 m with the ratio % 20 / < L d . The parameters of
the tank and liquid are shown in Table 5-2. The first natural frequency is calculated by
Eq. (5-1) according to the tank dimensions and filling condition. Resonant sloshing is
151
considered with external excitation frequency close to the first natural frequency of
the sloshing system. The excitation amplitude used in this study is 0.01 m.
Typical wave profiles at several time instants are shown in Figure 5-49 to Figure
5-56. Good agreement is observed between the CPM solutions and experimental
results.
From Figure 5-49 to Figure 5-56, it is observed that traveling wave does not
appear immediately after the tank starts to move. Instead, there is a transient period
during which the wave changes gradually from a standing wave to a traveling wave (i.
e. one wave crest traveling in the tank). After the initial formation of the traveling
wave (t=4.80 s), multi-crest traveling waves form (t=6.10 s). At about t=6.75 s the
wave crest hits the right wall of the tank and hydraulic run-up of the water appear.
The wave then falls back to the tank and results in the formation of bore (t=6.90 s).
Bores are formed when a wave travels towards an on-coming wave and superimposes
with the on-coming wave. The bore only has a short duration of existence (less than
0.1 s) and splits into multi-crested traveling waves or, in some cases, breaks in the
middle of the tank. As seen in the figures, at t=6.90 s, on the left of the bore, the water
surface is almost flat and the free surface elevation is very small; and on the right of
the bore, the wave elevation is significantly higher. The bore quickly splits into multi-
crested traveling waves (t=7.1 s). After that, multi-crested waves travel in the tank and
hit the left wall of the tank at t=7.85 s, where the water run-up along the wall of the
tank and breaking waves appear when the impact occurs. As can be seen from Figure
5-49 to Figure 5-56, the CPM solutions of water waves agree well with the
experimental results.
152
The long time simulation of this example is performed and the time history of the
free surface elevation at P
1
is shown in Figure 5-57. It is shown in the figure that the
CPM solution agrees very well with the experimental result. The comparison of
pressure history at P
2
is plotted in Figure 5-58. Again very good agreement between
numerical simulation and experimental results is observed. The animations of CPM
simulation and experiments of this case are recorded in the attached CD with file
name Slosh_low.avi.
The pressure contours at different time instants are shown in Figure 5-59. The
pressure field of the liquid is smooth and stable. In addition, large total pressure
occurs beneath the crest of the traveling waves and travels with the crest of the wave.
Furthermore, when the wave peaks reaches the wall, the pressure around the wall is as
expected larger that the hydrostatic pressure.
In order to illustrate the bore formation more clearly, the velocity fields of the
particles are shown in Figure 5-60 at the time when a bore is generated. It is shown in
the figure, at t=6.75 s the first crest of the traveling wave hits the right wall of the tank.
Breaking particles appear during the impact. After that the water particles falls back to
the main water body with velocities in the reverse direction (t=6.80 s). The wave then
travels towards the second wave crest and superimposes with this on-coming wave, as
shown in the figure (from t=6.80 s to t=6.90 s) The superimposed wave elevation
becomes highest at t=6.90 s, where bore formation is clearly observed. After that the
two wave crests split (at about t=6.95 s) and travel in opposite directions (from 6.95 s
to 7.10 s) until the second wave crest impacts with the tank wall and reverses its
direction of travel.
153
5.3.3 Experiments with sloshing wave impact on the tank ceiling
In the previous sections, liquid sloshing in deep and shallow water depths are
studied numerically and experimentally, with the main focus on resonant sloshing.
Breaking wave on the tank ceiling due to sloshing is not investigated since the two
tanks used in the experiments are open. In order to investigate the breaking waves, a
closed tank is built and breaking waves on the tank ceiling is studied.
A closed rectangular tank with dimension 0.57 m (L) x 0.3 m (H) x 0.3 m (B)
was build. The parameters of the tank and liquid are shown in Table 5-3. The filling
depth of the tank is 60% which is 0.18 m. Sloshing waves under resonant excitation
is considered with external excitation amplitude is 10 mm.
Figure 5-61 shows the free surface profiles of the liquid at different time instants
inside the tank. The comparison between numerical simulation and experimental
results shows good agreement. As shown in the figure, under resonant excitation
frequency, the sloshing waves become higher and higher. Finally the wave reaches the
top of the tank and breaking waves appear. It can be seen in the figure that the
numerical solution successfully captures the appearance of breaking waves through
breaking particles. The particle separation due to breaking and coalescence after
falling back to the fluid is well represented by the proposed CPM. The animations of
CPM simulation and experiments of this case are recorded in the attached CD with
file name Slosh_breaking.avi.
Since a closed tank was used in the experiment to study wave impacts onto the
tank top, no wave probe could be inserted to measure the free surface elevation. Video
camera was thus used in the experiment to record the wave profile. A transparent
template with regular grid of 1 cm interval was glued to the front face of the tank to
154
facilitate measurement of free surface elevation. Hence the free surface elevation at
the right wall of the tank could be obtained. The free surface elevation data is plotted
in Figure 5-62 and compared with the CPM solution. It is shown in the figure that the
CPM solution agrees very well with the experimental result. The free surface
elevation gets higher and higher and quickly reaches to a constant value 0.3, which
corresponds to the ceiling height of the tank.
Figure 5-63 plots the pressure contours at four representative time instants. It is
shown in the figure that, the numerical simulation results present satisfactory pressure
field. The pressure contours are smooth even when breaking waves occur. One can
observe in the figure that when sloshing waves reach the tank ceiling, the impact
pressure generated at the corner of the tank ceiling become high (at 4.56 s and 8.44 s).
Due to the sensor constraint, we did not measure the impact pressure at the tank
ceiling. The CPM solutions of pressure history at the two corners of the left wall
(marked in Figure 5-63) are shown in Figure 5-64 (left bottom corner P
2
) and Figure
5-65 (left ceiling corner P
1
). By comparing the two figures it can be seen that the
impact pressure at ceiling corner P
1
is very high (about 1700 Pa) and even larger than
the hydrodynamic pressure at bottom corner P
2
. More severely, the impact occurs in
very short time instant, as shown in Figure 5-65. Impact pressure due to sloshing is
another important aspect which should be addressed in future study.

In this chapter, sloshing waves in rectangular tanks are studied experimentally
and numerically.
155
Sloshing in deep water exhibits nonlinear standing waves. The highest free
surface elevation and largest hydrodynamic pressure are observed when resonance
occurs. Larger excitation amplitude produces higher free surface elevation and
hydrodynamic pressure. A filling depth of approximately 45%-70%H produces the
highest sloshing pressure and free surface elevation. Sloshing in low-filling tank
produces traveling waves which lead to the formation of bore when they superpose
one another. Breaking waves with high impact pressure at the tank ceiling appear due
to violent liquid sloshing in a closed tank.
The MPS solution of pressure history in liquid sloshing exhibits severe spurious
fluctuation. The proposed two ways to improve pressure fluctuation in MPS method,
i.e. the arc method for free surface particle recognition and incompressibility
adjustment of free surface particles, are shown effective to improve particle
distribution and help suppress pressure fluctuation amplitude.
The CPM simulation of liquid sloshing presents evenly distributed particle
distribution without hot or cold spot. The formation of the bore in low-filling sloshing
and breaking waves with smooth pressure contours are successfully demonstrated in
CPM simulation. The pressure contour of CPM solution is smooth even when
breaking occurs. Compared with MPS results, the CPM solutions agree much better
with the experimental results particularly in terms of particle distribution and pressure.

156

157

Table 5-1. Parameters of the tank and liquid for high-filling sloshing
L
d
H

g
0

Length of tank
Still water depth
Height of tank
Density of water
Dynamic viscosity of water
Gravity
First natural frequency

0.6 m
0.3 m
0.6 m
1000 kg/m
3

1x 10
-3
Pa.s
9.78 m/s
3

6.85 rad/s

Table 5-2. Parameters of the tank and liquid for low-filling sloshing
L
d
H

g
0
A
Length of tank
Still water depth
Height of tank
Density of water
Gravity
Excitation amplitude

0.6 m
0.03 m
0.6 m
1000 kg/m
3

1x 10
-3
Pa.s
9.78 m/s
3

2.825 rad/s
0.01 m

Table 5-3. Parameters of the tank and liquid for high-filling sloshing with breaking
L
d
H

g
0
A
Length of tank
Still water depth
Height of tank
Density of water
Gravity
Excitation amplitude

0.57 m
0.18 m
0.3 m
1000 kg/m
3

1x 10
-3
Pa.s
9.78 m/s
3

6.393 rad/s
0.01 m

158

159

Figure 5-1. The experimental setup

Figure 5-2. Tank with pressure sensor mounted

Figure 5-3. Experimental apparatus and working principle
Pressure
sensor
160

Figure 5-4. Definition of parameters for liquid sloshing in a rectangular tank

Figure 5-5. Fixing tools of rectangular tank on the shake table

Figure 5-6. Schematic Diagram and picture of wave probe
600mm

161

Figure 5-7. (a) Wave probe mounted to an adjustable stand (b) Calibration of the wave probes

Wave Probe 1 Calibration
y =-4.955x +25.989
R
2
=1.000
13
15
17
19
21
23
25
27
0 0.5 1 1.5 2 2.5 3
Vol tage (V)
W
a
t
e
r

L
e
v
e
l

(
c
m
)

Figure 5-8. Calibration results of wave probe 1
y =-4.747x +21.988
R
2
=1.000
7
9
11
13
15
17
19
21
23
0 0.5 1 1.5 2 2.5 3
Vol tage (V)
W
a
t
e
r

L
e
v
e
l

(
c
m
)

(a) (b)
Wave
probe
162
y =-4.985x +26.037
R
2
=1.000
13
15
17
19
21
23
25
27
0 0.5 1 1.5 2 2.5
Vol tage (V)
W
a
t
e
r

L
e
v
e
l

(
c
m
)


Figure 5-11. Experimental installation of pressure sensor

Figure 5-12. (a) Devices used for the calibration of pressure sensor (b) calibration of the pressure
sensor

(a) (b)
Bubble
level
Pressure
sensor
163
Pressure Sensor Calibration
y =10.387x - 2.708
R
2
=1.000
0
10
20
30
40
50
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
Vol tage (V)
W
a
t
e
r

L
e
v
e
l

(
c
m
)

Figure 5-13. Calibration result of the pressure sensor

Figure 5-14. Displacement transducer and its experimental installation

Figure 5-15. High speed camera set-up
Displacement
transducer
DC
Light
Data
processor
164

Figure 5-16. Securing technique

Figure 5-17. Displacement signal of shake table (5mm amplitude)

1
( 0 . 1 /
0
= )
165
(a)
(b)

(c)
Figure 5-19. Particle distribution simulated by MPS method and comparison with experimental results
(a) Experimental results; (b) MPS results; (c) zoom in of part in (b)

t=2.10 s t=4.35 s
t=2.10 s t=4.35 s
166

Figure 5-20. Pressure contours simulated by MPS method

t=2.10 s t=4.35 s
Cold spots area where loosely-spaced particles
are inaccurately recognized with zero pressure
Hot spots area with
accumulated particles
167

Figure 5-21. Particle distribution of MPS simulation with arc method and IA of free surface particles

Figure 5-22. Pressure contours of MPS solution with arc method and IA of free surface particles
t=2.10 s t=4.35 s
t=2.10 s t=4.35 s
168

2
( 0 . 1 /
0
= )
169

Figure 5-24. Particle distribution of CPM simulation

Figure 5-25. Pressure contours of CPM solution
t=2.10 s t=4.35 s
t=2.10 s t=4.35 s
170

2
( 0 . 1 /
0
= )

2
( 0 . 1 /
0
= )

1
( 0 . 1 /
0
= )

171

0
= ) by CPM

t=0.00 s t=0.00 s
t=2.10 s t=2.10 s
t=4.35 s t=4.35 s
t=5.775 s t=5.775 s
t=7.18 s t=7.18 s
t=8.55 s t=8.55 s
t=9.48 s t=9.48 s
t=9.975 s t=9.975 s
172

1
( 1 . 1 /
0
= )

2
( 1 . 1 /
0
= )

173

1
( 1 . 1 /
0
= )

2
( 1 . 1 /
0
= )

174

0
= ) by CPM
t=0.00 s t=0.00 s
t=3.70 s t=3.70 s
t=5.05 s t=5.05 s
t=6.40 s t=6.40 s
t=7.95 s t=7.95 s
t=9.50 s t=9.50 s
t=11.57 s t=11.57 s
t=17.32 s t=17.32 s
175

1
( 9 . 0 /
0
= )

2
( 9 . 0 /
0
= )

1
( 9 . 0 /
0
= )
176

2
( 9 . 0 /
0
= )

1
( 583 . 0 /
0
= )

2
( 583 . 0 /
0
= )

177

1
( 583 . 0 /
0
= )

2
( 583 . 0 /
0
= )

178

Figure 5-43. Maximum and minimum free surface elevations vs. excitation frequency

Figure 5-44. Maximum and minimum hydrodynamic pressure vs. excitation frequency

179

Figure 5-45. Maximum and minimum free surface elevations under different excitation amplitudes

Figure 5-46. Maximum and minimum hydrodynamic pressure under different excitation amplitudes

180

Figure 5-47. Maximum and minimum free surface elevations vs. filling depths

Figure 5-48. Maximum and minimum hydrodynamic pressure vs. filling depths

181

Figure 5-49. Standing Wave at the initial stage, t =2.75 s

Figure 5-50. Traveling wave starts to form, t =4.80 s

Figure 5-51. Multi-crested waves traveling, t =6.10 s

Figure 5-52. Hydraulic run-up (before formation of bore), t =6.75 s

CPM
CPM
CPM
CPM
182

Figure 5-53. Formation of bore, t =6.90 s

Figure 5-54. Bore splits into multi-crested traveling waves, t =7.10 s

Figure 5-55. Multi-crested traveling waves, t =7.50 s

Figure 5-56. Hydraulic run-up, t =7.85 s
CPM
CPM
CPM
CPM
183
0.00
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0 2 4 6 8 10 12 14 16 18
Ti me (s)
F
r
e
e

s
u
r
f
a
c
e

e
l
e
v
a
t
i
o
n

(
m
)
Experiment CPM

Figure 5-57. Comparison of the free surface elevation at Point P
1

0
200
400
600
0 2 4 6 8 10 12 14 16 18
Ti me (s)
P
r
e
s
s
u
r
e

(
P
a
)
Experiment CPM

Figure 5-58. Comparison of the pressure history at Point P
2

184

Figure 5-59. Pressure contours for different time instants by CPM

t=2.75 s t=4.80 s
t=6.10 s t=6.75 s
t=6.90 s t=7.10 s
t=7.50 s t=7.85 s
185

Figure 5-60. Velocity field of the fluid particles by CPM
t=6.80 s
t=6.85 s
t=6.90 s
t=6.95 s
t=7.00 s
t=6.75 s
t=6.60 s
t=7.10 s
First crest hits the wall
Traveling wave crests
Formation of bore
Highest elevation
Bore starts to split
Bore splits to traveling wave crests
2
nd
crest
1
st
crest
186

Figure 5-61. Free surface profiles at different time instants. Left: experiments; Right: CPM simulation.

t=0 s
t=2.19 s
t=2.66 s
t=3.66 s
t=4.6 s
t=5.08 s
187
0
0.1
0.2
0.3
0.4
0 2 4 6 8 10
Ti me (s)
F
r
e
e

s
u
r
f
a
c
e

e
l
e
v
a
t
i
o
n

(
m
)
Experiment CPM

Figure 5-62. Time history of free surface elevations at the right wall

Figure 5-63. Pressure contours for different time instants by CPM
t=4.56 s t=4.57 s
t=5.08 s t=8.44 s
P
1

P
2

188
0
500
1000
1500
2000
2500
3000
0 2 4 6 8 10
Ti me (s)
P
r
e
s
s
u
r
e

(
P
a
)
CPM

2

0
500
1000
1500
2000
0 2 4 6 8 10
Ti me (s)
P
r
e
s
s
u
r
e

(
P
a
)
CPM

1

Chapter 6 Conclusions and Future Research
189
6.1 Conclusions
In this thesis, a new particle method named CPM is developed and applied to
simulate incompressible free surface flows. The CPM is Lagrangian particle method
which solves the Navier-Stokes equations in a semi-implicit time stepping scheme.
There are three main features in CPM compared with other particle methods such as
SPH and MPS as follows.
(1) The discretization of Laplace operator needed in the Poisson equation of
pressure and gradient model are consistent with Taylor series expansion. As a result,
the approximation of Laplacian in CPM gives much better accuracy and convergence
compared with MPS and SPH methods, particularly when the particles are not
regularly spaced.
(2) A free surface particle recognition method (arc method) is applied to capture
the free surface more accurately. The arc method produces more accurate boundary
particle recognition results compared with the density threshold method in MPS.
(3) The incompressibility condition of the free surface particles is enforced by an
adjustment scheme. The scheme is shown to be able to improve the distribution of the
free surface particles and hence improves the stability of the algorithm.
These improvements are shown to be beneficial in improving numerical stability
of the algorithm and in particular resolve the problem of spurious pressure fluctuation.
A benchmark dam break example is first presented to demonstrate the performance of
CPM in comparison with the MPS method. The CPM solution shows remarkably
190
better particle distribution and smoother pressure contours than the MPS solution. The
CPM solution also agrees well with published experimental results both in terms of
free surface profiles and its development over time. A parametric study is performed
to investigate the effects of different parameters used in CPM. The inverse distance
weighting function is found to be suitable in the weighted least-square solution of
derivatives in CPM. The initial particle distance (
0
L ) is dependent on the geometry
and fluid motion type which may require a convergence study by starting with a
coarse distribution of particles. The influence radius
0
1 . 2 L r
e
= or
0
5 . 2 L r
e
= is
recommended as a good compromise between capturing sufficient particles for
accuracy of derivatives and computational efficiency in the simulation of both gentle
and violent fluid motion. In terms of computational efficiency, the time used in the
addition time required to solve derivatives based on local particles in CPM is not
significant particularly in problems with large particle numbers.
The performance of CPM is demonstrated by more numerical examples, i.e.
hydrostatic pressure computation in a static tank, free oscillation of liquid in a
container, dam break in large rectangular tank and dam break in a tank with obstacle
block. In all these examples, the CPM solution of free surface profiles agree very well
with published numerical and/or experimental results. For the dam break problems
without and with obstacle, the pressure solutions of CPM agree well with
experimental results. No severe pressure fluctuation appears even in violent free
surface motion with wave breaking. This is a great advantage of the CPM since
conventional particle methods (such as SPH and MPS) tend to give severe and
spurious pressure fluctuation particularly in long time simulation.
191
The proposed CPM is thereafter applied to study liquid sloshing in rectangular
tanks. Experiments were conducted on rectangular tanks partially filled with water,
making use of a 1.5 m x 1 m shake-table. Sloshing waves in high and low filling
levels are investigated. Parametric studies with respect to excitation frequency,
excitation amplitude and filling depths are performed experimentally and numerically.
Standing waves, traveling waves and breaking waves in different filling depths are
well captured by CPM. The CPM simulation of liquid sloshing presents evenly
distributed particle distribution without hot or cold spot and smooth pressure contour
even when breaking occurs. Compared with MPS results, the CPM solutions agree
much better with the experimental results particularly in terms of particle distribution
and pressure.
In summary, the proposed CPM is shown to be robust and accurate in simulation
of incompressible free surface flows in all the cases considered, in terms of free
surface profile, wave motion and breaking, particle distribution and pressure,
especially the good solution of the pressure filed.
6.2 Future work
In this study, incompressible free surface flow problems are considered and
simulated by CPM. The capability of CPM in other kinds of fluid motion, such as
two-phase flow, has not been addressed. Further work should be conducted to identify
the effectiveness of the proposed CPM in simulation of sloshing involving two phase
fluids (e.g. water and air).
While the proposed CPM has shown great potential in 2D simulation, more
works are needed to extend it to 3D free surface flow problems, such as strong waves
192
impacting on coastal and offshore objects. Flow motions with more complicated
boundaries and excitation conditions should also be tested so as to fully exploit the
research and commercial potential of the CPM.
In addition, due to the facility constraint of the shake-table in the water sloshing
experiments, only unidirectional regular excitation signal was used. Multi-degree of
freedom movements of the tank can be applied in future study. Furthermore, more
complicated or irregular external excitation signals may be investigated. Different
tank shape other than rectangular tank should be studied in future to better model the
real tank geometry.
Lastly, rigid tank wall has been assumed in the study of liquid sloshing, which
means there is no deformation of the wall under the liquid pressure during the
sloshing procedure. This may not be true in the real tank situation, e. g. in membrane
LNG tanks. Fluid-structure interaction should be considered in such a case to account
for the effects of tank flexibility.

References
193
References
Abramson, H. N. (1996). The dynamic behaviour of liquids in moving containers,
Report of NASA SP, 106.
Abramson, H. N. (2003). Dynamics of contained liquids: A personal odyssey, Applied
Mechanics Reviews, 56(1), R1-R7.
Ahmadi, A., Badiei, P. and Namin, M. M. (2007). An implicit two-dimensional non-
hydorstatic model for free-surface flows, International Journal for Numerical
Methods in Fluids, 54, 1055-1074.
Akyildiz, H. and Unal, E. (2005). Experimental investigation of pressure distribution
on a rectangular tank due to the liquid sloshing, Ocean Engineering, 32(11-12), 1503-
1516.
Alam, A., Kai, H. and Suzuki, K. (2007). Two-dimensional numerical simulation of
water splash phenomena with and without surface tension, Journal of Marine Science
and Technology, 12, 59-71.
Apsley, D. D. and Hu, W. (2003). CFD simulation of two- and three-dimensional
free-surface flow, International Journal for Numerical Methods in Fluids, 42, 465-
491.
Armenio, V. (1997). An improved MAC method (SIMAC) for unsteady high-
Reynolds free surface flows, International Journal for Numerical Methods in Fluids,
24, 185-214.
Armenio, V. and La Rocca, M. (1996). On the analysis of sloshing of water in
rectangular containers: Numerical study and experimental validation, Ocean
Engineering, 23(8), 705-739.
Ataie-Ashtiani, B. and Farhadi, L. (2006). A stable moving-particle semi-implicit
method for free surface flows, Fluid Dynamics Research, 38(4), 241-256.
Ataie-Ashtiani, B., Shobeyri, G. and Farhadi, L. (2008). Modified incompressible
SPH method for simulating free surface problems, Fluid Dynamics Research, 40, 637-
661.
Balendra, T., Ang, K. K., Paramasivam, P. and Lee, S. L. (1982a). Seismic design of
flexible cylindrical liquid storage tanks, Earthquake Engineering and Structural
Dynamics, 10, 477-496.
Balendra, T., Ang, K. K., Paramasivam, P. and Lee, S. L. (1982b). Free vibration
analysis of cylindrical liquid storage tanks, International Journal of Mechanical
Sciences, 24(1), 47-59.
Balendra, T., Wang, C. M. and Cheong, H. F. (1995). Effectiveness of tuned liquid
column dampers for vibration control of towers, Engineering Structures, 17(9), 668-
675.
References
194
Balendra, T., Wang, C. M. and Rakesh, G. (1999). Vibration control of various types
of buildings using TLCD, Journal of Wind Engineering and Industrial Aerodynamics,
83, 197-208.
Balendra, T., Wang, C. M. and Yan, N. (2001). Control of wind-excited towers by
active tuned liquid column damper, Engineering structures, 23 (9), 1054-1067.
Bass, R. L., Bowles, E. B., Trudell, R. W., Navickas, J ., Peck, J . C., Yoshimura, N.,
Endo, S. and Pots, B. F. M. (1985). Modeling criteria for scaled LNG sloshing
experiments, Transactions of the ASME, 107, 272280.
Bathe, K. J . (1996). Finite Element Procedures, Prentice Hall Englewood Cliffs, New
J ersey.
Batina, J . T. (1992). A gridless Euler/Navier-Stokes solution algorithm for complex
two dimensional applications, NASA Technical Memorandum, 107631.
Batina, J . T. (1993). A gridless Euler/Navier-Stokes solution algorithm for complex
aircraft applications, AIAA paper 93-0333, 31st Aerospace Sciences Meeting and
Exhibition, Reno, NV, USA.
Bermudez, A., Rodriguez, R., and Santamarina, D. (2003). Finite element
computation of sloshing modes in containers with elastic baffle plates, International
Journal for Numerical Methods in Engineering, 56(3), 447-467.
Bettess, R. (1981). Operation counts for boundary integral and finite element methods,
International Journal for Numerical Methods in Engineering, 15, 306- 308.
Biswal, K. C., Bhattacharyya, S. K., and Sinha, P. K. (2003). Free-vibration analysis
of liquid-filled tank with baffles, Journal of Sound and Vibration, 259(1), 177-192.
Brown, S. J . (1982). A survey of studies into the hydrodynamic response of fluid
coupled cylinders, ASME Journal of Pressure Vessel Technology, 104(1), 2-20.
Bucchignani, E. (2004). A numerical study of non-linear dynamics in a tank for
aerospace applications, Applied Numerical Mathematics, 49, 307-318.
Bucchignani, E., Stella, F. and Paglia, F. (2004). A partition method for the solution
of a coupled liquid-structure interaction problem, Applied Numerical Mathematics, 51,
463-475.
Chen, B. F. and Chiang, H. W. (2000). Complete two-dimensional analysis of sea-
wave-induced fully non-linear sloshing fluid in a rigid floating tank, Ocean
engineering, 27, 953-977.
Chen, B. F. and Nokes, R. (2005). Time-independent finite difference analysis of fully
non-linear and viscous fluid sloshing in a rectangular tank, Journal of Computational
Physics, 209, 47-81.
Chen, W., Haroun, M. A., and Liu, F. (1996). Large amplitude liquid sloshing in
seismically excited tanks, Earthquake Engineering and Structural Dynamics, 25(7),
653-669.
References
195
Chen, Y. G., Djidjeli, K. and Price, W. G. (2008). Numerical simulation of liquid
sloshing phenomena in partially filled containers, Computers and Fluids, 38, 830-842.
Chew, C. S., Yeo, K. S.,and Shu, C. (2006). A generalized finite-difference (GFD)
ALE scheme for incompressible flows around moving solid bodies on hybrid
meshfree-Cartesian grids, Journal of Computational Physics, 218, 510-548.
Chiu, C. H. (2006). Commercial and technical considerations in the developments of
offshore liquefaction plant, Proceedings of the 23rd World Gas Conference,
Amsterdam 2006.
Cho, J . R. and Lee, H. W. (2004). Non-linear finite element analysis of large
amplitude sloshing flow in two-dimensional tank, International Journal for
Numerical Methods in Engineering, 61(4), 514-531.
Colagrossi, A. and Landrini, M. (2003). Numerical simulation of interfacial flows by
smmothed particle hydrodynamics, Journal of Computational Physics, 191, 448-475.
Courant, R., Friedrichs, K. and Lewy, H. (1967). On the partial difference equations
of mathematical physics, IBM Journal, 215-234 (English translation of the 1928
German original).
Delorme, L., Colagrossi, A., Souto-Iglesias, A., Zamora-Rodrguez, R. and Botia-
Vera, E. (2009). A set of canonical problems in sloshing, Part I: Pressure field in
forced rollcomparison between experimental results and SPH, Ocean Engineering,
36, 168-178.
Dilts, G. A. (2000). Moving least square particle hydrodynamics II: Conservation and
boundaries, International Journal for Numerical Methods in Engineering, 48, 1503-
1524.
Dutta, S. and Laha, M. K. (2000). Analysis of the small amplitude sloshing of a liquid
in a rigid container of arbitrary shape using a low-order boundary element method,
International Journal for Numerical Methods in Engineering, 47(9), 1633-1648.
Eatock Taylor, R., Wu, G. X., Bai, W. and Hu, Z. Z. (2008). Numerical wave tanks
based on finite element and boundary element modeling, Journal of Offshore
Mechanics and Arctic Engineering, 130, 031001.
Edelsbrunner, H., Kirkpatrick, D. G. and Andseidel, R. (1983). On the shape of a set
of points in the plane, IEEE Transactions on Information Theory, IT-29(4), 551-559.
Edelsbrunner, H. and Mucke, E. P. (1994). Three-dimensional alpha shapes, ACM
Transactions on Graphics, 13(1), 43-72.
Energy Market Authority (2006). Integrated summary report for proposed Singapore
LNG terminal, Singapore, Aug 2006.
Faltinsen, O. M. (1978). A numerical non-linear method of sloshing in tanks with two
dimensional flow, Journal of Ship Research, 18(4), 224-241.
References
196
Faltinsen, O. M., Rognebakke, O. F. and Timokha, A. N. (2003). Resonant three-
dimensional nonlinear sloshing in a squre-base basin, Journal of Fluid Mechanics,
487, 1-22.
Faltinsen, O. M., Rognebakke, O. F. and Timokha, A. N. (2006). Transient and
steady-state amplitudes of resonant three-dimensional sloshing in a square base tank
with a finite fluid depth, Physics of Fluids, 18, 012103.
Faltinsen, O. M., Rognebakke, O. F., Lukovsky, I. A. and Timokha, A. N. (2000).
Multidimensional modal analysis of nonlinear sloshing in a rectangular tank with
finite water depth, Journal of Fluid Mechanics, 407, 201-234.
Fang, J ., Parriaux, A., Rentschler, M. and Ancey, C. (2009). Improved SPH methods
for simulating free surface flows of viscous fluids, Applied Numerical Mathematics,
59, 251-271.
Fekken, G. (1998). Numerical simulation of green water loading on the foredeck of a
ship, MSc-thesis, August 1998.
Fenner, R. (1983). The boundary integral equation (boundary element) method in
engineering stress analysis, Journal of Strain Analysis for Engineering Design, 18(4),
199-205.
Fletcher, C. A. J . (1991). Computational techniques for fluid dynamics, Springer:
Berlin, 1991.
Foss, M. M. (2007). Introduction to LNG, Center for Energy Economics, University
of Texas, Austin.
Frandsen J . B. (2004). Sloshing motions in excited tanks, Journal of Computational
Physics, 196, 53-87.
Frandsen, J . B. and Borthwick, A. G. L. (2003). Simulation of sloshing motions in
fixed and vertically excited containers using a 2-D inviscid sigma-transformed finite
difference solver, Journal of Fluids and Structures, 18(2), 197-214.
Gavete, L., Gavete, M. L. and Benito, J . J . (2003). Improvements of generalized finite
difference method and comparison with other meshless method, Applied
Mathematical Modelling, 27, 831-847.
Gavrilova, E. (2004). Coupled frequencies of a fluid-structure interaction cylindrical
system, Proceedings of the International Congress of Theoretical and Applied
Mechanics (XXI), 15-21 Aug, Warsaw, Poland.
Gedikli, A. and Erguven, M. E. (2003). Evaluation of sloshing problem by variational
boundary element method, Engineering Analysis with Boundary Elements, 27(9), 935-
943.
Godderidge, B., Turnock, S., Earl, C. and Tan, M. (2009). The effect of fluid
compressibility on the simulation of sloshing impacts, Ocean Engineering, 36(8),
578-587.
References
197
Gotoh, H., Ikari, H., Memita, T. and Sakai, T. (2005). Lagranginan particle method
for simulation of wave overtopping on a vertical seawall, Coastal Engineering
Journal, 47(2-3), 157-181.
Gotoh, H. and Sakai, T. (1999). Lagrangian simulation of breaking waves using
particle method, Coastal Engineering Journal, 41(3-4), 303-326.
Greaves, D. M. (2007). Viscous waves and wave-structure interaction in a tank using
adapting quadtree grids, Journal of Fluids and Structures, 23, 1149-1167.
Greco, M., Faltinsen, O. M. and Landrini, M. (2005). Shipping of water on a two-
dimensional structure, Journal of Fluid Mechanics, 525, 309-332.
Greco, M., Landrini, M. and Faltinsen, O. M. (2004). Impact flows and loads on ship-
deck structures, Journal of Fluids and Structures, 19, 251-275.
Griebel, M., Dornseifer, T. and Neunhoeffer T. (1998). Numerical Simulation in Fluid
Dynamics, Society for Industrial and Applied Mathematics.
Grilli, S. T. and Svendsen, I. A. (1990). Corner problems and global accuracy in the
boundary element solution of nonlinear wave flows, Engineering Analysis with
Boundary Elements, 7(4), 178-195.
Gu, H., Li, Y., Li, S. and Zhang, Q. (2005). The level set and particle level set method
for tracing interface, Journal of Hydrodynamics, Series A, 20(2), 152-160.
Guillot, M. J . (2006). Application of a discontinuous Galerkin finite element method
to liquid sloshing, Journal of Offshore Mechanics and Arctic Engineering, 128, 1-10.
Harlow, F. H. (1963). The particle-in-cell method for numerical solution of problems
in fluid dynamics, in: Proceedings of Symposia in Applied Mathematics, 1963.
Harlow, F. H. and Welch, J . E. (1965). Numerical Calculation of Time-Dependent
Viscous Incompressible Flow of Fluid with Free Surface, Physics of Fluids, 8(12),
2182-2189.
Haynes, W. M. ed. (2010). CRC Handbook of Chemistry and Physics, 91st Edition, CRC
Press/Taylor and Francis, Boca Raton, FL.
Hill, D. F. (2003). Transient and steady-state amplitudes of forced waves in
rectangular basins, Physics of Fluids, 15(6), 1576-1587.
Hirt, C. W. and Nichols, B. D. (1981). Volume of fluid (VOF) method for the
dynamics of free boundaries, Journal of Computational Physics, 39, 201-225.
Housner, G. W. (1957). Dynamic pressures on accelerated fluid containers, Bulletin of
the Seismological Society of America, 47, 15-35.
http://en.wikipedia.org/wiki/Properties_of_water
http://webbook.nist.gov/chemistry
References
198
http://www.businesswire.com/multimedia/home/20081217005080/en/1735488/Exxon
Mobil-Technology-Yields-World%E2%80%99s-Largest-LNG-Carrier
Hu, C. H. and Kashiwagi, M. (2004). A CIP method for numerical simulations of
violent free surface flows, Journal of Marine Science and Technology, 9(4), 143-157.
Huang, S., Duan, W. and Zhu, X. (2010). Time-domain simulation of tank sloshing
pressure and experimental validation, Journal of Hydrodynamics, 22(5), 556-563.
Huang, Z. J . and Hsiung, C. C. (1996). Nonlinear shallow-water flow on deck,
Journal of Ship Research, 40(4), 303-315.

Huijsmans, R. H. M., Tritschler, G., Gaillarde, G. and Dallinga, R. P. D. (2004).
Sloshing of partially filled LNG carriers, Proceedings of the Fourteenth International
Offshore and Polar Engineering Conference, Toulon, Fance, May 23-28, 2004.
Ibrahim, R. A. (2005). Liquid sloshing dynamics: theory and applications, New York,
Cambridge University Press.
Idelsohn, S. R, Onate, E. and Del Pin, F. (2004). The particle finite element method: a
powerful tool to solve incompressible flows with free-surfaces and breaking waves,
International Journal for Numerical Methods in Engineering, 61, 964-989.
Idelsohn, S. R. and Onate, E. (2006). To mesh or not to mesh. That is the question,
Computer Methods in Applied Mechanics and Engineering, 195, 4681-4696.
Idelsohn, S. R., Onate, E., Calvo, N. and Del Pin, F. (2003). The meshless finite
element method, International Journal for Numerical Methods in Engineering, 58(6),
893-912.
Idelsohn, S. R., Storti, M.A. and Onate, E. (2001). Lagrangian formulations to solve
free surface incompressible inviscid fluid flows, Computer Methods in Applied
Mechanics and Engineering, 191, 583-593.
Ikegawa, M. (1974). Finite element analysis of fluid motion in a container, In: Finite
Element Methods in Flow Problem. UAH Press, Huntsville, Alabama, 737738.
J ohnson, N. L. (1996). The legacy and future of CFD at Los Alamos, Proceedings of
the 1996 Canadian CFD Conference, Ottawa, Canada.
Kershaw, D. S. (1978). The incomplete Cholesky-conjugate gradient method for the
iterative solution of systems of linear equations, Journal of Computational Physics,
26, 43-65.
Khayyer, A. and Gotoh, H. (2008). Development of CMPS method for accurate
water-surface tracking in breaking waves, Coastal Engineering Journal, 50(2), 179-
207.
Khayyer, A. and Gotoh, H. (2009). Modified moving particle semi-implicit methods
for the prediction of 2D wave impact pressure, Coastal Engineering, 56, 419-440.
References
199
Kim, J . K., Koh, H. M., and Kwahk, I. J . (1996). Dynamic response of rectangular
flexible fluid containers, Journal of Engineering Mechanics-ASCE, 122(9), 807-817.
Kim, M. H., Lee, S. M., Lee, J . M., Noh, B. J . and Kim, W. S. (2010). Fatigue
strength assessment of Mark-III type LNG cargo containment system, Ocean
Engineering, 37(14-15), 1243-1252.
Kim, Y., Nam, B. W., Kim, D. W. and Kim, Y. S. (2007). Study on coupling effects
of ship motion and sloshing, Ocean Engineering, 34(16), 2176-2187.
Kim, Y., Shin, Y. and Lee, K. H. (2004). Numerical study on slosh-induced impact
pressures on three-dimensional prismatic tanks, Applied Ocean Research, 26, 213-226.
Kim, Y., Shin, Y. S., Lin, W. M. and Yue, D. K. P. (2003). Study on sloshing problem
coupled with ship motion in waves, Proceedings of the Eighth International
Conference on Numerical Ship Hydrology, Busan, Korea.
Kleefsman, K. M. T., Fekken, G., Veldman, A. E. P., Buchner, B., Bunnik, T. H. J .
and Iwanowski, B.(2002). Prediction of green water and wave loading using a Navier-
Stokes based simulation tool, Proceedings 21th ASME Offshore Mechanics and Arctic
Engineering, American Society of Mechanical Engineers, paper 28480.
Koh, C. G., Mahatma, S. and Wang, C. M. (1994). Theoretical and experimental
studies on rectangular liquid dampers under arbitrary excitations, Earthquake
Engineering and Structural Dynamics, 23, 17-31.
Koh, C. G., Mahatma, S. and Wang, C. M. (1995). Reduction of structural vibrations
by multiple-mode liquid dampers, Engineering Structures, 17, 122-128.
Koh, H. M., Kim, J . K. and Park, J . H. (1998). Fluid-sturcutre interaction analysis of
3-D rectangular tanks by a variationally coupled BEM-FEM and comparison with test
results, Earthquake Engineering and Structural Dynamics, 27, 109-124.
Koshizuka, S., Nobe, A. and Oka, Y. (1998). Numerical analysis of breaking waves
using the moving particle semi-implicit method, International Journal for Numerical
Methods in Fluids, 26(7), 751-769.
Koshizuka, S. and Oka, Y. (1996). Moving-particle semi-implicit method for
fragmentation of incompressible fluid, Nuclear Science and Engineering, 123(3), 421-
434.
Koshizuka, S., Tamako, H. and Oka, Y. (1995). A particle method for incompressible
viscous flow with fluid fragmentation, Computational Fluid Dynamic Journal, 4(1),
29-46.
Larese, A., Rossi, R., Onate, E. and Idelsohn, S.R. (2008).Validation of the particle
finite element method (PFEM) for simulation of free surface flows, Engineering
Computations, 25, 385-425.
La Rocca, M., Sciortino, G., Adduce, C. and Boniforti, M. A. (2005). Experimental
and theoretical investigation on the sloshing of a two-liquid system with free surface,
Physics of Fluids, 17, 062101.
References
200
Lee, C. J . K., Noguchi, H. and Koshizuka, S. (2007a). Fluid-shell structure interaction
analysis by coupled particle and finite element method, Computers and Structures,
85(11-14), 688-697.
Lee, D. H., Kim, M. H., Kwon, S. H., Kim, J . W. and Lee, Y. B. (2007b). A
parametric sensitivity study on LNG tank sloshing loads by numerical simulations,
Ocean Engineering, 34, 3-9.
Lee, E. S., Moulinec, C., Xu, R., Violeau, D., Laurence, D. and Stansby, P. (2008).
Comparisons of weakly compressible and truly incompressible algorithms for the
SPH mesh free particle method, Journal of Computational Physics, 227, 84178436.
Lee, S. J ., Kim, M. H., Lee, D. H., Kim, J . W. and Kim, Y. H. (2007c). The effects of
LNG-tank sloshing on the global motions of LNG carriers, Ocean Engineering, 34,
10-20.
Lee, T., Zhou, Z. and Cao, Y. (2002). Numerical simulations of hydraulic jumps in
water sloshing and water impacting, Journal of Fluids Engineering, 124, 215-226.
Liszka, T. and Orkisz, J . (1980). The finite difference method at arbitrary irregular
grid and its application in applied mechanics, Computers and Structures, 11, 83-95
Liu, D. and Lin, P. (2008). A numerical study of three-dimensional liquid sloshing in
tanks, Journal of Computational Physics, 227(8), 3921-3939.
Liu, M. B., Xie, W. P. and Liu, G. R. (2005). Modeling incompressible flows using a
finite particle method, Applied Mathematical Modelling, 29, 1252-1270.
Lloyds Register (2005). Comparative sloshing analysis of LNG ship containment
systems.
Lloyds Register (2008). Guidance on the Operation of Membrane LNG Ships to
Reduce the Risk of Damage due to Sloshing.
Lo, E. Y. M. and Shao, S. D. (2002). Simulation of near-shore solitary wave
mechanics by an incompressible SPH method, Applied Ocean Research, 24(5), 275-
286.
Lohner, R., Yang, C. and Onate, E. (2006). On the simulation of flows with violent
free surface motion, Computer Methods in Applied Mechanics and Engineering, 195,
5597-5620.
Losasso, F., Talton, J . O., Kwatra, N. and Fedkiw, R. (2008). Two-way coupled SPH
and particle Level Set fluid simulation, IEEE Transactions on Visualization and
Computer Graphics, 14(4), 797-804.
Luo, Y. and Haussler-Combe, U. (2002). A generalized finite-difference method
based on minimizing global residual, Computer Methods in Applied Mechanics and
Engineering, 191, 1421-1438.
References
201
Mendez, B. and Velazquez, A. (2004). Finite point solver for the simulation of 2-D
laminar incompressible unsteady flows, Computer Methods in Applied Mechanics and
Engineering, 193, 825-848.
Mikelis, N. E. and J ournee, J . M. J . (1984). Experimental and Numerical Simulations
of Sloshing Behaviour in Liquid Cargo Tanks and its Effect on Ship Motions,
National conference on numerical methods for Transient and coupled problems,
Venice, Italy.
Ming, P. and Duan, W. (2010). Numerical simulation of sloshing in rectangular tank
with VOF based on unstructured grids, Journal of Hydrodynamics, 22(6), 856-864.
Mitra, S. and Sinhamahapatra, K. P. (2005). Coupled slosh dynamics of liquid filled
containers using pressure based finite element method, Exploring Innovation in
Education and ResearchiCEER-2005 Tainan, Taiwan, 1-5 March 2005.
Modi, V. J ., Akinturk, A. and Tse, W. (2003). A Family of Efficient Sloshing Liquid
Dampers for Suppression of Wind-Induced Instabilities, Journal of Vibration and
Control, 9(3-4), 361-386.
Modi, V. J . and Seto, M. L. (1997). Suppression of flow-induced oscillations using
sloshing liquid dampers: analysis and experiments, Journal of Wind Engineering and
Industrial Aerodynamics, 67-8, 611-625.
Molteni, D. and Colagrossi, A. (2009). A simple procedure to improve the pressure
evaluation in hydrodynamic context using SPH, Computer Physics Communications,
180, 861-872.
Monaghan, J . J . (1988). An introduction to SPH, Computer Physics Communications,
48(1), 89-96.
Monaghan, J . J . (1994). Simulating free surface flows with SPH, Journal of
Computational Physics, 110, 399-406.
Nomura, T. (1994). ALE finite-element computations of fluid-structure interaction
problems, Computer Methods in Applied Mechanics and Engineering, 112(1-4), 291-
308.
Nakayama, T. and Washizu, K. (1980). Nonlinear analysis of liquid motion in a
container subjected to a forced pitching oscillation, International Journal for
Numerical Methods in Engineerings, 15, 12071220.
Nakayama, T. and Washizu, K. (1981). The boundary element method applied to the
analysis of two-dimensional nonlinear sloshing problems, International Journal for
Numerical Methods in Engineerings, 17, 6311646.
Ockendon, H. and Ockendon, J . R. (1973). Resonant surface waves, Journal of Fluid
Mechanics, 59, 397-413.
Oger, G., Doring, M., Alessandrini, B. and Ferrant, P. (2006). Two-dimensional SPH
simulations of wedge water entries, Journal of Computational Physics, 213, 803-822.
References
202
Onate, E., Idelsohn, S. R., Zienkiewicz, O. C., Taylor, R. L. and Sacco, C. (1996). A
stabilized finite point method for analysis of fluid mechanics problems, Computer
Methods in Applied Mechanics and Engineering, 139(1-4), 315-346.
Onate, E., Sacco, C. and Idelsohn, S. R. (2000). A finite point method for
incompressible flow problems, Computing and Visualization in Science, 3(1-2), 67-75.
Osher, S. and Sethian, J . A. (1988). Fronts propagating with curvature-dependent
speed: algorithms based on HamiltonJ acobi formulations, Journal of Computational
Physics, 79, 1249.
Pal, N. C., Bhattacharyya, S. K. and Sinha, P. K. (1999). Coupled slosh dynamics of
liquid-filled composite cylindrical tanks, Journal of Engineering Mechanics, 491-495.
Pal, N. C., Bhattacharyya, S. K. and Sinha, P. K. (2003). Non-linear slosh dynamics
of liquid-filled laminated composite container: a two-dimensional finite element
approach, Journal of Sound and Vibration, 261, 729-749.
Pawell, A. (1997). Free surface waves in a wave tank, International Series in
Numerical Mathematics, 124, 311-320.
Perrone, N. and Kao, R. (1975). A general finite difference method for arbitrary
meshes, Computers and Structures, 5, 45-58.
Qian, L., Causon, D. M., Mingham, C. G. and Ingram, D. M. (2006) A free-surface
capturing method for two fluid flows with moving bodies, Proceedings of the Royal
Society, A8, 462(2065), 21-42.
Radovitzky, R. and Ortiz, M. (1998). Lagrangian finite element analysis of Newtonian
fluid flows, International Journal for Numerical Methods in Engineering, 43, 607-
619.
Ramaswamy, B. (1990). Numerical simulation of unsteady viscous free surface flow,
Journal of Computational Physics, 90, 396-430.
Ramaswamy, B. and Kawahara, M. (1987). Lagrangian finite element analysis
applied to viscous free surface fluid flow, International Journal for Numerical
Methods in Fluids, 7, 953-984.
Ramaswamy, B., Kawahara, M. and Nakayama, T. (1986). Lagrangian finite element
method for the analysis of two-dimensional sloshing problems, International Journal
for Numerical Methods in Fluids, 6, 659-670.
Rammerstofer, F. G., Scharf, K. and Fischer, F. D. (1990). Storage tanks under
earthquake loading, ASME Applied Mechanics Reviews, 43(1l), 261-283.
Rognebakke, O. F., Hoff, J . R., Allers, J . M., Berget, K. and Berge, B. O. (2005).
Experimental approaches for determining sloshing loads in LNG tanks, Proceedings
of the SNAME Maritime Technology Conference and Expo and Ship Production
Symposium, Houston Tx ETATS-UNIS, 113, 384-401.
References
203
Romero, J . A., Ramirez, O., Fortanell, J . M., Maritinez, M. and Lozano, A. (2006).
Analysis of lateral sloshing forces within road containers with high fill levels,
Proceedings of IMechE Part D: Journal of Automobile Engineering, 220, 303-312.
Shao, S. D. and Lo, E. Y. M. (2003). Incompressible SPH method for simulating
Newtonian and non-Newtonian flows with a free surface, Advances in Water
Resources, 26(7), 787-800.
Shankar, K. and Balendra, T. (2002). Application of the energy flow method to
vibration control of buildings with multiple liquid dampers, Journal of Wind
Engineering and Industrial Aerodynamics, 90, 1893-1906.
Shibata, K. and Koshizuka, S. (2007). Numerical analysis of shipping water impact on
a deck using a particle method, Ocean Engineering, 34(3-4), 585-593.
Shibata, K., Koshizuka, S. and Tanizawa, K. (2009). Three-dimensional numerical
analysis of shipping water onto a moving ship using a particle method, Journal of
Marine Science and Technology, 14, 214-227.
Solaas, F. and Faltinsen, O. M. (1997). Combined numerical and analytical solution
for sloshing in two-dimensional tanks of general shape, Journal of Ship Research,
41(2), 118-129.
Soulaimani, A. and Saad, Y. (1998). An arbitrary Lagrangian-Eulerian finite element
method for solving three-dimensional free surface flows, Computer Methods in
Applied Mechanics and Engineering, 162, 79-106.
Souto-Iglesias, A., Perez-Rojas, L. and Zamora Rodriguez, R. (2004). Simulation of
anti-roll tanks and sloshing type problems with smoothed particle hydrodynamics,
Ocean Engineering, 31(8-9), 1169-1192.
Souto-Iglesias, A., Delorme, L., Perez-Rojas, L. and Abril-Perez, S. (2006). Liquid
moment amplitude assessment in sloshing type problems with smooth particle
hydrodynamics, Ocean Engineering, 33(11-12), 1462-1484.
Stolbetsov, V. I. (1967). On oscillations of a fluid in a container having the shape of a
rectangular parallelepiped, Mekh Zhidk Gaza (Fluid Dynamics) N1, 67-76 (in
Russian).
Strandberg, L. (1978). Lateral stability of road tanks, National Road and Traffic
research Institute, Report No. 138A.
Su, T. C., Lou, Y. K., Flipse, J . E. and Bridges, T. J . (1982). A Numerical Analysis of
Large Amplitude Liquid Sloshing in Baffled Containers, US Department of
Transportation, Final Report, MA-RD-940-82046.
Sudharsan, N. M., Ajaykumar, R., Murali, K. and Kumar, K. (2004). A comparative
study of dynamic mesh update methods used in the simulation of fluid-structure
interaction problems with a nonlinear free surface, Proceedings of the Institution of
Mechanical Engineers 218 part C: Journal of Mechanical Engineering Science, 283-
300.
References
204
Sueyoshi, M., Kashiwagi, M. and Naito, S. (2008). Numerical simulation of wave-
induced nonlinear motions of a two-dimensional floating body by the moving particle
semi-implicit method, Journal of Marine Science and Technology, 13(2), 85-94.
Taylor, Sir Geoffrey (1953). An experimental study of standing waves, Proceedings
of the Royal Society of London, Series A, Mathematical and Physical Sciences,
218(1132), 44-59.
Tiwari, S. and Kuhnert, J . (2007). Modeling of two-phase flows with surface tension
by finite pointset method, Journal of Computational and Applied Mathematics, 203,
376-386.
Turnbull, M. S., Borthwick, A. G. L. and Eatock Taylor, R. (2003). Numerical wave
tank based on a -transformed finite element inviscid flow solver, International
Journal for Numerical Methods in Fluids, 42, 641-663.
Tveitnes, T., Ostvold, T. K., Pastoor, L. W. and Sele, H. O. (2004). A sloshing design
load procedure for membrane LNG tankers, Proceedings of the Ninth Symposium on
Practical Design of Ships and Other Floating Structures, Luebeck-Travemuende,
Germany.
Valentine, D. T. and Frandsen, J . B. (2005). Nonliear free-surface and viscous-
internal sloshing, Journal of Offshore Mechanics and Arctic Engineering, 127, 141-
149.
Veletsos, A. S. and Tang, Y. (1986). Dynamics of vertically excited liquid storage
tanks, ASCE Journal of Structural Engineering, 112(6), 1228-1246.
Viccione, G., Bovolin, V. and Pugliese Carratelli, E. (2008). Defining and optimizing
algorithms for neighbouring particle identification in SPH fluid simulations,
International Journal for Numerical Methods in Fluids, 58, 625-638.
Wang, C. Z. and Khoo, B. C. (2005). Finite element analysis of two-dimensional
nonlinear sloshing problems in random excitations, Ocean Engineering, 32(2), 107-
133.
Wang, W., Li, J . F., and Wang, T. S. (2006). Damping computation of liquid sloshing
with small amplitude in rigid container using FEM, Acta Mechanica Sinica, 22(1), 93-
98.
Wang, W., Wang, X., Wang, J ., and Wei, R. (1996), Dynamical behavior of
parametrically excited solitary waves in Faraday's water trough experiment, Physics
Letters A, 219, 74-78.
Webber, D. M., Gant, S. E. and Ivings, M. J . (2008). LNG source term models for
hazard analysis: a review of the state of the art and an approach to model assessment,
Health and Safety Laboratory Report, MSU/2008/24.
Wu, C. H. and Chen, B. F. (2009). Sloshing waves and resonance modes of fluid in a
3D tank by a time-independent finite difference method, Ocean Engineering, 36(6-7),
500-510.
References
205
Wu, G. X., Ma, Q. W. and Eatock Taylor, R. (1998). Numerical simulation of
sloshing waves in a 3D tank based on a finite element method, Applied Ocean
Research, 20(6), 337-355.
Wu, G. X. and Eatock Taylor, R. (1994). Finite element analysis of two-dimensional
non-linear transient water waves, Applied Ocean Research, 16, 363-372.
Yan, G., Rakheja, S. and Siddiqui, K. (2009). Experimental study of liquid slosh
dynamics in a partially-filled tank, Journal of Fluids Engineering, 131(7), 071303.
Yi, W. and Natsiavas, S. (1990). Seismic response of anchored fluid-filled tanks using
finite elements, Proceedings of ASME Pressure Vessels and Piping Conference, PVP-
191, 25-30.
Yoon, H. Y., Koshizuka, S. and Oka, Y. (1999). A particle-gridless hybrid method for
incompressible flows, International Journal for Numerical Methods in Fluids, 30(4),
407-424.
Zhang, S., Morita, K., Fukuda, K. and Shirakawa, N. (2006). An improved MPS
method for numerical simulations of convective heat transfer problems, International
Journal for Numerical Methods in Fluids, 51(1), 31-47.
Zhang, X., Sudharsan, N. M., Ajaykumar, R. and Kumar, K. (2005). Simulation of
free-surface flow in a tank using the Navier-Stokes model and unstructured finite
volume method, Proceedings of IMechE Part C: Journal of Mechanical Engineering
Science, 219, 251-266.
Zhang, X., Teng, B. and Ning, D. (2004). Simulation of fully nonlinear 3-D numerical
wave tank, China Ocean Engineering, 18(1), 59-68.

References
206

Appendix A
207
Appendix A: CD for animation files and explanation note
The CD contains the CPM simulation animations for selected numerical
examples. The video files for sloshing experiments in Chapter 5 are also included.
Detailed summary of the files are listed as follows.
File name Case name Section
Corresponding
figure
Dam_marin_w.avi Dam break with 5 . 0 / =
w
L d Sect. 4.5.2 Figure 4-29
Dam_marin_p.avi Dam break with 5 . 0 / =
w
Dam_marin_v.avi Dam break with 5 . 0 / =
w
Dam_Obstacle.avi Dam break with obstacle Sect. 4.5.3 Figure 4-37
Slosh_1_0.avi
Resonant water sloshing with
1 /
0
=
Sect. 5.3.1.1 Figure 5-29
Slosh_1_1.avi
water sloshing with
1 . 1 /
0
=
Sect. 5.3.1.2 Figure 5-34
Slosh_low.avi
Experiments of sloshing waves
in low-filling tank
Sect. 5.3.2
Figure 5-49 to
Figure 5-56
Slosh_breaking.avi
Experiments with sloshing
wave impact on the tank ceiling
Sect. 5.3.3 Figure 5-61

GaoM PDF

Uploaded by

Copyright:

Available Formats

GaoM PDF

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

GaoM PDF

Uploaded by

Copyright:

Available Formats

NUMERICAL SIMULATION OF LIQUID SLOSHING IN

RECTANGULAR TANKS USING CONSISTENT PARTICLE

You might also like