A study of the stability properties in simulation of wave propagation with SPH method

When ordinary Smoothed Particle Hydrodynamics (SPH) method is used to simulate wave propagation in a wave tank, it is usually observed that the wave height decays and the wave length elongates along the direction of wave propagation. Accompanied with this phenomenon, the pressure under water decays either and shows a big oscillation simultaneously. The reason is the natural potential tensile instability of modeling water motion with ordinary SPH which is caused by particle negative stress in the computation. To deal with the problems, a new sextic kernel function is proposed to reduce this instability. An appropriate smooth length is given and its computation criterion is also suggested. At the same time, a new kind dynamic boundary condition is introduced. Based on these improvements, the new SPH method named stability improved SPH (SISPH) can simulate the wave propagation well. Both the water surface and pressure can be well expressed and the oscillation of pressure is nearly eliminated. Compared with other improved methods, SISPH can truly reveal the physical reality without bringing some new problems in a simple way.


Introduction
SPH was firstly created by Lucy (1977), Gingold and Monaghan (1977) to simulate nonaxisymmetric phenomena in astrophysics involving a compressible gas which always has a positive gas pressure. The governing equations which control the motion of the particles result in particles repelling each other with equal and opposite forces. Later on, it has been applied to fluid simulation. However, the pressure of fluid particles is not always positive. When the particle stress is negative, the particles no longer repel, but attract. This can easily cause particles to pair up and clump, which results in failure of the simulation. This SPH instability is a consequence of the particles being free to move under negative stress. It exists in the fluid model with any dimensional accuracy and time accuracy.
In order to improve the stability of SPH, Phillips and Monaghan (1985) added a constant stress tensor on the original one. It is a simple way to remove the instability, but it cannot be applied to all other models. Dyka and Ingel (1995) calculated the pressure at the stress points which are away from SPH nodes. Randles and Libersky (2000) expanded this method to multidimensional space. However, the tensile instability cannot be removed completely and some compressive instability is brought. From the stability analysis of the wave simulation with SPH, Morris (1994Morris ( , 1995 found that the stability properties of SPH are improved by using the kernels whose first derivatives go to zero with decreasing inter particle spacing more rapidly. In other words, using higher order spline kernel function can improve the stability of SPH. Chen et al. (1999) proposed corrective smoothed particle method (CSPM) through Taylor expansion. It increased the derivative accuracy of SPH to improve the stability. But serious splash of particles is observed at fluid free surface during the computation. Monaghan (2000), Gray et al. (2001) added an artificial stress in the momentum equation in order to repel the particles when a clustering of particles occurred in gases and MHD. Zhang et al. (2010) talked about the relationship between smooth length and the stability of SPH, and gave a suggestion about how to select an appropriate smooth length to maintain the error system of SPH stable. Guilcher et al. (2007) found wave height decaying and wave length stretch during wave propagation, and combined Godunov scheme and CSPM to simulate the wave propagation. The result demonstrated that this method can simulate the wave propagation well. However, because of the complexity of this method, the viscosity term is not considered. In addition, due to use of the HLLC approximate Riemann solver in Godunov scheme, some physical discontinuities are erased. In order to reduce this effect, smaller particle space is required (the space is proposed as 0.002 m in the article). Meanwhile for a certain particle space, the stable computation time is limited, because the correction in Riemann solver accumulated to a certain degree will affect the stability of SPH. Yang and Liu (2012) proposed an enhanced cubical spline kernel which improved a lot the ordinary SPH stability in hydrostatic problems, but it was not good enough in wave propagation simulation. Like the method used by Guilcher et al. (2007), Gao (2011) simulated wave propagation using SPH combined with CSPM and another rapid non-conservative Riemann solver, and the result was also good. The viscosity term was considered, but it was calculated individually, and only the pressure term was calculated using Riemann solver. The erasing of physical discontinuities such as water breaking surface still exists during the simulation.
In order to improve the stability of SPH, Phillips and Monaghan (1985) added a constant stress tensor on the original one. It is a simple way to remove the instability, but it cannot be applied to all other models. Dyka and Ingel (1995) calculated the pressure at the stress points which are away from SPH nodes. Randles and Libersky (2000) expanded this method to multidimensional space. However, the tensile instability cannot be removed completely and some compressive instability is brought. From the stability analysis of the wave simulation with SPH, Morris (1994Morris ( , 1995 found that the stability properties of SPH are improved by using the kernels whose first derivatives go to zero with decreasing inter particle spacing more rapidly. In other words, using higher order spline kernel function can improve the stability of SPH. Chen et al. (1999) proposed corrective smoothed particle method (CSPM) through Taylor expansion. It increased the derivative accuracy of SPH to improve the stability. But serious splash of particles is observed at fluid free surface during the computation. Monaghan (2000), Gray et al. (2001) added an artificial stress in the momentum equation in order to repel the particles when a clustering of particles occurred in gases and MHD. Zhang et al. (2010) talked about the relationship between smooth length and the stability of SPH, and gave a suggestion about how to select an appropriate smooth length to maintain the error system of SPH stable. Guilcher et al. (2007) found wave height decaying and wave length stretch during wave propagation, and combined Godunov scheme and CSPM to simulate the wave propagation. The result demonstrated that this method can simulate the wave propagation well. However, because of the complexity of this method, the viscosity term is not considered. In addition, due to use of the HLLC approximate Riemann solver in Godunov scheme, some physical discontinuities are erased. In order to reduce this effect, smaller particle space is required (the space is proposed as 0.002 m in the article). Meanwhile for a certain particle space, the stable computation time is limited, because the correction in Riemann solver accumulated to a certain degree will affect the stability of SPH. Yang and Liu (2012) proposed an enhanced cubical spline kernel which improved a lot the ordinary SPH stability in hydrostatic problems, but it was not good enough in wave propagation simulation. Like the method used by Guilcher et al. (2007), Gao (2011) simulated wave propagation using SPH combined with CSPM and another rapid non-conservative Riemann solver, and the result was also good. The viscosity term was considered, but it was calculated individually, and only the pressure term was calculated using Riemann solver. The erasing of physical discontinuities such as water breaking surface still exists during the simulation.
In the present paper, firstly the relevant theory of SPH and the discrete governing equations are introduced. Based on the ordinary SPH method with two kinds of boundary conditions, the wave propagation in a wave tank is simulated. The cause of the instability of ordinary SPH method and the reason of wave propagation failure are investigated in detail. Refer to the studies of Morris (1994Morris ( , 1995 and Zhang et al. (2010), some improvements are conducted to reduce the instability of ordinary SPH Method. A sextic kernel function is proposed first, and then the criterion of smooth length selection is given, because using a proper smooth length calculated from the criterion can also improve the stability of ordinary SPH method. At last, based on the dynamic and repulsive boundary conditions, a newly enhanced dynamic boundary condition is presented. Compared with the method that combined CSPM and Riemann solver, this improved SPH method can calculate the viscosity term exactly, and the physical characteristics at any location can be truly described. And then, the same example using the improved SPH method with the improvements is conducted. The results show that the stability of SPH is improved a lot and the wave propagation can be simulated accurately.

Basic theory of SPH
In SPH, the function at a particle can be obtained by a summation interpolate of the function at the particles in the support area, where i is the number of the particle where the function is calculated; j is the number of the support particle of particle i; n is the total number of support particles of particle i; f i is the function to be calculated; f j is the function at support particles; m j and ρ j are the quality and density of particle j, respectively; W ij (h) is the smooth function or kernel func-tion; and h is the smooth length. The gradient of the functions at a particle can be obtained as: ∇ ∇ i where is the Laplacian operator, and means to calculate the gradient of kernel with respect to particle i.
The equations governing the motion of fluid particles are Navier-Stokes (N-S) equations, which can be written in the form of SPH as: (3) and Eq. (4) are the continuity equation and momentum equation, respectively. is the density of particle, p is the particle pressure, is the particle position, r ij is the distance between two particles, F is the external field force, and refers to the viscous force. Here the viscous force can be calculated by, α in which the top right denotes the coordinate directions. τ ε μ For Newtonian fluids, the viscous shear stress should be proportional to the shear strain rate denoted by through the dynamic viscosity , where in two dimensions α β in which the top right and denote the different coordinate directions.
This approach of the viscous force is straightforward and can model the variable viscosity for different fluids.
There still needs another equation to calculate the pressure to let the N-S equations have solutions. Monaghan (1994) applied the following equation of the state for water to model free surface flows, where is a constant, and to simulate an incompressible fluid is usually taken 7. is the reference density. B is a problem dependent parameter, which sets a limit for the maximum change of the density.
Variable time step is applied in this paper. Time step control is dependent on the CFL condition, the forcing ∆t terms, the smooth length, the sound speed, and the viscous diffusion term (Monaghan, 1989). A variable time step is calculated according to the following formulation (Monaghan and Kos, 1999), where is a constant dependent on concrete situations. is the force per unit mass, is the current sound speed of water particle i and it follows the expression where c 0 is the reference sound speed. Symplectic time integration algorithm (Leimkuhler, 1997) is chosen as the time scheme in this paper. Firstly, the density, position and velocity are calculated at the middle of the time step as: p n+ 1 2 i where the superscript n denotes the time step, pressure is calculated by using the equation of state. At the second stage, these variables are calculated at the end of the time step by (12) Fig. 1. Schematic diagram of repulsive boundary condition.

Simulation of wave propagation with ordinary SPH method and the analysis of its instability
3.1 Kernel function and boundary conditions used in ordinary SPH method The performance of an SPH model is critically dependent on the choice of the kernel functions. Quintic kernel function (Wendland, 1995) is frequently used in ordinary SPH method since it is found that the tensile correction (Monaghan, 2000) is automatically activated. The function is written as: is the dimensional parameter as in 2D, and r is the distance between two particles. The first derivative of quintic kernel function is In most literatures, the smooth length of ordinary SPH method is usually taken as . Two boundary conditions are commonly used in ordinary SPH method: one is the repulsive boundary condition (RBC) (Monaghan, 1994;Rogers and Dalrymple, 2008), the other is the dynamic boundary condition (DBC) (Crespo et al., 2007;Dalrymple and Knio, 2001). RBC can ensure that a water particle can never cross a solid boundary. Analogous to inter-molecular forces, the particles that constitute the boundary exert normal central forces on the neighbor fluid particles in support area as shown in Fig. 1.
The force f experienced by a water particle acting normal to the wall is given by where n is the outward unit normal vector of the solid wall. The distance is the perpendicular distance of the particle from the wall, while is the projection of interpolation location joining the two adjacent boundary particles and is the velocity of the water particle projected onto the normal.
Here the terms in Eq. (15) are calculated by ∆b where is the distance between two adjacent boundary particles, z is the ordinate value of the boundary particle and the origin of the coordinate is at the bottom of water, and d is the water depth.
In ordinary SPH method with another boundary condition, the DBC, boundary particles are forced to satisfy the same equations as fluid particles such as the momentum equation (4), the continuity equation (3), and the state equation (8). However, they do not move according to Eq. (11) and Eq. (12) and remain fixed in the position (fixed boundaries) or move according to some externally imposed functions (moving objects like wavemakers etc.). Thus the boundary particles can naturally generate a repulsive force on the neighbor water particles, and no need to be set manually as the repulsive boundary condition. This boundary condition is composed of two layers, the particles are organized in a staggered manner as shown in Fig. 2, and the distance between two layers is half of the initial particle space. 3.2 Wave propagation simulation with ordinary SPH method As examples, the simulations of a regular wave propagation in a wave tank by using ordinary SPH method with both RBC and DBC mentioned above are presented in this section. The depth of the numerical wave tank is 0.4 m, and the wave height, wave period, wave number, wave length,   and the Ursell number are 0.12 m, 1.2 s, 3.25 m -1 , 1.94 m, and 7.03, respectively. The numerical wave tank is shown in Fig. 3. The length of the wave propagation area in this example is 10 times the wave length. At the left boundary, a piston wave paddle is arranged to generate waves. At the end of the tank, an absorbing region with one wave length is set to avoid the reflection of waves. The time duration of the simulation is 24 s, which is 20 times the wave period. The time cost for the simulation is 12 h 14 min on the platform κ with one core of Intel CPU I7 3770k. The water surface and pressure of the particles whose position is 0.1 m below the still water level obtained from the last step of simulation are presented. In this simulation, the initial particle space is 0.02 m, and the constant in Eq. (9) is 0.1 to ensure less variation of density and time step during computation.
The water surface at the last step in the model with RBC is shown in Fig. 4 and it is compared with the theoretical waves calculated by the second-order Stokes wave theory. The horizontal axis is the relative position , and the vertical axis is the ratio of wave elevation and wave amplitude .
Obviously, from Fig. 4, we can see that the wave height decays and the wave length elongates along the flume. Also, it can be seen that the average water level at the end of propagation area is lower than zero. The variation of the particle dynamic pressure at the position 0.01 m below the still water level along the flume at the last time step is shown in Fig. 5. Also the results calculated by the second-order Stokes wave theory are given for comparison. The figure shows that similar to the wave height, the particle pressure decays along the wave tank. The oscillation of the pressure around the theoretical result is serious, especially at the end of the propagation area, the wave pressure fluctuation almost disappears.
Similarly, the comparisons of the water surface and particle pressure simulated with DBC and theoretical results are shown in Fig. 6 and Fig. 7, respectively. It can be seen that although the wave simulated with DBC is a little better than that with RBC, the wave height decay is still serious and the wave length elongates. Also, the particle pressure is unsatisfactory.

Instability analysis in wave propagation simulation with ordinary SPH method
To resolve the problem in the wave simulation with ordinary SPH method as shown above, the reason causing the problem should be cleared. Analyzing the phenomena presented above, we can find that one reason may be the instability of ordinary SPH. The source of this instability with ordinary SPH method is the existence of negative particle pressure. For a particle whose density is smaller than the reference density, the pressure calculated from the state equation (8) is negative. At the beginning of the calculation, all particles density are larger than or equal to the reference density. However, when the particles depart from each other during the computation, the term in continuity equation (3) is positive. From Eq. (14), we can see the kernel derivative term is negative, so the change rate of density is negative now. At the next time step, the density of particle decreases. As long as the density is reduced to be smaller than the reference density, the pressure of this particle keeps negative. In contrast, as the particles approach each other, their pressure increases. So the particles, clustering reduces the distance between particles, and finally causes the water level dropping mentioned in previous examples. Typical example of the clustering of particles is shown in Fig. 8. Except the clustering of particles, negative particle pressure also impedes the wave propagation. The state of particles at a time during the wave propagation simulation is presented in Fig. 9. The computational domain is divided into several parts according to the water surface fluctuation. Water particles in Region A shown in Fig. 9 approach with each other due to the water surface drops. Thus the pressure of these particles increases, meanwhile the pressure of particles in Region B decreases because of the water surface rising. An extra force generated by the pressure gradient from the wave trough to the wave peak is formed, which is noted with the solid arrow in Fig. 9. This extra force is the source which elongates the wave length, and it also induces the wave amplitude in some degrees. Compressed particles in Region A and stretched particles in Region B generate other extra upward and downward forces respectively shown as the dotted arrow in Fig. 9. These extra forces are other sources to make the wave height decay during propagation. At the same time, extra increased pressure and decreased pressure make the original pressure under water surface changed. The particles compression and stretch degree under water surface vary with different boundary conditions applied in the model. Hence, a well designed boundary condition will be helpful for wave propagation simulation. These forces always come with the wave simulation by use of SPH method. In other words, they cannot be eliminated completely. However, the extra forces can be reduced dramatically when the properly improved method is used, which makes the result of the wave propagation simulation with SPH acceptable.

Development of stability improved SPH method (SIS-PH)
4.1 SISPH method As mentioned above, the reason causing the instability of the ordinary SPH method is the negative particle pressure during simulation. Morris (1995) studied this instability of ordinary SPH method using dispersion relation for linearity in both one and two dimensions systematically, and pointed that the instability cannot be removed under any circumstances. But he also stated that the stability properties of SPH are improved by using the kernels whose first derivatives go to zero with decreasing inter particle spacing more rapidly. To describe it informally, the stability properties of SPH method with higher order spline kernel function are better than the one with lower order, and the higher the order is, the better the stability will be. Zhang et al. (2010) studied the stability of SPH with small perturbation and linearization method. In their works, they studied the relationship between smooth length and the stability of SPH method. Appropriate smooth lengths for different smoothing functions are derived in the literature which can improve the stability of error system of SPH method. Monaghan (1994) and Morris et al. (1997) pointed out that if the density variation relative degree is smaller than 1%, the simulation can be thought stable, but this effect on the stability of SPH method is very small and in most researches, this criteria can be satisfied because under this criterion the water can be considered incompressible. So combined with the instability analysis presented above, the improvements of ordinary SPH method can be processed in these three aspects: a higher order spline kernel function, an appropriate smooth length, and a proper boundary condition. These will be introduced seriatim.
Quintic kernel function is applied in the wave propagation simulation with ordinary SPH method. It causes several problems as discussed above. Therefore, higher order kernel function should be used. As referring to the quintic kernel function form, a sextic spline kernel function is proposed in this paper. The formulation is where the dimensional parameter in 2D, and its first derivative function is The calculation method of the smooth length used before has no benefit on improving the stability of ordinary SPH method. Here, with the suggestion given by Zhang et al. (2010), the most appropriate smooth length can be calculated from This means that the most appropriate smooth length can be determined with the fact that the second derivative of the kernel function at is zero. For the sextic kernel function, the most appropriate smooth length . As mentioned above, the performance of wave propagation simulation with DBC is a little better than that with RBC, but is not yet good enough. Based on the DBC, a new DBC is proposed in this paper. The boundary can be treated as follows. The thickness of boundary particles is set as two times the smooth length, and thus the water particles near the boundary can acquire full support from boundary particles. These boundary particles are arranged neatly as shown in Fig. 10.
The same with DBC, boundary particles here are forced to satisfy the same equations as fluid particles and do not move either. But it is different from DBC, and the boundary particles only exert force on water particles at the normal direction of solid wall as RBC. This ensures that the water particle near the boundary can move along the tangential direction of the wall freely. As shown in Fig. 10, the neighbor boundary particles of water particle i exert the normal force on particle i. This new DBC combines the advantages of both DBC and RBC, therefore it is called the enhanced DBC (EDBC).
After applying these improvements in ordinary SPH method, it can reduce the SPH instability greatly. Thus the new SPH method is called stability improved SPH (SISPH) method.

Simulation of wave propagation with SISPH method
In order to verify the performance of SISPH method, the regular wave propagation simulated with ordinary SPH method above is computed with SISPH method in parallel.
In this case, the initial particle space is 0.02 m, the same as the one with ordinary SPH method, but it costs 36 hours 11 minutes which is almost 3 times the one with ordinary SPH method on the same platform. Similarly, the water surface at the last time step of the simulation is presented in Fig. 11, and it is also compared with the second-order Stokes wave theoretical results. We can see from Fig. 11 that the water surface generated by SISPH method agrees with the theoretical result very well, and the wave height does not decay and the wave length does not elongate along the flume. There is also no dropping of the water level at the end of the propagation area. This indicates that the space of particles does not decrease.
The dynamic pressure at the position 0.01 m below the still water level is also obtained from the simulation with SI-SPH, and the comparison between it and the second-order Stokes wave theoretical result is shown in Fig. 12.
The result calculated with SISPH also agrees with the theoretical results quit well, meanwhile the oscillation of pressure is nearly eliminated.
Moreover, the particle distribution at a time is also given in Fig. 13, and it shows that the clustering does not exist   anymore as expected.
From the result of this example, we can see that the stability in the simulation with SIPSH is improved significantly. The wave propagation is simulated accurately and both the water surface and wave dynamic pressure agree well with the theoretical results.
In order to further verify the affectivity of the SISPH method for other situations and to determine the scope of its application, three other cases with different wave parameters are also calculated with SISPH. As the limitation of wave generation with a piston wave paddle for strong nonlinear waves, the calculated results are compared with the second-order Stokes waves. So the case of waves with the largest wave number but no wave breaking, with the strongest nonlinearity but in the scope of the second-order Stokes theory (the Ursell number smaller than ), and with nonlinear strength between the previous two cases are chosen. The selected wave parameters are listed in Table 1 Similarly to the first case, Fig. 14-Fig. 16 give the comparison of the wave surface and the wave pressure between the results at the end step of the simulation and the corresponding theoretical results, respectively. From the figures, we can see that the waves simulated with SISPH agree very well with the theoretical ones for the three cases. So considering the previous case, four cases with SISPH are verified in this paper. It can be concluded that the propagation of all the waves in the scope of the second-order Stokes theory without breaking can be simulated very well with SISPH.
SISPH improves the stability of SPH from the basic SPH theory in a simple way, and no extra special processing technology which would bring some other new problems is applied in this method. So, the simulation with SISPH can reveal the physical reality, which is the most attractive property of SISPH compared with other improved method.
However, simulation with SISPH has some other disadvantages due to the long smooth length which is much longer than the one in ordinary SPH method. From the com-

180
CHANG Jiang et al. China Ocean Eng., 2017, Vol. 31, No. 2, P. 173-181 κ putation time comparison between SISPH and ordinary SPH method, we can see that the direct disadvantage is the high cost of computation as the number of contributing neighbor particles increases. In addition, according to Eq. (9), the time step becomes longer for a long smooth length. This means that the time accuracy is reduced with SISPH. Fortunately this can be offset by setting a smaller constant in Eq. (9), which also increases computation time of the simulation. Besides, using EDBC requires thicker boundary layer with more particles, and the computation consumption on these particles is increased. In a word, much more cost of computation is needed in the simulation with SISPH.

Conclusions and discussion
The reason of the instability in the wave propagation simulation using ordinary SPH method can be attributed to the negative particle stress which is caused by the fact that the particle density is smaller than the reference density. Because of the extra force generated in the simulation of wave propagation, the wave amplitude decays and the wave length elongates along the path. The appearance of extra force during wave propagation is the result of the alternation of pressure increasing and decreasing along the wave propagation. During the simulation of wave propagation using ordinary SPH method with RBC and DBC, similar problems occur although the results with DBC is a little better than that with RBC. Also, the water level drops during wave propagation with these two boundary conditions because of the particles clustering. The pressure under water surface is not satisfactorily simulated in the model with ordinary SPH method.
Applying the sextic kernel function constructed in this paper with the most appropriate smooth length in SPH is the key to improve the stability of SPH method. On the other hand, a proper boundary condition is very helpful to simulate wave propagation. This new SPH method with enhanced DBC (EDBC) is called SISPH method. It can reduce the instability significantly with no need of smaller particle spacing and can maintain stable for a very long time duration. The simulated example using SISPH shows that the particle clustering and oscillation of the particle pressure is nearly eliminated, and the water surface and the pressure agree with the theoretical results quite well. Another advantage of SISPH is that the process is simpler and can truly reveal the physical reality without bringing some new problems compared with other methods.
But the large cost of the computation limits the application of SISPH, and thus applying parallel algorithm and GPU acceleration technology in SISPH is necessary for a large-scale simulation. This is the future work about the SI-SPH.