Table of Contents Author Guidelines Submit a Manuscript
Journal of Applied Mathematics
Volume 2012, Article ID 139583, 24 pages
Research Article

A Stabilized Incompressible SPH Method by Relaxing the Density Invariance Condition

1Department of Civil Engineering, Kyushu University, 744 Motooka, Nishi–ku, Fukuoka 819-0395, Japan
2Faculty of Education and Human Science, Yokohama National University, 79-1 Tokiwadai, Hodogaya-ku, Yokohama 240-8501, Japan

Received 5 January 2012; Accepted 16 March 2012

Academic Editor: Hiroshi Kanayama

Copyright © 2012 Mitsuteru Asai et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


A stabilized Incompressible Smoothed Particle Hydrodynamics (ISPH) is proposed to simulate free surface flow problems. In the ISPH, pressure is evaluated by solving pressure Poisson equation using a semi-implicit algorithm based on the projection method. Even if the pressure is evaluated implicitly, the unrealistic pressure fluctuations cannot be eliminated. In order to overcome this problem, there are several improvements. One is small compressibility approach, and the other is introduction of two kinds of pressure Poisson equation related to velocity divergence-free and density invariance conditions, respectively. In this paper, a stabilized formulation, which was originally proposed in the framework of Moving Particle Semi-implicit (MPS) method, is applied to ISPH in order to relax the density invariance condition. This formulation leads to a new pressure Poisson equation with a relaxation coefficient, which can be estimated by a preanalysis calculation. The efficiency of the proposed formulation is tested by a couple of numerical examples of dam-breaking problem, and its effects are discussed by using several resolution models with different particle initial distances. Also, the effect of eddy viscosity is briefly discussed in this paper.

1. Introduction

The meshless particle methods have been applied in many engineering applications including the free-surface fluid flows. In the particle methods, the state of a system is represented by a set of discrete particles, without a fixed connectivity; hence, such methods are inherently well suited for the analysis of moving discontinuities and large deformations such as the free-surface fluid flows with breaking and fragmentation.

The SPH technique was originally proposed by Lucy [1] and further developed by Gingold and Monaghan [2] for treating astrophysical problems. Its main advantage is the absence of a computational grid or mesh since it is spatially discretized into Lagrangian moving particles. This allows the possibility of easily modeling flows with a complex geometry or flows where large deformations or the appearance of a free surface occurs. At the present time, it is being exploited for the solution of problems appearing in different physical processes. Monaghan [3] has provided a fairly extensive review of SPH methods. The SPH method had been applied into compressible and incompressible viscous flow problems [47]. The SPH is originally developed in compressible flow, and then some special treatment is required to satisfy the incompressible condition. One approach is to run the simulations in the quasi-incompressible limit, that is, by selecting the smallest possible speed of sound which still gives a very low Mach number ensuring density fluctuations within 1% [4, 5]. This method is known as the Weakly Compressible Smooth Particle Hydrodynamics (WCSPH). In the WCSPH, the artificial viscosity, which is originally developed by Monaghan, has been widely used not only for the energy dissipation but also for preventing unphysical penetration of particles. Recently, a proposal for constructing an incompressible SPH model has been introduced, whose pressure is implicitly calculated by solving a discretized pressure Poisson equation at every time step [817].

Lee et al. [13] presented comparisons of a semi-implicit and truly incompressible SPH (ISPH) algorithm with the classical WCSPH method, showing how some of the problems encountered in WCSPH have been resolved by using ISPH to simulate incompressible flows. They used the function of temporal velocity divergence for discretized source terms of Poisson equation of pressure to ensure truly incompressible flow. Khayyer et al. [14, 15] proposed a corrected incompressible SPH method (CISPH) derived based on a variational approach to ensure the angular momentum conservation of ISPH formulations to improve the pressure distribution by improvement of momentum conservation and the second improvement is achieved by deriving and employing a higher-order source term based on a more accurate differentiation.

The source term in pressure Poisson equation (PPE) for ISPH is not unique; it has several formulations in the literature; one of them as a function of density variation and the other utilizes velocity divergence condition to formulate the source term. The former formulation with the density variation can keep a uniform particle distribution, although evaluated pressure include high unrealistic fluctuation. On the other hand, the formulation of the divergence-free condition evaluates much smoother pressure distribution, but density errors may occur due to particle clustering. Then, modified schemes have been proposed to satisfy the above two conditions: density invariant and divergence-free condition. Pozorski and Wawrenczuk [9] proposed a modified scheme, in which both the PPEs are solved separately at two intermediate states in each time step. Hu and Adams [16, 17] introduced internal iterations to satisfy both conditions accurately at the same moment. These modified schemes need to solve multiple PPEs in each time step, and these computational costs become expense compared to the conventional ISPH.

Recently, in the framework of MPS, there is a trend to introduce a higher-order source term in the PPE. Kondo and Koshizuka [18] proposed a new formulation with a source term composed by three parts; one is main part and another two terms related to error-compensating parts. Tanaka and Masunaga [19] introduced a similar high-order source term with two components incorporated with quasi-compressibility. Note that the number of PPEs per time step in their higher-order source term formulations is just one and its numerical cost is almost same as the original scheme. In this study, we reformulate a source term of the PPE which contains both contributions from velocity-divergence-free and density invariance conditions. Only one PPE per time step should be solved as the recent development in the MPS, but our formulation with a relaxation coefficient is unique. Note that, the relaxation coefficient depends on the initial particle distance, and a suitable relaxation coefficient can be obtained from the hydrostatic pressure calculations as a preanalysis. The accuracy and the efficiency of the proposed model are investigated in a couple of examples which are previously selected in published papers.

The turbulence models in the SPH are also important issue and the effects in the WCSPH have been nicely investigated by Violeau and Issa [20]. Lee et al. have introduced the same turbulence model such as -ε model into ISPH. Gotoh et al. [21] and Shao and Gotoh [22] introduced the static Smagorinsky model into the ISPH, and he discussed the effect of additional eddy viscosity shortly. In this paper, we also discuss the effect of eddy viscosity from our simulation results.

2. Typical Incompressible Smoothed Particle Hydrodynamics (ISPHs) Formulation

In this section, typical ISPH formulation, which is similar procedure in moving particle semi-implicit method (MPS) proposed by Koshizuka and Oka [24], is summarized. The main feature is that semi-implicit integration scheme is applied into particle discretized equations for the incompressible flow problem. The original idea of the semi-implicit scheme is called by projection method, which has been widely used in the finite difference method and in the finite element method. After the basic application of projection method into SPH is described here, several similar schemes will be categorized by the difference of treatment of PPE in the next section.

2.1. The Governing Equations for Incompressible Flow

In the Lagrange description, the continuity equation and the Navier-Stokes equations can be written as where and are density and kinematic viscosity of fluid, and are a velocity vector and pressure of fluid, respectively, is gravity acceleration, and indicates time. The turbulence stress is necessary to represent the effects of turbulence with coarse spatial grids, and its application into the particle simulation has been initially developed by Gotoh et al. [21]. In the most general incompressible flow approach, the density is assumed by a constant value with its initial value . Then, the aforementioned governing equations lead to

2.2. Projection Method

In the projection method [25], the velocity-pressure-coupled problem has been solved separately for velocity and pressure. Here, all the state variables may update from a previous time step to current time step. Below, superscripts and () indicate previous and current time step, respectively. In the first predictor step, intermediate state without pressure gradient is assumed, and its velocity field is indicated by . The intermediate velocity field can be evaluated by solving the following equation: Then, the following corrector step introduces an effect of remaining “current” pressure gradient term as follows: where indicates the incremental velocity from the predicted velocity .

The key point here is the evaluation of “current” pressure value. By taking the divergence of correction step (2.7) as Then, the incompressible condition (2.3) leads to By substituting (2.10) into (2.9), this leads to the following pressure Poisson equation (PPE): The above corrector step can be implemented by substituting the pressure gradient with the solution of PPE.

2.3. The SPH Methodology

A spatial discretization using scattered particles, which is based on the SPH, is summarized. First, a physical scalar function at a sampling point can be represented by the following integral form: where is a weight function called by smoothing kernel function in the SPH literature. In the smoothing kernel function, and are the distance between neighbor particles and smoothing length, respectively. For SPH numerical analysis, the integral equation (2.12) is approximated by a summation of contributions from neighbor particles in the support domain. where the subscripts and indicate positions of labeled particle, and and mean density and representative mass related to particle , respectively. Note that the triangle bracket means SPH approximation of a function . The gradient of the scalar function can be assumed by using the above defined SPH approximation as follows: Also, the other expression for the gradient can be represented by In this paper, quintic spline function [26] is utilized as a kernel function. where is , respectively, in two- and three-dimensional space. It has been observed that a cubic spline produces fluctuations in the pressure and velocity fields for fluid dynamics simulation, and the quintic spline shown in (2.16) gives more stable solutions.

2.4. Projection-Based ISPH Formulations

Here, the projection method for incompressible fluid problem, which is summarized in Section 2.2, is descretized into particle quantities based on the SPH methodology. For this purpose, the gradient of pressure and the divergence of velocity are approximated as follows: Although the Laplacian could be derived directly from the original SPH approximation of a function in (2.17), this approach may lead to a loss of resolution. Then, the second derivative of velocity for the viscous force and the Laplacian of pressure have been proposed by Morris et al. [5] by an approximation expression as follows: where is a parameter to avoid a zero dominator, and its value is usually given by . For the case of and , the Laplacian term is simplified as Similarly, the Laplacian of pressure in pressure Poisson equation (PPE) is given by The PPE after SPH interpolation is solved by a preconditioned (diagonal scaling) Conjugate Gradient (PCG) method [27] with a convergence tolerance (= 1.0 × 10−9).

2.5. Modeling of the Turbulence Stress

When dealing with the turbulent flows, the turbulent stress in (2.2), which are called by sub-particle scale stress in the particle simulations, needs to be modeled. In this paper, a large eddy simulation approach [21, 22] is used for modeling the turbulent stress as where and are the turbulence eddy viscosity and the turbulence kinetic energy, respectively. indicates the strain rate of the mean flow, and is the Kronecker delta. It is assumed in this paper that the eddy viscosity is modeled by the static Smagorinsky model as , in which is the Smagorinsky constant (taken as the analytical value in this paper), the constant is taken as , in which is the smoothing length defined in (2.12). The local strain rate can be calculated in the SPH formulation as Violeau and Issa [20].

2.6. Treatment of No-Slip Boundary Condition

The boundary condition on the rigid bodies has an important role to prevent penetration and to reduce error related to truncation of the kernel function. Takeda et al. [28] and Morris et al. [5] have introduced a special wall particle which can satisfy imposed boundary conditions. Recently, Bierbrauer et al. [29] described a consistent treatment of boundary conditions, utilizing the momentum equation to obtain approximations to velocity of image particles.

In our research, dummy particles technique, in which dummy particles are regularly distributed at the initial state and have zero velocity through the whole simulation process, is utilized just for simplicity in the implementation. In the following simulation, the pressure Poisson equation is solved for all particles including these dummy particles to get an enough repulsive force preventing penetration.

3. Stabilizations of Pressure Evaluation in Pressure Poisson Equation

Here, the pressure Poisson equation is reconsidered to overcome the error of artificial pressure fluctuation in the ISPH. The key points are related to the accuracy of density representation in SPH formulation and the treatment of pressure Poisson equation.

3.1. Keeping Divergence Free Scheme

Divergence free condition in the projection-based ISPH has been initially proposed by Cummins and Rudman [8]. Lee et al. [13] has applied it into the Reynolds turbulence model which uses an averaging in time. They called by “truly” incompressible SPH since the initial density is assumed constant for each particle. Then the divergence of the intermediate velocity has been used to calculate the PPE as mentioned above in (2.11). The PPE can be written in SPH approximation by substituting (2.18) and (2.21) as follows:

3.2. Keeping Density Invariance Scheme

The alternative scheme can be derived by using a density invariance condition [14, 15]. Here, the “particle” density in the SPH is evaluated by The particle position updates after each predictor step in the density invariance scheme and the particle density is updated on the intermediate particle positions. The intermediate particle density is indicated by . By assuming incompressibility condition with , the mass conservation law (2.1) can be rewritten for each particle as follows: By substituting (3.3) into (2.11) and using the SPH form, the PPE for ISPH can be approximately redefined by The main difference between the keeping divergence-free and keeping density-invariance scheme appeared in the source term of the PPE. Note that this keeping density-invariance scheme is analogous to the formulation in the MPS, although the MPS utilizes a “particle number” density instead of the particle density. The above two schemes have a relationship. Ataie-Ashtiani and Shobeyri [30] has converted from a PPE in the keeping density invariance scheme to a PPE in the keeping divergence-free scheme.

3.3. Combination Scheme of Both Divergence-Free and Density-Invariance Condition

A notable scheme was proposed by Hu and Adams [16]. The divergence-free and density-invariance conditions are sufficiently satisfied in their scheme. As they discussed, the divergence-free scheme calculates a smoothed pressure field but a large density variation will appear. Hu and Adams’s scheme includes an internal iteration at each step, and two kinds of PPEs should be solved until both the divergence-free and density-invariance conditions are approximately satisfied. According to Xu et al. [31], this scheme shows accurate and robust solutions, but total calculation time shows 4-5 times higher than that of the above two scheme.

3.4. Relaxed Density Invariance Scheme Incorporated with Divergence-Free Condition

Here, we proposed an efficient and robust ISPH scheme using both conditions without internal iterations. In the sense of physical observation, physical density should keep its initial value for incompressible flow. However, during numerical simulation, the “particle” density may change slightly from the initial value because the particle density is strongly dependent on particle locations in the SPH method. If the particle distribution can keep almost uniformity, the difference between “physical” and “particle” density may be vanishingly small. In other words, accurate SPH results in incompressible flow need to keep the uniform particle distribution. For this purpose, the different source term in pressure Poisson equation can be derived using the “particle” density. The SPH interpolations are introduced into the original mass conservation law before the perfect compressibility condition is applied By substituting (3.5) into (2.9) and using SPH form, the PPE can be represented by

Here, it is assumed that the current particle density is “hopefully” closed to initial density and the incremental particle density are defined by where the above integration scheme is called by the method of coordinate descent with a relaxation coefficient , and the PPE is modified as follows:

The similar equation, in which the density is replaced by a particle number density, was proposed by Losasso et al. [32], but they did not introduce the relaxation coefficient. Note that our proposed scheme couples the divergence-free and a relaxed density-invariance condition, and a special case using leads to the original divergence-free scheme. The effect of the relaxation coefficient will be tested in the later examples. Figure 1 shows the flow charts of these schemes to show the difference between existing and our proposed scheme. Similar modifications in the source term of PPE have been proposed in the MPS by Tanaka and Masunaga [19]. Recently, Khayyer and Gotoh [33] proposed a different higher-order source term without the unknown coefficient like the relaxation coefficient in this paper. It is important that the relaxation coefficient is strongly dependent on the initial particle distance, and the optimum value can be calibrated by a simple hydrostatic pressure test with the same initial particle distance as the final simulation model. The hydrostatic pressure test is called by preanalysis in this paper.

Figure 1: Flow charts in the ISPH schemes.
3.5. Tracking the Free Surface Boundary

Detection of free surface has an important role in the ISPH for free surface flow, because the pressure values on free surface particles should be equal to zero as Dirichlet boundary conditions of PPE. The method how to track the free surface may differ in each ISPH scheme.

Usually in the keeping density-invariance scheme, surface particles have been detected by referring the current particle density . The details have been discussed by Gotoh et al. [21], Shao and Gotoh [22], and Khayyer et al. [14, 15]. On the other hand, in the keeping divergence-free scheme, Lee et al. [13] proposed a new treatment with the divergence of a particle position vector. If the particle density keeps around its initial value, the former free surface detection method can be utilized. In our simulation, surface particles are simply judged by the total number of neighboring particles.

G. R. Liu and M. B. Liu [34] have investigated the number of neighboring particles to estimate an efficient variable smoothing length for the adaptive analysis. In the case of a simply cubic patterned lattice, is usually chosen as larger than 1.2 times of the initial particle distance . They showed that the number of neighboring particle within the support domain with for cubic spline kernel function should be about 21 in two-dimensional simulations. We checked a threshold for judging free surface particles for quintic spline kernel function with , and the threshold should be about 28 and 190 in 2D and 3D, respectively.

4. Preanalysis to Determine an Efficient Relaxation Coefficient

In this section, hydrostatic pressure evaluations are performed to investigate the effects of relaxation coefficient and to determine a suitable range of its value with reference to an initial particle distance.

The three particle models have been generated with different initial particles distances , and m as shown in Figure 2. The theoretical hydrostatic pressure is given by a law: with water density and a water height . Figure 3 shows pressure histories with different relaxation coefficients for each model. From this figure, the proper ranges of relaxation coefficient with initial particle distances , and m are approximately about (0.1 : 0.2), (0.0005 : 0.002) and (0.00005 : 0.0002), respectively. In this paper, a constant time is chosen by corresponding to Khayyer et al. [15]. Here note that the optimum parameter calibrated from this preanalysis can use the same value in the later examples if the model has the same particle resolution.

Figure 2: Schematic diagram of hydrostatic pressure at point A for three values of particle sizes.
Figure 3: Time history for pressure distributions under the effect of relaxation coefficient at different particle size models , and m, respectively.

5. Numerical Examples

Here, several numerical examples are solved by the current scheme with an efficient relaxation coefficient, which are calibrated by the hydrostatic pressure evaluation in the previous section.

5.1. The Effect of Relaxation Coefficient during Dam Break Simulation

A two-dimensional dam break analysis is performed to compare the proper relaxation coefficient between hydrostatic pressure and dam break simulation with the same particle distance. The geometry of a 2D dam break is shown in Figure 4, where the particle distance , the water width , water height , and the wall width  L.

Figure 4: Dam break analysis (unit: cm).

At first, Figure 5 shows the results of free surface detection by using the number of neighbors. It is seen that this simple free surface detection scheme is sufficiently accurate to determine the Dirichlet boundary conditions of the pressure Poisson equation and it is also suitable for any formulation of the pressure Poisson equation. The effects of relaxation coefficient are investigated by the density errors. Two boundary particles A and B, which position is marked in Figure 4, are selected to output an evaluated numerical density. Figure 6 shows the time histories of the density at particles A and B. From this observation with a particle distance 0.01 m, it seems that, too low relaxation coefficient (below 0.01) the density errors are high. A proper range of relaxation coefficient 0.1~0.25 leads to a stable solution. In addition, the density error fluctuations become serious when the relaxation coefficient is larger than this proper range. In the same way, the suitable ranges of relaxation coefficient for different particles distances m are evaluated by (0.0005 : 0.0025) and (0.00005 : 0.0001), respectively. Note that these proper ranges for each initial particle distance are close to the preevaluated proper ranges calibrated with the hydrostatic pressure test. Finally, the optimum values of relaxation coefficient in 2D dam break analysis for particle sizes , and 0.0025 m are determined by 0.15, 0.001, and 0.00006, respectively.

Figure 5: Detect free surface numerically using the number of neighbouring particles for using and for quintic spline kernel function at initial particle size .
Figure 6: Effect of relaxation coefficient on the density evaluation error at (a) particles A and (b) Particle B, respectively.

Figure 7 shows the pressure distributions for three models with different initial particle distances  m. A suitable relaxation coefficient is utilized for each model. The first water impact at the right wall generates highest pressure, and it returns in the form of a jet. Then, it becomes a stable state after two more water impacts act on both side walls. The snapshots for water impact, after the first water impact, reversing jet and water stable state are captured from each model with different particle distance. These snapshots at the same time show similar water shape. The pressure histories at the right corner are plotted in Figure 8. Although unrealistic pressure fluctuation appears in the case of lowest resolution model with  m, a similar tendency of pressure history can get from different resolution models. Adjusting suitable relaxation coefficient can increase the pressure smoothness. The hydrostatic pressure after this dam break analysis is analytically evaluated by 1000 , and our evaluated pressure after 4 seconds looks to converge into the analytical hydrostatic value. In this example, water front speed is plotted in Figure 9. Our results shows a good agreement with the experimental results obtained by Koshizuka and Oka [24] and Martin and Moyce [35], moreover the numerical results obtained by Lee et al. [23].

Figure 7: Pressure distributions for dam break analysis from different particle size models  m, respectively.
Figure 8: Time history for pressure distribution at different particle sizes with proper relaxation coefficient.
Figure 9: Comparison on dam break fronts in the dam break analysis.
5.2. Comparison of Configurations and Pressure during Dam Break Simulation

Next, water configurations and pressure distribution are compared with an experimental data by Zhou et al. [36] and with a result by original incompressible SPH, which is the same as a special case of our proposed scheme with . The schematic diagram is the same as Zhou’s experiment is shown in Figure 10, and the pressure measuring point is located at a point on the right wall (0.16 m). The particle initial distance is selected as . A proper relaxation coefficient for this resolution is selected by which is the same optimum value in the hydrostatic pressure test, and then the numerical solution is compared to the truly incompressible scheme with . Figure 11 shows the comparison results of the snapshots with pressure distribution from the initial state to the final stable state. The snapshots from each scheme are captured at the first water impact, running up along the right wall reversing development of splash-up and the stable state. Although the wave configurations show similarities, the pressure value from the truly incompressible scheme is less than that from our proposed scheme. In addition, the total volume of the water at the final stable state is compared between the proposed scheme and truly incompressible scheme using the water height. It seems that the proposed scheme conserves the total volume compared to the theoretical value of height about 0.22 m, while the truly incompressible scheme cannot conserve the volume at the final stable state. Figure 12 shows the comparison of pressure history at the right corner among our proposed scheme result with a proper relaxation coefficient, result from the truly incompressible scheme, and experimental data by Zhou et al. [36]. Although the pressure level from the truly incompressible scheme is lower than the experimental data in the entire simulation period, the evaluated pressure from our proposed scheme shows a good agreement with the experimental data. In this figure, imaginary pressure peak is evaluated around in the results without turbulence model. The combination with the proposed stabilized ISPH and turbulence model generates smoothed and accurate pressure distribution.

Figure 10: Dam break simulation corresponding to experiment by Zhou et al. [36].
Figure 11: Pressure distribution for dam break analysis by (a) stabilized ISPH with eddy viscosity effect, and (b) original ISPH without eddy viscosity effect, respectively.
Figure 12: (a) Comparison between the current stabilized model including/excluding eddy viscosity, divergence-free scheme condition, and experimental data by Zhou et al. [36].
5.3. 3D Dam Break Flow with an Obstacle

The last application is one of the benchmark test suggested by SPH European Research Interest Community (SPHERIC). The experimental tests on a dam break flow with an obstacle was carried out at the Maritime Research Institute Netherlands (MARIN) as reported by Kleefsman et al. [37].

Figures 13 and 14 show geometry of the experimental test and locations of pressure sensor, respectively. While ps1 to ps8 sensors were used in the experimental test, here only odd numbers of pressure sensor are utilized for the comparison. In the numerical modeling, the initial particle distance is fixed at 0.01 m for both regions of water and wall. The total number of particles is about 1.4 millions, and 0.67 million particles are located in the water. In order to evaluate an efficient relaxation coefficient, the same procedure as two-dimensional cases is applied. First, the hydrostatic pressure test has been implemented by using the same initial particle distance and time increment in the 3D dam break problem. Then the optimum relaxation coefficient was fixed by 0.1.

Figure 13: Geometry of the 3D dam break experiment.
Figure 14: Locations of the pressure sensor on the obstacle.

The pressure time history on the front (ps1 and ps3) and top (ps5 and ps7) is shown in Figure 15. Figures 16 and 17 show the snapshots with particle pressure values and labels related to free surface, respectively. In Figure 16, the numerical solutions by our proposed relaxed density invariant scheme (stabilizes ISPH) are compared with Kleefsman’s experimental results and numerical results by the keeping divergence-free scheme with (original ISPH). The first impact occurred at about 0.42 s both in the numerical and experimental test, although the time of secondary hit has about 0.5 s difference (0.45 s and 0.50 s, resp.). That is, our solution shows small delay as the time goes. The pressure resulting from the keeping divergence-free scheme shows lower value during the simulation, although a smooth pressure distribution can be generated as in Figure 16. It seems that the keeping divergence scheme cannot keep the total volume of water. On the other hand, except for the local pressure oscillation especially at ps5 and ps7, the pressure histories by our proposed scheme show good agreement with the experimental results.

Figure 15: Time history of Pressure for 3D dam break at ps1, ps3, ps5, and ps7.
Figure 16: Time sequence of 3D dam break simulation by (a) stabilized ISPH method with eddy viscosity effect and (b) without eddy viscosity effect.
Figure 17: Time sequence for detection free surface in 3D dam break simulation by (a) stabilized ISPH method with eddy viscosity effect and (b) without eddy viscosity effect.

Lee et al. [38] has been simulated to the same problem, and they have discussed the difference between weakly compressible SPH and their proposed truly incompressible SPH that is one of the keeping divergence-free scheme. According to their result, the weakly compressible SPH shows a critical error in the pressure and their truly incompressible SPH solution has the similar tendency as our results. However, the original ISPH scheme cannot keep the total volume as far as we have checked.

6. Conclusion

A stabilized incompressible smoothed particle hydrodynamics is proposed to simulate free surface flow. The modification is appeared in the source term of pressure Poisson equation, and the idea is similar to the recent development in Moving Particles Semi-implicit method (MPS). Although only one set of linear equations should be solved to evaluate pressure at each particle, both the velocity divergence-free condition and the density invariance condition can be approximately satisfied. The additional parameter is the relaxation coefficient, and its value can be calibrated by a simple hydrostatic simulation with a regular initial particle distribution. It has a uniform tendency that the relaxation coefficient becomes smaller due to decrease in the initial particle distance. The efficiency and its accuracy have been tested by the dam break in two- and three-dimension simulations compared to their reference solutions. Our proposed scheme shows the clear advantage to keep the total volume by comparing the original ISPH, and it may contribute to have an accurate pressure value. However, it still has an artificial oscillation in the pressure value with original viscosity. The additional viscosity based on the Subparticle Scale turbulence model shows an important role to generate smoother pressure distribution and to decrease the number of isolated particles after the splash in the dam break problems.


The first author is partially supported by the Ministry of Education, Science, Sports and Culture of Japan (Grant-in-Aid for Young Scientist (B) 2176036).


  1. L. B. Lucy, “A numerical approach to the testing of the fusion process,” Astronomical Journal, vol. 88, pp. 1013–1024, 1977. View at Google Scholar
  2. R. A. Gingold and J. J. Monaghan, “Smoothed particle hydrodynamics: theory and application to non-spherical stars,” Monthly Notices of the Royal Astronomical Society, vol. 181, pp. 375–389, 1977. View at Google Scholar · View at Zentralblatt MATH
  3. J. J. Monaghan, “Smoothed particle hydrodynamics,” Annual Review of Astronomy and Astrophysics, vol. 30, no. 1, pp. 543–574, 1992. View at Publisher · View at Google Scholar · View at Scopus
  4. J. J. Monaghan, “Simulating free surface flows with SPH,” Journal of Computational Physics, vol. 110, no. 2, pp. 399–406, 1994. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  5. J. P. Morris, P. J. Fox, and Y. Zhu, “Modeling low Reynolds number incompressible flows using SPH,” Journal of Computational Physics, vol. 136, no. 1, pp. 214–226, 1997. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  6. J. J. Monaghan, “Heat conduction with discontinuous conductivity,” in Applied Mathematics Reports and Preprints, Monash University, 1995. View at Google Scholar
  7. N. Okahci, A. Hirota, S. Izawa, Y. Fukunishi, and H. Higuchi, “SPH simulation of pulsating pipe flow at a junction,” in Proceedings of the 1st International Symposium on Advanced Fluid Information, pp. 388–391, 2001.
  8. S. J. Cummins and M. Rudman, “An SPH projection method,” Journal of Computational Physics, vol. 152, no. 2, pp. 584–607, 1999. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  9. J. Pozorski and A. Wawrenczuk, “SPH computation of incompressible viscous flows,” Journal of Theoretical and Applied Mechanics, vol. 40, pp. 917–937, 2002. View at Google Scholar
  10. S. Shao and E. Y. M. Lo, “Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface,” Advances in Water Resources, vol. 26, no. 7, pp. 787–800, 2003. View at Publisher · View at Google Scholar · View at Scopus
  11. X. Y. Hu and N. A. Adams, “Angular-momentum conservative smoothed particle dynamics for incompressible viscous flows,” Physics of Fluids, vol. 18, no. 10, Article ID 101702, 4 pages, 2006. View at Publisher · View at Google Scholar · View at Scopus
  12. M. Ellero, M. Serrano, and P. Espanol, “Incompressible smoothed particle hydrodynamics,” Journal of Computational Physics, vol. 226, no. 2, pp. 1731–1752, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  13. E.-S. Lee, C. Moulinec, R. Xu, D. Violeau, D. Laurence, and P. Stansby, “Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method,” Journal of Computational Physics, vol. 227, no. 18, pp. 8417–8436, 2008. View at Publisher · View at Google Scholar
  14. A. Khayyer, H. Gotoh, and S. Shao, “Corrected incompressible SPH method for accurate water-surface tracking in breaking waves,” Coastal Engineering, vol. 55, no. 3, pp. 236–250, 2008. View at Publisher · View at Google Scholar · View at Scopus
  15. A. Khayyer, H. Gotoh, and S. Shao, “Enhanced predictions of wave impact pressure by improved incompressible SPH methods,” Applied Ocean Research, vol. 31, no. 2, pp. 111–131, 2009. View at Publisher · View at Google Scholar · View at Scopus
  16. X. Y. Hu and N. A. Adams, “An incompressible multi-phase SPH method,” Journal of Computational Physics, vol. 227, no. 1, pp. 264–278, 2007. View at Publisher · View at Google Scholar · View at Scopus
  17. X. Y. Hu and N. A. Adams, “A constant-density approach for incompressible multi-phase SPH,” Journal of Computational Physics, vol. 228, no. 6, pp. 2082–2091, 2009. View at Publisher · View at Google Scholar · View at Scopus
  18. M. Kondo and S. Koshizuka, “Improvement of stability in moving particle semi-implicit method,” International Journal for Numerical Methods in Fluids, vol. 65, no. 6, pp. 638–654, 2011. View at Publisher · View at Google Scholar
  19. M. Tanaka and T. Masunaga, “Stabilization and smoothing of pressure in MPS method by quasi-compressibility,” Journal of Computational Physics, vol. 229, no. 11, pp. 4279–4290, 2010. View at Publisher · View at Google Scholar · View at Scopus
  20. D. Violeau and R. Issa, “Numerical modelling of complex turbulent free-surface flows with the SPH method: an overview,” International Journal for Numerical Methods in Fluids, vol. 53, no. 2, pp. 277–304, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  21. H. Gotoh, S. Shao, and T. Memita, “SPH-LES model for numerical investigation of wave interaction with partially immersed breakwater,” Coastal Engineering Journal, vol. 46, no. 1, pp. 39–63, 2004. View at Publisher · View at Google Scholar · View at Scopus
  22. S. Shao and H. Gotoh, “Turbulence particle models for tracking free surfaces,” Journal of Hydraulic Research, vol. 43, no. 3, pp. 276–289, 2005. View at Publisher · View at Google Scholar · View at Scopus
  23. B. H. Lee, J. C. Park, M. H. Kim, S. J. Jung, M. C. Ryu, and Y. S. Kim, “Numerical simulation of impact loads using a particle method,” Ocean Engineering, vol. 37, no. 2-3, pp. 164–173, 2010. View at Publisher · View at Google Scholar · View at Scopus
  24. S. Koshizuka and Y. Oka, “Moving-particle semi-implicit method for fragmentation of incompressible fluid,” Nuclear Science and Engineering, vol. 123, no. 3, pp. 421–434, 1996. View at Google Scholar · View at Scopus
  25. A. J. Chorin, “Numerical solution of the Navier-Stokes equations,” Mathematics of Computation, vol. 22, pp. 745–762, 1968. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  26. I. J. Schoenberg, “Contributions to the problem of approximation of equidistant data by analytic functions,” Quarterly of Applied Mathematics, vol. 4, pp. 45–99, 1946. View at Google Scholar · View at Zentralblatt MATH
  27. J. A. Meijerink and H. A. van der Vorst, “An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix,” Mathematics of Computation, vol. 31, no. 137, pp. 148–162, 1977. View at Google Scholar
  28. H. Takeda, S. Miyama, and M. Sekiya, “Numerical simulation of viscous flow by smoothed particle hydrodynamics,” Progress of Theoretical Physics, vol. 92, no. 5, pp. 939–960, 1994. View at Publisher · View at Google Scholar
  29. F. Bierbrauer, P. C. Bollada, and T. N. Phillips, “A consistent reflected image particle approach to the treatment of boundary conditions in smoothed particle hydrodynamics,” Computer Methods in Applied Mechanics and Engineering, vol. 198, no. 41–44, pp. 3400–3410, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  30. B. Ataie-Ashtiani and G. Shobeyri, “Numerical simulation of landslide impulsive waves by incompressible smoothed particle hydrodynamics,” International Journal for Numerical Methods in Fluids, vol. 56, no. 2, pp. 209–232, 2008. View at Publisher · View at Google Scholar
  31. R. Xu, P. Stansby, and D. Laurence, “Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach,” Journal of Computational Physics, vol. 228, no. 18, pp. 6703–6725, 2009. View at Publisher · View at Google Scholar
  32. F. Losasso, J. O. Talton, N. Kwatra, and R. Fedkiw, “Two-way coupled SPH and particle level set fluid simulation,” IEEE Transactions on Visualization and Computer Graphics, vol. 14, no. 4, pp. 797–804, 2008. View at Publisher · View at Google Scholar · View at Scopus
  33. A. Khayyer and H. Gotoh, “Enhancement of stability and accuracy of the moving particle semi-implicit method,” Journal of Computational Physics, vol. 230, no. 8, pp. 3093–3118, 2011. View at Publisher · View at Google Scholar
  34. G. R. Liu and M. B. Liu, Smoothed Particle Hydrodynamic, a Meshfree Particle Method, World Scientific, 2003.
  35. J. C. Martin and W. J. Moyce, “An experimental study of the collapse of liquid columns on a rigid horizontal plane,” Philosophical Transactions of the Royal Society of London A, vol. 244, no. 882, pp. 312–324, 1952. View at Publisher · View at Google Scholar
  36. Z. Q. Zhou, J. O. De Kat, and B. Buchner, “A nonlinear 3-D approach to simulate green water dynamics on deck,” in Proceedings of the 7th International Conference on Numerical Ship Hydrodynamics, Nantes, France, July 1999.
  37. K. M. T. Kleefsman, G. Fekken, A. E. P. Veldman, B. Iwanowski, and B. Buchner, “A volume-of-fluid based simulation method for wave impact problems,” Journal of Computational Physics, vol. 206, no. 1, pp. 363–393, 2005. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  38. E. S. Lee, D. Violeau, R. Issa, and S. Ploix, “Application of weakly compressible and truly incompressible SPH to 3-D water collapse in waterworks,” Journal of Hydraulic Research, vol. 48, no. 1, pp. 50–60, 2010. View at Publisher · View at Google Scholar · View at Scopus