A Stabilized Incompressible SPH Method by Relaxing the Density Invariance Condition

Abstract

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 [4–7]. 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 [8–17].

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 𝐷𝜌𝐷𝑡+𝜌∇⋅𝐮=0, (2.1)𝐷𝐮1𝐷𝑡=−𝜌∇𝑃+𝜐∇²1𝐮+𝐠+𝜌∇⋅𝝉, (2.2) 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∇⋅𝐮=0, (2.3)𝐷𝐮1𝐷𝑡=−𝜌⁰∇𝑃+𝜐∇²1𝐮+𝐠+𝜌⁰∇⋅𝝉. (2.4)

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 (𝑛+1) 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:𝐮^∗−𝐮^𝑛Δ𝑡=𝜐∇²𝐮^𝑛1+𝐠+𝜌⁰∇⋅𝝉, (2.5)(Predictor)∶𝐮^∗=𝐮^𝑛+Δ𝑡𝜐∇²𝐮^𝑛1+𝐠+𝜌⁰∇⋅𝝉, (2.6) Then, the following corrector step introduces an effect of remaining “current” pressure gradient term as follows:𝐮^𝑛+1−𝐮^∗1Δ𝑡=−𝜌⁰∇𝑃^𝑛+1, (2.7)(Corrector)∶𝐮^𝑛+1=𝐮^∗+Δ𝐮^∗=𝐮^∗1−Δ𝑡𝜌⁰∇𝑃^𝑛+1, (2.8) 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𝐮∇⋅^𝑛+1−𝐮^∗1Δ𝑡=−∇⋅𝜌⁰∇𝑃^𝑛+1. (2.9) Then, the incompressible condition (2.3) leads to𝐮∇⋅^𝑛+1Δ𝑡=0. (2.10) By substituting (2.10) into (2.9), this leads to the following pressure Poisson equation (PPE):∇²𝑃^𝑛+1=𝜌⁰∇⋅𝐮^∗Δ𝑡. (2.11) 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:𝜙𝑥_𝑖,𝑡=𝑊𝑥_𝑖−𝑥_𝑗,ℎ𝜙𝑥_𝑗,𝑡𝑑𝑣=𝑊𝑟_𝑖𝑗,ℎ𝜙𝑥_𝑗,𝑡𝑑𝑣, (2.12) 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.𝜙𝑥_𝑖,𝑡≈𝜙_𝑖=_𝑗𝑚_𝑗𝜌_𝑗𝑊𝑟_𝑖𝑗,ℎ𝜙𝑥_𝑗,𝑡, (2.13) 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:∇𝜙𝑥_𝑖≈∇𝜙_𝑖=1𝜌_𝑖𝑗𝑚_𝑗𝜙_𝑗−𝜙_𝑖∇𝑊𝑟_𝑖𝑗,ℎ. (2.14) Also, the other expression for the gradient can be represented by∇𝜙𝑥_𝑖≈⟨∇𝜙_𝑖⟩=𝜌_𝑖𝑗𝑚_𝑗𝜙_𝑗𝜌²_𝑗+𝜙_𝑖𝜌²_𝑖∇𝑊𝑟_𝑖𝑗,ℎ. (2.15) In this paper, quintic spline function [26] is utilized as a kernel function.𝑊𝑟_𝑖𝑗,ℎ=𝛽_𝑑𝑟3−^𝑖𝑗ℎ⁵𝑟−62−^𝑖𝑗ℎ⁵𝑟+151−^𝑖𝑗ℎ⁵,0≤𝑟_𝑖𝑗𝑟<ℎ,3−^𝑖𝑗ℎ⁵𝑟−62−^𝑖𝑗ℎ⁵,ℎ≤𝑟_𝑖𝑗𝑟<2ℎ,3−^𝑖𝑗ℎ⁵,2ℎ≤𝑟_𝑖𝑗<3ℎ,0,𝑟_𝑖𝑗≥3ℎ, (2.16) where 𝛽_𝑑 is 7/478𝜋ℎ²and3/358𝜋ℎ³, 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:∇𝑝𝑥_𝑖≈∇𝑝_𝑖=𝜌_𝑖𝑗𝑚_𝑗𝑝_𝑗𝜌²_𝑗+𝑝_𝑖𝜌²_𝑖∇𝑊𝑟_𝑖𝑗,ℎ, (2.17)∇⋅𝐮𝑥_𝑖≈∇⋅𝐮_𝑖=1𝜌_𝑖𝑗𝑚_𝑗𝐮_𝑗−𝐮_𝑖⋅∇𝑊𝑟_𝑖𝑗,ℎ. (2.18) 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:∇⋅(𝜐∇𝐮)𝑥_𝑖≈∇⋅𝜐∇⋅𝐮_𝑖=_𝑗𝑚_𝑗𝜌_𝑖𝜐_𝑖+𝜌_𝑗𝜐_𝑗𝜌_𝑖𝜌_𝑗𝑟_𝑖𝑗∇𝑊𝑟_𝑖−𝑟_𝑗,ℎ𝑟²_𝑖𝑗+𝜂²𝐮_𝑖𝑗, (2.19) where 𝜂 is a parameter to avoid a zero dominator, and its value is usually given by 𝜂²=0.0001 ℎ². For the case of 𝜐_𝑖=𝜐_𝑗 and 𝜌_𝑖=𝜌_𝑗, the Laplacian term is simplified as∇⋅𝜐∇⋅𝐮_𝑖=2𝜐_𝑖𝜌_𝑖𝑗𝑚_𝑗𝑟_𝑖𝑗∇𝑊𝑟_𝑖−𝑟_𝑗,ℎ𝑟²_𝑖𝑗+𝜂²𝐮_𝑖𝑗. (2.20) Similarly, the Laplacian of pressure in pressure Poisson equation (PPE) is given by∇²𝑝𝑥_𝑖≈∇²𝑝_𝑖=2𝜌_𝑖𝑗𝑚_𝑗𝑝_𝑖𝑗𝑟_𝑖𝑗∇𝑊𝑟_𝑖−𝑟_𝑗,ℎ𝑟²_𝑖𝑗+𝜂². (2.21) The PPE after SPH interpolation is solved by a preconditioned (diagonal scaling) Conjugate Gradient (PCG) method [27] with a convergence tolerance (= 1.0 × 10⁻⁹).

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𝜏_𝐼𝐽𝜌₀=2𝜐_𝑇𝑆_𝐼𝐽−23𝑘𝛿_𝐼𝐽, (2.22) 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 𝐶_𝑠=0.2 is the Smagorinsky constant (taken as the analytical value in this paper), the constant Δ is taken as 2ℎ, in which ℎ is the smoothing length defined in (2.12). The local strain rate |𝑆|=(2𝑆_𝐼𝐽𝑆_𝐼𝐽)^1/2 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:∇²𝑝_𝑖^𝑛+1=𝜌⁰Δ𝑡∇⋅𝐮^∗_𝑖. (3.1)

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𝜌𝑥_𝑖≈𝜌_𝑖=_𝑗𝑚_𝑗𝑊𝑟_𝑖𝑗,ℎ=_𝑗𝑚_𝑗𝑊_𝑖𝑗. (3.2) 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 ⟨𝜌_𝑖^𝑛+1⟩=𝜌⁰, the mass conservation law (2.1) can be rewritten for each particle as follows:𝜌^∗_𝑖−𝜌⁰Δ𝑡+𝜌⁰∇⋅𝐮^∗_𝑖=0. (3.3) By substituting (3.3) into (2.11) and using the SPH form, the PPE for ISPH can be approximately redefined by⟨∇²𝑝_𝑖^𝑛+1𝜌⟩=⁰−⟨𝜌^∗_𝑖⟩Δ𝑡². (3.4) 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∇⋅𝐮_𝑖^𝑛+11=−𝜌⁰𝜌_𝑖^𝑛+1−𝜌^𝑛_𝑖Δ𝑡. (3.5) By substituting (3.5) into (2.9) and using SPH form, the PPE can be represented by∇²𝑝_𝑖^𝑛+1=𝜌⁰Δ𝑡∇⋅𝐮^∗_𝑖+𝜌_𝑖^𝑛+1−𝜌^𝑛_𝑖Δ𝑡². (3.6)

Here, it is assumed that the current particle density is “hopefully” closed to initial density and the incremental particle density ⟨Δ𝜌^𝑛_𝑖⟩ are defined by𝜌_𝑖^𝑛+1=𝜌^𝑛_𝑖+Δ𝜌^𝑛_𝑖≈𝜌⁰, (3.7)Δ𝜌^𝑛_𝑖∶=𝛼𝜌⁰−𝜌^𝑛_𝑖, (3.8) where the above integration scheme is called by the method of coordinate descent with a relaxation coefficient 𝛼(0≤𝛼≤1), and the PPE is modified as follows:∇²𝑝_𝑖^𝑛+1=𝜌⁰Δ𝑡∇⋅𝐮^∗_𝑖𝜌+𝛼⁰−𝜌^𝑛_𝑖Δ𝑡². (3.9)

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 𝛼=0 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 𝑘=2 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 𝑘=3, 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 𝑑₀=0.01,0.005, and 0.0025m as shown in Figure 2. The theoretical hydrostatic pressure is given by a law: 𝑝=𝜌𝑔ℎ(=980N/m²) with water density 𝜌=1000kg/m³ and a water height ℎ=0.1m. 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 𝑑₀=0.01,0.005, and 0.0025m 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 Δ𝑡=0.1𝑑₀ 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.

(a) At particle size 𝑑0=0.01 m

(b) At particle size 𝑑0=0.005 m

(c) At particle size 𝑑0=0.0025 m

Figure 3 

Time history for pressure distributions under the effect of relaxation coefficient at different particle size models 𝑑0=0.01,0.005, and 0.0025m, 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 𝑑₀=0.01m, the water width 𝐿=0.20m, water height 𝐻_𝑤=2.5L, and the wall width 𝑊_𝑙=5 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 𝑑₀=0.005,0.0025m 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 𝑑₀=0.01,0.005, 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 𝑘=3 and ℎ=1.2𝑑0 for quintic spline kernel function at initial particle size 𝑑0=0.005m.

(a) At particle A

(b) At particle B

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 𝑑₀=0.01,0.005,and0.0025 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 𝑑₀=0.01 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 N/m², 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 𝑑0=0.01,0.005,and0.0025 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 𝛼=0. 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 𝑑₀=0.005m. A proper relaxation coefficient for this resolution is selected by 𝛼=0.001 which is the same optimum value in the hydrostatic pressure test, and then the numerical solution is compared to the truly incompressible scheme with 𝛼=0. 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 (𝛼=0) 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 𝑡=2.04(𝑡(𝑔/ℎ)^0.5=8.24) 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 (𝛼=0) 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 ps₁ to ps₈ 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 𝑑₀=0.01m and time increment Δ𝑡=0.001s 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 (ps₁ and ps₃) and top (ps₅ and ps₇) 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 𝛼=0 (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 ps₅ and ps₇, the pressure histories by our proposed scheme show good agreement with the experimental results.

(a) Output value at ps1

(b) Output value at ps3

(c) Output value at ps5

(d) Output value at ps7

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.

Acknowledgments

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