A numerical study of the Peridynamic Differential Operator discretization of incompressible Navier-Stokes problems

. We study the incompressible Navier-Stokes equations using the Projection Method. The applications of interest are the classical channel ﬂow problems such as Couette, shear, and Poiseuille. In addition, we consider the Taylor-Green vortex and lid-driven cavity applications. For discretization, we use the Peridynamic Diﬀerential Operator (PDDO). The main emphasis of the paper is the performance of the PDDO as a discretization method under these ﬂow problems. We present a careful numerical study with quantiﬁcations and report convergence tables with convergence rates. We also study the approximation properties of the PDDO and prove that the N -th order PDDO approximates polynomials of degree at most N exactly. As a result, we prove that the PDDO discretization guarantees the zero row sum property of the arising system matrix.


Introduction
The Peridynamic Differential Operator (PDDO) is a nonlocal operator introduced by Madenci and his colleagues.Madenci et al. [15] derived the peridynamic equation of motion by replacing the derivatives in Navier's equation by the PDDO.In addition, peridynamic form of the strain energy density function, the force density vectors, and the deformation gradient tensor can all be derived by the PDDO.Hence, this differential operator's intimate connection to Peridynamics [21,22] led to its naming.The first book on PDDO by Madenci et al. on PDDO [14] contains abundant applications in various fields.The second book [17] covers the advances in the field of Peridynamics as well as the PDDO.
Several aspects of the PDDO have been studied: applications with Helmholtz-type problems [20], the Eikonal equation [7], the solution of linear and nonlinear partial differential equations (PDEs) [16], and physics informed deep learning for the discovery of PDEs whose solutions develop sharp gradients [12].The PDDO was used for classical flow problems [11,19] similar to the ones we consider.
The PDDO has several strengths.In particular, changing the PDE can be done by modifying a few lines of code.It has the ability to handle various boundary conditions (BCs) with ease.Its extension to arbitrary dimensions is straightforward.It enjoys meshless quadrature, thereby furnishing simple system matrix assembly.

Introduction of the PDDO
Let x be a point in Ω and δ > 0 be given.Define the family of x as For x ′ ∈ H x , define ξ := x ′ − x.Assume a sufficiently smooth function u : Ω → R. The PDDO utilizes the Taylor expansion with remainder term where for some c between x and x + ξ.Henceforth, we consider N > 0 to be fixed.
The so-called PD functions g p N for p = 0, . . ., N are essential to the construction of the PDDO.They are chosen to be linear combinations of monomials ξ q , q = 0, . . ., N, multiplied by the weight function w, and are defined by a p q w(|ξ|)ξ q , p = 0, . . ., N, where the coefficients a p q will be specified later.Although there are many possible choices of the weight function, throughout the paper it is chosen as (2.4) Using g p N , the PDDO constructs an approximation to the desired derivative of f .The first step is to test u(x + ξ) against g p N (ξ) over H x .Using (2.1), the first step yields the following expression: In (2.6), u (n) (x) has been moved outside the integral because the integration in (2.5) is with respect to x ′ .The PD functions are designed by requiring to satisfy the following orthogonality property: On incorporating (2.7) into (2.6), the sum collapses Neglecting the remainder term in (2.8) yields the PDDO approximation to the desired derivative of u: (2.9) Thus, we define the p-th order PDDO derivative, a linear operator, as If one resorts to an x and x ′ only representation, it is easy to see that the PDDO is an instance of a convolution operator: where g p N (x ′ − x) plays the role of a kernel function.We constructed similar type of operators using various kernel functions to enforce local boundary conditions (BCs) in nonlocal problems [1,2,4,5].
By substituting the expression of g p N given in (2.3) in the orthogonality property (2.7), one obtains The coefficients a p q are stored in the matrix a ∈ R (N +1)×(N +1) .By defining A, b ∈ R (N +1)×(N +1) , the coefficients are determined by solving the system Aa = b, where the shape matrix A and the right hand side matrix b are given by N q=0 A nq a p q = b p n , n, p = 0, . . ., N, with entries defined by Here we slightly abuse the traditional matrix notation by denoting the (q, p)-th entry of the matrix a by a p q instead of the standard notation a qp , and likewise b p n instead of b np .

System Matrix Assembly
The stiffness matrix corresponding to Π p N is assembled using the collocation method with piecewise constant basis functions.Define the domain Ω := (a, b) and the uniform grid with grid size h = b−a ndof−1 with nodes x i , i = 1, . . ., ndof.At DOF x i , the action of the PDDO is The piecewise constant discretization of u is After plugging the expression of u in (3.1), we obtain where φ j is the piecewise constant basis functions with support The integral in (3.2) is computed by a piecewise constant approximation Then, Consequently, the entries of the system matrix are Due to the convolution nature of the operator as given in (2.10), we numerically observe that the employment of collocation method gives rise to a centrosymmetric (symmetric with respect to the center) [3] system matrix.The proof of this interesting centrosymmetry property will be subject of future work.

Approximation Property of the PDDO
The quality of the PDDO approximation to u (p) (x) is determined by the term Hx R N +1 g p N (ξ) dx ′ .For polynomials, this approximation is exact as shown below.
Corollary 4.2.Let u(x) ≡ c be a constant function.Then In other words, constant functions are in the kernel of Π p N , p = 1, . . ., N .
Proof.Since u is a polynomial of degree zero, by (4.1) Corollary 4.2 gives the following important linear algebra implication.
The ZRS is a key property for discretization operators; for instance, see [5,Sec. 7.1] for collocation, see [6] for meshless methods.We carefully monitor that the ZRS property is maintained in machine precision for all assembled system matrices.

The PDDO Discretization of a PDE
The PDDO discretization of a PDE consists of replacing the derivatives appearing in the PDE with their PDDO counterparts.To this end, we generalize the discussion in Sec. 2 to multivariate functions.Since all of the test problems considered in subsequent sections are in two spatial dimensions and also for clarity of the presentation, we restrict ourselves to 2D.Let x = (x 1 , x 2 ) and x ′ = (x ′ 1 , x ′ 2 ) be points in Ω and ξ = x ′ − x.Suppose that δ = (δ 1 , δ 2 ) with δ i > 0 be given.Define the family of x as The bivariate PDDO takes the form where p = (p 1 , p 2 ) with p i ≥ 0 and p 1 + p 2 ≤ N .The PD functions g p N are chosen to be linear combinations of bivariate monomials 1 ξ q 2 2 , with q i ≥ 0 and q 1 + q 2 ≤ N , multiplied by the weight function w, and are defined by We denote the partial derivatives by Define the set of possible derivatives as P = {p} = {(0, 0), (1, 0), (2, 0), . . ., (N, 0), . . ., (0, N )}.
Since the highest order of derivative is N , in 2D, the number of possible derivatives is equal to the cardinality of P , which is Define the PDDO approximation to ∂ p u as Thus, we define the p-th order PDDO derivative as where n! = n 1 !n 2 ! and δ np = δ n 1 p 1 δ n 2 p 2 .By substituting the expression of g p N given in (5.6) in the orthogonality property (5.4), one obtains It is instructive to calculate the number of coefficients a p q needed to construct the PD functions given in (5.6).For each fixed p, there are (N +2)(N +1) 2 terms in the double sum appearing in (5.6), which is equal to N +2 2 .On the other hand, there are N +2 2 possible choices of p as pointed out in (5.7).Hence, the coefficients a p q are stored in a matrix of size Similarly, defining the other two matrices with the same size, A, b ∈ R d 2 ×d 2 , the coefficients are determined by solving the system where (5.5) 5.1.Extension to Higher Dimensions.The extension to higher dimensions is straightforward.
Let M be the spatial dimension.Let x = (x 1 , . . ., x M ) and x ′ = (x ′ 1 , . . ., x ′ M ) be points in Ω and ξ = x ′ − x.Suppose that δ = (δ 1 , . . ., δ M ) with δ i > 0 be given.Define the family of x as The multivariate PDDO takes the form where p = (p 1 , . . ., p M ) with p i ≥ 0 and The PD functions g p N are chosen to be linear combinations of multivariate monomials N , multiplied by the weight function w, and are defined by Define the set of possible derivatives as P = {p}.In M -dimensions, the number of possible derivatives is equal to (5.7) The p-th order PDDO derivative is defined in exactly the same way as in (5.3).Following the steps similar to those in 2D, the coefficients are stored in the matrix a p q ∈ R d M ×d M .Defining the matrices A, b ∈ R d M ×d M , the coefficients are determined by solving the system where A and b are identical to those in (5.5).

The Problem Description and the Projection Method
To describe the fluid dynamics, we consider the Navier-Stokes equations in the primitive variables, the fluid velocity u and the pressure p where f is the external force, ρ is the constant density, and µ is the fluid viscosity.With appropriate BC on u, we solve (6.1) by utilizing the projection method [8,9,10].The projection method was also used with the PDDO for similar flow problems in the Eulearian framework [19].
Let Γ denote the boundary of the domain Ω.Let ∆t be the size of time step, u n be the fluid velocity at time t = n∆t for n = 0, 1, . ... To maintain the incompressibility condition (6.2), one applies the projection method as follows: First compute an intermediate velocity field u * by subject to the same BC as that for u.Then, compute p n+1 solving the boundary value problem Finally, obtain the updated velocity by using the prescribed BC for u.
Remark 6.1.The formulation of the projection method described above guarantees the incompressibility condition ∇ • u n+1 = 0 for all n = 0, 1, . ... This guarantee is valid for the analytical solution, without any reference to a numerical method except for the time discretization in (6.3), which is not our emphasis in this study.As will be shown in Sec. 9, we observed that incompressibility cannot be ensured for the numerical solution in machine precision for all problems.Thus, further research is necessary to modify the PDDO discretization, and perhaps the projection method, to ensure the incompressibility condition.

The Numerical Implementation of the Projection Method and The Corresponding PDDO Discretization
We describe the numerical implementation of the projection method in three steps by resorting to a componentwise notation.The PDDO is utilized for each partial derivative involved in the steps of the projection method.For clarity, we explicitly state the corresponding PDDO discretization as well.
The corresponding PDDO steps are as follows: Step 1-PDDO.
using the prescribed BC for (u 1 , u 2 ).

Classical Fluid Flow Problems
We study a set of classical fluid flow problems, i.e., Couette flow, shear flow, Poiseuille flow, and Taylor Green vortex.A study of these flows in the Lagrangian and Eulerian descriptions was reported in [11] and [19], respectively.
In all of the problems, the domain is chosen as Ω = (L, R) × (B, T ) with L = B = 0 and R = T = 1.In addition, the external force f is always assumed to be zero.

Couette Flow. The classical Couette problem involves fluid flow between two plates with initial condition u
and u 2 (x 1 , x 2 , 0) = 0, (x 1 , x 2 ) ∈ Ω, where u 1 and u 2 denote the horizontal and vertical velocities, respectively.
The two plates are placed at South and North edges, x 2 = B and x 2 = T , respectively, with Dirichlet BCs.The plates have the following BCs.For x 1 ∈ (L, R), The BCs on the East and West edges are periodic, i.e., u i (L, The corresponding Reynolds number is where U max = 1 m/s is the maximum velocity.
For the Couette flow, the exact solution of the flow front is a function of x 2 only, given by a series [11,18] where the stationary term is given by 8.2.Shear Flow.Similar to the Couette flow, the classical shear problem involves fluid flow between two plates with initial condition u 1 (x 1 , x 2 , 0) = 0, (x 1 , x 2 ) ∈ Ω, and u 2 (x 1 , x 2 , 0) = 0, (x 1 , x 2 ) ∈ Ω.
The two plates are placed at South and North edges, x 2 = B and x 2 = T , respectively, with Dirichlet BCs.The plates have the following BCs.For x 1 ∈ (L, R), The BCs on the East and West edges are periodic, i.e., u i (L, The corresponding Reynolds number is Re = 1. The exact solution can be obtained from (8.1) by a simple scaling and is a function of x 2 only, given by The BCs on the East and West edges are periodic, i.e., u i (L, The corresponding Reynolds number is Re = 0.25.
For plane Poiseuille flow, the exact solution is a function of x 2 only, given by where G = 2, µ = 1, and

8.4.
Taylor-Green Vortex.The Taylor-Green vortex is an unsteady flow of a decaying vortex whose exact solution is given by and is depicted in Fig. 8.1.The BCs are periodic and are given by The three sides of the box are placed at East, West, South edges x 1 = L, x 1 = R, x 2 = B, respectively, and an open top at the North edge, x 2 = T .The following Dirichlet BCs are enforced:

Numerical Experiments and History of Convergence
Our numerical study concentrates on reporting history of convergence tables with respect to DOF.We also report the divergence of velocity field to the monitor the incompressibility of velocity at , and the corresponding exact solutions (dotted lines), which are almost invisible due to the perfect agreement the numerical level.We start with the details of the numerical implementation used to report the tables and figures.In the tables, the "Grid" column specifies the grids used in the computations, more explicitly, Grid = i for i = 1, . . ., 6 means that the corresponding grid has equally spaced 8i + 1 nodes on each axis of the domain.DOF are identified with the grid nodes.As a result, a uniform grid with spacing h = 1 8i is used.Since all of our examples are in 2D, the number of degrees of freedom (NDOF) in the grid is N i = (8i + 1) 2 .The grid contains DOF at the physical boundary of the domain; see Fig. 9.1.The remaining DOF are placed inside the domain in an equispaced manner.In the column labeled "Order" we display an approximate order of convergence as follows.We report the ℓ ∞ -norm of e i with Grid = i and the quantity This formula is derived by assuming e i = C N −r i , for some constant C independent of N i , where r is the order of convergence.
The horizon is chosen as δ = 3 h in alignment with the choice of δ = (N + 1) h as in [14,15] with N = 2 since we are solving a second-order PDE.The weight function is chosen to be the Gaussian According to Lemma 4.1, since N = 2, the PDDO method captures polynomials of degree up to 2 exactly.As a result, the method is expected to approximate u 1 (x 1 , x 2 , ∞) term of (8.1), (8.2), and (8.3) up to machine precision.Thus, the errors displayed in Tables 9.1, 9.2, 9.3 are due to approximating the transient terms represented by the infinite sums in (8.1), (8.2), and (8.3).This is substantiated by examining the errors corresponding to time t = 2 at which time the exact solution is quite close u 1 (x 1 , x 2 , ∞) due to the presence of the exponentially decaying term.
In our implementation unlike [11,19], we do not need to employ fictitious DOF.Instead we simply identify the DOF corresponding to the periodic parts of the domain boundary.Note that this is possible because there are DOF located precisely on the physical boundary.This would not be possible if DOF were placed at the cell centers as usually done in the literature [11,14,15].We think that this identification approach is in a better agreement with the physics of the underlying problem.In addition, we find it easier to implement because all it takes is mapping the DOF to a cylinder (channel flow) or a torus (the Taylor-Green vortex).
For all channel flow problems, the numerical solution is in such a perfect agreement with the exact solution (dotted line) that it is invisible in Figs.9.2, 9.3, and 9.4.For the Taylor-Green vortex problem, the numerical solution is also in a perfect agreement with the exact solution; see Figs. 9.5 and 9.6.For lid-driven cavity problem, the exact solution is not available.The velocity profile and The figures are limited in conveying the convergence of the numerical method.That is why, we extended our efforts to quantifications in the form of history of convergence tables.
For all channel flows under consideration, u 2 = 0, hence, we only report u 1 .The results indicate that the order of convergence of u 1 is at least 2. Since N is chosen to be 2 and the fact that the method captures second order polynomials exactly, one may anticipate a convergence of order up to 3, which would be optimal.In this sense, the PDDO method converges suboptimally.Since the exact solution for velocity is available, one can derive the expression of the gradient of pressure from (6.1).We computed the gradient using the PDDO.We observe that the error of the gradient 0. For the Taylor-Green vortex, we have the exact formula for the pressure.Our numerical approximation to pressure converges, but it is hard to determine a clear convergence order.In the last two rows of Table 9.4 at t = 1.0 and t = 2.0, we did not display the orders of convergence because we believe that the error gets stagnated at around 1.0E-7 and 1.0E-10, respectively.For all channel flow problems, we observe that the divergence of velocity field is zero in machine precision.Hence, the incompressibility condition (6.2) is numerically satisfied; see Tables 9.1, 9.2, and 9.3.However, in the Taylor-Green vortex, at early simulation times, the incompressibility  condition can be achieved only up to 1.0E-4, at later times it is achieved up to 1.0E-7.For liddriven cavity, the exact solution is not available.That is why, we cannot report the error.Also, the incompressible condition is not satisfied due to singularities at the upper corners of the domain [13].
We did not display the numerical results with N = 4 because the matrix A in (5.5) becomes severely ill-conditioned.It would be worthwhile to study the case N = 4 only with an effective preconditioner, which is a topic out of the scope of this paper.

Conclusion
We solved the incompressible Navier-Stokes equation utilizing the projection method for classical flow problems such as channel flow, the Taylor-Green vortex, and the lid-driven cavity.One of the main contributions is a detailed history of convergence study in which we exhibit the convergence order of the method with respect to the number of DOF.Another main contribution is a nodal DOF version of the PDDO as opposed to cell-centered approach commonly used in the literature.The nodal version eliminates the necessity of fictitious nodes for periodic BC.To the best of the authors' knowledge, this is a novel way of implementing the PDDO for the flow problems.We uncover an exactness property of the PDDO, namely, polynomials of degree N can be approximately exactly, thereby, leading to the zero row sum property.Also, in channel flow problems, the exactness property results in a high order of convergence.

Figure 8 . 1 .
Figure 8.1.The initial velocity contour map of the Taylor-Green vortex

8. 5 .
The Lid-Driven Cavity.The lid-driven cavity problem involves fluid flow in an open topped 2D box with initial condition u 1

Figure 9 . 3 .Figure 9 . 4 .
Figure 9.3.The evolution of the flow front of shear flow at the middle of the channel, i.e., x 1 = 1 2 , and the corresponding exact solutions (dotted lines), which are almost invisible due to the perfect agreement

2 Figure 9 . 5 . 6 Figure 9 . 6 .Figure 9 . 7 .
Figure 9.5.The comparison of the history of the velocity magnitude distribution of Taylor-Green vortex (left column) and the corresponding exact solutions (right column)

Table 9 .
1. Error report for the Couette flow

Table 9 .
3. Error report for the Poiseuille flow

Table 9 .
4.Error report for the Taylor-Green vortex