A Mesh-Free Particle Method for Simulation of Mobile-Bed Behavior Induced by Dam Break

The mesh-free particle method is the latest method in the computational fluid dynamics field. It has become well-known due to its application for simulating any kind of interface deformation and fragmentation, including large disturbance on the interface. In this study, the moving particle semi-implicit (MPS) method has been developed for multi-phase to simulate mobile (erodible)-bed under violent flows such as a dam break. In order to develop the MPS method, multi-phase particle interaction models including self-buoyancy and surface tension models were employed. All of the newly adopted particle interaction models were modified by the author. When fluids collide with the soil bed, the erodible bed can be liquefied. In this regard, polyethylene oxide (PEO) solution was used as a mobile bed. The water column captured by the gate collapsed after opening the gate, and it then propagated over a mobile bed. The mobile bed has various thicknesses based on corresponding experiments. Moreover, the dam break on an inclined tank was also simulated to investigate the swash uprush problem. The maximum velocity and height of mixed flows were measured at a certain point, which aligned with the corresponding experiment. All of the simulation results were compared with those obtained from some experimental work.


Introduction
Due to the highly concentrated energy, the flows induced by a dam break are very violent and include a large deformation of free surface and non-linearity.When it propagates along the river and channel, a severe flood and structural damage can follow.In this regard, the dam break problem has been theoretically, experimentally, and numerically studied by Chanson et al. and Liao et al. [1,2].Furthermore, Janosi et al. investigated the soil bed behavior due to the flows induced by the dam break [3][4][5].In the experimental and theoretical method, the dynamics of strong turbulence at the free surface and wave profiles were investigated by Brocchini et al. [6,7].Postacchini et al. outlined sediment transport due to dam break swash uprush numerically [5].In Reference [5], the collision of the collapsed water column to the soil bed in the inclined tank was investigated both experimentally and numerically by Postacchini et al.Despite many efforts to analyze it, its complication and non-linearity cannot be solved easily.The complexity increases dramatically especially when the flows collide with the structures or movable soil bed.Therefore, accurate analysis is very important for estimating consequences and preparing a counter measure.
Computational fluid dynamics (CFD) is an effective method for simulating the violent flows with the large deformation of a free surface including the dam break problem.Conventionally, CFD methods are based on the Eulerian approach, which uses mesh/grid to calculate flows and flux.In this method, flux across a fixed grid is used to calculate velocity and pressure of flows.However, several problems still exist in relation to expecting large deformation of the free surface.In this regard, the method Appl.Sci.2018, 8, 1070 2 of 15 of volume of fluid (VOF) and the level set are implemented into the programme [8].Despite using these special treatments, it has been difficult to simulate the complicated phenomena due to large deformation, fragmentation, and coalescence in order to solve those problems [9].Another problem of the conventional grid-based CFD method is a numerical diffusion induced by an advection term of the governing equation [10,11].Juez et al. detailed heterogeneous erodible beds under transient flows by using the Hirano model and the Exner equation [12,13].In References [12,13], two-dimensional (2D) morphodynamics was solved by a set of two equations under the assumption of influenced thickness and roughness of the bed due to time-varying surface textures by Juez et al.Those two investigations contribute and extend the possibility of having heterogeneous beds.The particle method has been used to investigate solid object transportation under large flows.The large floating body transport in the flash floods problem was investigated by Amicarelli et al. [14].Albano et al. applied the three-dimensional smoothed particle hydrodynamics (SPH) method on three-dimensional (3D) solid body transport in free surface flows [15].In their work, the weakly compressible-SPH model for the main flow implemented a free-slip boundary condition on the solid body.The moving particle semi-implicit (MPS) method was also employed to simulate solid object movement under large waves.Shakibaeinia and Jin investigated a mobile-sea bed problem by using the MPS method, and introducing the particle interaction model between two different types of fluids [16].
The particle method is based on the fully Lagrangian approach CFD method.In the particle method, all of the material, including the fluid, is comprised of particles.Each particle travels through the computational domain carrying physical phenomena including position, velocity, acceleration, and pressure.In this regard, the non-linearity, including the large disturbance of the free surface, sharp coalescent, and fragmentation, can be handled straightforwardly by tracking each particle.Another advantage of the particle method is the advection term in the Navier-Stokes equation, which is one of the governing equations that cannot be used.Therefore, numerical diffusion such as losing fluid volume can be resolved.In the particle method, SPH and the moving particle semi-implicit (MPS) method are well-known methods.SPH were introduced by Monaghan [17].Another well-known particle method is MPS, which was originally introduced by Koshizuka et al. [18].In the SPH, the pressure on the particle can be calculated by the equation of state.In this regard, the calculated time of simulation is faster when compared with MPS.In the MPS, the pressure can be implicitly estimated by the Poisson Pressure Equation (PPE).As a result, the fluctuation of pressure was suppressed, and it led to a more accurate pressure calculation while its calculation time is much more than SPH.Another difference between SPH and MPS is the measurement of influences of neighboring particles.SPH is based on continuous functions to approximate the physical quantities, while MPS uses the functionalized particle interaction model to solve gradient and diffusion.More comparative study can be found [19].In this regard, MPS may have a greater likelihood to simulate a sharp interface under large deformation and fragmentation.Tanaka et al. found that the MPS method was developed by applying a multi-source term of PPE to reduce pressure fluctuation [20].Lee et al. modified multisource terms of PPE by introducing a relaxation parameter [11].The surface tension and buoyancy adjustment models for multi-phase flows were investigated by Nomura et al. [21].The surface tension model was modified by Kim et al. [22] with more robust physics.
In this study, the newly developed MPS method for multi-phase flows by adding multi-phase particle interaction models was used to simulate the mobile-bed behavior induced by the dam break.All of the new particle interaction models, including buoyancy correction and surface tension models, were modified by the author with a robust method.The water column was closed at the gate, and it collapsed spontaneously by gravitational force when the gate was open.After the gate, a mobile-bed with an arbitrary thickness was located.In this study, the PEO solution was used instead of soil.According to study of Postacchini et al., it is possible to simulate PEO solution as a fluid particle due to liquefaction [5].Therefore, multi-phase MPS is well-matched for simulation.For the inclined tank, the maximum height and velocity of waves of water-PEO solution mixed flows were measured at a certain point.Furthermore, for the plat tank, which is not inclined, the snapshots were compared Appl.Sci.2018, 8, 1070 3 of 15 with experimental results.Both comparisons show that they align well between the experimental and numerical simulation.

Governing Equations
The continuity and Navier-Stokes equations are governing equations of the MPS method for incompressible viscous flows.
where ρ is the density of a fluid, t is time, → u is a velocity of a particle, P is pressure, ν is kinematic viscosity, → g is gravitational acceleration, σ is surface tension coefficient, κ is the curvature of an interface, → n is a unit normal vector for an interface, ∇ is a gradient, ∇ 2 is a Laplacian, and → F is an acceleration of external force.For the term of the surface tension, it automatically vanishes when the particle is not on the interface or free surface.
Since the MPS method follows the Lagrangian approach, the differential operator needs to be replaced by the particle interaction method.In the particle interaction method, the kernel function was introduced, which represents the effects of neighboring particles with respect to the distance between the center and neighboring particle.It was originally suggested by Koshizuka and Oka [18] and modified with a third-order function.
where subscript ij denotes the substitution of arbitrary physical quantities between the center particle and neighboring particle (r ij = r j − r i ), and r e is an effective range of the particle interaction.According to Equation (3), the value of the kernel function will be zero when the particle distance exceeds an effective range.In this study, the default value of the effective range is set to 2.1 × l 0 , where l 0 is the particle distance at an initial arrangement.

Particle Interaction Model
As stated above, all of the differential operators in the governing equation should be replaced by the particle interaction model.The gradient model calculating the averaged local weight of the gradient vectors is expressed below: where denotes the particle interaction model, φ is an arbitrary physical quantity, d is a number of dimension, and n 0 is a particle number density at the initial stage, which will be discussed later.
The divergence model represents the divergence of the physical quantity between the center and the neighboring particles.

∇ • φ
The Laplacian model represents the diffusion of a fluid below. where: λ is an adjustment parameter to equal the Laplacian value to an analytical value.
The particle number density corresponds to the density of fluid in the MPS method.It can be calculated by the summation of the kernel function.
MPS follows the Lagrangian approach.Therefore, the particle without physical volume is used for simulation, while the conventional grid method uses grid as a volume of fluid.Without a physical volume or a very small volume on the particle, the fluid density may increase to an infinite value.In order to avoid this issue, the particle number density has been used instead of the density of the fluid.Moreover, the particle number density can be calculated by adding the neighboring particle with respect to kernel function.Therefore, the MPS method can avert the violation of continuum and manage the effects of neighboring particles.
In order to simulate incompressible flows, the incompressibility model is employed, which has a similar algorithm to the simplified marker-and-cell (SMAC) method in the mesh-based CFD system.The incompressibility model consists of two stages, which include an explicit and an implicit stage.In the explicit stage, the intermediate position and velocity of each particle can be calculated from the viscosity and gravity directly.According to the continuity equation, which is one of governing equations, the particle number density at the intermediate stage, n * , may differ to it at the initial arrangement, n 0 .It might also lead to the dissatisfaction of the continuity equation.Therefore, the particle number density at an intermediate stage should be adjusted implicitly.
where n is the corrective particle number density, and it is related to corrective velocity, u , which can be obtained through the following mass conservation equation: The corrective velocity can be derived by the following implicit pressure gradient term: By recalling Equations ( 9)- (11), a Poisson equation for pressure calculation can be obtained: The left-hand side of Equation ( 12) can be discretized by the Laplacian model, which is given in Equation (7); its right side is expressed by the deviation of the particle number density with the constant value.These simultaneous linear symmetric matrix equations can be solved by the incomplete Cholesky conjugate gradient method.By substituting φ by P n+1 in Equation (5), which is the gradient model, the pressure gradient can be obtained.Equation ( 12) is used as a source term for PPE.However, an exaggerated fluctuation of pressure took place.Reference [14] introduced the multi-source term to solve this problem.By adding a divergence of velocity in Equation ( 12), the pressure fluctuation was suppressed.For this treatment, Lee et al. applied a relaxation parameter to distribute the weight of each terms.The source term of PPE used in this study is expressed below: where is the intermediate velocity, which is the temporal quantities obtained in the explicit stage, P n+1 is pressure at time n + 1, and γ 0 is the relaxation parameter.
The kernel function is a key factor of MPS, but it has an omni-spread function.It may occur under the distinct pressure head, and the buoyancy force can be underestimated at multi-phase flows.It can be solved by applying a self-buoyancy force to the center particle.Self-buoyancy is an adjustment of underestimated buoyant values.This term was applied to the velocity corrector term: where w B is the kernel function for the buoyancy-correction model, z is a vertical position of a particle, and C B is a velocity-adjustment parameter.w B is the same equation as Equation ( 3), but its effective range is set to 1.5, which was obtained from a numerical test.The reason for using C B is due to the self-buoyancy term remedy being an underestimated buoyancy force.Therefore, the buoyancy correction term should be adjusted by a proper parameter.The value of C B can be chosen by a numerical experiment.The buoyancy correction model was suggested by Shirakawa et al. [23], and it has been modified by the author using a physical formulation.
In general, the surface tension is not a strong element to simulate the flow.However, it is not an ignored component for the interface of multi-phase flows.The fourth term on the left side of Equation ( 2) is the surface tension term.The particle method for tracking each particle is advantageous for a very large deformation of the free surface and interface.Yet, the weakness of the specialty lead particle method makes it difficult to indicate the interface layer as well.Therefore, the special treatment to define a curvature of the interface is introduced by Nomura et al. [21].For the treatment, two additional particle number densities are employed. where where r st e is the effective range of the surface tension, and it is set to 1.5.With these newly added particle number densities, the curvature of the interface can be calculated by using Equation (16).
After calculating the interface curvature, a normal vector of the interface needs to be estimated.Nomura et al. suggested calculating the normal vector of interface by using particle number density [21].Although using particle number density follows the Lagrangian approach, it may have a problem near the wall region because, over the wall, there are no physical particles.In this regard, the author suggested the new method to calculate the normal vector of interface.Figure 1 shows a schematic model for the unit normal vector of the interface.For the center particle, the normal vector can be calculated by adding the distance between the center and the neighboring particles using the decomposed direction.Then, the value summation is divided by the absolute value of the normal vector.At this moment, all of the particles involving this calculation should be an interface particle.The unit normal vector can be calculated via the equations below: where x n and z n are the horizontal and a vertical distance between the center and neighboring particles, respectively.
After calculating the interface curvature, a normal vector of the interface needs to be estimated.Nomura et al. suggested calculating the normal vector of interface by using particle number density [21].Although using particle number density follows the Lagrangian approach, it may have a problem near the wall region because, over the wall, there are no physical particles.In this regard, the author suggested the new method to calculate the normal vector of interface.Figure 1 shows a schematic model for the unit normal vector of the interface.For the center particle, the normal vector can be calculated by adding the distance between the center and the neighboring particles using the decomposed direction.Then, the value summation is divided by the absolute value of the normal vector.At this moment, all of the particles involving this calculation should be an interface particle.The unit normal vector can be calculated via the equations below: where In the MPS method, the kinematic free surface and interface boundary condition can be automatically satisfied by tracing the particle directly.However, the dynamic boundary condition needs special treatment in order to be satisfied.For the free surface, it can be done by applying atmospheric pressure on the free surface particle.In order to apply atmospheric pressure on the free surface particle, it is vital to indicate the free surface particle accurately.A free surface particle can be identified under the following conditions: where i N is the number of particles within an effective range, and 0 N is the number of particles at the initial arrangement.Since the physical meaning of the interface between two fluids is similar to the free surface, the free surface search method in the MPS method can indicate an interface between two fluids.However, for the MPS method, the air region is regarded as a vacuum, which means that there is no particle for the air.The free surface searching method may not detect the particle, which disengaged from the In the MPS method, the kinematic free surface and interface boundary condition can be automatically satisfied by tracing the particle directly.However, the dynamic boundary condition needs special treatment in order to be satisfied.For the free surface, it can be done by applying atmospheric pressure on the free surface particle.In order to apply atmospheric pressure on the free surface particle, it is vital to indicate the free surface particle accurately.A free surface particle can be identified under the following conditions: where N i is the number of particles within an effective range, and N 0 is the number of particles at the initial arrangement.β 1 and β 2 are arbitrary numbers with a parameter that should be less than 1.0.In this study, β 1 and β 2 are set to 0.97 and 0.85, respectively.Since the physical meaning of the interface between two fluids is similar to the free surface, the free surface search method in the MPS method can indicate an interface between two fluids.However, for the MPS method, the air region is regarded as a vacuum, which means that there is no particle for the air.The free surface searching method may not detect the particle, which disengaged from the same fluid particle, because it may be surrounded by another fluid particle as the interface particle.In this study, the interface particle searching method, in which the free surface searching method is modified, can be calculated using the equation below: where subscription ξ denotes the kind of particle, and β 3 is a parameter for the interface searching method that is set to 0.35.All of the parameters in the particle searching method are selected by a numerical test.As shown in Equation ( 21), the particle number density can be calculated by using a similar particle.Moreover, wall and dummy particles excluded this calculation.All of the validations of particle interaction models were found by Kim [24].

Simulation of Dam Breaking
To demonstrate the performance of the MPS method on the aimed problem in this study, the classic example for validation of the Lagrangian formulation is simulated and compared with the experiment.The schematic model is shown in Figure 2. The tank partially contains the water column, which is blocked by the gate.The width of the tank is 0.6 m, and the height of the tank is 0.6 m.The water column is 0.15 m in width and 0.3 m in height.By opening the gate, the water column collapsed spontaneously, similar to the dam breaking.The pressure measuring point was 0.01 m above the tank bottom on the right wall.The initial distance of the particle was set to 0.0025 m.Therefore, a total of 11,628 particles were used for simulation.The density of water is 1000 kg/m 3 , and the kinematic viscosity of water is × 10 −6 m 2 /s.The gravitational acceleration is 9.81 m/s 2 .The computational time for 3-s simulation time was averaged by 183 s among 10 simulations.
Figure 3 shows the comparison between numerical simulation and the corresponding experiment by Koshizuka and Oka [18].According to snapshots in Figure 3, the collapsed water ran along the bottom to the right wall, and then collided.The collided water ran up along the wall, and then it fell by gravitational force at 0.8 s.The comparison of snapshots from numerical simulations and the experiment have good agreement.Figure 4 represents the pressure history measured at a reference point.According to the time history of pressure, the first peak was observed at 0.3 s.By using a contrasting study between snapshots and the time history of pressure, the first peak was due to the collision of the water column on the wall.The second peak was observed around 0.8 s, and it was related to the falling water after collision.Figure 5 shows the comparative study on the initial particle distance, which is based on the grid size of the conventional CFD method.As shown in Figure 5, the initial particle distance need to be less than 0.005 m.Therefore, the minimum particle distance at the initial arrangement was set to less than 0.005 m.
same fluid particle, because it may be surrounded by another fluid particle as the interface particle.In this study, the interface particle searching method, in which the free surface searching method is modified, can be calculated using the equation below: where subscription  denotes the kind of particle, and 3  is a parameter for the interface searching method that is set to 0.35.All of the parameters in the particle searching method are selected by a numerical test.As shown in Equation ( 21), the particle number density can be calculated by using a similar particle.Moreover, wall and dummy particles excluded this calculation.All of the validations of particle interaction models were found by Kim [24].

Simulation of Dam Breaking
To demonstrate the performance of the MPS method on the aimed problem in this study, the classic example for validation of the Lagrangian formulation is simulated and compared with the experiment.The schematic model is shown in Figure 2. The tank partially contains the water column, which is blocked by the gate.The width of the tank is 0.6 m, and the height of the tank is 0.6 m.The water column is 0.15 m in width and 0.3 m in height.By opening the gate, the water column collapsed spontaneously, similar to the dam breaking.The pressure measuring point was 0.01 m above the tank bottom on the right wall.The initial distance of the particle was set to 0.0025 m.Therefore, a total of 11,628 particles were used for simulation.The density of water is Figure 3 shows the comparison between numerical simulation and the corresponding experiment by Koshizuka and Oka [18].According to snapshots in Figure 3, the collapsed water ran along the bottom to the right wall, and then collided.The collided water ran up along the wall, and then it fell by gravitational force at 0.8 s.The comparison of snapshots from numerical simulations and the experiment have good agreement.Figure 4 represents the pressure history measured at a reference point.According to the time history of pressure, the first peak was observed at 0.3 s.By using a contrasting study between snapshots and the time history of pressure, the first peak was due to the collision of the water column on the wall.The second peak was observed around 0.8 s, and it was related to the falling water after collision.Figure 5 shows the comparative study on the initial particle distance, which is based on the grid size of the conventional CFD method.As shown in Figure 5, the initial particle distance need to be less than 0.005 m.Therefore, the minimum particle distance at the initial arrangement was set to less than 0.005 m.

Simulation of Dam Breaking with Soil Bed
Newly developed MPS for multiple liquids was applied to simulate the broken dam problem with the soil bed shown in Figure 6.The tank is divided by a gate, which is located 1 m away from the left side of the tank.The left side of the gate is filled with water particles with 0.15 m depth, while the soil bed with various thicknesses is on the right side of the tank.For the soil particle, polyethyleneoxide (PEO) solution was used, which corresponded to the experiment [3].In the simulation, the thickness of the PEO solution was set to 0.015 m, 0.030 m, 0.058 m, and 0.070 m.The kinematic viscosities of the water particles were

Simulation of Dam Breaking with Soil Bed
Newly developed MPS for multiple liquids was applied to simulate the broken dam problem with the soil bed shown in Figure 6.The tank is divided by a gate, which is located 1 m away from the left side of the tank.The left side of the gate is filled with water particles with 0.15 m depth, while the soil bed with various thicknesses is on the right side of the tank.For the soil particle, polyethyleneoxide (PEO) solution was used, which corresponded to the experiment [3].In the simulation, the thickness of the PEO solution was set to 0.015 m, 0.030 m, 0.058 m, and 0.070 m.The kinematic viscosities of the water particles were , respectively.The densities of water and PEO are 3 1000 kg/m and 3 960 kg/m , respectively.The gate opened upwards of 1.5 m/s.For this simulation, the total particle usage was 91,640 for 0.015 m of depth of PEO, 100,380 for 0.030 m of depth of PEO, 120,760 for 0.058 m of depth of PEO, and 141,420 for 0.070 m of depth of PEO.The initial particle distance is set to 0.001 m, and the time interval is 0.0001 s.All of the variances that use the initial particle distance and time interval were obtained by the numerical convergence test.

Simulation of Dam Breaking with Soil Bed
Newly developed MPS for multiple liquids was applied to simulate the broken dam problem with the soil bed shown in Figure 6.The tank is divided by a gate, which is located 1 m away from the left side of the tank.The left side of the gate is filled with water particles with 0.15 m depth, while the soil bed with various thicknesses is on the right side of the tank.For the soil particle, polyethyleneoxide (PEO) solution was used, which corresponded to the experiment [3].In the simulation, the thickness of the PEO solution was set to 0.015 m, 0.030 m, 0.058 m, and 0.070 m.The kinematic viscosities of the water particles were , respectively.The densities of water and PEO are 3 1000 kg/m and 3 960 kg/m , respectively.The gate opened upwards of 1.5 m/s.For this simulation, the total particle usage was 91,640 for 0.015 m of depth of PEO, 100,380 for 0.030 m of depth of PEO, 120,760 for 0.058 m of depth of PEO, and 141,420 for 0.070 m of depth of PEO.The initial particle distance is set to 0.001 m, and the time interval is 0.0001 s.All of the variances that use the initial particle distance and time interval were obtained by the numerical convergence test.

Simulation of Dam Breaking with Soil Bed
Newly developed MPS for multiple liquids was applied to simulate the broken dam problem with the soil bed shown in Figure 6.The tank is divided by a gate, which is located 1 m away from the left side of the tank.The left side of the gate is filled with water particles with 0.15 m depth, while the soil bed with various thicknesses is on the right side of the tank.For the soil particle, polyethylene-oxide (PEO) solution was used, which corresponded to the experiment [3].In the simulation, the thickness of the PEO solution was set to 0.015 m, 0.030 m, 0.058 m, and 0.070 m.The kinematic viscosities of the water particles were 1.0 × 10 −6 m 2 /s and 6.0 × 10 6 g/mol, respectively.The densities of water and PEO are 1000 kg/m 3 and 960 kg/m 3 , respectively.The gate opened upwards of 1.5 m/s.For this simulation, the total particle usage was 91,640 for 0.015 m of depth of PEO, 100,380 for 0.030 m of depth of PEO, 120,760 for 0.058 m of depth of PEO, and 141,420 for 0.070 m of depth of PEO.The initial particle distance is set to 0.001 m, and the time interval is 0.0001 s.All of the variances that use the initial particle distance and time interval were obtained by the numerical convergence test.The simulation ended at 1 s, and the computational time of each case was 11,720 s, 13,668 s, 19,271 s, and 22,518 s for cases of various thickness, which are 0.015 m, 0.03 m, 0.058 m, and 0.07 m, respectively.
The simulation ended at 1 s, and the computational time of each case was 11,720 s, 13,668 s, 19,271 s, and 22,518 s for cases of various thickness, which are 0.015 m, 0.03 m, 0.058 m, and 0.07 m, respectively.Figure 7 represents comparisons of results with 0.015-m depth of PEO solution between numerical simulation and experiment.By opening the gate, the water column collapsed spontaneously by gravitational force; then, the collided water particle collided with the PEO solution.According to Janosi et al., the PEO solution was used to simulate the solid particle [3].However, its behavior seems to be fluid.This is the reason that the multi-phase MPS method was employed to investigate it.The disturbed PEO solution and the collided water particle form mixed flows and propagated along the tank bed.The PEO solution was thrusted, and then it turned into waves.By recalling the dam break problem without erodible beds, the collapsed water column runs over the bottom of the tank rather than plunging.However, when the collapsed water column collided with the beds, it turned to the plunging waves and mixed, which is shown in both the experiment and the numerical simulation.Moreover, erodible beds were disturbed by an up thrust, which generated waves and propagated fluids.According to comparisons, the overturning wave features and the mixing process were almost the same as the corresponding experiment.There are slight differences between them, but this could be because of small 3D effects from the experiment and the small influence of the roughness of the gate.Since the discrepancy is very small, it can be negligible in any case.
Figures 8 and 9 shows the comparisons between the numerical simulations and corresponding experiments with various thicknesses of the PEO solution.The increased amount of PEO solution resisted the propagation of water waves.In this regard, the size of overturning was reduced, and the timing of the generation was delayed.Moreover, in the case of 0.15-m thickness of the PEO solution, the water particle overturned the PEO solution.However, an increased thickness of the PEO solution prevented it.Due to the weight of the soil bed and the collapsed water column, there was more upthrust.These phenomena were observed in both the experiments and the numerical simulation.In the cases of the thick soil layer, which were 0.058 m and 0.07 m, the collapsed water particle ran beneath the PEO solution rather than running over the soil bed.
In order to check the accuracy of the numerical simulation, the tip velocity of the plunged waves was measured and compared by the root mean square error method.The celerity of waves from the experiment was measured by image and captured time.The error rate was 94% for 0.015 m in depth and 97% for 0.03 m, 0.058 m, and 0.07 m in depth.Figure 7 represents comparisons of results with 0.015-m depth of PEO solution between numerical simulation and experiment.By opening the gate, the water column collapsed spontaneously by gravitational force; then, the collided water particle collided with the PEO solution.According to Janosi et al., the PEO solution was used to simulate the solid particle [3].However, its behavior seems to be fluid.This is the reason that the multi-phase MPS method was employed to investigate it.The disturbed PEO solution and the collided water particle form mixed flows and propagated along the tank bed.The PEO solution was thrusted, and then it turned into waves.By recalling the dam break problem without erodible beds, the collapsed water column runs over the bottom of the tank rather than plunging.However, when the collapsed water column collided with the beds, it turned to the plunging waves and mixed, which is shown in both the experiment and the numerical simulation.Moreover, erodible beds were disturbed by an up thrust, which generated waves and propagated fluids.According to comparisons, the overturning wave features and the mixing process were almost the same as the corresponding experiment.There are slight differences between them, but this could be because of small 3D effects from the experiment and the small influence of the roughness of the gate.Since the discrepancy is very small, it can be negligible in any case.
Figures 8 and 9 shows the comparisons between the numerical simulations and corresponding experiments with various thicknesses of the PEO solution.The increased amount of PEO solution resisted the propagation of water waves.In this regard, the size of overturning was reduced, and the timing of the generation was delayed.Moreover, in the case of 0.15-m thickness of the PEO solution, the water particle overturned the PEO solution.However, an increased thickness of the PEO solution prevented it.Due to the weight of the soil bed and the collapsed water column, there was more upthrust.These phenomena were observed in both the experiments and the numerical simulation.In the cases of the thick soil layer, which were 0.058 m and 0.07 m, the collapsed water particle ran beneath the PEO solution rather than running over the soil bed.
In order to check the accuracy of the numerical simulation, the tip velocity of the plunged waves was measured and compared by the root mean square error method.The celerity of waves from the experiment was measured by image and captured time.The error rate was 94% for 0.015 m in depth and 97% for 0.03 m, 0.058 m, and 0.07 m in depth.Figure 7 represents comparisons of results with 0.015-m depth of PEO solution between numerical simulation and experiment.By opening the gate, the water column collapsed spontaneously by gravitational force; then, the collided water particle collided with the PEO solution.According to Janosi et al., the PEO solution was used to simulate the solid particle [3].However, its behavior seems to be fluid.This is the reason that the multi-phase MPS method was employed to investigate it.The disturbed PEO solution and the collided water particle form mixed flows and propagated along the tank bed.The PEO solution was thrusted, and then it turned into waves.By recalling the dam break problem without erodible beds, the collapsed water column runs over the bottom of the tank rather than plunging.However, when the collapsed water column collided with the beds, it turned to the plunging waves and mixed, which is shown in both the experiment and the numerical simulation.Moreover, erodible beds were disturbed by an up thrust, which generated waves and propagated fluids.According to comparisons, the overturning wave features and the mixing process were almost the same as the corresponding experiment.There are slight differences between them, but this could be because of small 3D effects from the experiment and the small influence of the roughness of the gate.Since the discrepancy is very small, it can be negligible in any case.
Figures 8 and 9 shows the comparisons between the numerical simulations and corresponding experiments with various thicknesses of the PEO solution.The increased amount of PEO solution resisted the propagation of water waves.In this regard, the size of overturning was reduced, and the timing of the generation was delayed.Moreover, in the case of 0.15-m thickness of the PEO solution, the water particle overturned the PEO solution.However, an increased thickness of the PEO solution prevented it.Due to the weight of the soil bed and the collapsed water column, there was more upthrust.These phenomena were observed in both the experiments and the numerical simulation.In the cases of the thick soil layer, which were 0.058 m and 0.07 m, the collapsed water particle ran beneath the PEO solution rather than running over the soil bed.
In order to check the accuracy of the numerical simulation, the tip velocity of the plunged waves was measured and compared by the root mean square error method.The celerity of waves from the experiment was measured by image and captured time.The error rate was 94% for 0.015 m in depth and 97% for 0.03 m, 0.058 m, and 0.07 m in depth.

Swash Uprush Problem
In the slope bed such as the beach region, the uprush flows disturbed the sea bed, including small creatures and sediment.It affected the environment and piled structures.In this regard, the accurate expectation of soil behavior was due to large waves.For this simulation, the tilted tank was used.On the left side of the tank, water was contained and captured by the gate, and the waves could run over soil particles.All of the simulation conditions were followed by the corresponding experiment by Postacchini et al. [5].The 3.0-m water tank was inclined by 2.87 • , and the gate was located 1.0 m away from the left wall.The averaged water depth was 0.2 m, and the thickness of the soil bed was 0.2 m on the acrylic bed.The property of the soil bed was the PEO solution.The speed of opening the gate was 1.5 m/s upward.The schematic model of the simulation is shown in Figure 10.In this simulation, 94,320 particles were used in total.The simulation took near 100,000 s for 5 s of computational time.It was obtained by averaging 10 simulations.
When the gate opened, the water column spontaneously collapsed by gravitational acceleration and propagated toward the soil bed.The collapsed water column turned into very large waves, and it collided with the PEO solution.The PEO solution is not a rigid body.Therefore, it can be disturbed, and it was also propagated with water waves along the acrylic bed. Figure 11 shows the example of the numerical simulation.The elevation of the mixed flows and its maximum velocity were measured at the measuring point.The simulation results were compared to the corresponding experimental results shown in Figure 12.The left column of Figure 12 was captured by Postacchini et al. [5].In the comparison between the experimental and numerical processes, the results do not agree.The error rate from root mean square error (RMSE) was 82%.However, the tendency of it was matched.This is because the MPS method follows the fully Lagrangian approach, which means that all of the particles were tracked, and its physical quantities were observed on the individual particle.When there is fragmentation of the free surface, usually the conventional CFD treatment is ignored, while the particle method is not.As shown in Figure 11, the plunged and broken waves splashed water and the soil particle through the air region.This worked for measuring the maximum velocity and height.By ignoring it and taking only its tendency, both the experimental and numerical results can be regarded as appropriately observed.

Swash Uprush Problem
In the slope bed such as the beach region, the uprush flows disturbed the sea bed, including small creatures and sediment.It affected the environment and piled structures.In this regard, the accurate expectation of soil behavior was due to large waves.For this simulation, the tilted tank was used.On the left side of the tank, water was contained and captured by the gate, and the waves could run over soil particles.All of the simulation conditions were followed by the corresponding experiment by Postacchini et al. [5].The 3.0-m water tank was inclined by 2.87°, and the gate was located 1.0 m away from the left wall.The averaged water depth was 0.2 m, and the thickness of the soil bed was 0.2 m on the acrylic bed.The property of the soil bed was the PEO solution.The speed of opening the gate was 1.5 m/s upward.The schematic model of the simulation is shown in Figure 10.In this simulation, 94,320 particles were used in total.The simulation took near 100,000 s for 5 s of computational time.It was obtained by averaging 10 simulations.
When the gate opened, the water column spontaneously collapsed by gravitational acceleration and propagated toward the soil bed.The collapsed water column turned into very large waves, and it collided with the PEO solution.The PEO solution is not a rigid body.Therefore, it can be disturbed, and it was also propagated with water waves along the acrylic bed. Figure 11 shows the example of the numerical simulation.The elevation of the mixed flows and its maximum velocity were measured at the measuring point.The simulation results were compared to the corresponding experimental results shown in Figure 12.The left column of Figure 12 was captured by Postacchini et al. [5].In the comparison between the experimental and numerical processes, the results do not agree.The error rate from root mean square error (RMSE) was 82%.However, the tendency of it was matched.This is because the MPS method follows the fully Lagrangian approach, which means that all of the particles were tracked, and its physical quantities were observed on the individual particle.When there is fragmentation of the free surface, usually the conventional CFD treatment is ignored, while the particle method is not.As shown in Figure 11, the plunged and broken waves splashed water and the soil particle through the air region.This worked for measuring the maximum velocity and height.By ignoring it and taking only its tendency, both the experimental and numerical results can be regarded as appropriately observed.

Swash Uprush Problem
In the slope bed such as the beach region, the uprush flows disturbed the sea bed, including small creatures and sediment.It affected the environment and piled structures.In this regard, the accurate expectation of soil behavior was due to large waves.For this simulation, the tilted tank was used.On the left side of the tank, water was contained and captured by the gate, and the waves could run over soil particles.All of the simulation conditions were followed by the corresponding experiment by Postacchini et al. [5].The 3.0-m water tank was inclined by 2.87°, and the gate was located 1.0 m away from the left wall.The averaged water depth was 0.2 m, and the thickness of the soil bed was 0.2 m on the acrylic bed.The property of the soil bed was the PEO solution.The speed of opening the gate was 1.5 m/s upward.The schematic model of the simulation is shown in Figure 10.In this simulation, 94,320 particles were used in total.The simulation took near 100,000 s for 5 s of computational time.It was obtained by averaging 10 simulations.
When the gate opened, the water column spontaneously collapsed by gravitational acceleration and propagated toward the soil bed.The collapsed water column turned into very large waves, and it collided with the PEO solution.The PEO solution is not a rigid body.Therefore, it can be disturbed, and it was also propagated with water waves along the acrylic bed. Figure 11 shows the example of the numerical simulation.The elevation of the mixed flows and its maximum velocity were measured at the measuring point.The simulation results were compared to the corresponding experimental results shown in Figure 12.The left column of Figure 12 was captured by Postacchini et al. [5].In the comparison between the experimental and numerical processes, the results do not agree.The error rate from root mean square error (RMSE) was 82%.However, the tendency of it was matched.This is because the MPS method follows the fully Lagrangian approach, which means that all of the particles were tracked, and its physical quantities were observed on the individual particle.When there is fragmentation of the free surface, usually the conventional CFD treatment is ignored, while the particle method is not.As shown in Figure 11, the plunged and broken waves splashed water and the soil particle through the air region.This worked for measuring the maximum velocity and height.By ignoring it and taking only its tendency, both the experimental and numerical results can be regarded as appropriately observed.

Conclusions
Since the multi-phase flow has more interface than single-phase flows, multi-phase particle interaction models are required to be employed in order to achieve an accuracy of numerical simulation and avoid the numerical error.The newly employed multi-phase model, including selfbuoyancy and surface tension models, was introduced in this study.The dam break problem with erodible beds was simulated by developing the MPS method with the multi-phase particle interaction models, and then comparing the well to the corresponding experimental results including a singlephase dam break, dam break with erodible beds, and the swash uprush problem in the inclined tank.
In order to validate the developed MPS program, a single-phase dam break, which is a classical problem of the particle method, was simulated.The single-phase problem does not activate multiphase particle interaction models.However, it is necessary for being validated.By comparing snapshots between experiments and numerical simulation, both results show good agreement.Regarding the pressure calculation, the time history of pressure was well matched with visual results by taking pressure peaks.Moreover, the initial particle distance, which is the same as the grid size in a conventional CFD method, also investigated the initial particle distance, which was selected by a numerical test.
For the erodible beds problem, the PEO solution was used as a bed property.The PEO solution is not exactly a solid object.However, by considering the liquefaction of erodible beds, the PEO solution was an appropriate model.The water captured by the gate on the left side of the tank and erodible PEO solution beds was located beyond the gate with various thicknesses corresponding to experimental conditions.When the gate opened, the captured water collapsed and ran toward erodible beds.In the case of a relatively thin layer of the soil bed that was 0.015 m in thickness, the collapsed water runs over the beds, while in other cases such as 0.030 m, 0.058 m, and 0.070 m, the

Conclusions
Since the multi-phase flow has more interface than single-phase flows, multi-phase particle interaction models are required to be employed in order to achieve an accuracy of numerical simulation and avoid the numerical error.The newly employed multi-phase model, including self-buoyancy and surface tension models, was introduced in this study.The dam break problem with erodible beds was simulated by developing the MPS method with the multi-phase particle interaction models, and then comparing the well to the corresponding experimental results including a single-phase dam break, dam break with erodible beds, and the swash uprush problem in the inclined tank.
In order to validate the developed MPS program, a single-phase dam break, which is a classical problem of the particle method, was simulated.The single-phase problem does not activate multi-phase particle interaction models.However, it is necessary for being validated.By comparing snapshots between experiments and numerical simulation, both results show good agreement.Regarding the pressure calculation, the time history of pressure was well matched with visual results by taking pressure peaks.Moreover, the initial particle distance, which is the same as the grid size in a conventional CFD method, also investigated the initial particle distance, which was selected by a numerical test.
For the erodible beds problem, the PEO solution was used as a bed property.The PEO solution is not exactly a solid object.However, by considering the liquefaction of erodible beds, the PEO solution was an appropriate model.The water captured by the gate on the left side of the tank and erodible PEO solution beds was located beyond the gate with various thicknesses corresponding to experimental conditions.When the gate opened, the captured water collapsed and ran toward erodible beds.In the case of a relatively thin layer of the soil bed that was 0.015 m in thickness, the collapsed water runs over the beds, while in other cases such as 0.030 m, 0.058 m, and 0.070 m, the water dug into the erodible beds and lifted the beds.This is because the weight of the bed layer affected the flow development that was observed in both experiments and numerical simulations.
The swash uprush problem was also investigated by the newly developed multi-phase MPS method.For the swash uprush problem, the tilted tank, which contained the water and the PEO solution divided by the gate, was considered.The tilting angle was 2.87 • , which is the same as the corresponding experiment.Upon comparing the wave height and maximum velocity at the certain measuring point between the experiment and numerical simulation, it was found that have satisfactory agreement.There were several irregular data points in both the wave height and maximum velocity history, but these were related to the fragmented particle by a splash, which did not affect all of the results.
These observations cannot fully explain the effects of the soil bed downstream of the collapsed fluid, because PEO is also regarded as a fluid particle.Moreover, the physical characteristics of PEO solution are similar to a water particle.However, in the most realistic situation, the soil bed contains a fluid.Therefore, it can be assumed to be a fluid due to liquefaction.Despite stated flaws in this simulation, the agreement between numerical simulation and the experiment may help to understand the phenomena of soil behavior due to very violent forces induced by large-scale flows.It can be used to design wave prevention technology to avoid massive disasters.

zFigure 1 .
Figure 1.(a) Schematic of a curvature of interface and (b) finding method for a normal vector of interface.

1  and 2  1  and 2 
are arbitrary numbers with a parameter that should be less than 1.0.In this study, are set to 0.97 and 0.85, respectively.

Figure 1 .
Figure 1.(a) Schematic of a curvature of interface and (b) finding method for a normal vector of interface.
m/s .The computational time for 3-s simulation time was averaged by 183 s among 10 simulations.

Figure 2 .
Figure 2. Schematic model of the dam break problem.Figure 2. Schematic model of the dam break problem.

Figure 2 .Figure 3 .
Figure 2. Schematic model of the dam break problem.Figure 2. Schematic model of the dam break problem.

Figure 4 .
Figure 4. Time history of pressure at measuring point.

Figure 5 .
Figure 5. Convergence test for the initial particle distance.

Figure 4 .
Figure 4. Time history of pressure at measuring point.

Figure 5 .
Figure 5. Convergence test for the initial particle distance.

Figure 4 .
Figure 4. Time history of pressure at measuring point.

Figure 3 .
Figure 3. Evolution of the collapsed water column.

Figure 4 .
Figure 4. Time history of pressure at measuring point.

Figure 5 .
Figure 5. Convergence test for the initial particle distance.

Figure 5 .
Figure 5. Convergence test for the initial particle distance.

Figure 6 .
Figure 6.Schematic of the broken dam with soil bed.

Figure 6 .
Figure7represents comparisons of results with 0.015-m depth of PEO solution between numerical simulation and experiment.By opening the gate, the water column collapsed spontaneously by gravitational force; then, the collided water particle collided with the PEO solution.According to Janosi et al., the PEO solution was used to simulate the solid particle[3].However, its behavior seems to be fluid.This is the reason that the multi-phase MPS method was employed to investigate it.The disturbed PEO solution and the collided water particle form mixed flows and propagated along the tank bed.The PEO solution was thrusted, and then it turned into waves.By recalling the dam break problem without erodible beds, the collapsed water column runs over the bottom of the tank rather than plunging.However, when the collapsed water column collided with the beds, it turned to the plunging waves and mixed, which is shown in both the experiment and the numerical simulation.Moreover, erodible beds were disturbed by an up thrust, which generated waves and propagated fluids.According to comparisons, the overturning wave features and the mixing process were almost the same as the corresponding experiment.There are slight differences between them, but this could be because of small 3D effects from the experiment and the small influence of the roughness of the gate.Since the discrepancy is very small, it can be negligible in any case.Figures8 and 9shows the comparisons between the numerical simulations and corresponding experiments with various thicknesses of the PEO solution.The increased amount of PEO solution resisted the propagation of water waves.In this regard, the size of overturning was reduced, and the timing of the generation was delayed.Moreover, in the case of 0.15-m thickness of the PEO solution, the water particle overturned the PEO solution.However, an increased thickness of the PEO solution prevented it.Due to the weight of the soil bed and the collapsed water column, there was more upthrust.These phenomena were observed in both the experiments and the numerical simulation.In the cases of the thick soil layer, which were 0.058 m and 0.07 m, the collapsed water particle ran beneath the PEO solution rather than running over the soil bed.In order to check the accuracy of the numerical simulation, the tip velocity of the plunged waves was measured and compared by the root mean square error method.The celerity of waves from the experiment was measured by image and captured time.The error rate was 94% for 0.015 m in depth and 97% for 0.03 m, 0.058 m, and 0.07 m in depth.Time [s] Experiment Simulation Appl.Sci.2018, 8, x FOR PEER REVIEW 10 of 15 and 22,518 s for cases of various thickness, which are 0.015 m, 0.03 m, 0.058 m, and 0.07 m, respectively.

Figure 6 .
Figure 6.Schematic of the broken dam with soil bed.

Figure 9 .
Figure 9.Comparison between the experiment and numerical simulation with various thicknesses of the soil bed at 0.6 s.

Figure 7 .Figure 7 .Figure 8 .Figure 9 .
Figure 7. Result comparisons between the experiment and numerical simulation for a broken dam with 0.15-m depth of soil bed.

Figure 8 .Figure 7 .Figure 8 .Figure 9 .
Figure 8.Comparison between the experiment and numerical simulation with various thicknesses of the soil bed at 0.3 s.

Figure 9 .
Figure 9.Comparison between the experiment and numerical simulation with various thicknesses of the soil bed at 0.6 s.

Figure 10 .
Figure 10.Schematic numerical model of the swash uprush problem.

Figure 11 .
Figure 11.Snapshot of simulation for the swash uprush problem.

Figure 10 .
Figure 10.Schematic numerical model of the swash uprush problem.

Figure 10 .
Figure 10.Schematic numerical model of the swash uprush problem.

Figure 11 .
Figure 11.Snapshot of simulation for the swash uprush problem.

Figure 11 .
Figure 11.Snapshot of simulation for the swash uprush problem.

Figure 12 .
Figure 12.Time history of change of height and velocity from (a) experiment and (b) numerical simulation.

Figure 12 .
Figure 12.Time history of change of height and velocity from (a) experiment and (b) numerical simulation.