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

: The Lagrangian meshfree particle-based method has advantages in solving ﬂuid 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 ﬁnite di ﬀ erence (GFD) scheme is proposed to simulate time-dependent weakly compressible viscous ﬂow. The ﬂow is described with Lagrangian particles, and the partial di ﬀ erential 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 artiﬁcial 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 di ﬀ erent kernel functions are evaluated.


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],

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: 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 The equation of state for weakly compressible fluid flow is where c is the speed of sound.

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 Symmetry 2019, 11, 1086 4 of 21 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 Fi and Fj, 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: where Equation (5) is for two dimensions. x and y are the spatial coordinates of the particles, and hj = xj -xi, kj = yj -yi.
For the fourth-order FDPM, ignoring the high-order terms, the approximation of F is denoted by f : 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: where W(hj, kj, rs) is a kernel function in 2D and rs 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 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: 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 : 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: 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 According to Equation (7), the norm of G equals 0, so we obtain: where (10) Equation (9) gives us the following equation: where 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.

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, (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 where u i and v i are the velocity of particle i in two directions and where 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 v i t+ where v i (t+ 1 2 ∆ t) is the velocity of fluid particle i at time t+ 1 2 ∆t, and ∆t is the time step.

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 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 r i is the average distance between the neighboring particles of particle i: 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.

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

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 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. 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 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.
virtual particles. For no-slip walls, the particle-based boundary conditions are as follows 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.
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].

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 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: Symmetry 2019, 11, 1086 9 of 21

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 where u 0 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 1 = 30.0, C 2 = -1.0 × 10 6 , C 3 = 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 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. 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 rs is 3.2 times Δx in the computation.     (26), of FDPM simulation using different particle spacing ∆x at different time; (b) maximum error, see Equation (27), of FDPM simulation using different particle spacing ∆x at different time. 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 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.

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.
The numerical simulation for Poisseuille flow is obtained under weakly compressible (Ma =

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.
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 −6 m 2 s −1 , L = 10 −3 m, ρ = 10 3 kgm −3 , and F = 10 −4 ms −2 , so the maximum velocity is 1.25 × 10 −5 ms −1 and the Reynold number is Re = 1.25 × 10 −2 . 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 −4 s. The initial particle spacings ∆x and ∆y are both set as 5 × 10 −5 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, rs remains 3.2 times Δx. 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.
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, rs remains 3.2 times Δx.  (26), of FDPM simulation using different particle spacing Δx at different time; (b) maximum error, see Equation (27), of FDPM simulation using different particle spacing Δx at different time. 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 -3 , indicating that the computational results agree well with the theoretical solutions. For different error evaluation indexes, the RERMS and REMAX 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.  (26), of FDPM simulation using different particle spacing ∆x at different time; (b) maximum error, see Equation (27), of FDPM simulation using different particle spacing ∆x at different time. 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 −3 , 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 -3 , cubic spline, and quantic spline.

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  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 −3 , cubic spline, and quantic spline.

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 0 , the theoretical solution of the flow velocity over time in the direction perpendicular to the plate [48] is as follows: Couette flow is numerically simulated under weakly compressible (Ma = 0.0125) conditions. The parameters for Couette flow are v k = 10 6 m 2 s −1 , L = 10 −3 m, ρ= 10 3 kgm −3 , and u 0 = 1.25 × 10 −5 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 0 . 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 −5 s. The initial particle spacings ∆x and ∆y are both set as 2.5 × 10 −5 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.  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 −2 , 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.
velocity at different times along the y direction (from the stationary plate to the sliding plate).
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 -2 , indicating that the computational results agree well with the theoretical solution. When t = 1.0 s, the two errors result in an RERMS value of 1.7 and REMAX 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 rs/Δ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. 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.

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.

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.

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 vk = 10 -6 m 2 s -1 , the cylindrical radius R = 2 × 10 -2 m, the volume force F0 = 1.5 × 10 -7 ms -2 , and the speed of sound c = 5.77 × 10 -4 ms -1 . Δx and Δy are 0.003 m, rs 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]. The size of the computational domain L = 0.1 m, the kinematic viscosity v k = 10 −6 m 2 s −1 , the cylindrical radius R = 2 × 10 −2 m, the volume force F 0 = 1.5 × 10 −7 ms −2 , and the speed of sound c = 5.77 × 10 −4 ms −1 . ∆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]. 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.
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.

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.

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

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 0 . The size of the cavity L = 1.0 m, the kinematic viscosity v k = 0.01 m 2 s −1 , the sliding velocity u 0 = 1.0 m/s, and the speed of sound c = 10.0 ms −1 . ∆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. velocity u0 = 1.0 m/s, and the speed of sound c = 10.0 ms -1 . Δx and Δy are 0.025 m, rs 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. 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.

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,

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.