Proceedings of Coastal Engineering, JSCE, Vol.3, 2012, November

DEVELOPMENT OF A NOVEL 3D HIGHER
ORDER LAPLACIAN MODEL FOR ENHANCED
PREDICTION OF WAVE IMPACT PRESSURE
CALCULATION IN 3D MPS-BASEDSIMULATIONS
Abbas Khayyer1, Hitoshi Gotoh2 and Hiroyuki Ikari3
1Dept. of Civil and Earth Resources Engineering, Kyoto University, khayyer@particle.kuciv.kyoto-u.ac.jp
2Dept of Civil and Earth Resources Engineering, Kyoto University, gotoh@particle.kuciv.kyoto-u.ac.jp
3Dept of Civil and Earth Resources Engineering, Kyoto University, ikari@particle.kuciv.kyoto-u.ac.jp

A 3D higher order Laplacian model is proposed for enhanced wave impact pressure calculation by
the Moving Particle Semi-implicit (MPS) method. The Laplacian model is derived by meticulously taking
the divergence of a SPH (Smoothed Particle Hydrodynamics) gradient model and is then utilized fordiscretization of Laplacian of pressure corresponding to the Poisson Pressure Equation (PPE). The
enhancing and stabilizing effect of the 3D higher order Laplacian model is shown through simulations ofdesigned exponentially excited sinusoidal pressure oscillations and a schematic dam break with anobstacle.
Key Words : wave impact pressure, particle method, MPS, unphysical pressure fluctuations

1. INTRODUCTION
As a Lagrangian gridless method, the MPS(Moving Particle Semi-implicit; Koshizuka and Oka,1996) has been applied to a wide range ofengineering
applications,
includingcoastalhydrodynamic flows. Despite its wide range ofapplicability, the MPS method has a few major
drawbacks that may substantially affect itsperformance. Non-conservation of momentum(Khayyer and Gotoh, 2008), unphysical pressure
fluctuations (Khayyer and Gotoh, 2009, 2010) andnumerical instability (Khayyer and Gotoh, 2011) areamong the major drawbacks associated with theMPS method.
Through the past couple of years, the authors
have been working on enhancement of MPS methodby revisiting the derivation of differential operator
models and by proposing more accurate, consistentschemes, while trying to maintain the simplicity ofthe original method. Khayyer and Gotoh (2008)
proposed a Corrected version of MPS method,abbreviated as CMPS, characterized by ananti-symmetric pressure gradient model. Later, aHigher order Source term, abbreviated as HS, wasderived for enhancement of pressure calculation(Khayyer and Gotoh, 2009).
Another step towards enhancement andstabilization of pressure calculation by aprojection-based particle method is to apply a moreaccurate Laplacian model for discretization of
Laplacian of pressure in the Poisson Pressure
Equation (PPE). Khayyer and Gotoh (2010)highlighted the importance of the mathematicalconsistency of the Laplacian model and discretizedsource term of the PPE and derived a 2D Higher
order Laplacian model, abbreviated as HL, for theMPS method. In most cases, however, thementioned enhanced particle methods have been
applied to and verified by 2D calculations. On theother hand, most hydrodynamic flows areessentially three-dimensional. Hence, developmentof 3D accurate particle methods becomes
indispensable.
The main aim of this study is to develop a 3D
higher order Laplacian model for enhancement andstabilization of pressure field in 3D MPS-based
simulations. The 3D Laplacian model is derived bymeticulously taking the divergence of a commonly
applied SPH (Smoothed Particle Hydrodynamics)- 11 -

gradient model (Monaghan, 1992). It will be shownthat the Laplacian model derived in a 3D framework
differs to that corresponding to a 2D one and thatthe same approach can be applied for derivation of
consistent Laplacian models in the SPH context.
The enhancing effect of the 3D Laplacian model
will be shown by simulating designed exponentially
excited pressure oscillations (Khayyer and Gotoh,2010, 2011) and a schematic dam break with an
obstacle (Kleefsman et al., 2005).

 2 wij    wij (1)In 3D Cartesian coordinates, the gradients of ij andwij are expressed as:
ij wij ij rij
rij xij
wij rij
rij xij

ij riji

rij yiji

wij rij
rij yijwhere ij   j  i    ji
From Eq. (2): ij  wij jij rijjrij zijkwij rij
rij zij(2)krij x ij rij x ij



 ij rij wij rij
rij y ij rij y ij

 ij rij wij rij
rij z ij rij z ij
 ij wij x ij2
rij rij r2ij

 ij wij y ij2rij rij r2ij

 ij wij z ij2
rij rij rij2

 ij wij x ij2  y ij2  z ij2
rij rij
 ij wij
rij rij

rij

rij



rij

(5)

It should be noted that the only simplification inderivation of the above higher order Laplacianmodel (from Eq. 1) corresponds to the considered
first-order accurate finite difference approximationin Eq. (3), and all the other terms are exact
expressions. By considering the standard MPS
kernel, the higher order Laplacian in 3D would besimplified to:   i (3)On the other hand: 2 ji wij
 2 wij 2ij wij 1
 ij


n0 i j  rij rij
rij rij rij2

  2 wij 1   ij
n0 i j 
rij2 2 ji wij
rij

2Comparing Eqs (3) and (4) with theircorresponding equations in 2D framework (Eqs 11
and 12 in Khayyer and Gotoh, 2010), the only
difference corresponds to the expression of Eq. (4),where the effect of space dimensionality is appearedin the summation of the terms including second
partial derivatives of rij with respect to rij. This
would result in a different expression of Laplacianin 3D. From Eqs. (1), (3) and (4):rij2

 ji   ij wij

2wij  2 rij   i 

2(4)and ij   j  i .

 ij rij wij rij

 2 wij  rij  wij  2 rij  2 wij  rij  
rij2  xij rij xij2rij2  yij  2 wij  rij  wij  2 rij rij2  zij rij zij2
rij yij
 2 wij xij2  yij2  zij2 wij  3 xij2  yij2  zij2  rij  rij
rij2rij2rij3 2 wij 2 wij
rij2rij rij

By considering the Laplacian at a target particle
i as the divergence of the gradient calculated at thattarget particle and by applying the commonlyapplied SPH gradient model (Monaghan, 1992), theLaplacian at a target particle i would be formulated
as (Khayyer and Gotoh, 2010):
1 ij  wij  ij  2 wij
n0 i j  wij rij zij  rij zij 22. DERIVATION OF A 3D LAPLACIAN   i   wij rij    wij rij xij  rij xij  yij  rij yij  2 ij re 1
n0 i  j  rij3 (6)The 3D CMPS-HS with the above Higher order
Laplacian will be referred to as 3D CMPS-HS-HL.
Table 1 summarizes the derived higher order
Laplacian and its main comprising terms in 2D and3D.

- 12 -Table 1 Higher order Laplacian and its comprising terms

ij  wij

2D3D

2 ji wij

2 ji wij

rij rijrij rij

 wij
2 2 wij   irij2



1 wijrij rij  2 wij ij wij 1

 ij


n0 i  j rij rij 
rij2

 2 wijrij2

2 wijrij rij  2 wij 1

 ij

n0 i  j 
rij2 Fig. 1 Snapshots of water particles and pressure field

3. VERIFICATIONSThe enhancing effect of the 3D higher orderLaplacian model is demonstrated by performing two
numerical tests, namely, designed exponentially
excited sinusoidal pressure oscillations (Khayyerand Gotoh, 2010, 2011) and a schematic dam breakwith an obstacle (Kleefsman, 2005). To make thecomparisons more comprehensible, the fluid isconsidered to be inviscid so that the stabilizing
effect of the viscous forces will be entirely omitted.
Fig. 2 Time histories of pressure at measuring point A (Fig. 1)(1) Designed Exponentially Excited Sinusoidal
Pressure Oscillations
Hydrostatic pressure calculations with designedexponentially
excitedsinusoidalpressure
oscillations have been carried out to verify theenhanced performance of 3D CMPS-HS-HL method.
The calculation conditions are as follows: d0 =particle diameter = 1.0E-2 m and t = 1.0E-3 s. The
water tank is a cube of 0.18 m with the water depth
being 0.16 m. The modified gravitational
acceleration is defined as:
2 tgt) exp(0.03 )g d  g  Sin (2TT(7)where gd represents designed external accelerations
applied to fluid particles; g is the earth’sgravitational acceleration (= 9.81 m/s2); t stands
for the simulation time and T denotes the period ofsinusoidal term variations (= 0.02 s).
Fig. 1 shows the snapshots of water particles
together with the pressure field at t = 0.245 s. At thisinstant, the spatial distribution of pressure by the 3D
CMPS-HS-HL method appears to be almost
consistent with a hydrostatic pressure distribution.
In contrast, the 3D CMPS-HS method has depicted
an unphysical spurious pressure distribution
characterized by several zero-pressure inner fluidparticles.

Irregular Initial Distribution of Particles

Fig. 3 Time histories of pressure at measuring point A (Fig. 1)
for irregular initial distribution of particles

Fig. 2 depicts the time histories of pressure at a
fixed point (measuring point A in Fig. 1) by 3D
CMPS-HS and 3D CMPS-HS-HL methods. In spiteof relatively accurate approximations at thebeginning of the calculation, the pressure
calculation by the 3D CMPS-HS has become
destabilized after a few periods of designed
oscillations. The 3D CMPS-HS-HL method, on theother hand, has provided a notably more accurate
and stable pressure calculation.
Fig. 3 illustrates the time variations of calculatedpressure at point A when the initial distribution ofparticles is randomly altered and half of fluid
particles are displaced by  0.05d0 in both x and ydirections. As a result of an irregular distribution of- 13 -180.0P7 P5P4P8 P6
P3 P2m
99

P1m

0m

0.2

0.0water

H3H20.55 m

092

0.96.49m

68

0.16 m0.40 m1.00 m1.2m

0.30 m64

0.6m

1.128mm6m

0.1Fig. 4 Schematic sketch of domain - dam break with obstacle

particles (an initial compression or expansion andviolation of fluid incompressibility), the pressure
results by both methods appear to deviate from theanalytical solution at the beginning of thecalculation. This deviation tends to be smooth andrelatively stable in case of the 3D CMPS-HS-HLmethod. In contrast, the pressure trace by the 3D
CMPS-HS method is initiated by frequent
large-amplitude unphysical oscillations. After one
period of oscillation, the results by both methods
seem to become closer to the analytical solution.
However, the 3D CMPS-HS results still appear to
be inaccurate and inconsistent, particularly, in thevicinity of the peak and trough points. In contrast,
the 3D CMPS-HS-HL tends to provide a relatively
stable and accurate pressure calculation and thegrowth of instabilities tends to be less pronounced
when the HL scheme is incorporated.
(2) A Schematic Dam Break with an Obstacle
Schematic dam break flows have beencommonly used for the verification of numericalmethods. An experiment on a schematic dam breakflow and its impact against an obstacle, carried out
at the MAritime Research Institute Netherlands(MARIN), has been considered as a benchmark test
for the validation of both 3D particle methods and
3D grid-based methods. A schematic sketch of thecalculation domain including the positions of thewave height probes and pressure sensors installed in
the experiment is depicted in Fig. 4. Detailed
descriptions of the experiment have been provided
by Kleefsman et al. (2005). For all the calculations
performed in this section, d0 is set to be 2.0E-2 m,
resulting in a total number of 215940 particles.
Fig. 5 presents a qualitative comparison in
between 3D CMPS-HS and 3D CMPS-HS-HL byillustrating snapshots of water particles together
with the pressure field at t = 0.49 s. From this figure,it is evident that application of the new LaplacianFig. 5 Qualitative comparison in between 3D CMPS-HS and 3D
CMPS-HS-HL at t = 0.49 s - dam break with obstacle

model has significantly reduced the existing
numerical noise seen in the 3D CMPS-HS snapshot.
Further, the 3D CMPS-HS-HL has also resulted in amore integrated water jet.
Fig. 6(a-d) presents a qualitative comparison in
between simulation and experiment at t = 0.56 s.
Comparing Fig. 6(a-b) with (c), the 3D
CMPS-HS-HL has provided a more integrated jet.
Moreover, from Fig. 6(a) and (c), the height of thejet and other macroscopic features of flow, e.g.
backward curl of the jet towards the incoming flow,
appear to be in acceptable qualitative agreementwith the experiment. From Fig. 6(d), the high order
VOF method of Kleefsman (2005) has resulted in astraight fully integrated jet and the backward curl
has not been reproduced well.
Fig. 7(a-b) depicts time histories ofexperimental and calculated pressures at measuring
points P1 and P5 (Fig. 4). The figure shows thestep-by-step improvements in pressure calculation
by incorporating the higher order source term andthe higher order Laplacian model. Compared withthe 3D CMPS-HS, the 3D CMPS-HS-HL method
has resulted in a less-fluctuating, more-accuratepressure trace in a better agreement with theexperiment. In particular, a more consistentcalculation of the peak pressure rise and declination(Fig. 7(a)) is illustrated by the 3D CMPS-HS-HL
method. The high-order VOF method of Kleefsman
et al. (2005) has resulted in quite smooth variationsof pressure. Nevertheless, the pressure traces by thishigh-order grid-based method are characterized bysome instabilities particularly at the impact instants.

- 14 -Fig. 8 Time histories of water surface height at probe H3 (Fig.
4) - dam break with obstacle

Fig. 8 depicts the time variations of water
surface height at probe H3 and shows thestep-by-step enhancing effects of the higher ordersource term and the higher order Laplacian inproviding smoother and more accurate time
variations of water surface height.

4. CONCLUDING REMARKS

Fig. 6 Qualitative comparison in between 3D CMPS-HS, 3D
CMPS-HS-HL and experiment and VOF (Kleefsman etal., 2005) at t = 0.56 s - dam break with obstacle(a)

A 3D higher order Laplacian model is proposed
for further enhancement of wave impact pressure
calculation in 3D particle-based simulations. The
new Laplacian is derived by meticulously taking thedivergence of a particle-based gradient model. The
enhancing effect of the 3D Laplacian model isshown by simulating designed sinusoidal pressure
oscillations together with a schematic dam breakwith an obstacle.

P1

REFERENCES

(b)

P5

Fig. 7 Time histories of calculated and experimental(Kleefsman, 2005) pressure at measuring points P1and P5 (Fig. 4)

Koshizuka, S. and Y. Oka (1996): Moving particle
semi-implicit method for fragmentation of incompressible
fluid, Nuclear Science and Engineering, 123, 421-434.
Khayyer, A. and H. Gotoh (2008): Development of CMPS
method for accurate water-surface tracking in breaking
waves, Coastal Eng. J., 50(2), 179-207.
Khayyer, A. and H. Gotoh (2009): Modified Moving Particle
Semi-implicit methods for the prediction of 2D wave impact
pressure, Coastal Eng., 56(4), 419-440.
Khayyer, A. and H. Gotoh (2010): A Higher Order Laplacian
Model for Enhancement and Stabilization of Pressure
Calculation by the MPS Method, Applied Ocean Research,
32(1), 124-131.
Khayyer, A. and H. Gotoh (2011): Enhancement of stability andaccuracy of the Moving Particle Semi-implicit method,Journal of Computational Physics, 230(8), 3093-3118.
Monaghan, J.J. (1992): Smoothed particle hydrodynamics, Ann.
Rev. Astron. Astrophys, 30(1), 543-574.
Kleefsman, K.M.T., Fekken, G., Veldman, A.E.P., Iwanowski
B. and B. Buchner (2005): A Volume-of-Fluid based
simulation method for wave impact problems, Journal ofComputational Physics, 206, 363-393.
Kleefsman, K.M.T. (2005): Water Impact Loading on Offshore
Structures, A Numerical Study, Ph.D. thesis, University ofGroningen, The Netherlands, November.
(Received June 15, 2012)- 15 -