Extension of Smoothed Particle Hydrodynamics (SPH), Mathematical Background of Vortex Blob Method (VBM) and Moving Particle Semi-Implicit (MPS)

SPH has a reasonable mathematical background. Although VBM and MPS are similar to SPH, their mathematical backgrounds seem fragile. VBM has some problems in treating the viscous diffusion of vortices but is known as a practical method for calculating viscous flows. The mathematical background of MPS is also not sufficient. Not with standing, the numerical results seem reasonable in many cases. The problem common in both VBM and MPS is that the space derivatives necessary for calculating viscous diffusion are not estimated reasonably, although the treatment of advection is mathematically correct. This paper discusses a method to estimate the above mentioned problem of how to treat the space derivatives. The numerical results show the comparison among FDM (Finite Difference Method), SPH and MPS in detail. In some cases, there are big differences among them. An extension of SPH is also given.


Introduction
The author has once shown how to obtain the space derivative in an irregular mesh using the moving least square method [1].A similar problem is discussed from a different viewpoint in SPH.Gingold & Monaghan [2] and Lucy [3] have developed Smooth Particle Hydrodynamics method (SPH).In SPH, the continuous mass distribution is approximated by the finite number of particles.Namely, the continuous quantities are represented by the finite number of the discrete quantities.At the same time, the space derivatives of the distributed quantities are expressed algebraically by the discrete quantities.Then, they transform the continuous system into a discrete system convenient for the numerical solution of the initial and boundary value problem.
Vortex Blob Method (VBM) is one of the practical numerical tools for high Reynolds number flows [4].The abduction is approximated reasonably.However, the diffusion of vortices can't be approximated precisely.We have shown that this problem can be solved, if we apply the ideas developed by SPH.
Moving Particle Semi-implicit method (MPS) uses a similar discretization of the initial and boundary value problem as SPH and VBM.However, the mathematical background is not sufficient.Not with standing, the numerical results seem reasonable in many cases [5].This suggests us that we may find the mathematical background [1] [6].We have shown an interesting relationship between SPH and MPS as far as the gradient operator is concerned.
The numerical results show the comparison among FDM, SPH and MPS in detail.In some cases, there are big differences among them.
A classification of numerical solutions is shown in Figure 1.The weak solutions: FEM, FVM, IRM and GIRM are Finite Element Method, Finite Volume Method, Integral Representation Method and Generalized Integral Representation Method, respectively.The strong solutions: FDM, SPH, MPS, LBM and ColM are Finite Difference Method, Smooth Particle Hydrodynamics, Moving Particle Semi-Implicit Method, Voltex Blob Method, Lattice Boltzmann Method and Collocation Method, respectively.The general characteristics of the various numerical solutions are shown in Table 1.Since SPH, VBM and MPS belong to the strong solutions and are mesh-less methods, the low computational cost becomes very important in some cases.e e e e i j k be a position vector, where the suffix in Greek letter is used to refer to the component of the coordinates, and the summation convention is used for the Greek suffixes.We define a function ( ) η x as an integral of ( ) ( ) ( )

Discretization Used in SPH and the Extension
x with respect to a volume V, where ( ) ( ) ξ x and ( ) ρ x are a weight function, an auxiliary function and the density of a fluid, respectively: The function ( ) η x and ( ) ξ x can be scalars, vectors and tensors.If we introduce an discretization of the vo- lume V and the density ρ as , 0,1, , 1, where the suffix in Roman letter is used to discriminate the point, and j M is the mass of the volume element d j V .We have an approximation of Equation ( 1): Furthermore, if we assume for the weight w If σ tends to 0, ( ) w x tends to ( ) δ x , where ( ) δ x is Dirac's delta function, namely: Using Equation ( 6), we have a following approximation: From Equations ( 4) and ( 7), we obtain If we assume ( ) , then, we have from Equation ( 8) Hence, we have This means that Equation ( 9) is a plausible approximation.Substituting Equation ( 9) into Equation ( 8), we obtain Setting x to i x in Equation (12), we have using Equation (10) ( ) ( ) ( ) ( ) on both sides of Equation ( 8), we obtain If we assume ( ) ( ) We derive Substituting Equations ( 9), ( 16) and (17) into Equation ( 14), we obtain Setting x to x i in Equation ( 18), we have where we assume ( ) Now, we derive formulas for second derivatives.Firstly, from Equation (14), we have .
Setting x to i x in Equations ( 22), ( 23) and (24), we where we assume ( ) 0 w′′ is finite.From Equation (19), we obtain for scalar φ and vector a [ ] ( ) ( ) From Equation ( 27), we obtain for scalar φ ( ) where d is the number of the dimension.

Application to Vortex Blob Method (VBM)
The basic equations for an incompressible viscous flow are given by Navier-Stokes equation [7]: where µ is the coefficient of viscosity, and the continuity equation: The pressure p satisfy ( ) The viscous flow is determined by Equations ( 31) or (32), ( 33) and (34).
If we introduce vorticity ω, the basic equations can be rewritten as Equation ( 36) can be written as where ν µ ρ = is the kinematic viscosity.If ω is determined, the velocity u is obtained by an integral representation [7]: , respectively, and The second term on the right hand side of Equation ( 38) is usually called Bio-Savart law.
Vortex Blob Method (VBM) uses Equations ( 37) and ( 38) and discretizes the vortex field as an assembly of concentrated vortices.The position of the concentrated vortices are determined by where a is the Lagrangian coordinates.It becomes very important how to obtain ⋅ ∇u ω and ∆ω .In this problem, we use the vorticity ω corresponds to the density ρ in Section 2.Then, we consider instead of Equation ( 1) Equation ( 3) is replaced by Then, we have instead of Equations ( 8) and ( 9) In all equations in Section 2, if we replace ρ and M with ω and Ω, respectively, we obtain the differential formulas in Section 3.
Now, we consider a weight ( ) where r = x .The weight ( ) w r has an asymptotic form: ( ) This asymptotic form is similar to ( ) k r defined by Equation (4) in Ref. [5].If e r is big, ( ) w r satisfies Eq- uation (49).In this case, Equation (52) has similar forms as given by Equation (2) in Ref. [5], if we assume  .The author has once discussed this problem from a different angle [1] and had the conclusion that the discrete differential operators used in MPS, especially, Laplace operator don't have the strict background from the mathematical viewpoint.Those operators should be considered to be a kind of experimental formulas that is derived numerically.If we apply them to irregular grids, the results vary irregularly.We should say the operator estimate the derivatives statistically.Quite naturally, the results are not unique.It may give a good estimate in one time and a wrong result in the other time.Hence, if we need to verify the accuracy, we should obtain several numerical results of the same problem using the various discretizations of the region and take the statistical average.
Although we do not deny this kind of approach, we wish to ensure the reliability.For the purpose, we need much discussion.Recently, a new paper [6] discusses the accuracy of Laplace operator in MPS in detail.Hopefully, we wish to prove mathematically, if possible, that, if we decrease the grid size zero, then, the numerical error approaches stably to zero.

Numerical Verification of Discretized Differential Operators of SPH with Gaussian Kernel
Numerical calculations were conducted to verify the validity of the discretized differential operators such as Equations ( 28) and (30).For simplicity, one-dimensional (1D) cases were considered.The first and second derivates of a scalar function ( ) were calculated using Equations ( 28) and (30), respectively.
The region is divided into N intervals.First, we consider the uniform mesh.Let dx i and x i ( ) As the weight function ( ) w r , we use Gaussian function: ( ) The weight function ( ) The 1D version of Equations ( 28) and (30) are given by ( .

Uniform Density and Regular Mesh
The density i ρ and the length of element d i x are given by 1 0, 1, , where the computational region is defined as 0 x L < < .

Verification of Laplace Operator
From Equation (62), we have The numerical results are shown in Figure 3.If we extend the region beyond the boundaries, the estimations within the original region are improved.

Non-Uniform Density and Regular Mesh
The density i ρ and the length of element d i x are given by ( ) respectively.The other parameters are same as in Section 5.1.1.
The numerical results are shown in Figure 4.The numerical result agrees well with the exact result.

Uniform Density and Non-Uniform Mesh
The density i ρ and the length of element d i x are given by 1 0,1, , 1, respectively.The other parameters are same as in Section 5.1.1.
The numerical results are shown in Figure 5.The numerical result agrees well with the exact result.

Non-Uniform Density and Non-Uniform Mesh
The density i ρ and the length of element d i x are given by ( )   ( ) ( ) respectively.The other parameter are same as in Section 5.1.1.
The numerical results are shown in Figure 6.The numerical result agrees well with the exact result.

1D Formulas for Discrete Gradient and Laplacian Operator
For convenience, we summarize the 1D discrete differential operators used in SHP and MPS as follows.From Equations ( 61) and ( 62), we have for SHP From Refs.[5] and [8], we have for MPS ( ) ( ) where ( ) ,

Effects of Mesh and Weight on Estimation of Gradient and Laplacian of a Given Function
We summarize the meshes and weights used in the following in Table 2 and Table 3, respectively.The parameters used in Sections 6.2.1-6.2.5 are given below.In Sections 6.2.1-6.2.5, the computational results are shown in Figures 7-11.The results are summarized Table 4 and Table 5 in Section 6.2.6.

Regular Mesh
As shown in Figure 7, there are no big errors in all methods SPH0, SPH1, SPH2 and MPS.The accuracy of SPH0 is very high.( ) Table 5. Summary on d 2 φ/dx 2 (○: good, △: not good, ×: bad).

Algebraic Mesh
As shown in Figure 8, a big error occurs in 2 2 d dx φ of MPS.

Geometric Mesh
As shown in Figure 9, a big error occurs in 2 2 d dx φ of MPS.

Random Mesh
As shown in Figure 10, the errors are rather small in all methods SPH0, SPH1, SPH2 and MPS not only for d dx φ but also for 2

Summary of Results
From the above mentioned numerical results, the following summaries are obtained.

1D Fluid Motion without Pressure and Viscosity
The motion of vast number of particles distributed in space under the action of the gravitational force may be treated as a fluid motion without pressure and viscosity [9] [10]: where ρ , u and Π is the density, velocity and gravitational potential.G is the gravitational constant.The solution of the problem defined above by Eulerian method is given as follows: 1) At time t, assume ρ , u and Π are given.
2) t ρ ∂ ∂ and u t ∂ ∂ are obtained from Equations (72a) and (72b), and Π is obtained from Equation (72c).3) ρ , and u at time d t t + is calculated.4) Repeat this process.In the solution by Lagrangian method, first, the equations in Eulerian form are transformed into Laglangian form: The following procedures give the solution of the problem defined above by Lagrangian method: 1) At time t, assume ρ , u and Π are given.
3) x , ρ and u of the material point at time x is obtained, then, d i x and i ρ is obtained as shown below: ( ) ( ) ( ) Hence, from the theoretical viewpoint, this problem can be solved without using the gradient operator.However, we use the gradient operator to obtain i ρ using the continuity equation.

1) Trapezoidal Distribution of the Initial Density
The initial conditions are given by ( ) The boundary conditions are specified as The computational conditions are as follows:  The numerical results are shown in Figure 14.The FDM uses the central differences for the first and second space derivatives, and the precision of FDM is considered high, if we neglect the spurious oscillation.The distribution pattern of MPS is different from those of FDM and SPH0.

1D Fluid Motion without Pressure but with Viscosity
In the present section, the viscosity is introduced: , where ν is the kinematic viscosity.Since the discrete Laplacian operator of SPH0 generates a big error at the discontinuity, the initial density distribution was smoothened using the five point running average, and, in the second example below, an small artificial numerical viscosity μ = 0.0000005 in case of N = 41 was added at every time step: ( ) ( ) In the third examples below, a big difference has occurred between SPH and MPS solutions.

1) Exponential Distribution of the Initial Density
The initial conditions are given by ( ) The numerical results are shown in Figure 15.The FDM uses the central differences for the first and second space derivatives, and the precision of FDM is considered high.In this example, the FDM, SPH and MPS solutions becomes similar. 2

) The Trapezoidal Distribution of the Initial Density
The initial conditions are given by ( ) where the initial density distribution was smoothened using the five point running average: ( ) The initial conditions are given by ( ) 0.5 when 0.25 2 0.5 when 0.25 0.75 0.5 when 0.75 The boundary conditions are specified as The computational conditions are as follows: Since we assume G = 0 in this example, the equation becomes Burgers equation.

a) Comparison of FDM Solution with the Analytical One
The central differences were used for the first and second space derivatives.As shown in Figure 17.The FDM solution is a very good approximation of the analytical solution [11].
b) Comparison of FDM, SPH0 and MPS Solutions Figure 18 and Figure 19 show the comparisons of ρ and u among FDM, SPH0 and MPS solutions.SPH0 and MPS solutions do not give good approximations.With respect to the velocity u, the difference between SPH0 and MPS is big.

Modified Gaussian Weight
Gaussian-type weights of finite support with C1 continuity are given as follows: ( )  where α satisfies The first and the second derivatives of ( ) w x given by Equation (90a) are given by ( ) Numerical examples are shown in Figure 20.When h is small, the accuracy becomes low.As has already pointed out in Section 5.1.1,there exist large errors at the boundaries.The numerical results agree well with the exact ones on overall.In the examples, γ is equal to dx = L/N = 0.1 and the accuracy seems sufficient when h is bigger than or equal to 4dx.

Conclusions
The author has once shown how to obtain the partial derivative in an irregular mesh using the moving least The discrete gradient and Laplacian operators of SPH include the first, first and second derivatives of the weight function with respect to the space coordinates, respectively.On the other hand, those of MPS include only the weight function itself.This may be the biggest difference between the discrete differential operators in SPH and MPH.
In the present paper, solutions by FDM (Finite Difference Method), SPH and MPS are compared numerically in detail.
1) Effects of mesh on the discrete gradient and Laplacian operators were studied.MPS showed good results with respect to the random mesh.
2) The FDM, SPH and MPS were applied to the initial value problems, and the effects of the difference of the solution method were studied.In some cases, the solutions of initial value problem showed a difference between SPH and MPS.
3) The discrete Laplacian operator of SPH is sensitive to the spacial discontinuities of the solution function.Hence, in the case of the initial value problem, the discontinuity in the initial condition should be smoothened beforehand, and the small amount of the artificial viscosity should be introduced.
4) The author has once studied a mathematical background of MPS theoretically [1] and did some numerical calculations.In very limited cases, the discrete Laplacian operator of MPS can be obtained mathematically.However, the generalization to the general mesh was not obtained.Hence, we are obliged to apply the Laplacian operator as a bold approximation in case of the general mesh.
5) Recently, Ng, Hwang and Sheu [6] discussed the accuracy of the discrete Laplacian operator of MPS theoretically and numerically.They clarified an important aspect of the properties of the operator.As one of the properties, they pointed out in 2 in conclusion of Ref. [6] that MPS gave generally a favorable result in case of the irregular mesh.We also had a similar impression as expressed in (1) above.
The support of Gaussian weight in SPH is infinite.In the present paper, weights of a Gaussian-type of finite support with C1 continuity were also given.

(
r with a small parameter 0 e r ε < = and constant α : be the length of the interval and the midpoint of the interval i, respectively: results are shown in Figure 2.Although there exist large errors at the boundaries as shown in Figure 2(a), the numerical results agree very well with the exact except the neighborhood of the boundaries.If we extend the region beyond the boundaries, the estimations within the original region are improved as shown in Figure 2(b).

Figure 2 .
Figure 2. A comparison between calculated and exact values of dφ/dx ((a) Not using extension of region; (b) Using extension of region).

Figure 3 .
Figure 3.A comparison between calculated and exact values of d 2 φ/dx 2 ((a) Not using extension of region; (b) Using extension of region).

.
The errors in MPS are surprisingly small.6.2.5.Sinusoidal MeshAs shown in Figure11, the errors are rather big in all methods SPH0, SPH1, SPH2 and MPS not only for d φ in MPS.The errors in SHP0 are smaller than those in the other methods.

Figure 14 .
Figure 14.Rectangular distribution of the initial density with G = 0.0015.

ρ in x 0 ρ in x FDM SPH0 MPSFigure 15 .
Figure 15.Exponential distribution of the initial density.
for SPH0 and 2 for MPS end cor N = .

Figure 17 .
Figure 17.Comparison of FDM solution with the analytical one ((a) FDM solution; (b) Exact solution).

Table 2 .
Classification of mesh.