The Modification of Boundary Treatment in the Incompressible SPH for Pressure Calculation Accuracy on the Solid Boundary

The Incompressible Smoothed Particle Hydrodynamic (ISPH) is one of the particle methods and commonly used to solve some complicated physical problems including free surface flow problems. The study regarding the boundary treatment has become an active research area in the mesh-free or particle method recently for measuring the accurate and robust pressure near the boundary. The penetrations of fluid particles may be happened if the adequate pressure boundary condition on the solid boundary cannot be satisfied. In this paper, a simple boundary treatment, which can be satisfied the non-homogenous Neumann boundary condition on the solid boundary and Dirichlet condition on the water surface, is proposed. The key point of our proposed treatment is that these boundary conditions are automatically satisfied by solving a modified pressure Poisson equation. Lastly, the effectiveness and accuracy of boundary treatment proposed are then authenticated with couples of numerical analysis and compared with the experimental tests.


Introduction
Smoothed particle hydrodynamics (SPH) is a Lagrangian meshless particle method that has been popular and actively used by researchers especially in computational fluid dynamics nowadays.This SPH method was first developed by Lucy [1] and Gingold [2] which originated from simulation of the astrophysical problems about three decades ago.SPH method is a unique particle method which does not required any grid and has some advantages over grid based method such as simple in implementation and also eases of handling even for complex fluid and larger deformations.Then, the SPH method was extended to incompressible viscous flow.Nowadays, a stabilized Incompressible Smoothed Particle Hydrodynamic (ISPH) approach is one of the particle methods and commonly used to solve some of complicated physical problems including free surface flow problems [3].In addition, boundary treatment of the particle method is one of the important issues for accurate solutions.Many researchers tried to improving the boundary treatment from the ghost particle then fixed ghost particle by Marrone [4], and virtual marker by Asai [5] and until know this research topic still become a vigorous research area.Therefore, a modified boundary treatment is proposed, which resolves the problem of pressure accuracy near the boundary including the prevention the penetration.

A Stabilized Incompressible Smoothed Particle Hydrodynamics (ISPH)
In this paper, the basic of stabilized ISPH which introduced by Asai [3] was adopted.The following sections will be described the governing equations used and the concept of Smoothed Particle Hydrodynamics (SPH) with the modification for incompressible flow.

Governing equation
The continuity equation and the Navier-Stokes equation in the Lagrange description are given as where ρ and ν are density and kinematic viscosity of fluid, u and p are the velocity vector and pressure of fluid respectively.F is an external force, and t indicates time.The turbulence stress τ is necessary to represent the effects of the turbulence with coarse spatial grids.In the most general incompressible flow approach, the density is assumed by a constant value with its initial value.Then the above mentioned governing equations lead to

Stabilization of pressure in PPE
Generally, in SPH approach, there some difficulty occurred due to the accuracy of the density representation in the SPH formulation and the treatment of pressure Poisson equation (PPE).It is difficult for the numerical densityto keep its initial value.Therefore, here, the PPE is reconsidered to overcome the error of artificial pressure fluctuation in the ISPH.In density invariance scheme, the particle position updates after each predictor step and the particle density is updated on the intermediate particle positions.While in the correction step, the pressure term is used to update the particle velocity obtained from the intermediate step.Then, the final form of PPE in ISPH can be approximately redefined by the following equation whereby they are introducing the relaxation coefficient, α (0 ≤ α ≤ 1) and usually the small value is used around 1% or less.
3 Boundary Treatment

Slip and Non-Slip condition with virtual marker
The virtual marker proposed to be located in a position which is symmetrical to the wall particle across it solid boundary whereas the wall particles is placed on a grid like structure.Normally, the velocity and pressure on the marker are interpolated based on the concept of weighted average of neighboring particles, which is the fundamental equation of SPH as in the Equation ( 6).However, the portion of the particles within the compact support might be in the wall domain or unoccupied domain such as air phase.For that case, the interpolation approximation will used the modified weight function.

( ) ( ) ( )
The virtual marker is not directly associated with the SPH approximation but only a computational point to give the wall particle accurate physical properties.Then, the velocity on the virtual marker which is corresponding to the wall particle across its solid boundary is interpolated by the SPH approximation.Practically, it is necessary to discuss an optimized condition between slip and no-slip condition to refer the resolution of the simulation model.Therefore, the following equation is proposed which can be consider as mixed condition whereby β represents the ratio between slip and no-slip condition.Where M is a second order tensor to implement the mirroring processing, and it is symbolized by using inward normal vector of the wall and the kronecker delta and R is a mirror symmetric tensor

Pseudo Neumann boundary condition: Dirichlet conditions
Using the similar approach as introduced by Asai [5], the pressure on each virtual marker and the external force will evaluated by using the physical quantity in the particles neighborhood with the normalized weight function as in Equation (6).Thus, the pressure on the walls particles in the border processing technique does not directly evaluated by solving the modified pressure Poisson, Equation (5).By introduced solving the PPE will gave entire circumference of Dirichlet condition whereby approximation of the value of the wall particles pressure by referring the value of the virtual marker as to satisfy the Neumann boundary condition.

3.3.Neumann boundary treatment (completely satisfy the pressure Neumann conditions)
In this section, the Modified Virtual Marker with Regular grid (MVMRG) as new boundary treatment which satisfied the non-homogenous Neumann boundary condition will be described.The concept of this treatment is to give a wall particle accurate physical properties, velocity and pressure.As concern, the virtual marker used in pseudo Neumann condition is convenient to perform the analysis including the rate of real condition as slip or non-slip condition.However, by solving PPE with the giving the entire circumference Dirichlet condition causes a lower pressure value than the actual hydrostatic pressure distribution in the vicinity of wall particles.Therefore, the MVMRG is proposed in order to overcome this problem.The procedures will be described briefly.First, the wall particles placing similarly as in Pseudo Neumann condition method which placed equally spaced structural grid-like inside solid.Then, the location of the virtual markers are located on the boundary surface which now no longer in the fluid domain.The process of calculated the value of the inner speed on the virtual marker is similar with the previous section which using the SPH approximation as in Equation (6).Therefore, the modification of the coefficient matrix and the source term need to be conducted due to the flow velocity is not strictly mirror-symmetrical.Then, pressure Dirichlet conditions at the stage of solving the Poisson equation was applied only to the place given originally and the pressure approximation of wall particles is evaluated by using the following equation: By using the above techniques, the pressure solved by the PPE can be limited to the water particles.Hence, the special modification of the coefficient matrix and source terms when the particles were discretized.The relationship of the concept mentioned above between the virtual marker, water particles and the wall particles as shown in Figure 1.After solving the PPE abovementioned, the pressure needs to go through mapping process whereas the pressure and the external force will be analyzed by SPH formulation.
Hence, there are some verification and validation test was made using this new modification which is used MVMRG to see the performance and also in comparison with experimental test that will be discussed in the following sections.

Verification with hydrostatic test
The first test is the simplest hydrostatic problem which the rectangular tank with compatible surface.The comparison of the pressure distribution was made between different boundary treatments which are using pseudo Neumann condition and proposed MVMRG boundary treatment.This test was conducted to show the robustness of our modified boundary treatment.The tank was filled up with the water height of 0.2m and the particles size is 0.01m.Where, the point of pressure measurement is located in the middle of bottom tank.Figure 2(a) shows the different results on the contour pressure distribution especially the effects near the bottom tank which obtained using different treatment.It is clearly show that there are high pressure distribution using the proposed modified treatment, using MVMRG compared with the Pseudo Neumann boundary condition that give lower pressure distribution.Then, the numerical result of the hydrostatic pressure were used to compared with the theoretical of value of hydrostatic pressure as shown in Figure 2(b).Obviously, the result obtained by using modified boundary treatment of Neumann condition, MVMRG gives almost similar result with the theoretical.However, the result found from Pseudo Neumann treatment gives the lower value especially near the bottom surface due to the treatment used at the boundary surface not perfectly treat as non-homogenous of Neumann condition.Next, the same setup was used but with the additional of several geometrical object submerged into the water with size of the particles is reduced to 0.005m.Generally, Figure 3(a) shows the pressure distribution inside the water domain and also at the concave-convex shape boundary surface without any penetration occurred.Then, the numerical results of the pressure distribution were compared with the theoretical hydrostatic pressure as shown in

Validation with 3D dam break flow with the opening gate
The simulation of three dimensional dam break flow is used to validate the versatility of our proposed method and compared with the experimental results.The experimental test on the investigation of dynamic pressure loads during dam break with an opening gate that was carried out at Technical University of Madrid (UPM) as reported by Lobovsky [6] is used as the reference for simulate purposed and design the analysis model as shown in the schematic diagram (Figure 4).The clear comparison of free surface profile evolution snapshot between the numerical simulation and the experimental at certain times as shown in Figure 5 and the effect of the dam gate removal looks similar at every snapshot on the free surface shape is observed.Then, the observation of the location of the wave front also shows the same tendency with the experiment.The pressure contour is also including in this figure shows the pressure gradient from the static condition before an opening gate and after gate's removal and the impact after hit the wall which shows higher pressure.Furthermore, from the snapshot, it is clearly shows that the lowered sensors are exposed to the high pressure as the red color of contour is observed after the water smashed on the wall surface.Therefore the results obtained at the lowered sensor from the simulations are compared with the experimental data generally.It indicates the significant results, where the simulation results fall within the range of 100 tests of experimental results from sensor 1 with the simulation's result as shown in Figure 6.From the comparison, the arrival time from the numerical simulation is almost same with the arrival time in the experiment and the tendency of the pressure time history is similar.

Conclusions
The numerical simulation that been conducted using the proposed boundary treatment, MVMRG with perfectly Neumann condition shows the good tendency and more accurate pressure value near the boundary surface based on the hydrostatic tests and also the comparison with the experimental test for hydrodynamic pressure evaluation.The profile of the free surface also shows the similarity with experiment results.In addition, the comparison between both boundary condition show clear discrepancy in the pressure especially at the edge of bottom surface and this strongly proves that the proposed modified treatment gives the better treatment.

Figure 1 .
Figure 1.The relationship of the modified virtual marker on the boundary with the wall particles and water particles.

Figure 3
Figure 3(b).Based on the result obtained, the result using MVMRG gives higher value compared with the Pseudo Neumann condition.

Figure 2 .
Figure 2. (a) The pressure distribution by using different boundary condition (b) the comparison of pressure distribution with the depth by using hydrostatic pressure, perfectly Neumann boundary treatment and Pseudo Neumann boundary condition for compatible surface boundary.

Figure 3 .
Figure 3. (a) The pressure distribution on the water particles and wall boundary by using different boundary condition (b) the comparison of pressure distribution with the depth by using hydrostatic pressure, MVMRG with perfectly Neumann condition and Pseudo Neumann condition for compatible surface with complex geometrical objects.

Figure 4 .
Figure 4.The dam break simulation with an opening gate which corresponding to the experimental [6].

Figure 5 .
Figure 5.The comparison between simulation and experimental of the free surface evolution at times 0.0, 316.7, 413.4,463.3, and ±3.3ms.

Figure 6 .
Figure 6.The comparison of the non-dimensional pressure at sensor 1 between the simulation results with 100 tests of experiment data.