A Stabilized Incompressible SPH Method by Relaxing the Density Invariance Condition

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


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

Typical Incompressible Smoothed Particle Hydrodynamics (ISPHs) Formulation
In this section, typical ISPH formulation, which is similar procedure in moving particle semiimplicit 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.

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, u and P are a velocity vector and pressure of fluid, respectively, g is gravity acceleration, and t 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 ρ 0 .Then, the aforementioned governing equations lead to

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 n and n 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 u * .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 Δu * indicates the incremental velocity from the predicted velocity u * .
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.

The SPH Methodology
A spatial discretization using scattered particles, which is based on the SPH, is summarized.First, a physical scalar function φ x i , t at a sampling point x i can be represented by the following integral form: where W is a weight function called by smoothing kernel function in the SPH literature.In the smoothing kernel function, r ij |x i −x j | and h 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 i and j indicate positions of labeled particle, and ρ j and m j mean density and representative mass related to particle j, respectively.Note that the triangle bracket φ i 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

2.15
In this paper, quintic spline function 26 is utilized as a kernel function.
where β d is 7/478πh 2 and 3/358πh 3 , 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.

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 η 2 0.0001 h 2 .For the case of υ i υ j and ρ i ρ j , 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 .

Modeling of the Turbulence Stress
When dealing with the turbulent flows, the turbulent stress in 2.2 , which are called by subparticle 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 υ T and k are the turbulence eddy viscosity and the turbulence kinetic energy, respectively.S IJ indicates the strain rate of the mean flow, and δ IJ is the Kronecker delta.It is assumed in this paper that the eddy viscosity is modeled by the static Smagorinsky model as υ T C s Δ 2 |S|, in which C s 0.2 is the Smagorinsky constant taken as the analytical value in this paper , the constant Δ is taken as 2h, in which h is the smoothing length defined in 2.12 .The local strain rate |S| 2S IJ S IJ 1/2 can be calculated in the SPH formulation as Violeau and Issa 20 .

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

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.

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 ρ 0 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.1

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 Journal of Applied Mathematics particle density is indicated by ρ * i .By assuming incompressibility condition with ρ n 1 i ρ 0 , 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.

Combination Scheme of Both Divergence-Free and Density-Invariance Condition
A notable scheme was proposed by Hu and Adams 16 .The divergence-free and densityinvariance 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.

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 Δρ n i are defined by 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: 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.

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 ρ i .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.   a simply cubic patterned lattice, h is usually chosen as larger than 1.2 times of the initial particle distance d 0 .They showed that the number of neighboring particle within the support domain kh with k 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 k 3, and the threshold should be about 28 and 190 in 2D and 3D, respectively.

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 d 0 0.01, 0.005, and 0.0025 m as shown in Figure 2. The theoretical hydrostatic pressure is given by a law: p ρgh 980 N/m 2 with water density ρ 1000 kg/m 3 and a water height h 0.1 m. Figure 3 shows pressure histories with different relaxation coefficients for

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.

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 d 0 0.01 m, the water width L 0.20 m, water height H w 2.5 L, and the wall width W l 5 L.
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 d 0 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 d 0 0.005, 0.0025 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 d 0 0.01, 0.005, and 0.0025 m are determined by 0.15, 0.001, and 0.00006, respectively.Figure 7 shows the pressure distributions for three models with different initial particle distances d 0 0.01, 0.005, and 0.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 d 0 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 2 , 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 .

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 d 0 0.005 m.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 Journal of Applied Mathematics 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 t 2.04 t g/h 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.

3D Dam Break Flow with an Obstacle
The Figures 13 and 14 show geometry of the experimental test and locations of pressure sensor, respectively.While ps 1 to ps 8 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 d 0 0.01 m and time increment Δt 0.001 s in the 3D dam break problem.Then the optimum relaxation coefficient α was fixed by 0.1.
The pressure time history on the front ps 1 and ps 3 and top ps 5 and ps 7 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 5 and ps 7 , the pressure histories by our proposed scheme show good agreement with the experimental results.
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.

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

Figure 1 :
Figure 1: Flow charts in the ISPH schemes.

Figure 2 :
Figure 2: Schematic diagram of hydrostatic pressure at point A for three values of particle sizes.

cFigure 3 :
Figure 3: Time history for pressure distributions under the effect of relaxation coefficient at different particle size models d 0 0.01, 0.005, and 0.0025 m, respectively.

Figure 5 :Figure 6 :
Figure 5: Detect free surface numerically using the number of neighbouring particles for using k 3 and h 1.2 d 0 for quintic spline kernel function at initial particle size d 0 0.005 m.

d 0 Figure 8 :Figure 9 :Figure 10 :
Figure 8: Time history for pressure distribution at different particle sizes with proper relaxation coefficient.

Figure 12 :
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 .

Figure 14 :
Figure 14: Locations of the pressure sensor on the obstacle.

Figure 17 :
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.
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.
Liu and M. B.Liu 34have investigated the number of neighboring particles to estimate an efficient variable smoothing length for the adaptive analysis.In the case of Journal of Applied Mathematics 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 .
Geometry of the 3D dam break experiment.