A Particle Method Based on a Generalized Finite Difference Scheme to Solve Weakly Compressible Viscous Flow Problems

Abstract

The Lagrangian meshfree particle-based method has advantages in solving fluid dynamics problems with complex or time-evolving boundaries for a single phase or multiple phases. A pure Lagrangian meshfree particle method based on a generalized finite difference (GFD) scheme is proposed to simulate time-dependent weakly compressible viscous flow. The flow is described with Lagrangian particles, and the partial differential terms in the Navier-Stokes equations are represented as the solution of a symmetric system of linear equations through a GFD scheme. In solving the particle-based symmetric equations, the numerical method only needs the kernel function itself instead of using its gradient, i.e., the approach is a kernel gradient free (KGF) method, which avoids using artificial parameters in solving for the viscous term and reduces the limitations of using the kernel function. Moreover, the order of Taylor series expansion can be easily improved in the meshless algorithm. In this paper, the particle method is validated with several test cases, and the convergence, accuracy, and different kernel functions are evaluated.

1. Introduction

Problems of weakly compressible flows have attracted much attention in aerospace and oceanic applications, such as wind engineering problems, turbine flow, blood flow, and water wave motion. Accurate predictions of such flows are important in computational fluid dynamics. For fluids at low Mach numbers, the ratio between the speed of flow and the speed of sound is extremely small, and therefore, density fluctuations are not obvious. As a result, such a situation can be called weakly compressible flow. Generally, there are three numerical ways to model weakly compressible flow, namely, the Eulerian approach, the Lagrangian approach, and the hybrid approach. The Eulerian approach solves for quantities at fixed locations in space, and the Lagrangian approach uses individual particles that move through both space and time and have their own physical properties, such as density, velocity, and pressure, to represent the dynamically evolving fluid flow. The flow is described by recording the time history of each fluid particle. In the present work, we propose a pure Lagrangian meshfree particle-based method based on a meshless finite difference scheme to solve weakly compressible flow problems.

The Lagrangian meshfree method is rapidly advancing and has been widely used in recent years because it can be easily adapted to modeling problems with complex or time-evolving boundaries for single or multiple phases, such as numerical simulations of dam break flow [1], hydraulic jumps [2], rising bubbles [3] and coalescing [4]. The ability of this method to model non-Newtonian fluid and large scale diffuse fluids has been demonstrated in some recent works [5,6] by introducing different symmetric models. Moreover, because the computations are based on the support domain, which is much smaller than the complete computational region, the ill conditioned system problem is rarely encountered. Among all Lagrangian meshfree methods, the smoothed particle hydrodynamics (SPH) method was one of the earliest methods developed and has been widely applied in different fields. The SPH method was first pioneered independently by Lucy [7] and Gingold and Monaghan [8] to solve astrophysical problems in 1977. Details of the SPH method as a computational fluid dynamics method can be found in recent reviews [9,10,11,12] and Liu and Liu’s book [13]. Some successful applications of this method include coastal engineering, nuclear engineering, ocean engineering, and bioengineering. However, the accuracy of the conventional SPH method is unsatisfactory, and it is not easy to achieve an accurate high-order SPH approach.

As a meshfree Lagrangian method, the particle distribution generally tends to be irregular in the computations, which leads to inconsistency and low accuracy [14,15]. For that reason, in some cases, only the first-order term of the fluid dynamics equations, the Navier-Stokes equations, is solved, and the viscous term, which contains the second-order differential, is obtained through the artificial viscosity with artificial parameters in the SPH method. This issue can also occur in the incompressible SPH (ISPH) method [16]. To improve the consistency and accuracy of these methods, different modifications have been developed. After using Taylor series expansion to normalize the kernel function, the corrective smoothed particle method (CSPM) [17,18] and the modified smoothed particle method (MSPH) [19] were proposed. Both methods have better accuracy than the conventional SPH method. Nevertheless, it should be noted that both methods improve accuracy by improving the particle approximation of the kernel gradient term, which leads to more strict requirements on the kernel function. These requirements are related to the compact condition, normalization condition, and delta function behavior [20] and limit the selection of the kernel function, especially when the second-order gradient of the kernel function is required.

To avoid the solution of the gradient of the kernel function, a method with kernel gradient free (KGF) features can be developed, as discussed in detail in [21,22,23]; notably, a KGF-SPH method was proposed in 2015. When a particle method only involves the kernel function itself in kernel and particle approximation, the kernel gradient is not necessary in the computation, and this approach is thus referred to as a KGF method. The KGF-SPH method is used to solve for the viscous term directly without using the artificial viscosity, and the results are good for 2D models. Another KGF method is the consistent particle method (CPM) [24,25] for incompressible flow simulation. In the CPM, Poisson’s equation is used in the same way as the moving-particle semi-implicit (MPS) method based on the particle number density and the difference algorithm.

The purpose of our work is to combine a finite difference scheme and the particle method for solving weakly compressible viscous flow problems. In the method, the flow is described with Lagrangian particles, and the partial differential terms in the Navier-Stokes equations are represented as the solution of a symmetric system of linear equations through a generalized finite difference (GFD) scheme. It should be noted that this method is not a completely new method, but we will simply refer to it as the finite difference particle method (FDPM) to simplify the description in the subsequent sections.

Meshless finite difference approximation was first discussed for fully arbitrary meshes by Jensen [26] in 1972. Perrone and Kao [27] also contributed to the development of this method at that time. Subsequently, a variation using the moving least squares method was proposed by Lizska and Orkisz [28], and some recent works have been published [29,30,31]. The meshless finite difference scheme or GFD approximation we used came from Benito, Urena and Gavete [32], and they provided a discussion of the influence of several factors in the GFD scheme. A comparison between the GFD method and the element-free Galerkin method (EFGM) in solving the Laplace equation was presented in [33,34]. The GFD method was shown to be more accurate than the EFGM and the GFD scheme was used as a Eulerian meshfree method. In the present work, the GFD scheme is utilized to build a Lagrangian meshfree particle-based method, namely, the FDPM.

The FDPM has several advantages compared with the conventional SPH. First, the method is KGF. Only the kernel function itself and the positions of each particle are used to compute the spatial differential through a set of symmetric linear equations. Second, the method can be easily extended to high orders because it is based on Taylor series expansion. We show a fourth-order scheme for the FDPM. Additionally, only a few lines of code need to change to obtain a high-order FDPM, which is simple for users, especially when the users want to start with a low-order but fast computation. Third, the second-order differential term can be obtained without additional limitations on the kernel function. Thus, the viscous term in the Navier-Stokes equations can be computed directly without introducing any artificial parameters. Fourth, the FDPM is characterized by good compatibility. Most boundary conditions in the existing Lagrangian particle-based methods, such as the SPH and MPS methods, can be used directly. In the present work, we focus on the evaluation of the convergence, efficiency, and effects of the kernel function. The method is tested by modeling flow in a pipeline, Poisseuille flow, Couette flow and flow in porous media. These classical flows are used in different ways to solve fluid dynamics problems [35,36,37].

The present paper is organized as follows. In Section 2, the FDPM is given to solve the Navier-Stokes equations. In Section 3, applications of the particle method are shown. Section 4 summarizes the results of this work.

2. Finite Difference Particle Method for Weakly Compressible Flow

2.1. Lagrangian Form of the Governing Equations for Weakly Compressible Viscous Flow

The Lagrangian form of the Navier-Stokes equations, i.e., the continuity equation and the momentum equation, including viscous and external forces, are defined by Equations (1) and (2), respectively. The Lagrangian form of governing equations is as follows:

$$\frac{\mathrm{D}\rho }{\mathrm{D}t}=-\rho \nabla ·\mathit{u},$$

(1)

$$\frac{\mathrm{D}u}{\mathrm{D}t}=-\frac{1}{\rho }\nabla p+v_{\mathrm{k}}{\nabla }^2\mathit{u}+\mathit{F},$$

(2)

where ρ is the density, t is the time, u is the particle velocity, p is the pressure, v_k is the kinematic viscosity, and F is an external body force, such as gravity. All these variables are related to the physical properties of fluid particles that can move in both space and time, rather than remain at a fixed position.

The material derivative is written as follows

$$\frac{\mathrm{D}}{\mathrm{D}t}=\frac{\partial }{\partial t}+\mathit{u}·\nabla .$$

(3)

The equation of state for weakly compressible fluid flow is

$$\frac{\mathrm{D}p}{\mathrm{D}t}=c^2\frac{\mathrm{D}\rho }{\mathrm{D}t},$$

(4)

where c is the speed of sound.

2.2. Generalized Finite Difference Scheme

In the FDPM, flow is described with Lagrangian particles, and the GFD approximation [32] is utilized to solve for the spatial differential terms in the governing equations.

Consider a particle i surrounded by particles j = 1, 2, …, N, with all N + 1 of the particles in a compact support domain, as shown in Figure 1. For a circular support domain, r_s represents the radius of the support domain, which is called the smoothing length in the SPH method. Particles j are white circles around particle i, which is the orange circle, and Ω represents the computational domain. The closest nodes to particle i are selected as j particles, and these particles should be in the support domain at the same time.

The values of an infinitely differentiable function F at the positions of particles i and j are defined as F_i and F_j, respectively. This F can be the pressure, velocity or density of particles in the computation. Let us expand this term as the Taylor series of F for particle i:

$$\begin{array}{c}{F}_j=F_i+h_j\frac{\partial F_i}{\partial x}+k_j\frac{\partial F_i}{\partial y}+\frac{h_j^2}{2}\frac{{\partial }^2{F}_i}{\partial x^2}+\frac{k_j^2}{2}\frac{{\partial }^2{F}_i}{\partial y^2}+h_j{k}_j\frac{{\partial }^2{F}_i}{\partial x\partial y}+\frac{h_j^3}{6}\frac{{\partial }^3{F}_i}{\partial x^3}+\frac{k_j^3}{6}\frac{{\partial }^3{F}_i}{\partial y^3}+\frac{h_j^2{k}_j}{2}\frac{{\partial }^3{F}_i}{\partial x^2\partial y}+\frac{h_j{k}_j^2}{2}\frac{{\partial }^3{F}_i}{\partial x\partial y^2}\\ +\frac{h_j^4}{24}\frac{{\partial }^4{F}_i}{\partial x^4}+\frac{k_j^4}{24}\frac{{\partial }^4{F}_i}{\partial y^4}+\frac{h_j^3{k}_j}{6}\frac{{\partial }^4{F}_i}{\partial x^3\partial y}+\frac{h_j^2{k}_j^2}{4}\frac{{\partial }^4{F}_i}{\partial x^2\partial y^2}+\frac{h_j{k}_j^3}{6}\frac{{\partial }^4{F}_i}{\partial x\partial y^3}+\cdots ,\end{array}$$

(5)

where Equation (5) is for two dimensions. x and y are the spatial coordinates of the particles, and h_j = x_j - x_i, k_j = y_j − y_i.

For the fourth-order FDPM, ignoring the high-order terms, the approximation of F is denoted by f:

$$\begin{array}{c}{f}_j=f_i+h_j\frac{\partial f_i}{\partial x}+k_j\frac{\partial f_i}{\partial y}+\frac{h_j^2}{2}\frac{{\partial }^2{f}_i}{\partial x^2}+\frac{k_j^2}{2}\frac{{\partial }^2{f}_i}{\partial y^2}+h_j{k}_j\frac{{\partial }^2{f}_i}{\partial x\partial y}+\frac{h_j^3}{6}\frac{{\partial }^3{f}_i}{\partial x^3}+\frac{k_j^3}{6}\frac{{\partial }^3{f}_i}{\partial y^3}+\frac{h_j^2{k}_j}{2}\frac{{\partial }^3{f}_i}{\partial x^2\partial y}+\frac{h_j{k}_j^2}{2}\frac{{\partial }^3{f}_i}{\partial x\partial y^2}\\ +\frac{h_j^4}{24}\frac{{\partial }^4{f}_i}{\partial x^4}+\frac{k_j^4}{24}\frac{{\partial }^4{f}_i}{\partial y^4}+\frac{h_j^3{k}_j}{6}\frac{{\partial }^4{f}_i}{\partial x^3\partial y}+\frac{h_j^2{k}_j^2}{4}\frac{{\partial }^4{f}_i}{\partial x^2\partial y^2}+\frac{h_j{k}_j^3}{6}\frac{{\partial }^4{f}_i}{\partial x\partial y^3}.\end{array}$$

(6)

After rearranging these equations and multiplying by a kernel function W on both sides of the equation, the sum of these expressions for all particles j is obtained:

$${\sum }_{j=1}^N\left(\begin{array}{c}{f}_i-f_j+h_j\frac{\partial f_i}{\partial x}+k_j\frac{\partial f_i}{\partial y}+\frac{h_j^2}{2}\frac{{\partial }^2{f}_i}{\partial x^2}+\frac{k_j^2}{2}\frac{{\partial }^2{f}_i}{\partial y^2}+h_j{k}_j\frac{{\partial }^2{f}_i}{\partial x\partial y}\\ +\frac{h_j^3}{6}\frac{{\partial }^3{f}_i}{\partial x^3}+\frac{k_j^3}{6}\frac{{\partial }^3{f}_i}{\partial y^3}+\frac{h_j^2{k}_j}{2}\frac{{\partial }^3{f}_i}{\partial x^2\partial y}+\frac{h_j{k}_j^2}{2}\frac{{\partial }^3{f}_i}{\partial x\partial y^2}\\ +\frac{h_j^4}{24}\frac{{\partial }^4{f}_i}{\partial x^4}+\frac{k_j^4}{24}\frac{{\partial }^4{f}_i}{\partial y^4}+\frac{h_j^3{k}_j}{6}\frac{{\partial }^4{f}_i}{\partial x^3\partial y}+\frac{h_j^2{k}_j^2}{4}\frac{{\partial }^4{f}_i}{\partial x^2\partial y^2}+\frac{h_j{k}_j^3}{6}\frac{{\partial }^4{f}_i}{\partial x\partial y^3}\end{array}\right)W\left(h_j, k_j, r_{\mathrm{s}}\right)=0,$$

(7)

where W(h_j, k_j, r_s) is a kernel function in 2D and r_s represents the size of the support domain. For W, different kernel functions, including Gaussian, cubic spline, and quintic spline functions, can be found in [13]. In the following equations, we use W for the kernel function.

Function G can be defined in 2D as

$$G(f)={\sum }_{j=1}^N{\left[\right(\begin{array}{c}{f}_i-f_j+h_j\frac{\partial f_i}{\partial x}+k_j\frac{\partial f_i}{\partial y}+\frac{h_j^2}{2}\frac{{\partial }^2{f}_i}{\partial x^2}+\frac{k_j^2}{2}\frac{{\partial }^2{f}_i}{\partial y^2}+h_j{k}_j\frac{{\partial }^2{f}_i}{\partial x\partial y}\\ +\frac{h_j^3}{6}\frac{{\partial }^3{f}_i}{\partial x^3}+\frac{k_j^3}{6}\frac{{\partial }^3{f}_i}{\partial y^3}+\frac{h_j^2{k}_j}{2}\frac{{\partial }^3{f}_i}{\partial x^2\partial y}+\frac{h_j{k}_j^2}{2}\frac{{\partial }^3{f}_i}{\partial x\partial y^2}\\ +\frac{h_j^4}{24}\frac{{\partial }^4{f}_i}{\partial x^4}+\frac{k_j^4}{24}\frac{{\partial }^4{f}_i}{\partial y^4}+\frac{h_j^3{k}_j}{6}\frac{{\partial }^4{f}_i}{\partial x^3\partial y}+\frac{h_j^2{k}_j^2}{4}\frac{{\partial }^4{f}_i}{\partial x^2\partial y^2}+\frac{h_j{k}_j^3}{6}\frac{{\partial }^4{f}_i}{\partial x\partial y^3}\end{array}\left)W\right]}^2=0.$$

(8)

According to Equation (7), the norm of G equals 0, so we obtain:

$$\begin{array}{c}\frac{\partial G\left(f\right)}{\partial \left(\frac{\partial f_i}{\partial x}\right)}=2{\sum }_{j=1}^N{\Phi h}_j{W}^2=0,\frac{\partial G\left(f\right)}{\partial \left(\frac{\partial f_i}{\partial y}\right)}=2{\sum }_{j=1}^N{\Phi k}_j{W}^2=0,\frac{\partial G\left(f\right)}{\partial \left(\frac{{\partial }^2{f}_i}{\partial x^2}\right)}={\sum }_{j=1}^N{\Phi h}_j^2{W}^2=0,\\ \frac{\partial G\left(f\right)}{\partial \left(\frac{{\partial }^2{f}_i}{\partial y^2}\right)}={\sum }_{j=1}^N{\Phi k}_j^2{W}^2=0,\frac{\partial G\left(f\right)}{\partial \left(\frac{{\partial }^2{f}_i}{\partial x\partial y}\right)}=2{\sum }_{j=1}^N{\Phi h}_j{k}_j{W}^2=0,\frac{\partial G\left(f\right)}{\partial \left(\frac{{\partial }^3{f}_i}{\partial x^3}\right)}=\frac{1}{3}{\sum }_{j=1}^N{\Phi k}_j^3{W}^2=0,\\ \frac{\partial G\left(f\right)}{\partial \left(\frac{{\partial }^3{f}_i}{\partial y^3}\right)}=\frac{1}{3}{\sum }_{j=1}^N{\Phi k}_j^3{W}^2=0,\frac{\partial G\left(f\right)}{\partial \left(\frac{{\partial }^3{f}_i}{\partial x^2\partial y}\right)}={\sum }_{j=1}^N{\Phi h}_j^2{k}_j{W}^2=0,\frac{\partial G\left(f\right)}{\partial \left(\frac{{\partial }^3{f}_i}{\partial x\partial y^2}\right)}={\sum }_{j=1}^N{\Phi h}_j{k}_j^2{W}^2=0,\\ \frac{\partial G\left(f\right)}{\partial \left(\frac{{\partial }^4{f}_i}{\partial x^4}\right)}=\frac{1}{12}{\sum }_{j=1}^N{\Phi h}_j^4{W}^2=0,\frac{\partial G\left(f\right)}{\partial \left(\frac{{\partial }^4{f}_i}{\partial y^4}\right)}=\frac{1}{12}{\sum }_{j=1}^N{\Phi k}_j^4{W}^2=0,\frac{\partial G\left(f\right)}{\partial \left(\frac{{\partial }^4{f}_i}{\partial x^3\partial y}\right)}=\frac{1}{3}{\sum }_{j=1}^N{\Phi h}_j^3{k}_j{W}^2=0,\\ \frac{\partial G\left(f\right)}{\partial \left(\frac{{\partial }^4{f}_i}{\partial x^2\partial y^2}\right)}=\frac{1}{2}{\sum }_{j=1}^N{\Phi h}_j^2{k}_j^2{W}^2=0,\frac{\partial G\left(f\right)}{\partial \left(\frac{{\partial }^4{f}_i}{\partial x\partial y^3}\right)}=\frac{1}{3}{\sum }_{j=1}^N{\Phi h}_j{k}_j^3{W}^2=0,\end{array}$$

(9)

where

$$\begin{array}{c}\Phi =f_i-f_j+h_j\frac{\partial f_i}{\partial x}+k_j\frac{\partial f_i}{\partial y}+\frac{h_j^2}{2}\frac{{\partial }^2{f}_i}{\partial x^2}+\frac{k_j^2}{2}\frac{{\partial }^2{f}_i}{\partial y^2}+h_j{k}_j\frac{{\partial }^2{f}_i}{\partial x\partial y}+\frac{h_j^3}{6}\frac{{\partial }^3{f}_i}{\partial x^3}+\frac{k_j^3}{6}\frac{{\partial }^3{f}_i}{\partial y^3}+\frac{h_j^2{k}_j}{2}\frac{{\partial }^3{f}_i}{\partial x^2\partial y}+\frac{h_j{k}_j^2}{2}\frac{{\partial }^3{f}_i}{\partial x\partial y^2}\\ +\frac{h_j^4}{24}\frac{{\partial }^4{f}_i}{\partial x^4}+\frac{k_j^4}{24}\frac{{\partial }^4{f}_i}{\partial y^4}+\frac{h_j^3{k}_j}{6}\frac{{\partial }^4{f}_i}{\partial x^3\partial y}+\frac{h_j^2{k}_j^2}{4}\frac{{\partial }^4{f}_i}{\partial x^2\partial y^2}+\frac{h_j{k}_j^3}{6}\frac{{\partial }^4{f}_i}{\partial x\partial y^3}.\end{array}$$

(10)

Equation (9) gives us the following equation:

AD = B,

(11)

where

$$\mathbf{A}=\left[\begin{array}{cccccc}{\sum }_{j=1}^N{h}_j^2{W}^2& {\sum }_{j=1}^N{h}_j{k}_j{W}^2& {\sum }_{j=1}^N\frac{1}{2}{h}_j^3{W}^2& {\sum }_{j=1}^N\frac{1}{2}{h}_j{k}_j^2{W}^2& \cdots & {\sum }_{j=1}^N\frac{1}{6}{h}_j^2{k}_j^3{W}^2\\ & {\sum }_{j=1}^N{k}_j^2{W}^2& {\sum }_{j=1}^N\frac{1}{2}{h}_j^2{k}_j{W}^2& {\sum }_{j=1}^N\frac{1}{2}{k}_j^3{W}^2& \cdots & {\sum }_{j=1}^N\frac{1}{6}{h}_j{k}_j^4{W}^2\\ & & {\sum }_{j=1}^N\frac{1}{4}{h}_j^4{W}^2& {\sum }_{j=1}^N\frac{1}{4}{h}_j^2{k}_j^2{W}^2& \cdots & {\sum }_{j=1}^N\frac{1}{12}{h}_j^3{k}_j^3{W}^2\\ & & & {\sum }_{j=1}^N\frac{1}{4}{k}_j^4{W}^2& \cdots & {\sum }_{j=1}^N\frac{1}{12}{h}_j{k}_j^5{W}^2\\ & & & & \ddots & ⋮\\ \mathrm{symmetric}& & & & & {\sum }_{j=1}^N\frac{1}{36}{h}_j^2{k}_j^6{W}^2\end{array}\right],$$

(12)

$$\mathbf{D}={\left\{\begin{array}{cccccc}\frac{\partial f_i}{\partial x}& \frac{\partial f_i}{\partial y}& \frac{{\partial }^2{f}_i}{\partial x^2}& \frac{{\partial }^2{f}_i}{\partial y^2}& \cdots & \frac{{\partial }^4{f}_i}{\partial x\partial y^3}\end{array}\right\}}^{\mathrm{T}},$$

(13)

$$\mathbf{B}=\left\{\begin{array}{c}{\sum }_{j=1}^N{(f}_j{-f}_i{)h}_j{W}^2\\ {\sum }_{j=1}^N{(f}_j{-f}_i{)k}_j{W}^2\\ {\sum }_{j=1}^N\frac{1}{2}{(f}_j{-f}_i{)h}_j^2{W}^2\\ {\sum }_{j=1}^N\frac{1}{2}{(f}_j{-f}_i{)k}_j^2{W}^2\\ ⋮\\ {\sum }_{j=1}^N\frac{1}{6}{(f}_j{-f}_i{)h}_j^2{k}_j^3{W}^2\end{array}\right\}.$$

(14)

Since matrix A is symmetrical, Equation (11) can be solved, and the solution gives the values of the spatial derivatives in matrix D. Thus, the spatial derivatives in Equations (1) and (2) can be obtained by solving a set of symmetric linear equations, and the material derivatives in the equations can be integrated using a time integration scheme.

2.3. Particle Representation for Governing Equations

Taking particle i as an example, this section gives the particle representation of the governing equations and the solution to Equation (11). The solution includes the values of the spatial derivatives needed in the governing equations.

The coefficients of D and B are denoted by D_m(f_i) and B_m(f_i), respectively, with m = 1, 2, …, 5. For example, $D_2\left(f_i\right)=\frac{\partial f_i}{\partial y}$ (the second coefficient in Equation(13)), and $B_2\left(f_i\right)={\sum }_{j=1}^N{(f}_j{-f}_i{)k}_j{W}^2$ (the second coefficient in Equation (14)). In addition, the symmetric matrix A can be decomposed into the upper and lower triangular matrices A = LL^T. The coefficients of the matrix L are denoted by L(m, n), with m and n = 1, 2, 3, 4, 5.

By using the GFD scheme and Cholesky factorization to solve Equation (11), we obtain the solutions for the Lagrangian derivative terms in Equations (1) and (2) in two-dimensional form

$$\frac{{\mathrm{D}\rho }_i}{\mathrm{D}t}=-{\rho }_i{D}_1\left(u_i\right)-{\rho }_i{D}_2\left(v_i\right),$$

(15)

$$\frac{{\mathrm{D}u}_{\mathrm{i}}}{\mathrm{D}t}=-\frac{1}{{\rho }_i}{D}_1\left(p_i\right)+v_{\mathrm{k}}{D}_3\left(u_i\right),$$

(16)

$$\frac{{\mathrm{D}v}_i}{\mathrm{D}t}=-\frac{1}{{\rho }_i}{D}_2\left(p_i\right)+v_{\mathrm{k}}{D}_4\left(v_i\right)+\mathit{g},$$

(17)

where u_i and v_i are the velocity of particle i in two directions and

$$D_m\left(f\right)=\{\begin{array}{cc}\frac{1}{L(m,m)}\left[Y_m\left(f\right)-{\sum }_{n=m+1}^N{L(n,m)D}_n\left(f\right)\right]& m=1,2,3,4\\ \frac{Y\left(f\right)}{L(m,m)}& m=5\end{array},$$

(18)

where

$$Y_m\left(f\right)=\{\begin{array}{cc}\frac{b_m\left(f\right)}{L(m,m)}& m=1\\ \frac{1}{L(m,m)}\left[b_m\left(f\right)-{\sum }_{n=1}^{m-1}{L(m,n)Y}_n\left(f\right)\right]& m=2,3,4,5\end{array}.$$

(19)

This method considers changes in density and is able to simulate flow at low Mach numbers, so it is used to solve weakly compressible viscous flow problems. Calculations of particle motion and time integration are performed based on second-order leapfrog integration. The equations for updating the position and velocity of particles are

$${\mathit{v}}_i\left(t+\frac{1}{2}\Delta t\right)={\mathit{v}}_i\left(t-\frac{1}{2}\Delta t\right)+\Delta t\frac{{\mathrm{D}v}_i\left(t\right)}{\mathrm{D}t},$$

(20)

$${\mathit{r}}_i\left(t+\Delta t\right)={\mathit{r}}_i\left(t\right)+\Delta {t\mathit{v}}_i(t+\frac{1}{2}\Delta t),$$

(21)

where ${\mathit{v}}_i(t+\frac{1}{2}\Delta t)$ is the velocity of fluid particle i at time $t+\frac{1}{2}\Delta t$, and Δt is the time step.

2.4. Artificial Particle Displacement

In simulations of flow in porous media (Section 3.5), artificial particle displacement is suggested as a particle motion correction to avoid particles in the vicinity of the stagnation points of fluid flow [38] and to avoid poor particle distributions [39]. Artificial particle displacement can be expressed as

$${\delta \mathit{r}}_i=\alpha {\overline{{\mathrm{r}}_i}}^2{\mathit{v}}_{\max }\Delta t{\sum }_{j=1}^N\frac{{\mathit{r}}_{ij}}{r_{ij}^3},$$

(22)

where r_i is the position of particle i, α is a problem-dependent parameter that is usually set between 0.01 and 0.1, v_max is the maximum velocity of all particles in the computational domain, r_ij = r_i − r_j which is the distance between particles i and j, and $\overline{r_i}$ is the average distance between the neighboring particles of particle i:

$$\overline{r_i}=\frac{1}{N}{\sum }_{j=1}^N{r}_{ij}.$$

(23)

It is noted that the problem-dependent parameter α should be selected carefully. This value should be small enough not to affect the physics of the flow but also large enough to avoid the accumulation of particles to form groups. In the present work, the value of artificial particle displacement is less than 0.1% of the physical particle displacement for a given time step, which is consistent with the magnitude in [40].

After moving the particles, the pressure and velocity components should be corrected by Taylor expansion.

2.5. Boundary Conditions

Several layers of virtual particles are used to implement the boundary condition. Similar treatments can be observed in the SPH and MPS simulations. On a flat wall, virtual particles are obtained by extending the boundary particles to the outside of the computational region, and the distribution of virtual particles is regular. The number of layers can be chosen according to the scale of the support domain. Figure 2 is a sketch of the treatment of particles near the wall. i represents the particle number, and Δx is the particle spacing.

For a flat wall, both the no-slip and free-slip boundary conditions can be implemented using virtual particles. For no-slip walls, the particle-based boundary conditions are as follows

$$p_{i+3}=p_{i+2}=p_{i+1}=p_i, {\mathit{v}}_{i+3}={\mathit{v}}_{i+2}={\mathit{v}}_{i+1}=0,$$

(24)

For free-slip walls, the tangential velocity component of virtual particles is maintained the same as the boundary particles.

For a round surface, virtual particles are established based on a radial distribution inside the object domain with particle spacing Δx. The particle distribution is shown in Figure 3.

The boundary condition for a rigid wall satisfies the following equations.

$$p_{i+2}=p_{i+1}=p_i, {\mathit{v}}_{i+3}={\mathit{v}}_{i+2}={\mathit{v}}_{i+1}=0.$$

(25)

Since the FDPM simulation is still based on the local support domain, most boundary conditions in the existing Lagrangian particle-based methods, such as the SPH and MPS methods, can be used directly or implemented with minor changes. The particle representation of no-slip, free-slip and superhydrophobic surfaces [41,42,43] can be found in [44,45,46].

3. Applications of the Finite Difference Particle Method

3.1. Fundamental Definition

Several test cases are simulated with the FDPM based on second-order Taylor series expansion. The numerical accuracy is evaluated by the root mean square errors (ε_RMS) and the maximum errors (ε_MAX), which are defined as

$${\epsilon }_{\mathrm{RMS}}\left(S\right)=\frac{\sqrt{\frac{1}{NT}{\sum }_{\mathrm{k}=1}^{NT}{\left|S_{\mathrm{num}}\left(k\right){-S}_{\mathrm{ana}}\left(k\right)\right|}^2}}{\sqrt{\frac{1}{NT}{\sum }_{k=1}^{NT}{\left|S_{\mathrm{ana}}\left(k\right)\right|}^2}},$$

(26)

$${\epsilon }_{\mathrm{MAX}}\left(S\right)=\underset{1\ll k\ll NT}{\max }\left|S_{\mathrm{num}}\left(k\right){-S}_{\mathrm{ana}}\left(k\right)\right|,$$

(27)

where S_num(k) and S_ana(k) are the numerical and analytical results of variable k, respectively. k could be the velocity or pressure.

The convergence rate of the FDPM is evaluated based on the root mean square error convergence rate (R_ERMS) and the maximum error convergence rate (R_EMAX) as follows:

$$R_{\mathrm{ERMS}}=\left|\frac{\ln \left({\epsilon }_{\mathrm{RMS}}\left({NT}_{\max }\right)\right)-\ln \left({\epsilon }_{\mathrm{RMS}}\left({NT}_{\min }\right)\right)}{\ln \left({NT}_{\max }\right)-\ln \left({NT}_{\min }\right)}\right|,$$

(28)

$$R_{\mathrm{EMAX}}=\left|\frac{\ln \left({\epsilon }_{\mathrm{MAX}}\left({NT}_{\max }\right)\right)-\ln \left({\epsilon }_{\mathrm{MAX}}\left({NT}_{\min }\right)\right)}{\ln \left({NT}_{\max }\right)-\ln \left({NT}_{\min }\right)}\right|.$$

(29)

3.2. Unsteady Flow in a Pipeline

Unsteady flow field in the pipeline is simulated to verify the FDPM method. The theoretical solutions of unsteady flow in a 1D pipeline (chapter 3 in book [47]) are

$$u=\frac{2}{\gamma +1}\frac{x}{t}+C_1,$$

(30)

$$p=C_3{\left(\frac{\gamma -1}{C_3\gamma }\left(-\frac{1}{2}{\left(C_1+\frac{2x}{\left(\gamma +1\right)t}\right)}^2+\frac{C_1^2\left(\gamma +1\right)}{2\left(\gamma -1\right)}+\frac{C_2\left(\gamma -3\right)}{\gamma +1}{t}^{\frac{2-2\gamma }{\gamma +1}}+\frac{x^2}{{(\gamma +1)t}^2}\right)\right)}^{\gamma /(\gamma -1)},$$

(31)

$$\rho ={\left(\frac{\gamma -1}{C_3\gamma }\left(-\frac{1}{2}{\left(C_1+\frac{2x}{\left(\gamma +1\right)t}\right)}^2+\frac{C_1^2\left(\gamma +1\right)}{2\left(\gamma -1\right)}+\frac{C_2\left(\gamma -3\right)}{\gamma +1}{t}^{\frac{2-2\gamma }{\gamma +1}}+\frac{x^2}{{(\gamma +1)t}^2}\right)\right)}^{1/(\gamma -1)},$$

(32)

$$c^2={\left(\frac{\gamma -1}{\gamma +1}{-C}_1\right)}^2-\frac{3-\gamma }{\gamma +1}{(\gamma -1)C}_2{t}^{-2(\gamma -1)(\gamma +1)},$$

(33)

where u₀ is the flow velocity distribution and x is the coordinate along the length of the pipeline. The coefficients (the unit can be obtained from dimensional analysis) C₁ = 30.0, C₂ = –1.0 × 10⁶, C₃ = 82571.0, and γ = 1.4; moreover, the initial time is 12.5 s, and the pipe length x is 700 m.

The FDPM algorithm with a second-order Taylor truncation is used, and the time step (Δt) is 0.0029 s. The effect of viscosity in the process of fluid motion is not considered. Dirichlet boundary conditions are used at both ends of the boundary. The velocity, pressure and density of four particles at both ends are set based on theoretical values.

The space of the initial particle (Δx) is 5.0 m, and r_s is 3.2 times Δx. The cubic spline kernel function is used in the calculations.

Table 1. FDPM results and theoretical solutions of the position and velocity of a particle (particle number: 70) at different times.

Time (s)	Particle Method: x (m)	Theoretical Solution: u (m/s)	Particle Method: u (m/s)	Error (10−8)
0.25	2.48	30.16201462	30.16201544	2.74
0.50	10.08	30.64615858	30.64615942	2.75
0.75	17.80	31.11956025	31.11956110	2.75
1.00	25.64	31.58265402	31.58265470	2.13

provides data for comparing the numerical velocity and theoretical solution of a particular particle at different times and positions. Notably, although the particle moves from x = 2.48 m to x = 25.64 m, the FDPM results agree well with the theoretical values, and this result verifies the algorithm.

The convergence verification of the FDPM method for unsteady flow in a pipeline is shown in Figure 4. The numerical error curves at three different moments are given using different Δx values, and r_s is 3.2 times Δx in the computation.

Figure 4 shows that the FDPM method displays good convergence at different times. When t = 1.0 s, the convergence curve yields a R_ERMS value of 1.7 and R_EMAX value of 1.8.

Given that the FDPM is a KGF method, the effect of the type of kernel function on this method is evaluated. Four types of kernel functions, including $\frac{1}{r^3}$, Gaussian, cubic spline and quintic spline functions, are compared through two types of errors with different r_s conditions, as shown in Figure 5. The figure shows that the maximum error of the Gaussian kernel function is larger than that of the other methods. The errors of other types of functions are similar.

3.3. Poisseuille Flow

Steady, axisymmetric Poisseuille flow between two infinite plates is a classical test model in hydrodynamics. In this section, the model is used to verify the governing equations and the rigid wall boundaries. Assuming that the distance between two infinite plates is L, the volume force F is loaded on the fluid between plates in the x direction from time t = 0. The theoretical solution of the velocity distribution of the flow at a given time [48] is as follows.

$$u\left(y,t\right)=\frac{F}{v_{\mathrm{k}}}y\left(y-L\right)+{\sum }_{n=0}^{\infty }\frac{{4FL}^2}{v_{\mathrm{k}}{\mathsf{\pi }}^3{(2n+1)}^3}\sin \left[\frac{\mathsf{\pi }y}{L}(2n+1)\right]\exp \left[-\frac{{(2n+1)}^2{\mathsf{\pi }}^2{v}_{\mathrm{k}}}{L^2}t\right].$$

(34)

The numerical simulation for Poisseuille flow is obtained under weakly compressible (Ma = 0.0125) conditions. Based on reference [48], the parameters of the Poisseuille field are chosen as υ_k = 10⁻⁶ m²s⁻¹, L = 10⁻³ m, ρ = 10³ kgm⁻³, and F = 10⁻⁴ ms⁻², so the maximum velocity is 1.25 × 10⁻⁵ ms⁻¹ and the Reynold number is Re = 1.25 × 10⁻². The plate boundaries at the upper and lower ends are established using rigid walls. One layer of boundary particles and three layers of virtual particles are used. The FDPM with second-order Taylor truncation is utilized to perform the computation. The speed of sound c is taken as 0.001 m/s, as suggested in [49], and the time step Δt is 3.0 × 10⁻⁴ s. The initial particle spacings Δx and Δy are both set as 5 × 10⁻⁵ m, r_s is 3.2 times Δx, and the kernel function is selected as a cubic spline function.

A comparison of the numerical and theoretical solutions of the velocity of the flow field in the x direction at different times is shown in Figure 6. The particle velocity is obtained by bilinear interpolation. As time increases, the positions of particles gradually change until the uniformly distributed particles at the initial stage are completely mixed in disorder. At this time, the numerical solution of the particle velocity is still consistent with the theoretical solution.

During the computation, the particle velocity remains symmetrically distributed and gradually increases at different times before reaching a steady state. The velocity of particles in the middle of the two plates is the largest due to the viscous force, and the velocity is small near the plates. The FDPM solution is in good agreement with the theoretical solution.

Figure 7 shows the numerical error curves of the FDPM method at different times and is used to analyze the convergence of the FDPM method. During the computation, r_s remains 3.2 times Δx.

Figure 7 shows that ε_RMS and ε_MAX decrease as Δx decreases, which indicates that the numerical accuracy converges with the initial particle spacing at different times. ε_RMS is on the order of 10⁻³, indicating that the computational results agree well with the theoretical solutions. For different error evaluation indexes, the R_ERMS and R_EMAX values of the FDPM method are approximately 1.7 and 1.8, respectively, with good convergence at t = 1.0 s. Since the second-order Taylor expansion-based FDPM is implemented in the test, the convergence rate is reasonable.

An error analysis of the four different types of kernel functions is conducted to analyze the sensitivity of the FDPM method to the kernel function, as shown in Figure 8.

Figure 8 shows that different types of kernel functions can be used in the FDPM method and that the differences in the calculation errors are insignificant. The calculation errors of the four types of kernel functions from large to small exhibit the following order: Gaussian, r⁻³, cubic spline, and quantic spline.

3.4. Couette Flow

Couette flow considers the fluid flow between a stationary plate and a sliding plate. To accurately solve the flow distribution, the viscous term and boundary flow must be solved correctly. Initially, the two plates and the fluid between them remain stationary. At a constant speed, the upper plate begins to slide parallel to the lower plate. Assuming that the plate spacing is L and the sliding velocity is u₀, the theoretical solution of the flow velocity over time in the direction perpendicular to the plate [48] is as follows:

$$u\left(y,t\right)=\frac{u_0}{L}y+{\sum }_{n=1}^{\infty }\frac{{2u}_0}{n\mathsf{\pi }}{(-1)}^n\sin (\frac{n\mathsf{\pi }}{L}y\left)\exp \right(-\frac{n^2{\mathsf{\pi }}^2{v}_{\mathrm{k}}}{L^2}t).$$

(35)

Couette flow is numerically simulated under weakly compressible (Ma = 0.0125) conditions. The parameters for Couette flow are v_k = 10⁶ m²s⁻¹, L = 10⁻³ m, ρ = 10³ kgm⁻³, and u₀ = 1.25 × 10⁻⁵ m/s.

The plate boundaries at the upper and lower ends are obtained using rigid walls, and the upper plate is set with a constant velocity u₀. One layer of boundary particles and three layers of virtual particles are used. The FDPM with second-order Taylor truncation is utilized to perform the computation. The speed of sound c is taken as 0.001 m/s, and the time step Δt is 5.0 × 10⁻⁵ s. The initial particle spacings Δx and Δy are both set as 2.5 × 10⁻⁵ m, r_s is 3.2 times Δx, and the kernel function selected is the cubic spline function. Figure 9 shows a comparison between the FDPM method and the theoretical solution for the flow velocity at different times.

Before reaching a steady state, the particle velocity near the upper plate rapidly increases due to the viscous force, and that near the lower plate increases in a relatively slow manner. The velocity distribution of the particles between the two plates is nonlinear.

The velocity error (ε_RMS and ε_MAX) at different times and at different Δx values is used to evaluate the convergence of the FDPM, as shown in Figure 10.

From the numerical results, both error indexes gradually decrease with decreasing Δx, which suggests that the numerical accuracy converges with the initial particle spacing at different times. ε_RMS is on the order of 10⁻², indicating that the computational results agree well with the theoretical solution. When t = 1.0 s, the two errors result in an R_ERMS value of 1.7 and R_EMAX value of 1.8. Since the second-order Taylor expansion-based FDPM is implemented in the test, the convergence rate is reasonable.

To analyze the sensitivity of the FDPM method to the kernel function, four different types of kernel functions with different r_s/Δx values are applied, as shown in Figure 11. Different types of kernel functions can be used in the FDPM method, and the calculation error of the Gaussian kernel function is larger than that of the other methods.

3.5. Flow in Porous Media

In this section, the FDPM algorithm is used to simulate the flow in a simplified model of porous media [50]. The simplified model can be seen as flow around a circular cylinder, as shown in Figure 12, and four sides of the domain are periodic boundaries.

The size of the computational domain L = 0.1 m, the kinematic viscosity v_k = 10⁻⁶ m²s⁻¹, the cylindrical radius R = 2 × 10⁻² m, the volume force F₀ = 1.5 × 10⁻⁷ ms⁻², and the speed of sound c = 5.77 × 10⁻⁴ ms⁻¹. Δx and Δy are 0.003 m, r_s is 3.2 times Δx, Δt = 1.04 s with 2000 steps, and the coefficient of artificial particle displacement is 0.05. A rigid wall boundary is used for the cylindrical boundary, and a periodic boundary is used on the four sides of the computational domain. One layer of boundary particles and three layers of virtual particles are used. The FDPM with second-order Taylor truncation is utilized to perform the computations. The particle distribution and velocity contours at the initial time and the final steady state are shown in Figure 13.

At the initial time, the particle distribution is regular. Then, the fluid begins to flow, and particles are gradually scattered and evenly distributed in the computational domain.

The velocity distributions along lines 1 and 2 (dotted-dashed lines in Figure 12) are shown in Figure 14. Both the FDPM results and finite element method (FEM) results are given to evaluate the accuracy of the numerical method. The FEM results come from the data of figure 6 in [48].

When Δx > 0.004, the computation is not sufficiently stable and can easily collapse, so this condition is not shown in Figure 14. When Δx = 0.0038 m, the FDPM calculation results and the FEM results at x = 0 and 0.1 m produce significant differences. When Δx = 0.002 m, the FDPM results are similar to the FEM results. After convergence is obtained, the results of the FDPM method are consistent with the FEM results, which verifies the correctness of the numerical method and the boundary conditions.

A comparison of the FDPM results (Δx = 0.002 m) with different kernel functions and the FEM reference results is shown in Figure 15. The figure shows that the FDPM results with different types of kernel functions are comparable to the reference results and exhibit only minor differences.

3.6. Lid-Driven Cavity Flow

The lid-driven cavity flow is widely used as a benchmark test case and the model in [51] is used to verify the method. The case is in a square cavity with a sliding plate on the upper side and three fixed rigid walls around, as shown in Figure 16. Initially, the fluid in the cavity remain stationary. At a constant speed, the upper plate begins to slide horizontally and the simulation is at Re = 100.

Rigid wall boundary condition is used for four sides of the cavity, and the upper plate is set with a constant velocity u₀. The size of the cavity L = 1.0 m, the kinematic viscosity v_k = 0.01 m²s⁻¹, the sliding velocity u₀ = 1.0 m/s, and the speed of sound c = 10.0 ms⁻¹. Δx and Δy are 0.025 m, r_s is 2.7 times Δx, Δt = 0.001 s with 3000 steps, the coefficient of artificial particle displacement is 0.05, and the kernel function selected is the cubic spline function. One layer of boundary particles and three layers of virtual particles are used. The FDPM with second-order Taylor truncation is utilized to perform the computations. The particle distribution and velocity contours at the initial time and the final steady state are shown in Figure 17.

At the initial time, the particle distribution is regular. Then, the fluid begins to flow, and particles are gradually scattered and evenly distributed in the computational domain.

Horizontal velocity component profiles along horizontal and vertical geometric centerlines at t = 3.0 s, respectively, are shown in Figure 18. Although the present FDPM computation employed only 2/5 particles in the work [51], these profiles are in good agreement with the reference results. The particle distribution shows the method works well in geometries with corners.

4. Conclusions

In this paper, a particle method based on the GFD scheme is proposed to simulate weakly compressible viscous flow. This approach represents the partial differential terms in the Navier-Stokes equations as the solution of a symmetric system of linear equations. The convergence and accuracy of the symmetric particle-based method are tested by modeling flow in a pipeline, Poisseuille flow, Couette flow, flow in porous media, and lid-driven cavity flow. The numerical results exhibit close agreement with the theoretical solutions and finite element results. The particle method utilizes the kernel function itself instead of its gradient, which avoids using artificial parameters to solve for the viscous term and reduces the limitations on the choice of kernel function. Moreover, the order of the Taylor series expansion can easily be improved in the meshless algorithm. The convergence rate of the particle-based calculations with second-order Taylor truncation is approximately 1.7 in the tests, and four different kernel functions are tested and determined to be reliable.