## Abstract

Moving particle semi-implicit (MPS) method is one of the Lagrangian methods widely used in engineering issues. This method, however, suffers from unphysical oscillations in its original form. In the present study, a modified incompressible MPS method is proposed to suppress these oscillations and is used for simulating free surface problems. To demonstrate the stability of the presented method, different kernel functions are used in the case of numerical dam break modeling as a benchmark simulation. A simple form of definition of curved wall boundaries is suggested which eliminates dummy particles and subsequently saves CPU time. Flow over an ogee spillway is simulated for the first time with the I-MPS method and as a new test case which has several curved lines in its geometry. The comparisons between theoretical solutions/experimental data and simulation results in terms of free surface and pressure show the accuracy of the method.

## INTRODUCTION

The moving particle semi-implicit (MPS) method proposed by Koshizuka & Oka (1996) is one of the fully Lagrangian particle methods. The MPS method is used in a wide range of engineering problems, such as sloshing flows (Khayyer & Gotoh 2010), multiphase flows (Shakibaeinia & Jin 2012), and fluid–structure interaction problems (Hwang *et al.* 2014). The MPS method in its original form suffers from pressure oscillations and fluctuations. Various modifications of the MPS method have been suggested in previous studies in order to suppress these spurious oscillations. Kondo *et al.* (2008) proposed to add an artificial pressure to pressure obtained from the Poisson equation. This method caused suppression of the pressure oscillations. Khayyer & Gotoh (2009) suggested pressure Poisson equations (PPEs) and a conservative form of gradient model. Their studies show that their method introduces accurate results. Tanaka & Masunaga (2010) proposed a composed PPE with two components: one is formulated as a function of density variation with a relaxation coefficient and the other utilizes velocity divergence condition. Khayyer & Gotoh (2010) proposed a higher order Laplacian model for enhancement and stabilization of pressure calculation by the MPS method. Application of higher order Laplacian model stabilized the method and introduced smooth pressure fields in their study. Kondo & Koshizuka (2011) introduced a PPE term composed of three parts with a coefficient for each part. Choice of proper coefficients in their method suppresses unwanted pressure fluctuations. Lee *et al.* (2011) proposed a modified PPE based on Tanaka & Masunaga's (2010) study. Shakibaeinia & Jin (2010) developed weakly compressible MPS. In this manner, an equation of state (Tait's equation) is implemented to calculate pressure of particles instead of solving the Poisson equation. This method is more stable in comparison with the original MPS. The above-mentioned studies generate reasonable and smooth pressure fields. However, almost all of these methods have more than one term in PPE. This feature causes the use of pre-analysis of a problem to reach a proper coefficient for each term.

The present study utilizes a simple and new form of particle number density assignment which causes slight changes in the PPE, gradient, and Laplacian models. These changes or techniques in particle number density calculation generate smooth and accurate computational pressure fields. Moreover, use of this technique enables us to eliminate dummy particles, which allows saving CPU time through the reduction of numbers of particles. In this technique, particle number density of each particle updates at the beginning of each time step and this value is considered as initial particle number density in comparison with the original MPS. Therefore, the source term of the Poisson equation has only one term with a relaxation factor as used in the original MPS method.

## GOVERNING EQUATIONS

*t*is the time,

*P*is the pressure, is the dynamic viscosity, and is the gravitational acceleration. In Lagrangian coordinates, the convection terms are directly calculated by the motion of particles. Since Lagrangian derivatives involve advection terms, these terms are eliminated from the left side of the momentum equation. The right side terms expressed by differential operators should be replaced by particle interactions.

## PARTICLE INTERACTION MODELS

In the MPS method, particle *i* interacts with other particles, *j*, in its vicinity covered with a kernel function , where is the distance between particles *i* and *j* and is the support domain. In the original MPS developed by Koshizuka & Oka (1996) standard kernel (KF-1) (Figure 1) has been applied in the majority of calculations.

Implying of standard kernel in computations has some defects and privileges. An important advantage of implying standard kernel is the feature of infinity at , which prevents particle clustering. Nevertheless, the singularity at the mentioned point which is adopted for the weight function is one of the reasons for instability in the original MPS method (Kondo & Koshizuka 2011). However, it is possible to use other weight functions to stabilize the method and these have been applied in the literature. Comparative analysis has been done on the performance of six kernel functions, including the standard kernel by simulating a dam break flow, and B-spline kernel was found to be more stable than five other employed ones (Ataie-Ashtiani & Farhadi 2006). Similar studies have been done on the effect of kernel functions (Pan *et al.* 2008; Shakibaeinia & Jin 2010; Lee *et al.* 2011). In this study, different kernel functions, which have the kernel conditions such as positivity, compact, unity, and decay and delta function behavior with major focus on only showing the method stability (no numerical explosion during simulations) are employed in the dam break problem in order to evaluate the stability of the presented method and the effect of kernels on the results. The following kernels are implemented in this study (Figure 2):

KF-2 Wendland (1995) | (3b) | |

KF-3 Morris et al. (1997) | ||

KF-4 Kondo & Koshizuka (2011) | ||

KF-5 Shakibaeinia & Jin (2010) | ||

KF-6 Cummins & Rudman (1999) | ||

KF-7 Morris et al. (1997) |

KF-2 Wendland (1995) | (3b) | |

KF-3 Morris et al. (1997) | ||

KF-4 Kondo & Koshizuka (2011) | ||

KF-5 Shakibaeinia & Jin (2010) | ||

KF-6 Cummins & Rudman (1999) | ||

KF-7 Morris et al. (1997) |

where is the radius of the support domain which restricts the interaction area around each particle. Since the area that is covered by these kernel functions is bounded, a particle interacts with a finite number of neighboring particles. Compared with a function covering an infinite area, such as a Gaussian function, all of the used functions need less memory and computational time.

*n*. In Equation (4), the contribution from particle

*i*itself is not considered. The operator represents the kernel approximation. Assuming that the particles have the same mass

*m*, the fluid density is proportional to the particle number density and may be computed using the following equation: Thus, the continuity equation is satisfied if the particle number density is constant. This constant value is denoted by (Koshizuka & Oka 1996). In original MPS discretization, differential operators, such as gradient and Laplacian, are represented by the following particle interaction models using the weight function: where

*d*is the number of spatial dimensions and is an arbitrary scalar. The Laplacian model parameter is a correction parameter which causes the variance increase equal to the analytical solution and is defined as: The current model of Laplacian is conservative since the quantity lost by particle

*i*is just obtained by particles

*j*(Koshizuka & Oka 1996). However, the gradient model is non-conservative (Khayyer & Gotoh 2009).

## I-MPS SOLUTION ALGORITHM

A projection method (Cummins & Rudman 1999) is used for time integration. In this method each time step is split into two basic steps of prediction and correction. The method uses a fractional step with the velocity field integrated forward in time without enforcing incompressibility. According to the resulting intermediate velocity field, particles move and their provisional positions are obtained. Based on the predicted particle positions, particle number density in the prediction step is computed. The value of the predicted particle number density is not equal to that of the initial value. Therefore, a correction step is required. In this step the pressure term, obtained from the mass conservation, is used to enforce incompressibility in the calculation. Modification of the particle number density is related to the modification of the velocity through the mass conservation equation. The modification of the velocity is derived from the implicit pressure term in the momentum conservation equation. The original MPS algorithm is depicted in Figure 3.

## MODIFICATIONS ON I-MPS

Using the proposed technique for particle number density assignment, for each particle, the particle number density is updated at the beginning of each time step, which causes smooth pressure fields (Arami-Fadafan 2014). As mentioned in the MPS method, density is approximated by the particle number density. Density error calculations are described according to the following procedure based on the study of Ataie-Ashtiani *et al.* (2008), and there is a similar procedure for particle number density error:

*a*at time had been denoted by , the particle number density of this particle at time

*t*will be: where are particle number density of particle

*a*and particle number density error until time

*t*, respectively. At time

*t*and the prediction step, the provisional particle number density of the particle is shown with . At the correction step and enforcing the incompressibility, particle number density must be modified: By combining Equations (10) and (11): where is the deviation of particle number density due to ignoring the pressure term in Equation (2) at the prediction step. In the incompressible MPS method, the source term of the Poisson equation for the particle is computed by: that may be written in the following form: In the original MPS method, both parts of and are considered in calculations whereas only should contribute in the source term. Clearly after the correction step, particle number density has been changed due to change in particles' position. Therefore, considering this value as initial particle number density at the beginning of each time step for each particle, it causes the elimination of the extra part in Equation (13b). In this investigation, this technique is used in calculations and therefore only is considered in source term of the Poisson equation. Using the mentioned technique, it is possible to eliminate dummy particles because of ignoring the same value in both and in the numerator of the Poisson equation (in the presented method, is updated each time step for each particle). This enables us to use only one row of wall particles in simulations. In MPS methods, since the number densities of particles near the surface are small, usually negative values of pressure after solving Poisson's equation of the pressure are eliminated, and these negative values are set to zero (Kondo

*et al.*2008). However, due to update in particle number density at the beginning of each time step for each particle, the method is capable of considering negative values of pressure in calculations (this is shown in the pressure field of ogee spillway simulation). Utilizing of one row wall boundary causes a reduction in particle number density values. Therefore, the value of particle number density is corrected mathematically for each particle (wall particles, inner fluid particles, and surface particles). This corrected value is used in the denominator of the Poisson equation and in gradient and Laplacian models. Thus, the following equations are used in our method: As mentioned before, the operator represents the kernel approximation. Due to mathematical correction, particle number density error is almost zero. The obtained pressure term has some similarity in apparent form to the second term of Equation (18) proposed by Kondo & Koshizuka (2011): where and are dimensionless parameters. Khayyer & Gotoh (2011) modified Equation (18) as the following form: Despite similarity in apparent form of Equation (14) with the second term of Equations (18) and (19), there are significant differences between the proposed method and the above-mentioned methods:

- (a)
Proposed PPE is very similar to the original PPE form of MPS (only one term in PPE).

- (b)
The denominator of the proposed equation utilizes the corrected form of particle number density almost equal to (not exactly ) which considers the particle positions in evaluating particle number density, . Moreover, use of the proposed form of PPE (because of suggested form of the denominator) has the ability of defining one row wall boundary.

- (c)
Use of new form of Laplacian and gradient operator in modified MPS.

The proposed I-MPS method algorithm is shown in Figure 4.

## BOUNDARY CONDITIONS

### Wall

Solid boundaries are represented by one line of particles. The Poisson equation of pressure is solved on these particles. In the absence of any contrivances, by approaching solid boundaries, the density of particles decreases, which causes recognition of the wall particles as free surface particles (Koshizuka & Oka 1996). Thus, several lines of dummy particles are also placed outside of solid walls in order to keep the fluid density at the wall particles in the standard MPS method. The thickness of dummy particles depends on the kernel range. Moreover, the MPS method does not calculate pressure gradient between dummy particles and inner fluid particles (Shao & Lo 2003). Thus, if only mathematical effect of dummy particles be considered in calculations, it is possible to eliminate dummy particles. For fluid particles near the boundaries, the treatment is similar to wall particles. Therefore, Equation (17) has been used for all wall particles, inner particles, and free-surface particles in order to keep close to .

### Free surface

It is clear that for an inner fluid particle, cut-off radius of kernel covers a finite number of particles, namely neighbor particles (Figure 5). For an inner particle, the complete number of neighbor particles is achieved. Since, no fluid particle exists in the outer region of free surface line (or near wall boundaries in absence of dummy/ghost particles or any contrivances), by approaching the free surface, the number of particles in the vicinity of the central particle in the support domain deduces and, consequently, the density of the central particle decreases. If particle number density of the central particle be less than distinct criteria (i.e., ), the particle is recognized as free surface particle. Known pressure (zero pressure) is prescribed to the surface particles.

### Inflow

A new form of inlet boundary is proposed and used in the present study. In this manner, in order to reduce the pressure fluctuations at the inlet boundary, some layers of particles (whose role is as dummy particles for inlet particles) with prescribed velocity equal to inlet velocity are defined at the inflow boundary to compensate for density deficiency. The significant difference between the used boundary in this study and other studies is related to moving extra layers of particles, which have been defined just before inlet particles. The particles on the inner first line of these particles are involved in the pressure calculation. Inlet particles are added to the distance between these layers and the fluid particles. For inlet particles, velocity in the flow direction is set to known velocity inflow until these particles move far from the inlet position, at least equal to a determined distance such as . is the initial distance between neighboring particles in the initial configuration. Inflow particles are added to inflow boundary each K time step, related to inflowing velocity, time step size, and distance of particles.

### Outflow

Each particle leaves the computational domain, eliminated from the computations.

## MODEL APPLICATIONS

In order to show the accuracy and validity of the proposed method, different case studies of dynamic problems are simulated: dam break problems as common test problems, flow over an ogee spillway as an open channel flow, and water bubble problem. Solid boundaries of the latter problem consist of several curve lines. Simulations of problems with curved boundaries are an interesting issue in particle methods.

### Dam break

Dam break problem, experimentally implemented by Hu & Kashiwagi (2004), is simulated using the proposed I-MPS to show the accuracy of the method. The geometry of the experimental set-up is shown in Figure 6.

Experimental history time of pressure at point A in Figure 6 is shown in Figure 7, and compared with the results of the used I-MPS method to evaluate the accuracy of the proposed method.

Computational pressure fields at different times (t = 0.02, 0.3, 0.74, and 1.2 s) are shown in Figure 8 using the present method.

Simulations of the dam break problem of Koshizuka & Oka (1996) were done in order to show the stability and to calculate the convergence of the proposed method. Initial geometry of this problem is shown in Figure 9.

In Figure 10, the comparison between pressure distributions of the original MPS and present method at t = 0.7 s are depicted implementing KF-1 in simulations. The obtained pressure distributions utilizing the present method are generally much more reasonable compared with expected theoretical pressure values.

Time histories of pressure at points B and C in Figure 9 are shown for these simulations in the following figure.

*CR*represents the convergence rate,

*X*is the computed value of the leading edge. Subscript represents the configuration, and the superscript

*t*denotes the sample point at different times (Shobeyri & Afshar 2012). The convergence rate of this method is calculated to be 1.596. Pressure fields at different times are depicted in Figure 12 for different particle sizes.

The proposed method introduces smooth pressure fields using different sizes of particles. In simulations when the particle size of 0.004 m is used, results of pressure fields and the geometry of fluid have small differences in comparison with those of utilizing a particle size of 0.008 m. Therefore, the latter size of particle is used to show the stability of the method implementing the above-mentioned kernels. As depicted in Figure 13, almost similar trajectories of the number of free surface particles are introduced utilizing each mentioned kernel. The proposed method is more stable than the original MPS method since no numerical explosion occurred using various types of kernels in our simulations.

### Ogee spillway

Despite capabilities of the MPS method to simulate flows with large deformations and fragmentations, this method is rarely used for modeling of open channel flows (Shakibaeinia & Jin 2009; Arami-Fadafan & Hessami-Kermani 2015; Jafari-Nodoushan *et al.* 2016). Flow over a spillway is an example of open channel flows. The experimental conditions and measurements performed by Chatila & Tabbara (2004) are used to simulate this study. The geometry of the spillway crest is based on the hydraulic design charts of USACE-WES profiles. The spillway has been designed by a vertical upstream face, and the radius of the upstream curved surface defined by and . The downstream profile of the crest centerline is defined by and . The slope of the straight portion of the spillway face is 60° (or a slope ). Prescribed inflow velocity of corresponding to the water head of over the spillway crest is used in this simulation. KF-2 with cutoff radius of and with initial particle distance of 0.005 m is used for this discharge. Initial configuration of particles is depicted in Figure 14.

*t*, respectively. Therefore, for the prediction step, this discrepancy is retrieved mathematically. An auxiliary function, as described in the wall boundary section, determines which wall particle must contribute in the calculations. Figure 15 shows the particle positions together with pressure fields at different times.

Free surface particles colored by velocity magnitude and experimental measurements of free surface elevations are shown in Figure 15(d). Computed pressure field in the reservoir is smooth and almost equal to hydrostatic value . There is no flow separation or overestimate in free surface elevation using this method, which shows the true prescription of the solid boundary condition in the modeling.

### Water bubble problem

*a*) × semi-minor axis (

*b*) should remain constant, and the outer boundary surface remain smooth (Bonet & Lok 1999). The condition that the drop remains elliptic with time varying axes

*a*and

*b*is that: where and is the initial value of , and are the x and y velocity components, respectively. The condition that the fluid remains incompressible is that

*a*

*×*

*b*is constant (Monaghan 1994). This means that the volume of the bubble does not change. Equations (25) and (26) can be solved semi-analytically with a higher order of accuracy (utilizing finite difference method). The results obtained were compared with the MPS simulation. In this study, the initial geometry of the bubble is represented by 1,313 particles with an initial particle distance of 0.05 m. Results are shown in Figures 16 and 17.

### Flow over curved bed

In order to show the robustness of the proposed method to simulate hydraulic problems, another new test case, namely, flow over a curved bed experimentally performed by Sivakumaran *et al.* (1983) is modeled numerically utilizing one row wall boundary condition. Laboratory experiments were carried out in a 915 × 75 × 44.5 cm flume with a bed of 1.5 cm thick plywood. The bed was elevated 10 cm above the base of the flume and the flume width was vertically partitioned along the entire channel into two compartments with the larger of 30 cm wide. Curved bed shape profile follows , x, y in meters. Several measurements were conducted on two discharge values and corresponding to heads of 15 cm and 7 cm over the crest, respectively. In this study, unit discharge of is simulated. Initial configuration of particles is depicted in Figure 18.

In this simulation, at both the upstream and downstream of the crest, the shape of the bed has curvature; initial particle size of 0.015 cm was used for numerical modeling of the current problem. In Figure 19, pressure and velocity fields at t = 1 s are shown. Pressure at the upstream reservoir is relatively hydrostatic.

Experimental and computational pressure head along the curved bed is depicted in Figure 20. As flow ascends the curved bed, flow depth decreases and velocity increases. Pressure drops from a hydrostatic pressure at upstream to a minimum pressure downstream of the crest.

Finally, computed free surface was compared with experimental measurements in Figure 21, which shows the accuracy of the MPS method evaluating free surface line.

## CONCLUSION

A stabilized I-MPS method was proposed and used to simulate free surface flows. The model of the used I-MPS method considering used boundary conditions is capable of presenting smooth pressure fields. A simple form of PPE was suggested which is similar to the original PPE in the MPS method. It consists of one term in PPE. The method is stable since the test cases have been simulated successfully applying various kernels. Based on proposed wall boundary and PPE, it is possible to use one row of wall particles and modeling of problems with curved boundaries such as ogee spillways. The proposed procedure saves CPU time because of eliminating dummy particles.