# Development of a Particle Interaction Kernel Function in MPS Method for Simulating Incompressible Free Surface Flow

Authors: Tony W. H. Sheu; Chenpeng Chiao; Chinlong Huang
Journal: Journal of Applied Mathematics
 (2011)

---

## Abstract

We  aimed to derive a kernel function that accounts for the interaction among moving particles within the framework of particle method. To predict a computationally more accurate moving particle solution for the Navier-Stokes equations, kernel function is a key to success in the development of interaction model. Since the smoothed quantity of a scalar or a vector at a spatial location is mathematically identical to its collocated value provided that the kernel function is chosen to be the Dirac delta function, our guideline is to derive the kernel function that is closer to the delta function as much as possible. The proposed particle interaction model using the newly developed kernel function will be validated through the two investigated Navier-Stokes problems which have either the semianalytical or the benchmark solutions.


---

## Body

## 1. Introduction

Numerical methods that are now available for performing flow simulations can be divided into the Eulerian and Lagrangian two classes. In the Eulerian method, flow equations are solved at a fixed mesh system. On the contrary, in the Lagrangian method, which can be subdivided into the mesh-based and meshfree two subgroups, either meshes or particles, on the other hand, are advected with the fluid flow. The consequence is that in the transport equations we no longer have a need to approximate the convection terms. Avoidance of the approximation of convective terms is now known as one of the apparent advantages of employing the Lagrangian approach. Numerical errors of the cross-wind diffusion type can, as a result, be well eliminated. There exists another mathematically rigorous arbitrary Lagrangian-Eulerian (ALE) method [1] which involves the use of Lagrangian and Eulerian solution steps in the course of moving mesh.In the literature, several well-known grid-based methods, such as the marker-and-cell (MAC) [2], volume-of-fluid (VOF) [3], and level-set [4] methods, have been developed to predict free surface flows. There existed also a major class of meshfree methods, which are known as the particle methods. Particle methods deal with the prescribed particles moving in the Lagrangian sense. As a result, the convection terms can be directly calculated without incurring any numerical diffusion error. Particle methods can be also separated into the Eulerian particle method, such as the particle-in-cell (PIC) method [5] and the Lagrangian particle methods, which include the smoothed particle hydrodynamics (SPHs) and moving particle semi-implicit (MPS) two classes of meshless methods. SPH method, introduced firstly by Lucy [6] and Gingold and Monaghan [7] in 1977, was developed mainly for the simulation of compressible fluid flows based on the interpolation theory through the introduced kernel function (or smoothing kernel). SPH method was later extended to simulate also the incompressible free surface flow [8].Another gridless particle method, called the MPS, was more recently developed to simulate the incompressible Navier-Stokes fluid flows [9]. The motion of each particle in this method is calculated through the interaction of its neighboring particles by means of the kernel (or weight) function. This means that all the spatial derivatives can be calculated by the deterministic particle interaction without the need of generating meshes in the flow. This explains why MPS method is gradually becoming effective for use in simulating many practical problems involving the complicated geometry and complex physics. For the flow problems with an inflow or an outflow boundary, this method was however found to have difficulty of tracing a particle easily. In addition, MPS method requires an extra CPU time to search for all its neighboring points.The rest of this paper will be organized as follows. In Section2 we present the Navier-Stokes and the passive scalar equations, which involve convection and diffusion flux terms. This is followed by presenting the particle interaction models used for the approximation of the gradient and Laplacian differential operators. In Section 4 emphasis will be placed on subject to free surfaces the kernel function used in the particle method. In Section 5, we will validate the proposed particle model by solving two Navier-Stokes problems. In Section 6, conclusions will be drawn based on the predicted results.


## 2. Working Equations

In this study the following Navier-Stokes equations for the motion of an incompressible fluid flow will be considered in a simply connected domainΩ:(2.1)Du̲Dt=-1ρ∇p+ν∇2u̲+f̲,(2.2)1ρDρDt=-∇⋅u̲,
where Du̲/Dt=∂u̲/∂t+(u̲·∇)u̲. In (2.1), which is the momentum equation, ρ stands for the density of the investigated flow, ν the kinematic viscosity, and f̲ the gravitational acceleration vector.The following model equation will be considered in the development of kernel function since it is the key equation in the simulation of momentum equations for the incompressible fluid flow:(2.3)u∂ϕ∂x+v∂ϕ∂y=μ(∂2ϕ∂x2+∂2ϕ∂y2).
Both of the velocity components u and v and the fluid viscosity μ in (2.3) are assumed to be the constant values to facilitate the development of the proposed kernel function in Section 4.


## 3. Particle Interaction Models for Two Differential Operators

Due to the difficulties of applying the grid-based methods to simulate complex flow physics, particle-based methods become gradually popular in the simulation of incompressible Navier-Stokes equations subject either to a free surface or an interface. In particle methods, the differential operators for the mass and momentum conservation equations need to be replaced by their corresponding particle interaction operators. In other words, partial differential equations (2.1)-(2.2) written in the continuous context will be transformed to their corresponding discrete particle interaction equations so that the transport equations under investigation can be approximated by the chosen moving particles and their interaction. The degree of success depends on the kernel (or weighting) function for the particles that are distanced from each other by the user’s prescribed finite length for the approximation of velocity vectors and pressure in the incompressible fluid flows around these particles.Consider a particle at which a physical quantityfi is defined at a position r̲i. One can represent fi approximately at each location as follows by virtue of the following kernel function w(r̲)(3.1)〈f(r̲)〉⋅∑iw(|r̲i-r̲|)=∑ifiw(|r̲i-r̲|).
It is important to note that the smoothed quantity 〈f〉i for f at r̲i turns out to be exactly identical to the local value of fi if the kernel function in (3.1) is the Dirac delta function. This implies that the chosen kernel function w(r̲), which should be constrained by the integral constraint equation ∫Vw(r)dV=1 shown in (3.1), determines the prediction quality of the particle methods. One can find different kernel functions in [10].For a scalarϕ at a location r̲j, its value can be expanded in Taylor series with respect to the value of ϕ at r̲i as follows:(3.2)ϕj=ϕi+∇ϕ∣ij⋅(r̲j-r̲i)+H.O.T.
By dropping the higher-order terms shown above, we are led to get the first-order approximated equation ϕj-ϕi=∇ϕ∣ij·(r̲j-r̲i). By multiplying the term (r̲j-r̲i)-1 on both hand sides of the resulting simplified equation, we can get (3.3)∇ϕ∣ij=(ϕj-ϕi)(r̲j-r̲i)|r̲j-r̲i||r̲j-r̲i|.
Let ∇ϕ∣ij be f. We then substitute it into (3.1) to get the following smoothed representation of ∇ϕ, which is denoted by 〈∇ϕ〉 at the node ij, in a domain with d dimensions (3.4)〈∇ϕ〉∣ij=dn0∑j≠iϕj-ϕi|r̲j-r̲i|(r̲j-r̲i)w(r̲j-r̲i).
In the above, n0 denotes the particle number density and it is defined as ∑j≠iw(|r̲j-r̲i|). Note that ni=n0 is true only under the incompressible flow condition.One can similarly derive the Laplacian operator for a scalarϕ, which has been derived in [9] as follows(3.5)〈∇2ϕ〉∣i=2dλn0∑j≠i(ϕj-ϕi)w(|r̲j-r̲i|),
where      (3.6)λ=∫V′|r̲j-r̲i|2w(|r̲j-r̲i|)dV∫V′w(|r̲j-r̲i|)dV.
Note that V′ shown above denotes the volume excluding of a small interval that contains the point at i. One can also employ the Laplacian operator given below [11](3.7)〈∇2ϕ〉∣ij=2dn0∑j≠iϕj-ϕi|r̲j-r̲i|w(|r̲j-r̲i|).The advantage of employing particle methods becomes clear in that all the spatial derivatives can be calculated from the chosen kernel function without the need of generating meshes in the physical domain. Much of the computational effort in the generation of a good quality mesh for performing an accurate numerical simulation in a domain with moving boundaries is therefore avoided.Unlike the SPH particle method, calculations of the values for∇ϕ and ∇2ϕ in (3.4), (3.5), and (3.7) involve only the kernel function w(r̲). Since the derivative of kernel function needs not to be calculated in the approximation of differential operators ∇ and ∇2, numerical oscillations, which may be generated in the fixed grid Eulerian approach for the cases with high solution gradients, can be completely avoided. As a result, one has the flexibility to choose the kernel function that has a slope as steep as the Dirac delta function.


## 4. Development of the Kernel Functions

In view of (3.4) and (3.5)–(3.7), we know that the quality of approximating the operators 〈∇ϕ〉ij and 〈∇2ϕ〉i depends entirely on the chosen kernel function w(r̲) and the number of the prescribed particles. Moreover, the chosen particle locations and the employed kernel function will determine the subsequent particle location in moving particle methods. For this reason, the chosen kernel function will be derived below in detail.In the light of (3.1), we know that the Dirac delta function δ(r̲), which is constrained by ∫-∞∞δ(r̲)dr=1, is an ideal kernel function. Such a function is, unfortunately, not implementable in computational practice. We have therefore to resort to its corresponding smoothed delta function when carrying out the simulations based on the particle method. The smoothed function is sometimes called as the nascent delta function δε(r̲), which is defined as limε→0δε(r̲)=δ(r̲). In the literature, one can find other nascent delta functions such as the Gaussian function, Lorentz line function, impulse function, and sinc function. There exists also the other class of kernel functions which were developed irrelevantly to the nascent delta functions. Typical examples include the exponential, cubic spline, and quadratic spline functions proposed by Belytschko et al. [12] and the kernel functions proposed by Koshizuka and Oka in 1996 [13] and Koshizuka et al. in 1998 [11].


### 4.1. Development of Kernel Function for the Pure Diffusion Equation

In this study a new kernel function will be rigorously developed and it will be used to simulate the incompressible Navier-Stokes equations using the particle methods. The guidelines of developing the proposed kernel function will be given below. The kernel functionw(r̲) under current development falls into the category of nascent delta function. This implies limre→0w(r̲,re)=δ(r̲),, where re is the radius of a small circle. The weight between two particles that have a distance r apart will be diminished as r≥re. For the sake of accuracy, the kernel function will be developed to retain the following constraint condition in the Dirac delta function(4.1)∫-∞∞w(r̲,re)dr̲=1.Development of the current kernel function starts with representing it in terms ofr/re as (4.2)w(r)={are+bre(rre)+cre(rre)2+dre(rre)3+ere(rre)4;0≤r≤re,0;re<r.
Derivation of w(r̲) is followed by imposing the four constraint conditions given by w(r=re)=(∂w/∂r)∣r=re=(∂w/∂r)∣r=0=0 and ∫0r/re=1w(r)dr=1/2. These imposed conditions enable us to derive the algebraic equations given below(4.3)a+b+c+d+e=0,b+2c+3d+4e=0,a+b2+c3+d4+e5=12,b=0.
Given the above four equations, one can then easily express a, b, c, d in terms of e as a=1-(1/15)e, b=0, c=-3+(6/5)e, d=2-(32/15)e.By substituting the resulting kernel function into (3.7) for ∇2ϕ, one can get the corresponding discrete equation for the Laplace equation ∂2ϕ/∂x2+∂2ϕ/∂y2=0 based on the particle method(4.4)1[w(h)+2w(2h)]h2{2h(ϕi+1,j+ϕi-1,j+ϕi,j+1+ϕi,j-1-4ϕi)w(h)+(ϕi+1,j+1+ϕi-1,j+1+ϕi+1,j-1+ϕi-1,j-1-4ϕi)w(2h)}=0.
By performing the modified equation analysis on (4.4), we are led to get the following modified equation:(4.5)∂2ϕ∂x2+∂2ϕ∂y2=-[w(2h)w(h)+2w(2h)∂4ϕ∂x2∂y2+112∂4ϕ∂x4+112∂4ϕ∂y4]h2-[w(2h)12(w(h)+2w(2h))(∂6ϕ∂x4∂y2+∂6ϕ∂x2∂y4)+1360∂6ϕ∂x6+1360∂6ϕ∂y6]h4+⋯.
Let w(2h)=(1/4)w(h) the first term on the right hand side of (4.5) turns out to be zero due to ∂2ϕ/∂x2+∂2ϕ/∂y2=0.Provided thatre=2h, where h denotes the grid size, we can easily know from the following nine-point equation that the resulting approximated equation for (2.3) investigated at the condition of μ=∞ has accuracy of the fourth order.(4.6)ϕi+1,j+1+ϕi-1,j+1+ϕi-1,j+1+ϕi+1,j-1+4ϕi+1,j+4ϕi-1,j+4ϕi,j+1+4ϕi,j-1-20ϕi,j6h2+O(h4)=0.
Note that (4.6) is derived under the condition of (4.7)3a+(22-12)b+(2-14)c+(2-18)d+(1-116)e=0.
The resulting five free parameters can then be uniquely determined from (4.3)–(17) and (4.7) as a=(480√2-705)/(512√2-745), b=0, c=(-9602+1515)/(5122-745), d=-210/(5122-745), e=480(2-600)/(5122-745).In summary, the kernel function developed for simulating the pure diffusion equation in this study is given as follows:(4.8)w(r)={4802-7055122-7451re+-9602+15155122-7451re(rre)2+-2105122-7451re(rre)3+4802-6005122-7451re(rre)4;0≤r≤re,0;re<r
under the condition given below (4.9)w(2h)=14w(h).
Note that the kernel function w(r)=re/r-1 given in [13] has an infinitely large magnitude as the value of r approaches zero. The newly developed kernel function shown in Figure 1, on the other hand, has its local maximum value at r=0. It is also worthy to mention that the currently developed kernel function applied to predict the pure diffusion equation satisfies the following constraint condition, which is also embedded in the delta function (4.10)∫r/re=-1r/re=1w(r)dr=1.Plot of the developed kernel function. (a) cast inr/re coordinate; (b) plotted in xy plane.
(a)(b)


## 4.1. Development of Kernel Function for the Pure Diffusion Equation

In this study a new kernel function will be rigorously developed and it will be used to simulate the incompressible Navier-Stokes equations using the particle methods. The guidelines of developing the proposed kernel function will be given below. The kernel functionw(r̲) under current development falls into the category of nascent delta function. This implies limre→0w(r̲,re)=δ(r̲),, where re is the radius of a small circle. The weight between two particles that have a distance r apart will be diminished as r≥re. For the sake of accuracy, the kernel function will be developed to retain the following constraint condition in the Dirac delta function(4.1)∫-∞∞w(r̲,re)dr̲=1.Development of the current kernel function starts with representing it in terms ofr/re as (4.2)w(r)={are+bre(rre)+cre(rre)2+dre(rre)3+ere(rre)4;0≤r≤re,0;re<r.
Derivation of w(r̲) is followed by imposing the four constraint conditions given by w(r=re)=(∂w/∂r)∣r=re=(∂w/∂r)∣r=0=0 and ∫0r/re=1w(r)dr=1/2. These imposed conditions enable us to derive the algebraic equations given below(4.3)a+b+c+d+e=0,b+2c+3d+4e=0,a+b2+c3+d4+e5=12,b=0.
Given the above four equations, one can then easily express a, b, c, d in terms of e as a=1-(1/15)e, b=0, c=-3+(6/5)e, d=2-(32/15)e.By substituting the resulting kernel function into (3.7) for ∇2ϕ, one can get the corresponding discrete equation for the Laplace equation ∂2ϕ/∂x2+∂2ϕ/∂y2=0 based on the particle method(4.4)1[w(h)+2w(2h)]h2{2h(ϕi+1,j+ϕi-1,j+ϕi,j+1+ϕi,j-1-4ϕi)w(h)+(ϕi+1,j+1+ϕi-1,j+1+ϕi+1,j-1+ϕi-1,j-1-4ϕi)w(2h)}=0.
By performing the modified equation analysis on (4.4), we are led to get the following modified equation:(4.5)∂2ϕ∂x2+∂2ϕ∂y2=-[w(2h)w(h)+2w(2h)∂4ϕ∂x2∂y2+112∂4ϕ∂x4+112∂4ϕ∂y4]h2-[w(2h)12(w(h)+2w(2h))(∂6ϕ∂x4∂y2+∂6ϕ∂x2∂y4)+1360∂6ϕ∂x6+1360∂6ϕ∂y6]h4+⋯.
Let w(2h)=(1/4)w(h) the first term on the right hand side of (4.5) turns out to be zero due to ∂2ϕ/∂x2+∂2ϕ/∂y2=0.Provided thatre=2h, where h denotes the grid size, we can easily know from the following nine-point equation that the resulting approximated equation for (2.3) investigated at the condition of μ=∞ has accuracy of the fourth order.(4.6)ϕi+1,j+1+ϕi-1,j+1+ϕi-1,j+1+ϕi+1,j-1+4ϕi+1,j+4ϕi-1,j+4ϕi,j+1+4ϕi,j-1-20ϕi,j6h2+O(h4)=0.
Note that (4.6) is derived under the condition of (4.7)3a+(22-12)b+(2-14)c+(2-18)d+(1-116)e=0.
The resulting five free parameters can then be uniquely determined from (4.3)–(17) and (4.7) as a=(480√2-705)/(512√2-745), b=0, c=(-9602+1515)/(5122-745), d=-210/(5122-745), e=480(2-600)/(5122-745).In summary, the kernel function developed for simulating the pure diffusion equation in this study is given as follows:(4.8)w(r)={4802-7055122-7451re+-9602+15155122-7451re(rre)2+-2105122-7451re(rre)3+4802-6005122-7451re(rre)4;0≤r≤re,0;re<r
under the condition given below (4.9)w(2h)=14w(h).
Note that the kernel function w(r)=re/r-1 given in [13] has an infinitely large magnitude as the value of r approaches zero. The newly developed kernel function shown in Figure 1, on the other hand, has its local maximum value at r=0. It is also worthy to mention that the currently developed kernel function applied to predict the pure diffusion equation satisfies the following constraint condition, which is also embedded in the delta function (4.10)∫r/re=-1r/re=1w(r)dr=1.Plot of the developed kernel function. (a) cast inr/re coordinate; (b) plotted in xy plane.
(a)(b)


## 5. Numerical Results

### 5.1. Solitary Wave Propagation

Study of the solitary wave propagation in a rectangular channel has received a considerable attention since the resulting analysis can provide engineers with useful information about the wave loading on the structures [14]. The problem under consideration is schematic in Figure 2 where the undisturbed water height is d and the channel width is 16d. Above the undisturbed water level, the crest of the investigated solitary wave is given by [15](5.1)y=d+Hsech2[3H4d3x].Figure 2
Schematic of the solitary wave propagation problem.Here,H(=2) and d(=10) denote the initial wave height and the still water depth, respectively. Initially, the velocities are given by(5.2)u=gdHdsech2[3H4d3x],v=3gd(Hd)3/2(yd)sech2[3H4d3x]tanh[3H4d3x],
where g(=9.81) is the gravitational acceleration. In the current study, the initial distance between the neighbor particles is 0.8, and 3482 particles were used in the computation. For the wave propagating towards the right wall, the predicted run-up height on the right wall in Figure 3 is compared with the result given in [15]. The computed wave pressures at different times are shown in Figure 4.Figure 3
Comparison of the predicted run-up heights of the solitary wave on the right wall.The predicted free surfaces and pressures at different times. (a)t=0.5; (b) t=2.0; (c) t=4.0; (d) t=6.0; (e) t=8.0; (f) t=10.0.
(a)(b)(c)(d)(e)(f)


### 5.2. Dam Break Flow

The investigated dambreak problem is schematic Figure5, where L(=1) denotes the height and the width of the stillwater in the tank. The computed results were compared with the results in [16]. Figures 6(a)–6(f) show the results for the fluid flow with Re=500.  Figures 7(a)–7(f) are the results predicted in the case with Re=1000. The red lines stand for the simulated free surfaces, and they are compared with those shown in [16].Figure 5
Schematic of the dam break problem.Comparison of the predicted free surfaces by the MPS and the level set methods at different times for the case withRe=500. (a) t=0.0; (b) t=0.5; (c) t=1.0; (d) t=1.5; (e) t=2.0; (f) t=2.5.
(a)(b)(c)(d)(e)(f)Comparison of the predicted free surfaces by the MPS and the level set methods at different times for the case withRe=500. (a) t=0.0; (b) t=0.5; (c) t=1.0; (d) t=1.5; (e) t=2.0; (lf) t=2.5.
(a)(b)(c)(d)(e)(f)


## 5.1. Solitary Wave Propagation

Study of the solitary wave propagation in a rectangular channel has received a considerable attention since the resulting analysis can provide engineers with useful information about the wave loading on the structures [14]. The problem under consideration is schematic in Figure 2 where the undisturbed water height is d and the channel width is 16d. Above the undisturbed water level, the crest of the investigated solitary wave is given by [15](5.1)y=d+Hsech2[3H4d3x].Figure 2
Schematic of the solitary wave propagation problem.Here,H(=2) and d(=10) denote the initial wave height and the still water depth, respectively. Initially, the velocities are given by(5.2)u=gdHdsech2[3H4d3x],v=3gd(Hd)3/2(yd)sech2[3H4d3x]tanh[3H4d3x],
where g(=9.81) is the gravitational acceleration. In the current study, the initial distance between the neighbor particles is 0.8, and 3482 particles were used in the computation. For the wave propagating towards the right wall, the predicted run-up height on the right wall in Figure 3 is compared with the result given in [15]. The computed wave pressures at different times are shown in Figure 4.Figure 3
Comparison of the predicted run-up heights of the solitary wave on the right wall.The predicted free surfaces and pressures at different times. (a)t=0.5; (b) t=2.0; (c) t=4.0; (d) t=6.0; (e) t=8.0; (f) t=10.0.
(a)(b)(c)(d)(e)(f)


## 5.2. Dam Break Flow

The investigated dambreak problem is schematic Figure5, where L(=1) denotes the height and the width of the stillwater in the tank. The computed results were compared with the results in [16]. Figures 6(a)–6(f) show the results for the fluid flow with Re=500.  Figures 7(a)–7(f) are the results predicted in the case with Re=1000. The red lines stand for the simulated free surfaces, and they are compared with those shown in [16].Figure 5
Schematic of the dam break problem.Comparison of the predicted free surfaces by the MPS and the level set methods at different times for the case withRe=500. (a) t=0.0; (b) t=0.5; (c) t=1.0; (d) t=1.5; (e) t=2.0; (f) t=2.5.
(a)(b)(c)(d)(e)(f)Comparison of the predicted free surfaces by the MPS and the level set methods at different times for the case withRe=500. (a) t=0.0; (b) t=0.5; (c) t=1.0; (d) t=1.5; (e) t=2.0; (lf) t=2.5.
(a)(b)(c)(d)(e)(f)


## 6. Concluding Remarks

The idea of developing the proposed particle interaction model is to develop the kernel function which accommodates the property∫r/re=-1r/re=1w(r)dr=1 embedded in the smoothed Dirac delta function. To justify the proposed kernel functions for simulating the realistic flow, we investigate two well-known benchmark problems for the validation sake. The predicted results of the dam break problem and the solitary wave problem clearly show that the solution predicted by the particle interaction model using the proposed kernel function can capture well the motion of the fluid flow subject to free surface.



---

*Source: 793653-2011-02-16.xml*