Radial Basis Function Finite Difference Method Based on Oseen Iteration for Solving Two-Dimensional Navier–Stokes Equations

In this paper, the radial basis function finite difference method is used to solve two-dimensional steady incompressible Navier–Stokes equations. First, the radial basis function finite difference method with polynomial is used to discretize the spatial operator. Then, the Oseen iterative scheme is used to deal with the nonlinear term, constructing the discrete scheme for Navier–Stokes equation based on the finite difference method of the radial basis function. This method does not require complete matrix reorganization in each nonlinear iteration, which simplifies the calculation process and obtains high-precision numerical solutions. Finally, several numerical examples are obtained to verify the convergence and effectiveness of the radial basis function finite difference method based on Oseen Iteration.


Introduction
The Navier-Stokes equations are a set of equations used to describe fluid substance such as liquid and air, and can be used to simulate weather, ocean currents, water flows in pipes, the motion of stars in galaxies, etc. Therefore, they are of significant research value; this paper discusses the numerical solution method of steady incompressible Navier-Stokes equations.
Scholars have studied many numerical methods for Navier-Stokes equations [1]. For unsteady Navier-Stokes problems [2] and steady Navier-Stokes problems [3], a variety of traditional numerical methods have been proposed, including the finite element method, finite difference method, and finite volume method. In response to incompressible conditions, nonlinearity, long-time integration, and other difficulties encountered in numerical solutions of three-dimensional unsteady Navier-Stokes equations, this paper discusses the research status and latest research results of highly efficient and fully discrete finite element methods, which can be used to overcome these difficulties. In addition, this paper illustrates stability and error estimation of finite element-space discrete solutions and optimal error estimation of efficient fully-discrete finite element solutions for solving three-dimensional unsteady Navier-Stokes equations [4]. These methods are mesh-based, and the generation of meshes increases computational cost; moreover, the computational domain and the technology used for mesh quality and stabilization affect the computational accuracy. Especially for high-dimensional problems, it is difficult and costly to generate good meshes. Therefore, in recent years many scholars have introduced meshless methods [5,6] for solving partial differential equations, among which the radial basis function method [7] is increasingly popular.
The radial basis function method has been widely used for solving various partial differential equations, with the radial basis function finite difference method (RBF-FD) proposed considering ill-conditioned problems with a high occurrence rate in the matrix generated by the global radial basis function method. In this way, it is easy to realize RBF-FD discretization without needing to use a mesh at all. The RBF-FD method can be regarded as an extension of the standard finite difference and local refinement of global radial basis function method. This method uses the radial basis function to create weights for the RBF-FD formula, and the differential matrix obtained in this way is as sparse as that obtained by the standard difference method, which combines many advantages of radial basis function and traditional finite difference approximation. Thus, RBF-FD approximation [8] is an attractive substitute for the global radial basis function method.
In order to control the stagnation error, Natasha Flyer et al. introduced higher-order polynomials under the background of Gaussian (GA) and Polyharmonic Spline (PHS) interpolation [9,10]. PHS eliminates the stagnation error in combination with the polynomial to capture the basic physical properties of the problem, and achieves high-order accuracy without the need to adjust the shape parameters. The present paper chooses the PHS basis function and applies the radial basis function finite difference method with the polynomial to discretize the spatial operator of the equation. It is well known that the iterative method [11] is effective for Navier-Stokes equations with strong nonlinear terms, e.g., the Newton iterative method, Stokes iterative method, and Oseen iterative method. This paper focuses on the Oseen iterative method to deal with the nonlinear terms of the equation. The radial basis function finite difference method based on Oseen iteration proposed in this paper does not need a complete set of matrices in each nonlinear iteration, which is more suitable for nonlinear problems, and obtains high-precision numerical solutions.

Problem Setup
This paper studies the following two-dimensional steady incompressible Navier-Stokes equations [12].
where Ω is the bounded region on R 2 , u = (u, v) : Ω → R 2 is the velocity vector, p : Ω → R is the pressure, ν > 0 is the viscosity coefficient, and f = ( f 1 , f 2 ) is the external force acting on a unit volume of fluid, assuming that the uniqueness of p can be achieved by imposing the condition, Ω p dx = 0. We define the 2-norm and ∞-norm of the vector x as

Numerical Method
In this part, we focus on the radial basis function finite difference method with polynomial.

Radial Basis Function Finite Difference Method with Polynomial
First, this part provides a brief introduction to the radial basis function (RBF). The RBF, expressed as Φ(r) : R d → R, is a d-dimensional radially symmetric function, and is only dependent on r = x − x k 2 , where · 2 represents the Euclidean distance between two points, x is the point to be solved, and x k is the center position of the RBF. Common radial basis functions include Gauss (GA) (e −(εr) 2 ), multiquadratic function ( 1 + (εr) 2 ), and Polyharmonic Spline (PHS) (r m , m = 1, 3, 5, · · · ), where ε is the shape parameter that determines the radial basis function. Next, we introduce the radial basis function finite difference method with polynomial. First, L is defined as a linear operator, which can be ∆, ∂ ∂x , ∂ ∂y , etc. We can approximate the value of the operator L at node x c by linear combination of the function value {u k } n k=1 . For this, we need to select n nodes near the central node x c to form a node template, which is recorded as X = {x k } n k=1 , including x c , as follows: where n is the size of the node template and w k is the differential weight. In order to calculate the weight w k , we can specify that the linear combination of function values be accurate for interpolation s(x): The constraint condition is where p l (x) is degree l of the binary polynomial and Φ( x − x k 2 ) is the radial basis function. When d = 2 and the degree of the polynomial is l = 1, Equation (3) is with the constraint condition (4) Suppose that the function value at x k is f k , k = 1, · · · , n; then, the matrix form of Equations (5) and (6) is In this formula,Â is used to represent the matrix at the left end (n + 3) × (n + 3); thus, The linear operator L is used to calculate the result at x = x c : If f i = 1 and f j = 0, then j = i; thus, Ls(x)| x=x c = w i and i = 1, · · · , n, according to Formula (2). Therefore, where D is the identity matrix of (n × n) and O is the null matrix of (3 × n). An identity matrix is formed by adding three columns to the right matrix, which can be ignored. Thus, the left-hand matrix is added with [w n+1 , w n+2 , w n+3 ]. The right-hand matrix is multiplied byÂ and then transposed, obtaining the following linear equations: Therefore, the weights w 1 · · · w n are obtained from the above equations. In the following numerical examples, the resulting linear system changes with the varying degree of the polynomial. In this paper, the PHS radial basis function is mainly considered; Φ(r) = r m , m = 1, 3, 5, · · · , PHS is a piecewise smooth function, where r = x − x k 2 , and it is unnecessary to select a shape parameter ε .

Radial Basis Function Finite Difference Method Based on Oseen Iteration
In this part, we introduce discretization of the Navier-Stokes equations using the radial basis function finite difference method based on Oseen iteration. First, we use the radial basis function finite difference method with polynomial to discretize the component equations. Then, we linearize the equations based on Oseen iteration starting from this scheme, and finally obtain the fully discrete scheme of the equations.

Discrete Scheme of Equation
First, the radial basis function finite difference method with polynomial is used to discretize the spatial operator of the equation. Equation (1) is expressed as the following components: −ν( ∂u ∂x For Equation (9), the differential value at a point x c is where x cj n j is the n points around x c and w j n j is the weight coefficient corresponding to these points obtained by applying the interpolation condition f x cj = f j .
First, we assume that there are N nodes in a zone. At a certain point x c , Formula (9) is discretized according to RBF-FD interpolation and the following form is obtained: where n 1 represents the number of nearest nodes at the central node x c , used for approximating the weight coefficients of ∆u and ∇u, and n 2 represents the number of nodes at the central node x c , used to form a template of nearest nodes for approximating the weight coefficients of ∇p; moreover, w ∆ j is a weight coefficient for approximating ∆u, w ∇ x j is a weight coefficient for approximating operator ∂ ∂x , and w ∇ y j is a weight coefficient for approximating operator ∂ ∂y . Next, we use the Oseen iterative scheme to deal with the nonlinear terms of the steady Navier-Stokes equation. Iterating according to Discrete Equation (13), we obtain the following form: then −ν where k denotes the number of iterations, u(x k c ) denotes the value of u at the kth iteration, u(x k+1 c ) denotes the value of u at the (k + 1)th iteration, and the initial value of the velocity u 0 is obtained from the corresponding Stokes equation. In Equation (15), the value of u at the kth level is known; thus, the nonlinear term in Equation (9) is linearized. Likewise, there are similar discrete forms for Equations (10) and (11): We can now write (18) in the form of a block matrix: The stabilization parameter is added to Equation (19) to obtain a system of equations in the form of BU = F:  where B is the left-end weight matrix, W ∆+∇ u and W ∆+∇ v are (N × N) matrices composed of weight coefficients and boundary condition coefficients used to approximate u and v in Equations (15) and (16), respectively, W ∇ x u and W ∇ y v are a matrix of (N × N) composed of weighted coefficients for approximating ∂u ∂x and ∂v ∂y , W ∇ x p and W ∇ y p are a matrix of (N × N) composed of weighted coefficients for approximating ∂p ∂x and ∂p ∂y , F 1 and F 2 are known right-end terms, U is the solution of velocity u on point set X, V is the solution of velocity v on point set X, P is the solution of pressure p on point set X, D is the unit matrix of (N × N), and is the stabilization parameter, = 1.0 × 10 −6 .

Selection of Degrees of Added Polynomial and Number of Template Nodes
According to [13,14], without depending on the dimension, when using radial basis function PHS: Φ(r) = r p , when p is an odd number and the l-degree polynomials are supplemented to approximate the k-order derivative the convergence rate is not determined by the order of PHS, and is instead determined by the highest degree of the polynomial used. The solution accuracy is O(h l−k+1 ). For the PDE with only the first-order spatial derivative, the convergence speed can reach O(h l ), while for the PDE with the secondorder spatial derivative, the convergence speed can reach O(h l−1 ). The number of terms of the l-degree polynomial is . When approximating the value of an operator acting on the center point of the template, symmetry of the node templates is beneficial, as symmetrical node templates provide uniform information. In this paper, the numerical examples are considered in two layouts, namely, a right-angle node layout and a hexagonal node layout, with the node layouts shown in Figure 1. The right-angled node layout has completely symmetrical templates on numbers 5, 9, 13, 21, 25, 29, 37, 45, etc., while the hexagonal node layout has completely symmetrical templates on numbers 7, 13, 19, 31, 37, 43, 55, 61, etc. In this paper, the minimum distance between two points under the right-angled node layout is h = , while N is the total number of nodes.

Numerical Method
In this section, several numerical examples are used to demonstrate the effectiveness of the proposed methods; the relative error L 2 and error L ∞ are applied for comparison of the numerical solution and true solution. The relative error L 2 is defined as error L 2 = u−u h 2 u 2 , while the error L ∞ is defined as is the obtained numerical solution and u is the true solution x ∈ X.

Convergence Test
This example studies the two-dimensional Navier-Strokes equations under the layout of right-angled nodes. We assume that Equation (1) has the true solutions where Ω = [0, 1] × [0, 1], the viscosity coefficient ν = 1, and the nodes are uniformly distributed on Ω as shown in Figure 1 left. For purpose of the convergence test, the basis function Φ(r) = r 7 is applied to add third-order and second-order polynomials for approximation of the velocity u = (u, v), while the basis function Φ(r) = r 5 is applied to add second-order and first-order polynomials for approximation of the pressure p. For the number of node templates approximated at a central node of the velocity and pressure in Equations (15) and (16), we select n 1 = 21 and n 2 = 9, respectively. Figure 2 left shows the addition of the third-order polynomial for approximation of the velocity and the error results for pressure approximation when adding the secondorder polynomial. The figure shows the relative error L 2 when the number of nodes N is 25, 81, 289, 1089, and 4225. The error decreases with increasing number of nodes, and the convergence order of velocity and pressure reaches third-order and second-order. When the number of nodes N is 4225, the velocity error is 4.4583 × 10 −5 and the pressure error is 2.2254 × 10 −2 . The right0hand side of Figure 2 shows the error results of the approximation when adding the second-order polynomial to the velocity and adding the first-order polynomial to the pressure. When the number of nodes N is 4225, the relative errors L 2 of the velocity and pressure are 4.2393 × 10 −5 and 2.8000 × 10 −2 , respectively. By comparing the two figures, the error between the numerical solution and the true solution obtained from approximation by adding third-order and second-order polynomials to the velocity and pressure is far less than that obtained from approximation by adding second-order and first-order polynomials to the velocity and pressure. When the number of nodes is less than 1089, the former reaches the corresponding convergence order. As the number of nodes approaches 4225, the order of velocity is less than 2 due to the reduction in point density. This example illustrates that under the right-angled node layout, the radial basis function difference method based on Oseen iteration is convergent and effective for solving incompressible Navier-Stokes equations.

Small Viscosity Problem
In order to test the suitability of the method proposed in this paper for small viscosity problems, we begin by assuming that Equation (1) has a true solution (21). By setting the viscosity coefficient ν = 0.00001, the nodes are distributed in a right-angled node layout. In this example, the basis function Φ(r) = r 7 is applied to add the third-order polynomial for approximation of the velocity u = (u, v), and the basis function Φ(r) = r 5 is applied to add the second-order polynomial for approximation of the pressure p. Table 1 shows the error results when ν is 0.00001 and number of nodes N is 25, 81, 289, 1089, and 4225. When the number of nodes N is 4225, the relative error L 2 of the velocity is 6.7100 × 10 −5 and the relative error L 2 of the pressure is 2.0271 × 10 −2 . Figure 3 shows the corresponding error order.

Hexagonal Node Layout
In this section, in order to test the suitability of the method proposed in this paper for different node layouts, we consider the two-dimensional equation under the hexagonal node layout, as shown in Figure 1 right, assuming that Equation (1) has a true solution (21) where Ω = [0, 1] × [0, 1] and the viscosity coefficient ν = 1. In this example, the basis function Φ(r) = r 7 is applied to add the third-order polynomial for approximation of the velocity u and the basis function Φ(r) = r 5 is applied to add the second-order polynomial for approximation of the pressure p. For the velocity and pressure in a discrete equation, the number of node templates approximated at a central node (n 1 and n 2 ) are separately selected as n 1 = 31 and n 2 = 13. When number of nodes is 4699, the relative error L 2 of the velocity is 2.8638 × 10 −5 and the relative error L 2 of the pressure is 1.3900 × 10 −2 . Figure 4 shows the relative error L 2 between the numerical solution and the true solution when the number of nodes N is 68, 279, 1166, and 4699. In summation, these numerical experiments show that the radial basis function finite difference method based on Oseen iteration has good performance for solving Navier-Strokes equations under the hexagonal node layout.

Conclusions
This paper proposes a radial basis function finite difference method based on Oseen iteration for solving two-dimensional steady Navier-Stokes equations with discretization of the spatial operators of the Navier-Stokes equations using the radial basis function finite difference method with polynomial, then linearizing the equation based on Oseen iteration. In this paper, we provide numerical solutions for the Navier-Stokes equations under right-angled and hexagonal node layouts, analyze and compare the influence of the number of polynomial additions on the accuracy and convergence of the solution, and verify the effectiveness of the proposed method. Under the two node layouts studied in this paper, the proposed method offers high-precision numerical solutions for solving Navier-Stokes equations while demonstrating good performance. Its extension to unsteady problems will be the focus of our future research work.