Amiot Baptiste
Ecole Polytechnique del’Université de Nantes &Okayama University

EXPLICIT-MPS METHOD
Development of accuracy and computational time

Table of contents
Overview ............................................................................................................ 2I.

E-MPS for tsunami inquiries ......................................................................... 3
What is the E-MPS method?.................................................................................................................................................. 3
Features and applications ................................................................................................................................................. 3
Brief state of art (bibliographic researches) ..................................................................................................................... 4
Equations and discretization ................................................................................................................................................. 5
Transport equations ......................................................................................................................................................... 5
Kernel function ................................................................................................................................................................. 5
Density .............................................................................................................................................................................. 6
Discretization .................................................................................................................................................................... 6
Programming workflow ......................................................................................................................................................... 7
Intermediary time step ..................................................................................................................................................... 7
Pressure calculation .......................................................................................................................................................... 8II.

Ways of improvements ............................................................................. 9
Accuracy of E-MPS ................................................................................................................................................................. 9
Pressure Gradient ............................................................................................................................................................. 9
Velocity Laplacian ............................................................................................................................................................. 9
Study on pressure function ............................................................................................................................................. 10
Computation time ............................................................................................................................................................... 10
Density linearization ....................................................................................................................................................... 11
Velocity Laplacian modified ............................................................................................................................................ 13
Pressure Gradient modified ............................................................................................................................................ 13

III.

Experiments and Results ......................................................................... 15

Hydrostatic .......................................................................................................................................................................... 15
Accuracy ......................................................................................................................................................................... 17
Computation time .......................................................................................................................................................... 19
Dam break ........................................................................................................................................................................... 20
Accuracy ......................................................................................................................................................................... 21
Computation time .......................................................................................................................................................... 23

Conclusion........................................................................................................ 24
Bibliography ..................................................................................................... 25
Annexe ............................................................................................................. 26

1

Overview
How is it doable to anticipate facility damages due toheavy rains or tsunamis? Japan is perhaps the most
famous country where this question needs to beinvestigated to save inhabitant lives. To illustrate thisstatement, one relevant example is the Tohoku tsunami
which destroyed the Fukushima power plant. Powerrealized was around 2 megatons of TNT, and sevenyears later we can still observe its effects on the japan
east coast. Lately, Okayama prefecture was affected byheavy rains on July 2018 which destroyed many facilities Picture 1 : Japan railroad destroyed due to riveras buildings, or bridges. In that case, overflow is less overflow in Hiroshima ©japantimes newspapers
powerful but mixed with gravity action and soil erosion, water begins to be hazardous.
To anticipate damages, the most convenient way is to make experiments. But conditions of disasters
are quite complicated to build, even in reduced scales. Values that we want to know, in terms of
mechanics view, are data pressures on layout buildings (or apparatus), subjected at water flows. So
that, scientists and engineers are trying now to create computing simulations which are perfectly
fitting with that kind of horrific events in a way to conceive safer equipment.
The major difficulty for this kind of simulations concerning boundaries.
In the case of floods or waves, we need to deal with the interface
between liquid and gas. CFD is commonly used in fluid mechanics
analysis, but free-boundaries are not well represented. Moreover, thesimple fact of using a Eulerian approach involves digital noises and thenefficiency is reduced. Here we choose to work on the E-MPS method,
Picture 3 : Dam-break model
because it’s using the Lagrangian approach,means that fluid is represented as a lot of particles. So that, simulation is
not curbed by digital noises and does not include any meshes.
Furthermore, computational time is reduced compares to grid-mesheulerian solutions.
Thus, the aim of this internship at Okayama University was to develop a
previously created E-MPS code. It has been testing on two models called:hydrostatic and dam-break problem as shown in the picture below.
These models are currently used, and experimental data are available tocheck simulation reliability. In the future, we would like to simulate
phenomenon in 3 dimensions so that we need to work on two sides.
Picture 2 : Dam breakFirst, we need to look for the efficiency of pressure representation, and experiment by Martin andthen, we want to obtain fast simulations by reducing computation time. Moyce (1952)

2I.

E-MPS for tsunami inquiries

What is the E-MPS method?
Features and applications
E-MPS acronym means Explicit-Moving Particle Semi-implicit, which is a direct link with theprevious method called MPS: Moving Particle Semi-implicit. These two methods do not introduce any
meshes, any nodes or that we call weak equations. Also, they use Lagrangian description by
discretization of strong form partials equations as Navier-Stokes equations. So we directly solve themby using particle interactions and introducing “Kernel function” as a weight function that allows
particles to only consider their neighbours in the calculation of interactions. It’s quite different fromanother solver as FEM, and one of the best argument to prefer Lagrangian MPS concerns freeboundaries.
Finite Difference or Finite Elements are efficient to solve partial equation problems whenmodel boundaries are fixed because mesh must not be recalculated at each calculation steps. Whenboundaries are moving, recalculating mesh is mandatory to bear efficiency, so computation time isincreasing and restricted our simulations. Applications of free boundaries are prevalent for twophasic experiments, such as phase shifts or atmosphere studies. To show how waves are running,
this kind of simulations are enough because they allow studies about fluid velocity and pressure field.
Tsunamis are bigger then shore waves, but they share the same movement equations so we knowthat E-MPS or MPS (SPH too) will work well.
What are the differences between MPS and E-MPS? This answer is giving by workflows
comparisons. Even they are sharing the same flow; they solve density, location and velocity, thenpressure for each particle; the last parameter is not calculated in the same manner. We also knowthat fluid movement is linked with pressure field, so any modifications in that processing will modify
the experiment. Instead of getting the pressure field by using the Poisson equation as MPS approach,the explicit way only uses the prior calculation of density to obtain pressure field. Easier and lesscomputation time is required but E-MPS must introduce the concept of weakly compressing fluid.
Density needs to be consistent from start to end.
In most cases, this assumption is true for waves studies and using E-MPS or MPS show
accurate results. The key parameter is, at this stage, the computational time required. And due toexplicit workflow, E-MPS wins. It is more comfortable to run an E-MPS study then MPS because theaccuracy in terms of pressure will be at least reach same results, also depending on the experiment,
E-MPS will be 3 or 5 times faster. But, E-MPS needs reflection to ensure compressing inquiries.
To conclude this part, we just want to precise what kind of experiments are available with EMPS approach. Because it is developed to simulate water flow, major experiment concerns
oceanographic development as fluid-structure interaction. But also, it is commonly used when waterflow is mixed with sediment, and some applications as river flow digging or landside prediction are
made.

3

Brief state of art (bibliographic researches)
Basically, papers about E-MPS, MPS or SPH have been writing in Japanese. As a beginner inthis fabulous language, it was a burdensome work to understand these papers. Nevertheless, few
ones are written in English and so I mainly used them.
As we previously said, simulations were studied in coastal engineering and two major
methods are still existing: SPH and MPS. SPH is used for compressible fluid problems and MPS for
incompressible ones. But E-MPS introduce few aspects of compressible fluid, so research must becompleted by inspecting the SPH model and use advances in E-MPS.
In October 2017, GOTOH and OKAYASU have published a tremendous paper named:”
Computational wave dynamics for innovative design of coastal structures” which is the best paper
referring to SPH functions applicated on MPS models at this moment. In this one, a lot of informationis shown and most of them were tried in our codes. Main problems of MPS were due to pressure
fluctuation, themselves caused by the Lagrangian nature of moving calculation. They designed someanswers based on a higher order scheme of interpolation, but also by reducing tensile instabilities.
Two of their developments provide interesting results -and are now used in our code: CMPS and HL.
CMPS is a slightly different way of programming pressure gradient and allows the conservation of
momentum. HL is another technic to programming Laplacian in a higher scheme and so it reduces
instabilities. Then the paper is written for MPS scheme and a lot of enhancements are made to the
Poisson pressure equation which is not used in E-MPS. But they are proposing some keys to solve
instabilities caused by boundaries like walls, which were also tried to improve the final accuracy of
our codes.
Our work on boundaries lead us to another paper wrote by MITSUME, YOSHIMURA,
MUROTANI and YAMADA, from the University of Tokyo. “Explicitly represented polygon wallboundary model for the explicit MPS method” (2015), is a research about programming boundaries
and especially how to deal with wall conditions. The major challenge concerning walls is to avoid
using unnatural repulsive forces to coerce particles in the tank. They are solving this aim bydeveloping virtual particles called “mirror particles” and are minimizing unnatural repulsive forces by
only using it when it’s mandatory, most of the time when particles reach a huge velocity. Mirror
particles technic (also called ERP for Explicit Polygon Representation), was tried but due tocalculation time increases, we don’t use it in our new codes. But it was inspirational, and I developed
some different solutions based on mirror particle programming.
To validate our codes, it is necessary to compare our results to experimental ones. First
studies were done by MARTIN and MOYCE in 1952 and they proposed dimensionless results
concerning dam-break model. Even if it appears to be a long time from now, these results are quite
accurate. With new sensors and equipment, KOSHIZUKA, TAMAKO and OKA proposed new data
which are slightly different from previous ones. Also, last research on this experimental subject wasmade by MATSUDA from Nagoya Institute of Technology, with a new model and bigger apparatus.
His thesis focused on geotechnics and complex situation where the ground is porous and includes soil
etc. It was too difficult to use this result because our goal was to obtain major improvements of timeconsuming and accuracy, but in the future, it could be useful to build our codes on models provided
by this thesis. In our case, I decided to compare dam-break results with experimental data of MARTIN
4and MOYCE work because they are not too different from KOSHIZUKA and al. ones, and our previous
codes checked with their results.
Last paper I need to quote is the one proposed by KOSHIZUKA, SHIBATA and MUROTANI.
“Introduction to particle methods” (2014) is written in full Japanese, except the formulas! Sometricks are proposed to reduce time-consuming, mainly by decreasing step of calculation, modifying
the workflow to avoid unnecessary loop etc. They are suggesting reducing the accuracy as long as it
stays in suitable error area in order to reach good results for less time. This includes streamlining
mathematical functions as linearized functions.

Equations and discretization
With all these statements said, we are entering now in some mathematical features, byintroducing movement equations, interaction formulas and discretized operators. This part will justremind the purpose of E-MPS, all the improvements made during my internship at Okayama
University are explained in detail after.

Transport equations
We can write conservative momentum on the fluid area, this is Navier-Stokes equation𝐷𝑢1= − ∇𝑃 + ν. Δu + f
𝐷𝑡
𝜌
With

(1)

𝐷𝑢
the differential of velocity, ∇ the gradient operator, ν is cinematic viscosity, Δ the Laplacian
𝐷𝑡

operator and f all the outside strength, here we only consider gravity force.
Then we can say that fluid keep is mass constant so that we write𝜕𝜌= −∇𝑢 = 0𝜕𝑡(2)With ∇ the gradient operator on velocity.
According to meshless features, Navier-Stokes equations are written in particles viewpoint as𝐷𝑢𝑖
1= − ∇𝑃𝑖 + ν. Δui + fi
𝐷𝑡
𝜌𝑖

(3)

It will be computed in that way, solving each Navier-stokes equation for each particle and thenswitch to another particle. By the way, to consider the action of j particle on the I one, it ismandatory to introduce an interacting particle model as Kernel functions.

Kernel function
We need to define an interaction model because fluid particle velocities and pressure arelinked with each other. But we understand that if particles are too far, it will be meaningless5

interactions between them. Kernel functions are introduced to weight these actions and formulaused here is:𝑟𝑒− 1 (0 ≤ 𝑟 < 𝑟𝑒 )𝑤(𝑟) = { 𝑟(4)0 (𝑟𝑒 < 𝑟)
From which 𝑟𝑒 is the lowest distance to share physicalproperties, its value is taken as 𝑟𝑒 = 3.5𝑙0 with l the space
between particles. 𝑟 is the current distance between
particles. It would be easier to show the purpose of theselater, but one thing that we must kept in mind is thefunction shape of w.
Lots of kernel functions are still existing, but in E-MPS weuse one which is equal to zero for faraway particles andFigure 1: x-axis is distance r[m] and y-axis is w(r) [UA]equal to infinity when they are closer. It is also a simple
equation because it appears in every computing things, so the simplest shape is also the fastest
processing.

Density
The lagragian approach means that density has to be defined for every particle, so equations
are defined as a summation of nearby neighbours.
𝑛𝑖 = ∑ 𝑤(|𝑥⃗⃗⃗⃗⃗𝑖𝑗 |)
𝑗≠𝑖

(5)

With ⃗⃗⃗⃗⃗
𝑥𝑖𝑗 = ⃗⃗⃗
𝑥𝑗 − ⃗⃗⃗
𝑥𝑖 , in all the report we simplify this script by 𝑥𝑖𝑗 .
Also, to show the incompressibility, we should write
𝜌𝑖 = 𝑚. 𝑁𝑖 =𝑚. 𝑛𝑖
∫𝑉 𝑤(𝑟)𝑑𝑉(6)

Assuming that each particle has the same mass, 𝑛𝑖 steady values ensure incompressibility aims. In
the following papers, we are using 𝑛0 as the first computation of density, and we will adjust current
density with 𝑛0 to validate incompressibility.

Discretization
Gradient

For the particle i, the gradient vector is defined regarding all the particles in its neighbourhood. A
weight is given by the kernel function, so for the particle i surrounded by particles j, the gradientvector of φ is discretized as:

6

〈∇φ〉𝑖 =φi + φj𝑑∑[𝑥 𝑤(|𝑥𝑖𝑗 |)]0𝑛|𝑥𝑖𝑗 |² 𝑖𝑗

(7)

𝑗≠𝑖

with 𝑑 the space dimension, 𝑛0 the fluid density.
As shown before, the gradient is only used to represent pressure. So φ is the pressure function in EMPS method.
LaplacianFor a time-dependent diffusion problem, it can be shown that the variance of the distribution
increases by 2𝑑𝜈𝛥𝑡 per time step 𝛥𝑡, so it leads to the following discretization:〈∇φ〉𝑖 =

2𝑑∑[φij . 𝑤(|𝑥𝑖𝑗 |)]
𝑛0 𝜆0
𝑗≠𝑖

(8)

With
𝜆0 =

∫𝑉 𝑤(|𝑥𝑖𝑗 |)|𝑥𝑖𝑗 |2 . 𝑑𝑉
∫𝑉 𝑤(|𝑥𝑖𝑗 |). 𝑑𝑉(9)Note that V is the area close to the studied particle.
Laplacian is only used with velocity calculation, so φ is the particle velocity in the E-MPS method.

Programming workflow
From a computing pointof view, the figure below shows
how the functions are running,
and we could appreciate all thestructures previously shown. We
now can introduce some otherphysical features that aremandatory to compile thesimulation.

Intermediary time step
To ensure stability, EMPS is working by mid-step
calculation in which first only
gravity is considered, thenpressure and viscosity are added.
In that way, step k goes to step7

k* thanks to𝑢∗ = 𝑢𝑘 + 𝑔∆𝑡(10)
For velocity, and for locations:𝑥 ∗ = 𝑥 𝑘 + 𝑢∗ ∆𝑡(11)Note that this mid-step cut the incompressibility off, particles are moving in any direction without
paying attention to the neighbourhood. Calculating density as 𝑛∗ shows a slight difference to 𝑛0 , thatis the weakly compressible fluid hypothesis. To keep up the pressure field, we compare the newdensity as the previous one which was supposed constant.

Pressure calculation
Another point is that pressure calculation involves mid-step calculation density 𝑛∗. As shownin the figure, the formula is𝑛∗
𝑃𝑖 (𝑛∗ ) = 𝜌𝑐 2 ( 0 − 1)𝑛(12)With 𝑐 2 a constant parameter which allowed consistency. Nevertheless, in this paper w,e will workon another experimental pressure calculation because the previous one has an unphysical behaviour
in some cases.

8

II.

Ways of improvements

Thanks to the bibliography, I was able to develop some tricks to increase accuracy anddecrease the computational time needed to solve hydrostatic and dam break systems. This part also
shows few advances that I’ve made, it concerns accuracy by research of better pressure fieldcalculation but also the computation time by switching from wall particle type as only one-layer
particle wall.

Accuracy of E-MPS
Pressure GradientDue to instabilities, the pressure gradient is always
modified. Even if our model works, sometimes pressures
are fluctuating, and so gradient discretization is notadequate. I’m proposing here to develop CMPS model toobtain a larger gradient and improve stability. As an
example, the figure shows how we can deal with thisproblem. Firstly, we need to determine for each particle,
which particle in the neighbouring area has the smallerpressure. According to the broad movements, smaller
values are close to each other. Added in discretizedgradient, they erase volatility by reducing the gap whenneighbour values are too big.
Consequently, the discretized gradient has a new shape:〈∇P〉𝑖 =

Min(i)

Min(j)
Reij

Figure 2: Difference between the two minus
particles will erase the pressure gap between I
and j when gradient will be calculated(Pi + Pj ) − (𝑃𝑖𝑚𝑖𝑛 + 𝑃𝑗𝑚𝑖𝑛 )𝑑∑ [(𝑥𝑖𝑗 𝑤(|𝑥𝑖𝑗 |))]
2
𝑛0
|𝑥 |
𝑗≠𝑖

(13)

𝑖𝑗

This is working well in the case of dam breaks but not in the case of hydrostatic. In this second one,
we are preferring to keep model previously built.

Velocity Laplacian
Khayyer and Gotoh(2010) proposed to consider the mathematical definition of
Laplacian, as a divergence of gradient in MPS.
〈∇2 𝑣〉𝑖 =𝜕𝑣𝑖𝑗 𝜕𝑤𝑖𝑗
𝜕 2 𝑤𝑖𝑗 𝐷𝑠 − 1 𝜕𝑤𝑖𝑗
1∑{
+ 𝑣𝑖𝑗 ( 2 +)}𝑛0𝜕𝑟𝑖𝑗 𝜕𝑟𝑖𝑗𝜕 𝑟𝑖𝑗
𝑟𝑖𝑗 𝜕𝑟𝑖𝑗

(14)

𝑗≠𝑖

Due to kernel function9𝑟𝑒− 1 (0 ≤ 𝑟 < 𝑟𝑒 )𝑤(𝑟) = { 𝑟
0 (𝑟𝑒 < 𝑟)(15)

(5 − 𝐷𝑠 )𝑟𝑒
1
∑
𝑣𝑖𝑗
3
𝑛0
|𝑟 |

(16)

(5 − 𝐷𝑠 )𝑟𝑒𝑣𝑗 − 𝑣𝑖
∑
3
𝑛0
|𝑟 |

(17)

Velocity Laplacian is derived as〈∇2 𝑣〉𝑖 =𝑗≠𝑖

𝑖𝑗

Thus, the calculation begins〈∇2 𝑣〉𝑖 =𝑗≠𝑖𝑖𝑗

Study on pressure function
To reduce instabilities due to pressure calculation, previous study worked on a new formula. I
choose to go on that way and find out the best parameters for each problem.
𝑃0𝑛∞ − 𝑛 𝑓𝑠∗− 1) (0 ≤ 𝑛𝑖 < 𝑛0 )(∞−𝑛𝑛𝑓𝑠 1 + 𝛥𝑛𝑖
(1 − 𝑛∞ ) ( 𝛥 ) − 1
𝑃(𝑛𝑖 ) =𝑃0
𝑛0 (𝑛∞ − 𝑛 𝑓𝑠 ) 𝑛𝑖∗
∗ ( 0 − 1) + 𝑃0 (𝑛𝑖 ≥ 𝑛0 )𝑛𝑓𝑠 1 + 𝛥
(𝑛∞ − 𝑛0 )2𝑛{(1 − 𝑛∞ ) ( 𝛥 ) − 1
With 𝑛∞ = 𝑛0 (1 + 𝛥) and 𝛥 is chosen by user.

(18)

This function is based on 3 experimental shapes (or points) we consider. Near the surface,pressure is reset to 0 and due to kernel function, density near the surface has a hyperbolic
progressing. Then for density equal to 𝑛0 , pressure must be set at 𝑃0 , it is primordial point. Finally,pressure is proportional to density when density is bigger than 𝑛0 .

Delta parameter must be chosen wisely because it could modify wave shapes and create
instabilities by removing impermeability. For example, when density is bigger than basic density,
choose a little value of delta increasing pressure and pressure gradient will be tremendous. If it waschosen to little, pressure will be continuous, and movement will be quite different from theexperiment. I decided to figure out the best parameter by comparison with experimental data.

Computation time
Even if we can decrease the computation time by using parallelization codes, we consider inthis paper that the algorithm could be improved by reducing useless calculation as wall particles
10calculation. In fact, wall particles have a key role because they are managing density andimpermeability, they have also an impact on the pressure of fluid particles. As we saw before, eachparticle is concerned by its neighbouring, so we can’t just remove wall particle, build a “one-particlewall” and then running the program like that. We must replace each role of wall particles and thencreate new particles that fit with the previous model.

Density linearization
Wall

In the first place, wall particles allow calculating
density for particle near walls. As shown on the figure,each particle along the wall side are weighted by kernel
functions and by summation we obtain the “wall” density.
In the case in which we remove all the particles, except
the first line, the summation will not be representative of
the wall weight. Overall density will be less important forparticles near walls instead of being the same.

Re
Fluid density

If we consider particles by their locations, weWall density
show that we can create a link between wall density anddistance to walls. Walls and not wall, due to edges in 2D
Figure 3: behavior of particle near walls(and 3D), it’s necessary to introduce different function in
the case of particles are constrained by 1,2 or 3 walls. In
Wallthis paper we will just talk about function for 2
Dimensions, to extend in 3 Dimensions, procedures are atleast equals.
Using previous simulations, I created a list inwhich we find: wall densities, distance to the nearest
particle for each orientation. Means that for edges, 2
distances are saved, one for wall and the other for thefloor. Then, particles were sorted by their distance and ifthey were considered as being near edges. Finally, I
created a lot of matrix system to figure out the most

Re

riw
?

rifloor?

Figure 4 : determination of near wall particle
and particle lack is pointed out.

convenient function to join all particles.
For particles near only one wall, I
decided to determine a function as
3
2
𝑍(𝑟𝑖𝑤 ) = 𝑎 × 𝑟𝑖𝑤+ 𝑏 × 𝑟𝑖𝑤
+ 𝑐 × 𝑟𝑖𝑤 + 𝑑(19)

Figure 5 : New formula for particle near single wall11

Which is fitting well with experimental points. The most important error on the value is near the wall
because our model is based on a discrete wall, not a continuous one. Also, to show this accuracy weobtain the following equation
𝐸𝑟𝑟𝑜𝑟(%) = max (𝑍 𝑒𝑥𝑝 (𝑟𝑖𝑤 ) − 𝑍 𝑝𝑜𝑙𝑦 (𝑟𝑖𝑤 )) = 5%
𝑍 𝑒𝑥𝑝 (𝑟𝑖𝑤 )(20)

Concerning edges, I wrote this matrix equation :𝐴×𝑥 =𝑏In which b is the wall density, x is a vector including constant parameter and A is a matrix picturing
the system. Given that we have two parameters called 𝑟𝑖,𝑤𝑎𝑙𝑙 & 𝑟𝑖,𝑓𝑙𝑜𝑜𝑟 , we want this kind of
equations𝑍(𝑟𝑖,𝑤𝑎𝑙𝑙 = 𝑟𝑊 , 𝑟𝑖,𝑓𝑙𝑜𝑜𝑟 = 𝑟𝐹 ) = 𝑥1 × 𝑟𝑊 2 + 𝑥2 × 𝑟𝑊 1 × 𝑟𝐹 1 + 𝑥3 × 𝑟𝐹 2 + 𝑥4 × 𝑟𝑊 + 𝑥5 × 𝑟𝐹 + 𝑥6
Our system is composed of m=15 000 equations, x has only 6 lines so we changed the A shape as :𝐴 = 𝑖𝑛𝑣(𝐴′ × 𝐴) × 𝐴′
And thus𝑥 = 𝑖𝑛𝑣(𝐴′ × 𝐴) × 𝐴′ × 𝑏
The matrix system is solved by
OCTAVE software.
I did some experiments about the mostconvenient shape for A system, but in theend, I figure out that the best function has
Figure 6 : Prior particle densities compares to new ones, calculated by linearized
to be a third polynomial degree.
functions

By comparison, the polynomial form
appears to be a good compromise and we should keep it to decrease calculation time broadly. Hereis the function used :3
2
𝑍(𝑟𝑊 , 𝑟𝐹 ) = 𝑥1 𝑟𝑊+ 𝑥2 𝑟𝑊
𝑟𝐹 + 𝑥3 𝑟𝑊 𝑟𝐹2 + 𝑥4 𝑟𝐹3 + 𝑥5 𝑟𝑊 2 + 𝑥6 𝑟𝑊1 𝑟𝐹 1 + 𝑥7 𝑟𝐹 2 + 𝑥8 𝑟𝑊 + 𝑥9 𝑟𝐹
+ 𝑥10

(21)

And x vector is⃗ = (−𝟐𝟕𝟓𝟐𝟖𝟏; 𝟏𝟓𝟓𝟖𝟏𝟖; 𝟏𝟐𝟗𝟑𝟏𝟖; −𝟏𝟏𝟓𝟏𝟎𝟓; 𝟑𝟎𝟎𝟕𝟏; −𝟏𝟗𝟕𝟖𝟖; 𝟏𝟕𝟕𝟎𝟎; −𝟕𝟑𝟗; −𝟏𝟎𝟒𝟕; 𝟑𝟑, 𝟔)𝒙

12

Velocity Laplacian modified
In addition to the High-Laplacian improvement which would concern all the fluid particle, wedeveloped a new wall formula, based on previous Japanese papers. Considering that wall particles
are still removed, it was mandatory to keep the effect of displaced particles.
〈∇2 𝑣〉𝑖 = 〈∇2 𝑣〉𝑖,𝑓𝑙𝑢𝑖𝑑 + 〈∇2 𝑣〉𝑖,𝑤𝑎𝑙𝑙(22)

Wall terms were firstly〈∇2 𝑣〉𝑖,𝑤𝑎𝑙𝑙 =2𝑑∑ [𝑣𝑖𝑗 . 𝑤(|𝑥𝑖𝑗 |)]
𝑛0 𝜆0

(23)

𝑗∈𝑤𝑎𝑙𝑙Considering velocity between particle i and j as the difference between velocity of i and velocity of j,
we showed that〈∇2 𝑣〉𝑖,𝑤𝑎𝑙𝑙 =2𝑑∑ [(𝑣𝑗 − 𝑣𝑖 )𝑤(|𝑥𝑖𝑗 |)]
𝑛0 𝜆0

(24)

𝑗∈𝑤𝑎𝑙𝑙Owing to the general movement of particle j, and without angular moment, velocity Laplacian due tothe wall is〈∇2 𝑣〉𝑖,𝑤𝑎𝑙𝑙 =2𝑑(𝑣− 𝑣𝑖 ) ∑ [𝑤(|𝑥𝑖𝑗 |)]𝑛0 𝜆0 𝑤𝑎𝑙𝑙

(25)

2𝑑
(𝑣
− 𝑣𝑖 ) × 𝑍𝑤𝑎𝑙𝑙 (𝑟𝑖𝑤 )𝑛0 𝜆0 𝑤𝑎𝑙𝑙

(26)

𝑗∈𝑤𝑎𝑙𝑙We reduced this equation as〈∇2 𝑣〉𝑖,𝑤𝑎𝑙𝑙 =Instead of calculating a big summation which requires to loop around the neighbouring area, thisfunction only operates one time. Calculation time is inevitably decreased, and this equation fits wellwith our new model of density.

Pressure Gradient modified
Pressure Gradient was quite difficult to deal with. As shown in the bibliography, manymethods have been tried and according to our previous models some of them are not available.
Mirror particles, for example, gave us interesting results but matrix processing asked too much time
and we must found a better solution.
Based on the accurate polygon wall representation, I chose to use the same approach basedon fictive wall particles, but without building them. Given that our model allocates pressure to allwall particles, which constrains fluid particles and allow them to stay behind the wall, I decided tomake “pressure points” at every location in grey on the figure. Accordingly, I allocate pressure basedon existing wall particles. Means that pressure points have not the same values as removed particles.
13

Wall

Due to this new pressure field, it is mandatory
to adjust pressure points by another weighted function.
This function has to take wall geometry into account.
Flat walls are easy to imagine in terms of internal
forces, so the pressure will be more important on the
ground and less important near the fluid boundary. In
the case of fast velocities, sometimes impermeability is
not safe, and I added external forces rely on particle
bound.
〈∇𝑃〉𝑖 = 〈∇𝑃〉𝑖𝑓𝑙𝑢𝑖𝑑 + 〈∇𝑃〉𝑤𝑎𝑙𝑙+ 𝑓𝑖𝑚𝑝
𝑖

Re
Fluid

(27)
Figure 7 : virtual calculation of pressure points

Previously, pressure gradient was〈∇P〉𝑖 =Pi + Pj𝑑∑[𝑥 𝑤(|𝑥𝑖𝑗 |)]0𝑛|𝑥𝑖𝑗 |² 𝑖𝑗

(28)

𝑗≠𝑖

And wall pressure gradient is now modified as:〈∇P〉𝑤𝑎𝑙𝑙=𝑖

Pi + Pj𝑑(1 + 𝑤1 + 𝑤2 + 𝑤3 )]∑ [(2 𝑥𝑖𝑗 𝑤(|𝑥𝑖𝑗 |)) ×
0
𝑛
|𝑥
|
𝑖𝑗
𝑗∈𝑤𝑎𝑙𝑙(29)

Where weight functions on the right part are less important than 1.

14

III.

Experiments and Results

In this part, we will compare the previous model and our new one. In addition, an
experimental study named “An Experimental Study of the Collapse of Liquid Columns on a Rigid
Horizontal Plane” (MARTIN and MOYCE, 1952) is used to show the accuracy in the case of a dambreak problem.

Hydrostatic
Hydrostatic problem is one way to show the accuracy of our simulation. Indeed, it allows
validating the shape of pressure representation because the viscosity term in Navier-Stokes formula
would be quite little in front of the pressure term.

50 particles

10 particles

Number of fluid particles:

648

Total number of particles:

1828

Initial spacing:

0.012 mDensity of the fluid:

1000 kg.m-3

Kinematic viscosity:

0.001 m2.s-1

Gravity:

9.8 m.s-2

A particle

38 particles
Figure 8 : Initial configuration of hydrostatic experiment in E-MPS previous model

Figure 8 shows the first shape of our experiment concerning hydrostatic. Green particles
have properties of “particle wall type”. Which means they have the same process of calculation thenpink fluid particles, but they are not allowed to move and their pressures are only calculated bysummation of fluid particles in their neighbourhoods. They have been built to ensure impermeability
and to constrain pink fluid particles in the tank.
Note two major things: firstly, kinematic viscosity is 1 000 times better than real viscosity.
This is to solve consistent enquiries in our models. Secondly, green particles are two times more
important than fluid particles. This is just to show the action of wall particle type on the computation
15

time. Also, the appropriate number of green particles depends on kernel function distances 𝑟𝑒 , andspacing between particles. Here 4 particles wall would be better in term of computation.

46 particles

Now our new model is quite little compared to the previous one, but we have the same number of
fluid particles and we deal with the black particles with properties of “wall particle”. Difference
between green particles and black ones concerning density calculation. Black wall particle does nottake time to calculate its density, it is only focused on pressure calculation.

Number of fluid particles:

648

Total number of particles:

759

Initial spacing:

0.012 mDensity of the fluid:

1000 kg.m-3

Kinematic viscosity:

0.001 m2.s-1

Gravity:

9.8 m.s-2

A particle
20 particles
Figure 9 : Initial configuration of hydrostatic problem in E-MPS, new model

To analyze the accuracy, we proposed to visualize pressure which occurs on the ground. As a current
hydrostatic problem, we know the theoretical shape of pressure which is :𝑃𝐴 = 𝜌𝑔 × (ℎ𝑓𝑟𝑒𝑒𝑆𝑢𝑟𝑓𝑎𝑐𝑒 − ℎ𝐴 )(30)

Figure 10 presents the priormodel at 2 different time
steps. While the beginningof the simulation shows agood pressure scattering,
the end of the simulation
appears to be more
ambivalent. Pressure fieldon the ground is dispersed
and also edge pressures areasymmetric. At the top of
the column, we observe
fluid particles repelled by
Figure 10 : Left : time is 0.1s ; Right : time is 1.5s

16the wall. Note that in this model, CMPS is not working.
Figure 11 is the new model with thesame steps than figure 10.
Beginnings are at least slightly
different from previous experiments,
we notice that the pressure field in
the centre appears to be more
important than nearby walls.
At 1.5s, the scattering of pressureshows relevant points. Firstly, it iswell-cut step by step, with highpressures on the ground and low
pressures at the fluid/gas interface.
Then edge pressure seems to be
more important than pressure forthe same height, but it is lessFigure 11 : Left: time is 0.1s; Right: time is 1.5s
significant compares to figure 10.
Only one or two particles are concerned. Finally, it appears that wall particle repels more fluid
particle than a prior experiment, this is still a problem due to interface treatment and pressure wall
representation.

Accuracy
The shape of the hydrostatic problem seems interesting, but it surely depends on set up
parameters. Pressure calculation needs a delta value to be launched, and it relies on empirical results.
As we know the theoretical pressure for the A-point, it is possible to look for the best empirical delta
which fit with the theory.

COMPARISON FOR GROUND PRESSURE(PARTICLE A) DEPENDING ON DELTA
Δ = 0,4

Δ = 0,15

Δ=0,1

Δ = 0,25

6500
6000

PRESSURE (PA)

5500
5000
4500
4000
3500
3000
2500
0

0,2

0,4

0,6

0,8

1

1,2

1,4

1,6

TIME (S)Graphic 1: Pressure of A-point depends on delta values.

17Graphic 1 shows pressure volatility, furthermore for some delta values the stability is not enough.
Delta values less than 0.15 are fluctuating with high frequency in addition of lower frequency. Over
0.15, high frequency is less important, and we must choose the best parameter by referring to thetheoretical value.

PRESSURE ON THE GROUND (PARTICLE A)
VERSUS TIME, FOR DIFFERENT VALUES OF
DELTA
Δ=1

Δ = 0,15

Δ = 0,25

Δ = 0,4

Δ=0,21

4300

PRESSURE (PA)

4200
4100
4000
3900
3800
3700
3

3,2

3,4

3,6

3,8

4

4,2

4,4

4,6

4,8

5

TIME (S)Graphic 2: another scale for A pressure analysis

After 3 seconds some new features are relevant, for example with delta value around 0.2, pressure
begins to fluctuate. Something interesting is the mean value according to Delta. When delta
increases, mean pressure is decreasing. However even oscillations are less important, mean values
are drifting away from the theoretical value which is supposed to be 4,2 kPa.
I chose to keep delta equals to 0.25 because it seems to be a good middle ground between
oscillations and mean pressure value. To validate this analysis, we propose to study the whole field ofpressure.

18

COMPARISON BETWEEN THEORICAL
PRESSURE AND PRESSURE SIMULATEDWITH DELTA = 0,25
Theorical pressure

Δ = 0,25

PRESSURE (PA)

5000
4000
3000
2000
1000
0
0

0,05

0,1

0,15

0,2

0,25

0,3

0,35

0,4

0,45

0,5

DEPTH (M)Graphic 3: Analysis of the whole pressure field in blue, compared to the theoretical value in green

Pressure field with delta equals to 0.25 is similar than theoretical pressure. We point theshape out which is nearby the real one. Nevertheless, the beginning of the function does not follow
this same proportionality. That means we can still improve the pressure formula.
With that being said, the next step of our work is to show computation time improvements.

Computation time
We remind that calculation time depends on previously developed clustering of neighbours,
so the figure submits this number. We must bear in mind that previous model was not well updated,so that time consuming will be irrelevant for the hydrostatic approach. The purpose of this study is toshow that a lot of wall particles are requesting a big amount of time.

Number ofcandidates
17
12Time for the Previous
model (s)Time for the New
model (s)Iteration

95,06
90,23

50,62
45,09

7500
7500

As we expected, the new model reduces from 45% to 50% processing time. Also, neighbouring
parameters do not modify this result tremendously.
We just have shown that our new one wall model is faster and give accurate results. The pressure
gradient is also well represented, we can now work on the dam break problem.

19

Dam break

45 particles

This time, previous work is well built in term of optimization, walls are 4 particles thick.
72 particles

Particle B

Number of fluid particles:
Total number of particles:
Initial spacing:
Density of the fluid:
Kinematic viscosity:
Gravity:648
1458
0.012 m1000 kg.m-3
0.001 m2.s-1
9.8 m.s-2

Particle C
Figure 12 : Initial experiment for dam-break problem

In comparison with our new model, the number of particles is 45% less important.

41 particles

72 particles

Particle B

Number of fluid particles:
Total number of particles:
Initial spacing:
Density of the fluid:
Kinematic viscosity:
Gravity:648
8020.012 m1000 kg.m-3
0.001 m2.s-1
9.8 m.s-2

Particle C

Figure 13: New model with one wall particle layout

As it was too big to stay in the main document, you could find the results and the comparison ofpressure field in annexe 1. To sum up, due to the initial pressure field, some particles are exploding
at the beginning, which creates some unnatural movements in the left tank area. Then wave
development seems to be fitting with the reality, even when the wave crash happened. After thistime, new wave shape is merging as expected. But, it is quite complicated to know if the shape is real
because particle behaviour nearby left edge is strange. We propose in the next part to verify theaccuracy of all these steps.

20

Accuracy
MARTIN and MOYCE provide experimental results concerning wave development for a dambreak problem. They carry out different shape experiments in order to define dimensionlessoutcomes.
So that, to contrast simulated results and know the better delta parameter, I made the same
dimensionless operations which are𝑡′ = 𝑡 ∗ √

𝐿′ =

2𝑔
𝐿0

(31)

𝐿
𝐿0

(32)

With initial water column height equals to 𝐿0

WAVE FRONT COMPARISON
n²= 1; 2,25 inch Z théol' (Δ=0,06)

l' (Δ=0,6)l' (Δ=0,3)l' (Δ=0,15)

4,5

DIMENSIONLESS LENGTH

4
3,5
3
2,5
2
1,5
1
0,5

1

1,5

2

2,5

3

DIMENSIONLESS TIME

Figure 14: Comparison between experimental values of way front advance and simulated ones which depend on delta
parameter

It appears that the best delta value is around 0.3. Underneath, waves are not running well; andabove 0.3 major inquiries are to deal with exploding particles in the near beginning. In these respects,
I decided to use 0.3 value for delta parameter in the following experiments.

21

C PRESSURE (PA)C pressure (Pa)
6000

PRESSURE (PA)

5000
4000
3000
2000
1000
0
0

0,2

0,4

0,6

0,8

1

1,2

1,4

1,6

TIME (S)Figure 15: Pressure of C-point versus time

Similarly, other results are available concerning pressure on the right side of the tank area. It is very
important for geotechnics studies and shape obtains should be the same as experimental ones. We
compare our outcomes with another experiment made by HU and KASHIWAGI (annexe 2) and it
appears that our simulation has some difficulties to deal with the first wave impact. The pressure
value is 3 times less important than expected. I think that in a one-particle wall, calculation ofpressure should be different, and it needs to be reviewed. By the way, the broad shape is quite thesame, with another second wave impact slightly different from the experiment.

B PRESSURE (PA)B pressure (Pa)
7000

PRESSURE (PA)

6000
5000
4000
3000
2000
1000
0
0

0,2

0,4

0,6

0,8

1

1,2

1,4

1,6

TIME (S)Figure 16: Pressure at B-point versus time22This other graphic is used to validate pressure exploding, in the near beginning. The maximum value
is reached after explosions when particles are repelled by gravity. Minimal one is arriving just after,
and it is lower than expected. Pressure near 0 means that it is a free surface particle, which is notentirely true after a one wave development. It is like water particles disappear from this area,pressure should be a little more important in that case. Then flux development is quite similar to theone we can imagine. It should be interesting to find some experimental values to compare with oursights.

Computation time
Now concerning simulation time, comparison with the previous model is interesting because it was
efficient in terms of time cost. The number of neighbours was optimized to keep accuracy and it wasproviding a good time-consuming contrast to previous models.
ITERATION OF CANDIDATE
LIST
8
12
16

NEW METHOD
(S)
49,13
49,45
49,47

PREVIOUS METHOD
(S)
76,22
76,50
80,05

NUMBER OF
ITERATION
7500
7500
7500

Outcomes are unexpected to be at this level. Gain is near 35% and consistency is saved. This one
particle wall provides better results without using hardware or soft parallelization which should be
used soon.

23

Conclusion
To sum up this paper, it appears that Lagrangian procedures are useful in the case of fluid
movements with free surface boundaries. I was able to add some enhancements by using physicallaws in discretized ways, but major problem was to keep consistency by using computing tricks as
unnatural forces etc... Modelling is still composed by the compromise between accuracy andcalculation speed, and this paper points out that we can improve speed without losing a lot in termsof accuracy by using linearization and moreover by defining the thin wall. This last enhancement wasdifficult to build because it should use computing tricks and put distance with “real” laws. After a lotof research, it was mandatory to imagine another way of calculation to provide better outcomes. I
tried some experimental model but the last one is for the moment the best one. Actually, it appears
to be limited and a new way of pressure test needs to be done to be used by comparison withexperimental results.
3 dimensions problems were not discussed here because it demands too much time to run. Thanks
to my new research, I hope that time consuming could be less important, but it would be needed tocreate linearized functions. Also, 3D models have not any experimental results which could be used
to validate the accuracy.
This internship was also a great opportunity to figure out what is a research work in the field of
numerical analysis for fluid mechanics. Wide guidelines are complicated to deal with, but I tried tofind some ways to answer them. The main problem was to understand how physical expectations arereached by computing, which was sometimes hard because MPS papers are most of the time writingin the Japanese language. Despite this complex environment, I learnt a lot from all the people I met
and I’m finally grateful to have lived this research internship in Okayama University Lab.

24

Bibliography
Gotoh, H., & Okayasu, A. (2017). Computational wave dynamics for innovative design of coastal
structures. Proceedings of the Japan Academy, Série B Physical and Biological Sciences, 525546.
Koshizuka, Shibata, & Murotani. (2014). Introduction to particle methods.
Ma, Q., Zhou, Y., & Yan, S. (2016). A review on approaches to solving Poisson’s equation inprojection-based meshless methods for modelling strongly nonlinear water waves. Journal of
Ocean Engineering and Marine Energy, 279-299.
Martin, J., & Moyce, W. (1952). An Experimental Study of the Collapse of Liquid Columns on a Rigid
Horizontal Plane.
Mitsume, N., Yoshimura, S., Murotani, K., & Yamada , T. (2015). Explicitly represented polygon wallboundary model for the explicit MPS method.

25

Annexe
ANNEXE 1: DAM-BREAK COMPARISON

Figure 17 : Left: Previous model; Right: New Model; time equals to 0.25 s

Figure 18 : Left: Previous model; Right: New model; time is 0.5s

Figure 19 : Left: Previous model; Right: New model; time is 0.75s

26Figure 20 : Left: Previous mode; Right: New model; time is 1s

Figure 21 : Left: Previous model; Right: New model; time is 1,25s

Figure 22: Left: Previous model; Right: New model; time is 1,5s

27

ANNEXE 2: HU AND KASHIWAGI EXPERIMENT

28