The meshless local Petrov–Galerkin method based on moving Kriging interpolation for solving the time fractional Navier–Stokes equations

In this paper, we present a numerical scheme used to solve the nonlinear time fractional Navier–Stokes equations in two dimensions. We first employ the meshless local Petrov–Galerkin (MLPG) method based on a local weak formulation to form the system of discretized equations and then we will approximate the time fractional derivative interpreted in the sense of Caputo by a simple quadrature formula. The moving Kriging interpolation which possesses the Kronecker delta property is applied to construct shape functions. This research aims to extend and develop further the applicability of the truly MLPG method to the generalized incompressible Navier–Stokes equations. Two numerical examples are provided to illustrate the accuracy and efficiency of the proposed algorithm. Very good agreement between the numerically and analytically computed solutions can be observed in the verification. The present MLPG method has proved its efficiency and reliability for solving the two-dimensional time fractional Navier–Stokes equations arising in fluid dynamics as well as several other problems in science and engineering.

to the generalization of the notion of differentiation to fractional orders, including the Grünwald-Letnikov's definition, Riemann-Liouville's definition, Caputo's definition, Jumarie's definition, Chen's definition, and local fractional (fractal) Yang derivative. Each of these fractional derivatives presents both advantages and disadvantages. The two most frequently used generalizations of the derivative are the Riemann-Liouville's and Caputo's definition. The Riemann-Liouville fractional derivative is not suitable to deal with the physical problems in the real world since it requires the definition of fractional order initial conditions, which have no physically meaningful interpretation yet. Caputo introduced an alternative definition in which the initial conditions for the fractional order differential equations can be given in the same manner as for the integer order differential equations with a known physical interpretation. Another main advantage of the Caputo fractional derivative over the definition of Riemann-Liouville is that the Caputo derivative of a constant function is zero whereas the same does not hold for the Riemann-Liouville derivative. It is reasonable that the fractional derivative of a constant should be zero from a physical point of view. We refer the interested reader to the review article of De Oliveira and Tenreiro Machado (2014) in order to learn more details about a brief overview of proposed definitions of fractional integrals and derivatives.
The classical Navier-Stokes equations (NSE), developed by Claude-Louis Navier and George Gabriel Stokes in 1822, are the fundamental equations of fluid dynamics that extensively studied both theoretically and numerically. Many numerical researches have reported that the NSE are reasonably suitable for describing the mathematical modeling and numerical simulation relating to the flow behavior of fluids. Mathematically, the NSE are a time dependent system of nonlinear second-order partial differential equations (PDE) which can be derived from the basic conservation laws for mass and momentum. The time fractional NSE having been recently proposed by El- Shahed and Salem (2004) can be obtained by simply replacing the first-order time derivative term by derivative of non-integer order but retaining the first-and second-order space derivatives. To obtain the solution of the time fractional NSE has never been an easy task due to the nonlinearity which makes them quite difficult or impossible to solve. There are very few cases in which the solution of fractional NSE can be solved analytically. We have to make certain simplifications and assumptions about the state of fluid and configuration for the pattern of flow is to be considered. Several analytical methods have been proposed and developed for the solution of the time-fractional NSE. El- Shahed and Salem (2004) employed the Laplace, Fourier sine and finite Hankel transforms to obtain the exact solution for three different special cases of the time fractional NSE including the flow in a long circular pipe, the flow due to a plane wall suddenly set in motion, and the flow in plane Coutte motion. Momani and Odibat (2006) constructed the explicit solutions of the time-fractional NSE in cylindrical coordinates for an unsteady onedimensional motion of a viscous fluid by Adomian decomposition method (ADM). In this scheme, the analytic solution is formed in terms of series with easily computable term. The effective Adomian's method was extended to determine the analytical approximate solution of the two-dimensional time fractional NSE in Cartesian coordinates by Wang et al. (2009). The pressure gradient term was also assumed to be a constant which make the fractional model easy to implement. Shortly after that, the investigation with regard to solving the one-dimensional time fractional NSE in cylindrical coordinates has emerged continuously. Ragab et al. (2012) employed the homotopy analysis method (HAM) to obtain an approximate solution of the NSE with fractional order in the form of power series. They also pointed out that both the homotopy perturbation method (HPM), ADM and variational iteration method (VIM) are special cases of the HAM. Kumar et al. (2014) introduced a new analytical and approximate technique based on ADM and Laplace transform method (LTM) to obtain the solution of time-fractional NSE of a viscous fluid in a tube. Kumar et al. (2015) presented a reliable approach based on the new homotopy perturbation transform method (HPTM) to seek the analytical solution of NSE with time fractional derivative in a tube. Almost all the accomplishments on the theoretical analysis to acquire the solution of such a problem in the literature are only limited to one dimension. None of the above schemes is completely satisfactory when dealing with the multi-dimensional problems of the time fractional NSE. To fill this gap, particular attention in this work is devoted to the development of computational techniques and numerical strategies to the primitive variables formulation (velocity and pressure) of the time fractional NSE in two dimensions.
On the other hand, so far, only few researchers have attempted to extend and develop the so-called meshless methods to the solution of the time fractional NSE. The meshless or meshfree methods are the recently developed numerical technique which can be alternatively used to overcome the difficulties and limitations of mesh generation. Academic works regarding the applications of meshless methods have received considerable interest and have been published both theoretical and numerical aspects (Buhmann 2003;Kaufmann et al. 2010Kaufmann et al. , 2012Li et al. 2010;Abbasbandy et al. 2012Abbasbandy et al. , 2014Shirzadi et al. 2012;Francomano 2013, 2015;Dai et al. 2013a, b;Yimnak and Luadsong 2014;Khankham et al. 2015;Phaochoo et al. 2016a, b;Tu et al. 2016). The meshless methods are used to establish the system of algebraic equations for the whole problem domain without the use of a predefined mesh or domain discretization. Generally speaking, in accordance with the formulation procedure, meshfree methods fall into three categories: meshfree weak form methods, meshfree strong form methods, and meshfree weak-strong form methods. This work is particularly focused on the meshfree weak form methods. The common characteristic of weak form methods is that the original PDE of a problem is converted into an integral equation globally or locally based on a principle of weighted residual methods. There are a great number of meshfree methods that use local nodes for approximating the field variable, for example, the element free Galerkin (EFG) method (Belystchko et al. 1994), the meshless local Petrov-Galerkin (MLPG) method (Atluri and Zhu 1998a, b), the reproducing kernel particle method (RKPM) (Liu et al. 1995), the point interpolation method (PIM) (Liu and Gu 2001a), the radial point interpolation method (RPIM) (Liu and Gu 2001b;Liu 2000, 2002) and so forth. Some typical meshless methods based on global weak forms such as the EFG, RKPM, and RPIM, being "meshless" only in terms of the interpolation of the field variables, have to use background cells to evaluate integrals appearing in the local weak formulation. This is one reason why the aforementioned methods are not truly meshless method. One of the most popular meshless methods is the meshless local Petrov-Galerkin (MLPG) method which was first proposed by Atluri and Zhu (1998a, b) for solving linear potential problems. The MLPG approach is referred to as a one of the truly meshless methods which is used much more widely than other existing methods. In this method, a nodal sub-domain is used in place of background integration cells in order that all integrations are carried out locally over small sub-domains of regular shapes. In addition, the MLPG method is also different from the truly meshless method based on the local boundary integral equation (LBIE) method (Zhu et al. 1998a, b) in the fact that there is no singular integral in the present MLPG method, while some kinds of singular integrals have to be tackled in the meshless LBIE method (Atluri and Zhu 2000). The concept of shape function construction is one of the central and most important issues that significantly effect on the performance of meshfree methods. A number of ways to efficiently create shape functions have been proposed including the moving least squares (MLS) approximation, radial point interpolation method (RPIM), and moving Kriging (MK) interpolation. The MLS shape functions do not have the Kronecker delta property thereby making the imposition of essential boundary conditions complicated. Meanwhile, Dai et al. (2003) has proven that the Kriging interpolation is essentially the same as the RPIM as long as the same basis functions are used. That is properties found in RPIM should apply to Kriging shape function as well.
In this study, we employ the moving Kriging (MK) interpolation technique, which was first introduced in computational mechanics by Gu (2003). One notable feature of shape function constructed using the MK interpolation is that it possesses the Kronecker delta property which satisfies the essential boundary conditions automatically. The essential boundary conditions can be easily implemented without any special treatments.
As mentioned previously, there still lacks research efforts that have utilized the MLPG method to investigate the time fractional model of NSE in Cartesian coordinates. To the best of our knowledge, such a task has not yet been carried out while this work is being reported. The objective of this work is to extend further the application of the MLPG method to the NSE of fractional order. We organize the rest of this paper as follows. "The moving Kriging interpolation method" section gives some basic concepts of the moving Kriging interpolation for constructing shape functions in meshfree procedure. In "Research methodology", we introduce the governing two-dimensional time fractional NSE in Cartesian coordinate system and then we describe how to formulate the standard weak formulation and establish the discretized system equations. In "Numerical experiments and results" section, the numerical experiments are presented and discussed in detail to demonstrate and confirm the accuracy and efficiency of the proposed scheme. Finally, we complete the paper with some conclusion given in "Conclusion" section and some references are listed at the end.

The moving Kriging interpolation method
Kriging or Gaussian process regression was originally applied in geostatistics for spatial interpolation. Subsequently, Kriging interpolation was employed to construct shape functions for enhancement of the meshfree methods, intimately related to generalization of finite element methods (FEM). The procedure of constructing shape functions for meshfree methods using the MK interpolation is detailed in this section. Let us consider a subdomain, s ⊂ , the neighborhood of a point x and denoted as the domain of definition of the MK interpolation for the trial function at x. To approximate a distribution function u in s over a number of randomly located nodes {x i , i = 1, 2, 3, . . . N } where N is the total number of nodes in the sub-domain, the formulation of the MK interpolation u h (x), ∀x ∈ � s can be expressed in the form of linear combination of the shape functions: The matrices A and B are determined by where I is the unit matrix of size N × N and p(x) is a vector of the polynomial with m basis functions given by The commonly used linear basis in two-dimensional space is given by the quadratic polynomial basis is and the cubic polynomial basis is The matrix P has a size N × m and represents the collected values of (5) as and r(x) in Eq.
(2) has the form of where γ x, x j is the correlation between any pair of nodes located at x and x j , and it belongs to the covariance of the field value u( x, y, x 2 , xy, y 2 , m = 6, p T (x) = 1, x, y, x 2 , xy, y 2 , x 3 , x 2 y, xy 2 , y 3 , m = 10.
Many different correlation functions can be used for the correlation matrix. Gaussian function with a correlation parameter θ is often used to best fit the model due to its simplicity.
where r ij = x i − x j and θ > 0 is a parameter controlling the shape of the correlation function.
The partial derivatives of the shape function (x) with respect to x i are obtained as where (·) ,i and (·) ,ii denote the first-and second-order spatial derivatives, respectively.

Research methodology
In this section we introduce the governing time fractional NSE in Cartesian coordinate system and then the MLPG formulation and numerical implementation including local weak form and discretization techniques are described in detail. Let us now consider the time fractional NSE in two dimensions for an incompressible fluid in the following form: where u, v, and p are the velocity components in the x and y directions and pressure, respectively, Re represents the Reynolds number, f x and f y stand for externally imposed forces that act throughout the body of fluid along the x and y directions, respectively, α is the parameter indicating the order of the fractional, 0 < α < 1, and the time fractional derivative presented herein from Caputo's viewpoint is defined by where Ŵ(·) denotes the gamma function. In the case of α = 1, Eqs. (10), (11) and (12) reduce to the classical NSE.

Spatial discretization
Instead of giving the global weak form, the MLPG method constructs the weak form over local sub-domains, s , which is a small region taken for each node in the global domain . The local sub-domains overlap each other and cover the whole global domain. The local sub-domains (support domain) could be of any geometric shape and size, such as open or closed intervals in one dimension, circles or squares in two dimensions and spheres or cubes in three dimensions. For simplicity they are usually taken as a circle. The local weak form for Eqs. (10), (11) and (12) at each x i ∈ i s can be weighted by test functions and integrated over a local sub-domain. The following equations are obtained: where i s is a local sub-domain associated with the point i, i.e. a bounded circle centered at x i of radius r 0 , and w is a test function.
Substituting the trial function, u h , v h and p h , for u, v and p into the local weak forms except the nonlinear term gives In order to avoid the evaluation of any numerical integration in the weak form, the Kronecker delta function is chosen as the test function in each sub-domain. The Kronecker delta function maintains a unit value at the node and gives zero value at all other nodes. For this choice the support domain can be arbitrary. The local integral equations can be simplified to the following expression: Equations (20), (21) and (22)

Temporal discretization
For some positive integer M 1 , let t = T M 1 be the step size of time variable. The nodal points in the time interval [0, T ] are given by t n = n�t, n = 0, 1, 2, . . . , M 1 . The approximate solutions at the point u(x i , t n ) are abbreviated by u n i . In the current work, we apply a simple quadrature formula to obtain the discrete approximation of the time fractional derivative in Caputo's sense (Murio 2008 where ω (α) k = k 1−α − (k − 1) 1−α and σ α,�t = 1 Ŵ(1−α) 1 1−α 1 �t α . As shown in Eq. (23), the coefficient matrix B is itself a function of unknown u and v. Traditional iterative technique such as Newton-Raphson method might be applied to treat the nonlinearity, but often that is a very time-consuming process. Alternatively, to balance sufficient accuracy and acceptable computational expense, the linearization method by Taylor series expansion of a function can be adopted to approximate the nonlinear term. Since U(x, t n ) has the first-order continuous derivative, it follows that Substituting Eq. (24) into Eq. (23), making apply Eq. (25) to the unknown values contained within the matrix B, and omitting higher-order terms lead to the 3N × 3N discretized system of linear algebraic equations or equivalently For n = 1, we get and for n ≥ 2, A formula such as this which expresses one unknown value directly in terms of known values is called an explicit formula. This can be easily solved as a linear algebraic system of equations.

Numerical experiments and results
In this section, we provide two numerical examples to corroborate the accuracy and efficiency of the proposed method. The results of numerical experiments are compared with analytical solution by performing the root mean square (RMS) error: where U i and u i are the analytical and approximate solutions at points x i , respectively, and N is the total number of nodal points. The rate of convergence can be computed approximately by in which E 1 and E 2 are errors corresponding to grids with step size of time variable t 1 and t 2 , respectively. In the MK procedure, the cubic basis function p T (x) = 1, x, y, x 2 , xy, y 2 , x 3 , x 2 y, xy 2 , y 3 , m = 10, is used for all numerical computations because in general a cubic polynomial basis will yield a better result than quadratic and linear basis. The correlation parameter is another factor that has a significant effect on the solution. As a matter of fact, no exact rules can be derived appropriately to determine the optimal value of this parameter. However, it is often suggested to be θ = ω/d 2 c in order to smooth out small features in the data. The numerator ω = 0.2 is used according to the recommendation of Zheng and Dai (2011). The denominator d c is the nodal spacing near the point of interest. If the nodes are uniformly distributed, d c is simply the distance between two neighboring nodes in that the distance between two consecutive nodes in each direction is constant. When nodes are non-uniform, d c is the average distance of the nodes in the support domain. Throughout the experiment, the nodal points are assumed to be both regular (uniform) and irregular (nonuniform)   The initial and boundary conditions can be obtained from the exact solution. In this test problem, the nodal distribution is defined on the square domain = [0, 1] × [0, 1]. In Table 1 we show RMS errors and convergence rate of the velocity and pressure obtained in solving test problem 1 for different values of t with regular points. Figure 2 shows the graphs of approximate and exact solutions after 10 iterations with α = 0.99, �x = �y = 0.1, �t = 0.1, Re = 100 for the horizontal and vertical components of velocity and pressure. The numerical results are shown by dots while the exact solutions are generated by meshes. Also the velocity and pressure fields at Re = 100 are plotted in Fig. 3 with 11 × 11 points. To observe the behavior of the numerical solutions for long time (50 iterations), the RMS errors are plotted against the number of iterations, indicated in Fig. 4. The applicability of the proposed scheme to irregularly scattered nodes is also given in Table 2, Figs. 5, 6 and 7.

Test problem 2
The most widely used benchmark tests in incompressible unsteady flow simulation is the Taylor-Green vortex problem. The Taylor-Green vortex is an outstandingly canonical problem in computational fluid dynamics developed to study decaying vortex and turbulent transition. This problem serves as a well-known benchmark case study for testing and validating the reliability of a numerical scheme and is also used to perform significant test on behavior of the velocity and pressure properties. The exact closed form solution for decaying vortex problem is given by Re , Table 1 The RMS errors and convergence rate of the velocity and pressure obtained in solving test problem 1 for different values of t with regular points t �x = �y = 1/10 Convergence rate u v p u v p 1/10 9.9638 × 10 −5 9.9638 × 10 −5 6.4379 × 10 −4 ---  The initial and boundary conditions can be obtained from the exact solution. In test problem 2, we consider a nonlinear model of problem (10)-(12) in the unit square domain x, y ∈ [0, 1] × [0, 1], t ∈ [0, 1]. Table 3 gives RMS errors and convergence rate (36) v x, y, t = −cos(x)sin y e of the velocity and pressure obtained in solving test problem 2 for different values of t with regular points. The graphs of approximate and exact solutions after 10 iterations with α = 0.99, �x = �y = 0.1, �t = 0.1, Re = 100 for the horizontal and vertical components of velocity and pressure are depicted in Fig. 8. Under the same illustration as above, the numerical results are shown by dots while the analytical solutions are connected together by meshes. The velocity and pressure fields at Re = 100 are also shown in Fig. 9 with 11 × 11 points. The long-time behavior of the numerical solutions is displayed in Fig. 10. The distribution with random points is reported in Table 4, Figs. 11, 12 and 13 as well.    We can verify the temporal numerical accuracy with a fixed and sufficiently small spacing of x = y = 0.1 and different temporal step sizes t. As can be seen in Tables 1, 2, 3 and 4, the RMS errors decrease as the number of iteration is increased, that is to say we can get higher accuracy when the temporal step length is decreased. The performance of numerical accuracy is very satisfactory for the present choice of x, y and t. In the case of randomly scattered nodes, ω = 0.2 in the empirical formula for the correlation parameter is kept unchanged in the computation. The step size of time variable (�t) is chosen the same as regularly arranged nodes. Comparing the obtained results in Tables 1 and 2 also indicates that the irregular distribution of nodal points leads to better results than the regular distribution for pressure solutions only. The pressure errors for case of irregular nodal distribution are obviously much smaller than those of regular nodal distribution whilst the velocity solutions give slightly different numerical results. What is more, we notice from Fig. 4 that the velocity errors gradually reduce when the process is iterated until reaching to 25th calculation whilst the pressure errors behave like an exponentially decreasing function for the first few time steps and then seem to be reduced linearly, approximately near n = 15. The behavior of the solutions in Fig. 7c Fig. 12 The velocity and pressure fields at Re = 10 is also observed. The errors decrease exponentially after only a few time steps. It evidently dramatically tends to zero. In test problem 2, it is noteworthy that when the uneven nodal points are considered, accurate results are obtained at low Reynolds numbers. At Re = 100, the numerical error is not good enough. Herein, the Reynolds number is selected to be 10. Other Reynolds number values are possible but the error increases and an undesirable result is unavoidable. Presumably this is because the exact solution of the Taylor-Green vortex problem not only decays exponentially with time but also depends on the Reynolds number which has a significant effect on the solution. More investigation is needed. Nevertheless, the results presented through these figures and tables show the performance and validity of the proposed method to approximate the exact solution with high accuracy. All in all the results obtained from both of the test problem 1 and 2 are congruous with what we expected. In special case, when α → 1, it can be seen that the numerical results depicted in Figs. 2, 5, 8 and 11 tend to the analytical solution of classical NSE. Moreover, it is worth pointing out that the behavior of the solution by fractional model can be observed as the fractional derivative parameter is changed.

Conclusion
In this paper, we have presented a numerical scheme used to obtain the approximate solution of the time fractional Navier-Stokes equations. The truly meshless local Petrov-Galerkin approach based on a local weak formulation is employed to approximate the solution of field variables. The Kronecker delta function is chosen as the test function in each sub-domain to avoid the evaluation of any numerical integration in weak formulation. We apply the effective moving Kriging interpolation which possesses the delta function property for constructing shape functions at scattered points. A quadrature formula is used to obtain the discrete formulation of the time fractional derivative interpreted in the sense of Caputo. Besides, the nonlinear parts of the time fractional NSE can be treated by a simple linearization method based on Taylor series expansion. The capability and accuracy of the proposed scheme is demonstrated through solving the body force and Taylor-Green vortex problems. Very good agreement between the analytical and numerical results can be found. It is apparent that the present algorithm based on MLPG method can readily be extended to solve the unsteady incompressible Navier-Stokes equations involving non-integer order time derivative and is also found to be a computationally efficient and reliable method to deal with FDE.