A corrected solid boundary treatment method for Smoothed Particle Hydrodynamics

Smoothed Particle Hydrodynamics method (SPH) has a good adaptability for simulating of free surface flow problems. However, there are some shortcomings of SPH which are still in open discussion. This paper presents a corrected solid boundary handling method for weakly compressible SPH. This improved method is very helpful for numerical stability and pressure distribution. Compared with other solid boundary handling methods, this corrected method is simpler for virtual ghost particle interpolation and the ghost particle evaluation relationship is clearer. Several numerical tests are given, like dam breaking, solitary wave impact and sloshing tank waves. The results show that the corrected solid boundary processing method can recover the spurious oscillations of pressure distribution when simulating the problems with complex geometry boundary.


Introduction
The Smoothed Particle Hydrodynamics (SPH) method is a meshless, purely Lgrangian technique which was originally developed by Lucy (1977), and Gingold and Monaghan (1977). SPH was successfully used to simulate complex problems, like fluid mechanics (Monaghan, 1994;Zheng et al., 2014), explosion mechanics , fluid structure interaction (Liang et al., 2010), and so on. Although there are some shortcomings of SPH which are still in open discussion, it is a robust and vigorous technique for violent free surface flow simulation.
In SPH, the particles are scattered by moving nodes and carrying field variables such as pressure, density and velocity. The smoothing kernels are used to approximate a continuous field. Boundary condition is one of the key aspects for numerical simulation and special attention should be paid to complex boundary, so the boundary treatment method plays an important role in SPH. Many researchers have investigated the SPH boundary processing method and different researchers may use different handling methods, like Monaghan (1994) used a group of boundary virtual particles which can add force to fluid particles near the boundary to prevent fluid particles moving through boundary. Shao and Lo (2003) adopted the ghost particle to overcome the truncation by solid boundary. Lee et al. (2008) introduced the uniform distribution ghost particle to handling the solid boundary. The generalized boundary processing scheme sets up virtual particles along the geometry boundary and the physical properties of virtual particles are assigned by the fluid particles (Randles and Libersky, 1996).
Generally, solid boundary treatment for SPH follows two basic concepts. The first one is to fill the wall with boundary particles, in order to ensure that the support of kernel function near a wall is completely covered with the neighbor particles, e.g. Morris et al. (1997) fixed boundary particles to model the curved surfaces and treat them as real particles. Density and pressure of boundary particles are evaluated as time stepping, and they are considered in the continuity equation for fluid phase. Consequently, the pressure field will increase or decrease when particles move towards or away from the wall in order to prevent penetration. A more appropriate ghost-fluid method has been used to enforce the boundary condition on the body surface (Marrone et al., 2011). Specifically, the solid domain is modeled through a set of virtual ghost particles and all the fluid fields are extended on these virtual particles through appropriate interpolation techniques. In order to enhance the boundary treatment, Bouscasse et al. (2013) discretized the solid surface by using a set of equispaced body nodes and a layer of ghost particles is also disposed in the solid region as before. On the second step, either non-vanishing surface integral when smoothing the flow quantities close to the boundary is accounted for, or artificial repulsion forces are introduced to prevent the particles cross the interface. It is just similar to Liu and Liu (2010) who introduced a Lennard-Jones-like potential between fluid and wall particles to add a repulsion force normal to the boundary.
In this work we present a corrected solid boundary treatment method, and the main process is similar to Marrone et al. (2011) but some details are simpler. With the help of the wall treatment method, it can easily handle arbitrarily geometries. Firstly, it gives the process of this corrected boundary treatment method. Then a classic dam-breaking simulation is presented to validate the accuracy and efficiency of the wall boundary handling. Finally, solitary wave impact and violent sloshing wave are given. With the help of the corrected boundary treatment method, this corrected SPH method can improve the pressure distribution and recover the spurious oscillation on pressure time history effectively.

δ-SPH
In this section, the standard SPH method is briefly introduced. More details were provided by Monaghan (1994) or . In the SPH method, the Lagrangian form of the Navier-Stokes equations are written as follows: where ρ is the fluid density, u is the particle velocity, t is time, p is the particle pressure, f is the body, and ν is the laminar kinematic viscosity. It should be noted that Eq. (1) is written in the form of a weakly compressible flow. In the traditional SPH method, the pressure of each particle can be obtained by an explicit relationship, which can be shown as: ½ 0 where C 0 is the artificial sound speed, which can be taken as 10 times the maximum velocity, is the constant water density, and γ=7 for liquid water problem simulation.
The δ-SPH method is one of the improved SPH schemes, which is based on the assumption of the weak compressibility. With the help of density correction, this method can effectively improve the stability of pressure calculation (Marrone et al., 2011). The particle summation forms of the δ-SPH equation are written as follows: where, h is the smooth length, r ji =r j -r i , W ij is the kernel function, f i is the value of f, associated with the generic particle i, V j is the volume of particle j, is the renormalized density gradient defined as: where δ and α are the coefficients of density and velocity diffusion, respectively, α=0.01 and δ=0.1 in all simulation tests.
2.2 Solid boundary handling Generally, ghost particle is adopted to describe the wall boundary and strengthen the boundary conditions. The ghost particles are produced from fluid particles according to the mirror relationship. Consequently, the physical properties of ghost particles, like position, velocity and density, are based on fluid particles. It is clear that the positions of ghost particles are not in uniform distribution when the wall boundary becomes complex. It should introduce some special treatments. When the solid boundary is an arc in the 2D case, the velocity of fluid particle should be transformed into local coordinate firstly, and then the velocity of the ghost particle is transformed into the global coordinate system after the mirror operation. Fig. 1 is the sketch of the arc wall mirror operation. It is obvious that ghost particles are not in uniform distribution, the density of the ghost particle is smaller than that of the corresponding fluid particle in Fig. 1.
The irregular distribution of ghost particle may cause the penetration and nonphysical oscillations of pressure distri- CHEN Yun-sai et al. China Ocean Eng., 2017, Vol. 31, No. 2, P. 238-247 239 bution, and some new boundary treatments come out to recover solid boundary accuracy. The boundary handling method developed by Marrone et al. (2011) can depict the wall boundary accurately, while it is difficult to obtain the properties of the ghost particles which do not have corresponding mirror particles in the fluid domain when the boundary shape becomes complex. A set of equispaced body nodes are used to discretize the solid surface and enhance the boundary condition (Bouscasse et al., 2013). It is clear that the volume of the body nodes may be overlapped with the first layer ghost particles in the boundary domain. The accumulation of these particles will increase the density near the solid wall. In order to overcome these shortcomings, this paper presents a corrected boundary handling method. The scheme of this paper can be divided into four steps: (1) arranging the fixed virtual ghost particles along the boundary side; (2) setting the mirror particles for the corresponding fixed virtual particles; (3) assigning the physical properties of the fixed virtual mirror particles by the simplified high order interpolation method; (4) evaluating the virtual particles according to the mirror relationship.
The first step, third step and fourth step are similar to Randles and Libersky (1996), but the second step is different which is shown in Fig. 2. It should be noted that only the fluid particles are included when using high order interpolation method, such as the moving least square (MLS) method, to assign the physical properties for the fixed virtual ghost particles. The effective particles for this interpolation are indicated by the blue circle in Fig. 3. For example, the pressure of a normal ghost particle P G is evaluated as follows: where, d is the distance between ghost particle and solid boundary, g is the gravity acceleration, n is the normal unit vector to the solid boundary, and W is the kernel function. Marrone et al. (2011) used the moving least square (MLS) interpolation as the kernel function. Different from the former option, the simplified finite different interpolation (SFDI) is adopted in this paper. As shown in Zheng (2011), SFDI can obtain similar interpolation precision as the MLS interpolation, but the time consumption is much smaller than that of MLS. This is a big improvement of this corrected solid boundary method.
Another important processing of this solid boundary is to assign the physical properties of ghost particles (hereinafter denoted as "special ghost particles") which do not have the corresponding mirror particles in the fluid domain. Fig. 4 shows some special ghost particles (particles above the pink line) near the free surface in a sink boundary. In previous researches, such as Marrone et al. (2011) and Bouscasse et al. (2013), the properties of these special ghost particles are determined by their nearest ghost particles. It means that the properties of boundary ghost particles should be determined before assigning the physical properties of special ghost particles. In order to ensure the special ghost particles getting right physical properties, it is difficult to satisfy the precision of evaluation when the boundary shapes become more complex. In order to solve this problem, the present method does not consider the mirror process of those special ghost particles in the fluid domain. As a key problem, the pressure of special ghost particle is obtained through SF-DI interpolation of its near fluid particles. The density of special ghost particle could be obtained by inversing the state equation Eq. (3). The velocity of ghost particles is set up according to the wall boundary. The details of assignment procedure for these special particles are shown in Fig. 5.
In order to show the improvement of boundary handling method clearly, Table 1 shows some differences among diff-   erent improved methods mentioned in Marrone et al. (2011) and Bouscasse et al. (2013). The new boundary handling method combined with δ-SPH is named as SPH-NB. SPH-TB is for the traditional ghost boundary handling method combined with δ-SPH.
Lastly, in terms of the rules of ghost particle arrangement, many papers such as Monaghan (1994) and  had shown how to generate mirror particles according the regular ghost. Marrone et al. (2011) also described how to generate mirror particles at some special conditions in detail. This paper follows the method of mirror particle generation of Marrone et al. (2011), but the interpolation method of ghost particles and assignment relation of some special boundary particles are different.

Particle arrangement at initial time
When starting the SPH computation, the initial fluid particles are not uniformly arranged near the complex boundary. It should adopt a moving action to set the fluid domain to equilibrium state. A modified SPH governing equation was used to solve the problem (Colagrossi et al., 2012). This particle arrangement method is important for the pressure distribution at the initial time. As this paper is mainly focused on the solid boundary improvement, this part is just simply introduced.
The first step is to close the domain boundaries. As a consequence, this implies that the free surface has to be treated as a solid boundary. The domain boundary is modeled through fixed solid particles. The second step assumes the density, the pressure and the volume constant all over the fluid domain. We indicate them through ρ 0 , p 0 and V 0 , respectively. Then the new governing equation Eq. (7) is used to compute the fluid particle. The last step is to use the position of particles which are in the equilibrium state to ini-tialize the SPH simulation.
where β=2p 0 /ρ 0 and . is just used to ensure the convergence of the computing scheme, and we can choose a liner damping term: where d is the coefficient of spatial dimension, α 0 is a free dimensionless parameter ranging from 1.0×10 -3 to 5.0×10 -3 . Fig. 6 shows the fluid particles near the arc sharp wall and corner sharp wall. After the process of initial particle arrangement, the fluid particles fit the boundaries perfectly.

Dam breaking
Dam breaking is a benchmark for violent free surface flow. Although it is not new for SPH to simulate this problem, it can validate the accuracy and efficiency of the wall boundary handling method. In this example, a rectangular column of water is confined between two vertical walls. The width of the water column is L and the height of this water column is H. At the beginning of computation, the dam is instantaneously removed and the water is allowed to flow out along a dry horizontal bed. L w is the distance between the two vertical walls. There are three pressure sensors P 1 , P 2 , and P 3 on the right vertical wall. The model of dam breaking is given in Fig. 7.

Pressure distribution
During the simulation H=0.6 m, L/H=2.0, L w =5.366H, the particle spacing dx=4.0×10 -3 m, the total particle number is 45000 (150×300) and the space size is 1.6×10 -5 . Time step length dt=1.0×10 -4 s. Three layers particles are distributed along the wall boundary. The pressure is evaluated as follows: where j is particle number of the effect fluid particle and N is the total number of the neighbor particle number. In order to improve the calculation accuracy, the initial particle pressure is given according to ISPH results. As t=0 baffle opening, the right side of this water column becomes free surface. Fig. 7 includes the initial pressure distribution of the water column. The initial pressure is not the same as the hydrostatic pressure distribution. Fig. 8 gives snapshots of the pressure distribution at different time for dam breaking simulation by SPH-NB. It can be seen in Fig. 8 that the pressure field by SPH-NB is smooth, even during the process of wave impact. Fig. 9 shows the comparison of the pressure history on P 1 and P 2 with experimental data (Buchner, 2002), SPH-NB and SPH-TB. According to the comparison of Fig. 9a, the time of pressure raise for dam breaking flow of two numerical results is later than that of experiment data. Pressure time history of Point P 2 is given by Fig. 9b, and the numerical results have a good agreement with the experimental data. The results of SPH-NB are more smooth and closer compared with that of SPH-TB. The results of SPH-TB have some slight oscillations in the pressure time history.

Solitary wave simulation
Solitary wave is one kind of extreme wave phenomena in complicated ocean environment. During the earthquake, the sea wave generated by the intense crustal movement will gradually become a solitary wave as the time goes on. It needs more attention to study the interaction between structure and wave, many researchers have simulated that problems e.g. Maxworthy (1976) did the experiment about the collision between two solitary waves, and Su and Mirie (1980) deduced an analytical solution for practical collision

242
CHEN Yun-sai et al. China Ocean Eng., 2017, Vol. 31, No. 2, P. 238-247 of two solitary waves. Fig. 10 gives the model of solitary wave simulation. There are two pressure sensors P 0 and P 1 on the right wall. The angle of the right side wall can be changed according to different cases. Fig. 10 shows the case when θ=120°. There are two wave gauges fixed which are 2 m and 8 m away from the left side of this tank, respectively. The length of the tank is 10 m and the depth of water is 0.25 m.
In this part, wave height λ=0.15 m, particle spacing dx=0.01 m, fluid particle number is 25000, stepping length dt=0.0001 s. A piston wave maker is set up on the left side of the tank, the motion of this plate can refer to Ma and Zhou (2009): , , and . The analytical solution for the wave profile can be derived from the equation where η is the water surface elevation, λ is the wave amplitude, d is the water depth and is the solitary wave velocity. The comparison of solitary wave profiles at different time is shown in Fig. 11. According to the results, the numerical wave elevations can obtain a good agreement with analytical results.

Solitary wave impact
Generally, the traditional ghost particle method can deal with simple shape boundary very well, but for some complex geometry, the traditional ghost particle method may meet some difficulties and the simulation may crash. In order to show the advantage of this corrected boundary method, this section has shown some comparison results of SPH-NB and SPH-TB for solidary wave impact simulation.
The model of solitary wave impact is shown in Fig. 11. The only difference is the slope angle at the end of the tank. The slope angle can be chosen 60° and 120°, respectively in this section. The solitary wave is generated the same as in Section 3.2.1. The particle size and stepping length are not changed. As the slope angle is not 90°, the total fluid particle number is not equal to 25000, the particle numbers are 24805 and 25166 for the cases θ=60° and 120° respectively, and the particle size is the same as that in Section 3.2.1.
In order to illustrate the advantages of SPH-NB clearly, Fig. 9. Comparison of pressure time history on P 1 (a) and P 2 (b) with experimental, SPH-NB and SPH-TB results.   Fig. 12 gives some snapshots of the pressure distribution of different slopes at t=6.0 s. Fig. 13 gives the comparison of experimental results and numerical results by SPH-NB at the end of the tank when t=6.0 s and θ=120°.The work of experiment for solitary wave impact is finished in HEU and the physical model is set up following the numerical test. Fig. 14 shows the comparison of the velocity and pressure distribution by SPH-NB and SPH-TB near the slope when t=3.0 s and θ=120°. From the comparison, the distributions of pressure and velocity of SPH-NB are more smooth and reliable than those of SPH-TB. In order to justify the improvement of SPH-NB, Fig. 15 gives the comparison of the pressure time history at different sensors, which can refer to Fig. 11. According to the results of Fig. 15, the results of SPH-NB are more smoothed. The pressure time history can get closer with experimental data. But it has obviously spurious oscillations for the results of SPH-TB.    Fig. 15. Comparison of the numerical pressure at P 0 and P 1 by experimental results, SPH-TB and SPH-NB.

Sloshing tank waves simulation
In order to show the advantage of this corrected boundary treatment method, it gives some results of more complex geometries. Firstly, it gives the results of small amplitude sloshing. It is very helpful to the verification of this corrected boundary treatment method for small amplitude wave simulation. Then this part considers some results of violent sloshing tank waves, even for complex geometry.

Small amplitude sloshing
The geometry of the sloshing tank is illustrated in Fig.  16. The length of the rectangular tank is L, its height is h, and the water depth is d. There is a pressure sensor fixed on the left wall. The height between the sensor and bottom is h 1 . The sway displacement of the tank is given by X s =a 0 [1-cos(Ωt)], where a 0 and Ω are the amplitude and the frequency of this motion, respectively.
When the motion amplitude of the tank is small, like , the frequency is not close to the natural frequency of the tank, and the free surface elevation can be analytically expressed by . More details about the solution can be found in Faltinsen (1976). In this simulation, d=1.0 m, L=2d, and a 0 =0.002d are selected.
The particle size dx=0.02 m, fluid particle number is 5000, and stepping time length dt=0.0001 s. Fig. 17 gives the comparison of the free surface elevations on the left tank wall. The excited frequency is Ω/ω 1 =0.8 and Ω/ω 1 =0.9. The results of different frequencies can obtain a good agreement with analytical results.

Violent sloshing
In order to show some results of violent sloshing flow, this section chooses the cases of Kishev et al. (2006). The geometry of this tank is a rectangle with L=0.6 m, h=0.5L, d= 0.4h. The parameters for this case are taken as a 0 =0.05 m and T 0 =1.5 s, which are the same as those in Kishev et al. (2006). Two pressure sensors fixed at point h 1 /L=0.1667 and h 2 /L= 0.0833 on the left wall. As the experimental pressure time history given in Kishev et al. (2006) did not show the transient period, it is just to do this comparison when it reaches the first pressure raise and the tank top impact occurs.
In order to test the convergence of this simulation, it does the comparison by using different particle numbers when h 1 /L=0.1667, which is shown in Fig. 18. From the comparison, the results of SPH-NB method can obtain a good convergence. The pressure time history can become smoother as the particle number increasing. Fig. 19 gives the comparison of the pressure time histories at h 2 /L=0.0833 on the left wall with the experiment, ISPH_R (Zheng et al.,   2014) and SPH-NB. According to the comparison of Fig.  19, the results of SPH-NB can obtain a good agreement with the experimental results and ISPH_R. It can give some helpful information, even weakly compressible SPH can get smooth and reliable pressure time history with the help of the improved boundary handling method. Fig. 20 shows some snapshots of the pressure distribution at different time by using SPH-NB.

Violent sloshing in a complex shape tank
In order to show the properties of the corrected boundary treatment method, this section gives the results of violent sloshing in a complex shape tank. The sway displacement of the tank is given by X s =a 0 [1-cos(Ωt)], which is the same as the previous simulation. The geometry of this sloshing tank and particle arrangement at the initial time are illustrated in Fig. 21. The particle size is 0.02 m and there are 7482 fluid particles in this tank.
In order to obtain complex wave profile and violent flow during this simulation, the excited amplitude is a 0 =0.1 m and T 0 =1.5 s. The snapshots of the pressure distribution of violent sloshing in this tank are shown in Fig. 22. From these results, the pressure distribution of this complex shape   tank is still smooth and reliable. As it is lack of corresponding experiment data, the verification of this case needs to do further research.

Conclusions
This paper presents a reliable boundary treatment method for weakly compressible SPH method. According to the algorithm of δ-SPH, the results of the pressure distribution by SPH-TB are more reliable, and some spurious oscillation on pressure time history can be corrected efficiently. This thought of virtual ghost particle handling method is not new completely, but some details of numerical techniques are quite different from previous processes, which is simpler and more reliable for solid boundary condition. It also gives some numerical tests, such as dam breaking, solitary wave impact and sloshing tank wave. With the help of more uniform distribution of ghost particles, especially for the cases of curve or complex shape, this corrected boundary handling method can satisfy the boundary exactly. According to the numerical results, SPH-NB can obtain more reliable results for velocity and pressure distribution. So the SPH method based on the weakly compressible assumption is still a good choice to solve violent free surface problems. The shortcoming of spurious oscillation on the pressure distribution can be the improved by improved boundary handling method and δ-SPH. SPH-NB method can also be easily extended to more complex models and practical cases.