Moving particle semi-implicit method with improved pressures stability properties

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.

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):

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.

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.

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)

  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.

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 ⁠.

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.

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.

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

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 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.

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.

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.

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.

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.

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.