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

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.

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 Section 2 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 Ω:𝐷𝑢1𝐷𝑡=−𝜌∇𝑝+𝜈∇²𝑢+𝑓1,𝜌𝐷𝜌𝐷𝑡=−∇⋅𝑢, (2.1)(2.2) where 𝐷𝑢/𝐷𝑡=𝜕𝑢/𝜕𝑡+(𝑢⋅∇)𝑢. In (2.1), which is the momentum equation, 𝜌 stands for the density of the investigated flow, 𝜈 the kinematic viscosity, and 𝑓 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) Both of the velocity components 𝑢 and 𝑣 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 quantity 𝑓_𝑖 is defined at a position 𝑟_𝑖. One can represent 𝑓_𝑖 approximately at each location as follows by virtue of the following kernel function 𝑤(𝑟)𝑓𝑟⋅_𝑖𝑤𝑟_𝑖−𝑟=_𝑖𝑓_𝑖𝑤𝑟_𝑖−𝑟. (3.1) It is important to note that the smoothed quantity ⟨𝑓⟩_𝑖 for 𝑓 at 𝑟_𝑖 turns out to be exactly identical to the local value of 𝑓_𝑖 if the kernel function in (3.1) is the Dirac delta function. This implies that the chosen kernel function 𝑤(𝑟), which should be constrained by the integral constraint equation ∫_𝑉𝑤(𝑟)𝑑𝑉=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 𝑟_𝑗, its value can be expanded in Taylor series with respect to the value of 𝜙 at 𝑟_𝑖 as follows:𝜙_𝑗=𝜙_𝑖+∇𝜙∣_𝑖𝑗⋅𝑟_𝑗−𝑟_𝑖+𝐻.𝑂.𝑇. (3.2) By dropping the higher-order terms shown above, we are led to get the first-order approximated equation 𝜙_𝑗−𝜙_𝑖=∇𝜙∣_𝑖𝑗⋅(𝑟_𝑗−𝑟_𝑖). By multiplying the term (𝑟_𝑗−𝑟_𝑖)⁻¹ on both hand sides of the resulting simplified equation, we can get ∇𝜙∣_𝑖𝑗=𝜙_𝑗−𝜙_𝑖𝑟_𝑗−𝑟_𝑖𝑟_𝑗−𝑟_𝑖𝑟_𝑗−𝑟_𝑖. (3.3) Let ∇𝜙∣_𝑖𝑗 be 𝑓. We then substitute it into (3.1) to get the following smoothed representation of ∇𝜙, which is denoted by ⟨∇𝜙⟩ at the node 𝑖𝑗, in a domain with 𝑑 dimensions ∇𝜙∣_𝑖𝑗=𝑑𝑛⁰_𝑗≠𝑖𝜙_𝑗−𝜙_𝑖𝑟_𝑗−𝑟_𝑖𝑟_𝑗−𝑟_𝑖𝑤𝑟_𝑗−𝑟_𝑖. (3.4) In the above, 𝑛⁰ denotes the particle number density and it is defined as ∑_𝑗≠𝑖𝑤(|𝑟_𝑗−𝑟_𝑖|). Note that 𝑛_𝑖=𝑛⁰ 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∇²𝜙∣_𝑖=2𝑑𝜆𝑛⁰_𝑗≠𝑖𝜙_𝑗−𝜙_𝑖𝑤𝑟_𝑗−𝑟_𝑖, (3.5) where ∫𝜆=^𝑉′𝑟_𝑗−𝑟_𝑖²𝑤𝑟_𝑗−𝑟_𝑖𝑑𝑉∫_𝑉′𝑤𝑟_𝑗−𝑟_𝑖𝑑𝑉. (3.6) Note that 𝑉^′ 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]∇²𝜙∣_𝑖𝑗=2𝑑𝑛⁰_𝑗≠𝑖𝜙_𝑗−𝜙_𝑖𝑟_𝑗−𝑟_𝑖𝑤𝑟_𝑗−𝑟_𝑖. (3.7)

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 ∇²𝜙 in (3.4), (3.5), and (3.7) involve only the kernel function 𝑤(𝑟). Since the derivative of kernel function needs not to be calculated in the approximation of differential operators ∇ and ∇², 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 ⟨∇𝜙⟩_𝑖𝑗 and ⟨∇²𝜙⟩_𝑖 depends entirely on the chosen kernel function 𝑤(𝑟) 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 𝛿(𝑟), which is constrained by ∫^∞_−∞𝛿(𝑟)𝑑𝑟=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 𝛿_𝜀(𝑟), which is defined as lim_𝜀→0𝛿_𝜀(𝑟)=𝛿(𝑟). 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 function 𝑤(𝑟) under current development falls into the category of nascent delta function. This implies lim_𝑟𝑒→0𝑤(𝑟,𝑟_𝑒)=𝛿(𝑟),, where 𝑟_𝑒 is the radius of a small circle. The weight between two particles that have a distance 𝑟 apart will be diminished as 𝑟≥𝑟_𝑒. For the sake of accuracy, the kernel function will be developed to retain the following constraint condition in the Dirac delta function^∞_−∞𝑤𝑟,𝑟_𝑒𝑑𝑟=1. (4.1)

Development of the current kernel function starts with representing it in terms of 𝑟/𝑟_𝑒 as 𝑤(𝑟)=𝑎𝑟_𝑒+𝑏𝑟_𝑒𝑟𝑟_𝑒+𝑐𝑟_𝑒𝑟𝑟_𝑒²+𝑑𝑟_𝑒𝑟𝑟_𝑒³+𝑒𝑟_𝑒𝑟𝑟_𝑒⁴;0≤𝑟≤𝑟_𝑒,0;𝑟_𝑒<𝑟. (4.2) Derivation of 𝑤(𝑟) is followed by imposing the four constraint conditions given by 𝑤(𝑟=𝑟_𝑒)=(𝜕𝑤/𝜕𝑟)∣_{𝑟=𝑟𝑒}=(𝜕𝑤/𝜕𝑟)∣_𝑟=0=0 and ∫^{𝑟/𝑟𝑒}₀⁼¹𝑤(𝑟)𝑑𝑟=1/2. These imposed conditions enable us to derive the algebraic equations given below𝑏𝑎+𝑏+𝑐+𝑑+𝑒=0,𝑏+2𝑐+3𝑑+4𝑒=0,𝑎+2+𝑐3+𝑑4+𝑒5=12,𝑏=0. (4.3) Given the above four equations, one can then easily express 𝑎, 𝑏, 𝑐, 𝑑 in terms of 𝑒 as 𝑎=1−(1/15)𝑒, 𝑏=0, 𝑐=−3+(6/5)𝑒, 𝑑=2−(32/15)𝑒.

By substituting the resulting kernel function into (3.7) for ∇²𝜙, one can get the corresponding discrete equation for the Laplace equation 𝜕²𝜙/𝜕x²+𝜕²𝜙/𝜕y²=0 based on the particle method1𝑤(ℎ)+2𝑤√2ℎℎ²𝜙_𝑖+1,𝑗+𝜙_{𝑖−1,𝑗}+𝜙_𝑖,𝑗+1+𝜙_{𝑖,𝑗−1}−4𝜙_𝑖𝑤(ℎ)+𝜙_{𝑖+1,𝑗+1}+𝜙_{𝑖−1,𝑗+1}+𝜙_{𝑖+1,𝑗−1}+𝜙_{𝑖−1,𝑗−1}−4𝜙_𝑖𝑤√2ℎ=0. (4.4) By performing the modified equation analysis on (4.4), we are led to get the following modified equation:𝜕²𝜙𝜕𝑥²+𝜕²𝜙𝜕𝑦²𝑤√=−2ℎ𝑤(ℎ)+2𝑤√𝜕2ℎ⁴𝜙𝜕𝑥²𝜕𝑦²+1𝜕12⁴𝜙𝜕𝑥⁴+1𝜕12⁴𝜙𝜕𝑦⁴ℎ²−𝑤√2ℎ12𝑤(ℎ)+2𝑤√𝜕2ℎ⁶𝜙𝜕𝑥⁴𝜕𝑦²+𝜕⁶𝜙𝜕𝑥²𝜕𝑦⁴+1𝜕360⁶𝜙𝜕𝑥⁶+1𝜕360⁶𝜙𝜕𝑦⁶ℎ⁴+⋯. (4.5) Let w(√2ℎ)=(1/4)𝑤(ℎ) the first term on the right hand side of (4.5) turns out to be zero due to 𝜕²𝜙/𝜕𝑥²+𝜕²𝜙/𝜕𝑦²=0.

Provided that 𝑟_𝑒=2ℎ, where ℎ 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.𝜙_{𝑖+1,𝑗+1}+𝜙_{𝑖−1,𝑗+1}+𝜙_{𝑖−1,𝑗+1}+𝜙_{𝑖+1,𝑗−1}+4𝜙_𝑖+1,𝑗+4𝜙_{𝑖−1,𝑗}+4𝜙_𝑖,𝑗+1+4𝜙_{𝑖,𝑗−1}−20𝜙_𝑖,𝑗6ℎ²+𝑂ℎ⁴=0. (4.6) Note that (4.6) is derived under the condition of 3𝑎+2√12−21𝑏+2−4𝑐+√2−181𝑑+1−16𝑒=0. (4.7) The resulting five free parameters can then be uniquely determined from (4.3)–(17) and (4.7) as 𝑎=(480√2−705)/(512√2−745), 𝑏=0, 𝑐=(−960√2+1515)/(512√2−745), 𝑑=−210/(512√2−745), 𝑒=480(√2−600)/(512√2−745).

In summary, the kernel function developed for simulating the pure diffusion equation in this study is given as follows:𝑤(𝑟)=480√2−705512√12−745𝑟_𝑒+−960√2+1515512√12−745𝑟_𝑒𝑟𝑟_𝑒²+−210512√12−745𝑟_𝑒𝑟𝑟_𝑒³+480√2−600512√12−745𝑟_𝑒𝑟𝑟_𝑒⁴;0≤𝑟≤𝑟_𝑒,0;𝑟_𝑒<𝑟 (4.8) under the condition given below 𝑤√=12ℎ4𝑤(ℎ). (4.9) Note that the kernel function 𝑤(𝑟)=𝑟_𝑒/𝑟−1 given in [13] has an infinitely large magnitude as the value of 𝑟 approaches zero. The newly developed kernel function shown in Figure 1, on the other hand, has its local maximum value at 𝑟=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 ^{𝑟/𝑟𝑒=1}_{𝑟/𝑟𝑒=−1}𝑤(𝑟)𝑑𝑟=1. (4.10)

(a)

(b)

Figure 1 

Plot of the developed kernel function. (a) cast in 𝑟/𝑟𝑒 coordinate; (b) plotted in 𝑥𝑦 plane.

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 𝑑 and the channel width is 16𝑑. Above the undisturbed water level, the crest of the investigated solitary wave is given by [15]𝑦=𝑑+𝐻sech²3𝐻4𝑑³𝑥. (5.1)

Figure 2 

Schematic of the solitary wave propagation problem.

Here, 𝐻(=2) and 𝑑(=10) denote the initial wave height and the still water depth, respectively. Initially, the velocities are given by𝑢=𝐻𝑔𝑑𝑑sech²3𝐻4𝑑³𝑥,𝑣=𝐻3𝑔𝑑𝑑^3/2𝑦𝑑sech²3𝐻4𝑑³𝑥tanh3𝐻4𝑑³𝑥, (5.2) where 𝑔(=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.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4 

The predicted free surfaces and pressures at different times. (a) 𝑡=0.5; (b) 𝑡=2.0; (c) 𝑡=4.0; (d) 𝑡=6.0; (e) 𝑡=8.0; (f) 𝑡=10.0.

5.2. Dam Break Flow

The investigated dambreak problem is schematic Figure 5, where 𝐿(=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.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6 

Comparison of the predicted free surfaces by the MPS and the level set methods at different times for the case with Re=500. (a) 𝑡=0.0; (b) 𝑡=0.5; (c) 𝑡=1.0; (d) 𝑡=1.5; (e) 𝑡=2.0; (f) 𝑡=2.5.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 7 

Comparison of the predicted free surfaces by the MPS and the level set methods at different times for the case with Re=500. (a) 𝑡=0.0; (b) 𝑡=0.5; (c) 𝑡=1.0; (d) 𝑡=1.5; (e) 𝑡=2.0; (lf) 𝑡=2.5.

6. Concluding Remarks

The idea of developing the proposed particle interaction model is to develop the kernel function which accommodates the property ∫^{𝑟/𝑟𝑒=1}_{𝑟/𝑟𝑒=−1}𝑤(𝑟)𝑑𝑟=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.

Acknowledgments

The financial support by the National Science Council under Grants NSC96-2221-E-002-293-MY2 and NSC97-2628-M-002-022 and CQSE project 97R0066-69 are gratefully acknowledged.