Received 12 October 2014; revised 13 November 2014; accepted 24 November 2014

ABSTRACT

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.

Keywords:

SPH, MPS, VBM, Concentration of Mass or Vorticity, Estimation of Space Derivatives

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

2. Discretization Used in SPH and the Extension

Let 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 as an integral of with respect to a volume V, where, and are a weight function, an auxiliary function and the density of a fluid, respectively:

. (1)

The function and can be scalars, vectors and tensors. If we introduce an discretization of the volume V and the density as

(2)

(3)

where the suffix in Roman letter is used to discriminate the point, and is the mass of the volume element. We have an approximation of Equation (1):

. (4)

Furthermore, if we assume for the weight w

(5a)

and

(5b)

If tends to 0, tends to, where is Dirac’s delta function, namely:

(6)

Using Equation (6), we have a following approximation:

. (7)

From Equations (4) and (7), we obtain

. (8)

If we assume, then, we have from Equation (8)

. (9)

Substituting, we obtain

. (10)

Hence, we have

. (11)

This means that Equation (9) is a plausible approximation.

Substituting Equation (9) into Equation (8), we obtain

. (12)

Setting to in Equation (12), we have using Equation (10)

. (13)

Operating on both sides of Equation (8), we obtain

. (14)

If we assume

. (15)

We derive

, (16)

. (17)

Substituting Equations (9), (16) and (17) into Equation (14), we obtain

. (18)

Setting x to x_i in Equation (18), we have

. (19)

, (20)

where we assume.

Now, we derive formulas for second derivatives. Firstly, from Equation (14), we have

(21)

From Equations (16) and (17), we obtain

(22)

(23)

where is Chronecker’s delta: and when and, respectively. Substitut- ing Equations (9), (17), (22) and (23) into Equation (21), we derive

(24)

Setting to in Equations (22), (23) and (24), we have

, (25)

, (26)

(27)

where we assume is finite. From Equation (19), we obtain for scalar and vector a

, (28)

. (29)

From Equation (27), we obtain for scalar

(30)

where d is the number of the dimension.

3. Application to Vortex Blob Method (VBM)

The basic equations for an incompressible viscous flow are given by Navier-Stokes equation [7] :

, (31)

or

, (32)

where is the coefficient of viscosity, and the continuity equation:

. (33)

The pressure p satisfy

. (34)

The viscous flow is determined by Equations (31) or (32), (33) and (34).

If we introduce vorticity ω, the basic equations can be rewritten as

, (35)

. (36)

Equation (36) can be written as

, (37)

where is the kinematic viscosity.

If ω is determined, the velocity u is obtained by an integral representation [7] :

(38)

where when, , , respectively, and

, (39a)

(39b)

(39c)

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

, (40)

where a is the Lagrangian coordinates. It becomes very important how to obtain and.

In this problem, we use the vorticity ω corresponds to the density ρ in Section 2. Then, we consider instead of Equation (1)

. (41)

Equation (3) is replaced by

. (42)

Then, we have instead of Equations (8) and (9)

. (43)

. (44)

In all equations in Section 2, if we replace ρ and M with ω and Ω, respectively, we obtain the differential formulas in Section 3.

4. Mathematical Background of Moving Particle Semi-Implicit Method (MPS)

From Equations (19), (28) and (30), we have

, (45)

, (46)

(47)

If the weight satisfies

, (48)

then, we obtain

, (49)

where is a constant. Then, we have

(50)

Substituting Equations (49) and (50) into Equations (45), (46) and (47), we obtain

, (51)

, (52)

(53)

Now, we consider a weight with a small parameter and constant:

(54)

where. The weight has an asymptotic form:

(55)

This asymptotic form is similar to defined by Equation (4) in Ref. [5] . If is big, satisfies Equation (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.

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

(56)

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 be the length of the interval and the midpoint of the interval i, respectively:

(57), (58)

As the weight function, we use Gaussian function:

. (59)

The weight function satisfies

. (60)

The 1D version of Equations (28) and (30) are given by

, (61)

(62)

5.1. Uniform Density and Regular Mesh

The density and the length of element are given by

(63a)

(63b)

where the computational region is defined as.

5.1.1. Verification of Gradient Operator

From Equations (61), (63) and (45), we have

. (64)

We used parameters:, and. The numerical 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).

5.1.2. Verification of Laplace Operator

From Equation (62), we have

(65)

The numerical results are shown in Figure 3. If we extend the region beyond the boundaries, the estimations within the original region are improved.

5.2. Non-Uniform Density and Regular Mesh

The density and the length of element are given by

(66a)

(66b)

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.

5.3. Uniform Density and Non-Uniform Mesh

The density and the length of element are given by

(67a)

(67b)

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.

5.4. Non-Uniform Density and Non-Uniform Mesh

The density and the length of element are given by

(68a)

(68b)

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.

6. Comparison between SPH and MPS

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

, (69a)

(69b)

. (69c)

From Refs. [5] and [8] , we have for MPS

, (70a)

, (70b)

where

(71a)

(71b)

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

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

6.2.2. Algebraic Mesh

As shown in Figure 8, a big error occurs in of MPS.

6.2.3. Geometric Mesh

As shown in Figure 9, a big error occurs in of MPS.

6.2.4. Random Mesh

As shown in Figure 10, the errors are rather small in all methods SPH0, SPH1, SPH2 and MPS not only for but also for. The errors in MPS are surprisingly small.

6.2.5. Sinusoidal Mesh

As shown in Figure 11, the errors are rather big in all methods SPH0, SPH1, SPH2 and MPS not only for but also for except in MPS. The errors in SHP0 are smaller than those in the other methods.

6.2.6. Summary of Results

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

6.3. Application to Solution of Initial Value Problem

6.3.1. 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] :

, (72a)

, (72b)

, (72c)

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) and are obtained from Equations (72a) and (72b), and Π is obtained from Equation (72c).

Figure 7. Regular mesh.

Figure 8. Algebraic mesh.

Figure 9. Geometric mesh.

Figure 10. Random mesh.

Figure 11. Sinusoidal mesh.

3), and u at time is calculated.

4) Repeat this process.

In the solution by Lagrangian method, first, the equations in Eulerian form are transformed into Laglangian form:

, (73a)

, (73b)

, (73c)

. (73d)

The following procedures give the solution of the problem defined above by Lagrangian method:

1) At time t, assume, u and Π are given.

2), and is obtained from Equations (73a), (73b) and (73c), and Π is obtained from Equation (73d).

3), and u of the material point at time is calculated.

4) Repeat this process.

If is obtained, then, and is obtained as shown below:

(74a)

and

. (74b)

Hence, from the theoretical viewpoint, this problem can be solved without using the gradient operator. However, we use the gradient operator to obtain using the continuity equation.

1) Trapezoidal Distribution of the Initial Density

The initial conditions are given by

, (75a)

(75b)

. (75c)

The boundary conditions are specified as

(76a)

(76b)

The computational conditions are as follows:

, , , , , ,

, ,.

The numerical results are shown in Figure 12. The FDM (Finite Difference Method) uses the central differ-

Figure 12. Trapezoidal distribution of the initial density.

ences for the first and second space derivatives, and the precision of FDM is considered high. The Eulerian solution is used for FDM, and the Lagrangian solution is used for SPH0 and MPS. The distribution pattern of MPS is slightly different from those of FDM and SPH0.

2) Rectangular Distribution of the Initial Density with G = 0.001

The initial conditions are given by

, (77a)

(77b)

. (77c)

The boundary conditions are specified as

(78a)

(78b)

The computational conditions are as follows:

, , , , , ,

, ,.

The numerical results are shown in Figure 13. 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 slightly different from those of FDM and SPH0.

Figure 13. Rectangular distribution of the initial density with G = 0.001.

3) Rectangular Distribution of the Initial Density with G = 0.0015

The initial conditions are given by

(79a)

(79b)

(79c)

The boundary condition are specified as

(80a)

(80b)

Computational condition

, , , , , ,

, ,.

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.

6.3.2. 1D Fluid Motion without Pressure but with Viscosity

In the present section, the viscosity is introduced:

(81a)

(81b)

(81c)

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:

(82a)

(82b)

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

(83a)

(83b)

(83c)

The boundary conditions are specified as

(84a)

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

(84b)

The computational conditions are as follows:

, , , , , ,

, ,.

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

, (85a)

Figure 15. Exponential distribution of the initial density.

(85b)

, (85c)

where the initial density distribution was smoothened using the five point running average:

. (86)

The boundary conditions are specified as

(87a)

(87b)

The computational conditions are as follows:

, , , , , ,

, ,.

The numerical results are shown in Figure 16. The FDM uses the central differences for the first and second space derivatives, and the precision of FDM is considered high. In this example, small difference is observed in the FDM, SPH and MPS solutions.

3) The Flat-Slope-Flat Distribution of the Initial Velocity

Figure 16. Trapezoidal distribution of the initial density.

The initial conditions are given by

(88a)

(88b)

(88c)

The boundary conditions are specified as

(89a)

(89b)

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.

7. Modified Gaussian Weight

Gaussian-type weights of finite support with C1 continuity are given as follows:

(90a), (90b)

Figure 18. Comparison of density ρ.

where satisfies

(91)

The first and the second derivatives of given by Equation (90a) are given by

, (92a)

Figure 19. Comparison of velocity u.

. (92b)

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.

8. Conclusions

The author has once shown how to obtain the partial derivative in an irregular mesh using the moving least

Figure 20. Effect of support 2h (; L = 4, N = 40, γ = 0.1, h = 0.2, 0.4, 0.8).

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). We have extended SPH theoretically in the present paper.

In Vortex Blob Method (VBM), the abduction is approximated reasonably. However, the diffusion of vortices can’t be approximated precisely. We have shown in the present paper that this problem can be solved, if we apply the ideas developed by SPH.

Moving Particle Semi-implicit method (MPS) uses a similar dscretization of the initial and boundary value problem as SPH and VBM. However, the mathematical background of MPS is not sufficient. We have shown in the present paper that the mathematical background of the discrete gradient operator is strengthened by applying the ideas developed by SPH. However, that of the discrete Laplacian operator could not given.

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.

Acknowledgements

The author thanks sincerely Prof. G. Yagawa, Toyo Univ., Japan for his support and encouragement.

References