Abstract
The differential quadrature (DQ) method is able to obtain quite accurate numerical solutions of differential equations with few grid points and less computational effort. However, the traditional DQ method is convenient only for regular regions and lacks upwind mechanism to characterize the convection of the fluid flow. In this paper, an upwind local radial basis functionbased DQ (RBFDQ) method is applied to solve the NavierStokes equations, instead of using an iterative algorithm for the primitive variables. The nonlinear collocated equations are solved using the LevenbergMarquardt method. The irregular regions of 2D channel flow with different obstructions situations are considered. Finally, the approach is validated by comparing the results with those obtained using the wellvalidated Fluent commercial package.
1. Introduction
Many computational methods have been applied to solve the NavierStokes equation, prominent among which are the finite differential method (FDM), the finite element method (FEM) and the finite volume method (FVM). Although these methods have been successfully applied to many engineering situations, their accuracy highly depends on mesh quality. Specifically, appropriate mesh generation is typically a major challenge in computational fluid dynamics (CFD). Moreover, the mesh generation, modification and remeshing become difficult and timeconsuming when the geometry of the problem is complex. Thus, numerical simulations and solution procedures using these traditional methods are often limited by slow mesh generation.
In the past decade, the socalled meshfree methods have became one of the attractive research areas in computational mechanics. A number of meshless methods have been proposed to date. These methods include the smoothed particle hydrodynamics method (SPH) [1], the diffuse element method (DEM) [2], the elementfree Galerkin method (EFGM) [3], the reproducing kernel particle method (RKPM) [4], the partition of unity method (PUM) [5], the hpclouds method (HPCM) [6], the finite point method (FPM) [7], the meshless local PetrovGalerkin method (MLPGM) [8], the movingparticle semiimplicit method (MPSM) [9] and the general finite difference method (GFDM) [10]. Recently, a new meshless method was proposed based on the socalled radial basis functions (RBFs), have become attractive for solving Partial differential equations (PDEs). Initially, RBFs were developed for multivariate data and function interpolation, especially for higher dimension problems. However, their “truly” meshfree nature motivated researchers to use them to deal with PDEs. Kansa [11, 12] initiated the concept of solving PDEs by using RBFs. Since the resultant coefficient matrix of Kansa’s method is not symmetric, Fasshauer [13] proposed the Hermite type method, which can generate symmetric coefficient matrix that guarantees the related linear equations being solvable. Other great contributions to solve PDEs by using RBFs include the work presented by Demirkaya [14], Fornberg [15], Hon and Wu [16], Chen and Tanaka [17], Chen [18], Xiang [19] et al. The advantages of RBFsbased methods have been summarized by Wu [20] as follows: (1) it is a truly meshfree algorithm; (2) it is space dimension independent in the sense that the convergence order is of $O\left(hd+1\right)$ where $h$ is the density of the collocation points and $d$ is the spatial dimension.
Unfortunately, the numerical schemes based on RBFs to solve PDEs proposed up to date have one feature in common, i.e., they are actually based on the function approximation. In other words, these methods directly substitute the expression of function approximation by RBFs into a PDE, and then change the dependent variables into the coefficients of function approximation. The process is very complicated, especially for nonlinear problems. Thus, the method has not been extensively applied to solve practical problems. The objective of this paper is to present a numerical scheme, which is simple to be implemented for linear and nonlinear problems. The method preserves the properties of high accuracy and mesh freeness of RBFs and the derivative approximation directly through the differential quadrature (DQ) methodology [21]. The DQ technique originates from the integral quadrature. The basic idea of the DQ method is that the derivatives of unknown function can be approximated in terms of the function values at a set of points, either uniformly or randomly distributed. Numerical experiments showed that the combination of DQ technique and RBFs provides a novel and effective algorithm for the use of RBFs to solve PDEs.
Furthermore, as observed by Shu et al. [22] and Sun et al. [23], the DQ method is limited in application to regular regions by its quadrature rules, such as rectangular and parallelogram regions, or circular, concentric, and sectorial regions. When an irregular region must be considered, the implementation of the DQ method becomes difficult. In addition, the upwind mechanism is necessary in numerical computation of fluid flow, but the DQ method lacks the upwind mechanism to characterize the convection of the fluid flow.
To overcome the difficulties and drawbacks mentioned above, a local RBFbased differential quadrature method (RBFDQ) is developed in this paper. It can handle irregular regions conveniently and possesses the upwind mechanism. Specifically, compared with the classical iterative approach, a local RBFDQ method is proposed to directly approximate the derivatives, in which RBFs are formed by using local supporting knots. This method is very flexible as it does not requires extra effort on the division of the computational domain. We collocate the fully nonlinear NavierStokes equations and solve the resulting nonlinear system directly using the LevenbergMarquardt method. The approach present a novel approach for accurate and stable solution of the NavierStokes equations using upwind local RBFDQ method. It seems to produce accurate and stable results compared to the conventional iterative methods. The proposed numerical scheme is applied to determine the velocity and pressure for the flow through a channel with the obstructions in different patterns. The predicted channel flow results mirror the analytical solutions, while the 2D results compare favorably with those obtained using the wellvalidated CFDFluent commercial package.
2. Numerical methodology
2.1. Systems
Three types of geometry are considered as shown in Fig. 1. Fig. 1a shows a 2D flow between two flat plates with one rectangular obstacle placed at the center of the bottom plate, Fig. 1b shows a 2D flow between two flat plates with two rectangular obstacles, Fig. 1c shows a 2D flow between two flat plates with three rectangular obstacles in staggered pattern. The above models all belong to irregular step flow, which is classical CFD model.
The meshless method requires a complete set of boundary conditions in order to obtain a solution. In this study, the noslip condition is imposed on all walls, implying zero velocity on the walls for the 2D channel flow. The inflow velocity for the 2D channel flow is specified by a mean velocity ${U}_{\infty}$. The $x$component of velocity is assumed to be equal to ${U}_{\infty}$ at the inlet, the $y$component of velocity is assumed to be zero, and the pressure is a reference gage value of zero. The outlet is located sufficiently far from the inlet ($L=10\text{H}$) to ensure fullydeveloped flow in the channel. At the outlet, we assume that the flow exits with straight streamlines, the gradients of all variables are zero in the flow direction and the continuity equation is satisfied. The wall and outlet pressure boundary equations are obtained by taking the scalar product of the momentum equation with the outward normal to the boundary. These boundary conditions are summarized in Table 1.
Fig. 1Different models of the system used for computational experiments
Table 1Boundary conditions
Velocity  pressure  
${\Gamma}_{wall}$  ${u}^{*}=0,\mathrm{}{v}^{*}=0$  $\frac{\partial {p}^{*}}{\partial {y}^{*}}+\frac{1}{\mathrm{R}\mathrm{e}}\left(\frac{{\partial}^{2}{v}^{*}}{\partial {x}^{*2}}+\frac{{\partial}^{2}{v}^{*}}{\partial {y}^{*2}}\right)=0,(\mathrm{f}\mathrm{o}\mathrm{r}ny\ne 0\mathrm{a}\mathrm{n}\mathrm{d}nx=0)$ $\frac{\partial {p}^{*}}{\partial {y}^{*}}+\frac{1}{\mathrm{R}\mathrm{e}}\left(\frac{{\partial}^{2}{u}^{*}}{\partial {x}^{*2}}+\frac{{\partial}^{2}{u}^{*}}{\partial {y}^{*2}}\right)=0,(\mathrm{f}\mathrm{o}\mathrm{r}nx\ne 0\mathrm{a}\mathrm{n}\mathrm{d}ny=0)$ 
${\Gamma}_{inlet}$  ${u}^{*}=1,\mathrm{}{v}^{*}=0$  ${p}^{*}=0$ 
${\Gamma}_{outlet}$  $\frac{\partial {u}^{*}}{\partial x}=0,\mathrm{}\frac{\partial {v}^{*}}{\partial x}=0$  $\frac{\partial {p}^{*}}{\partial x}=0$ 
2.2. Governing equation
In this section, the local RBFDQ method with upwind mechanism is described. It is then used to predict the velocity and pressure distribution in a twodimensional channel flow. The flow is assumed to be steady, viscous and incompressible. The NavierStokes equations for the flow in Cartesian coordinates are:
where all variables are nondimensional, ${u}^{*}$ and ${v}^{*}$ represent the velocity components in the Cartesian coordinate directions ${x}^{*}$ and ${y}^{*}$, respectively. ${p}^{*}$ is the pressure, $Re$ is the Reynolds number ($\rho {U}_{\infty}D/{\mu}_{\infty}$) where $D$ is the characteristic length of the domain. ${\rho}_{\infty}$ and ${U}_{\infty}$ are the reference density and mean velocity, respectively. The Cartesian coordinates x* and y* are nondimensional with respect to $D$, ${u}^{*}$ and ${v}^{*}$ are nondimensional with respect to ${U}_{\infty}$, ${p}^{*}$ is nondimensional with respect to the dynamic pressure $\rho {U}_{\infty}^{2}$.
2.3. Upwind local RBFDQ method
The DQ method is a numerical discretization technique for approximation of derivatives. The essence of the DQ method is that the partial derivative of an unknown function with respect to an independent variable is approximated by a weighted linear sum of function values at all discrete points. Suppose that a function $f\left(x\right)$ is sufficiently smooth. Then its $m$th order derivative with respect to $x$ at a point ${x}_{i}$ can be approximated by DQ as:
where ${x}_{j}$ are the discrete points in the domain, $f\left({x}_{j}\right)$ and ${w}_{ij}^{\left(m\right)}$ are the function values at these points and the related weighting coefficients. Obviously, the key procedure in the DQ method is the determination of the weighting coefficients ${w}_{i,j}^{\left(m\right)}$.
Due to excellent performance of RBFs in the function approximation, they are chosen as the base functions to determine the weighting coefficients ${w}_{i,j}^{\left(m\right)}$ in this work. There are many RBFs (expression of $\phi $) available. The most commonly used RBFs are:
(a) Multiquadratics (MQs): $\phi \left(r\right)=\sqrt{{r}^{2}+{c}^{2}}$, $c>0$;
(b) Thinplate splines (TPS): $\phi \left(r\right)={r}^{2n}\mathrm{l}\mathrm{o}\mathrm{g}\left(r\right)$, $n$ is an integer;
(c) Gaussians: $\phi \left(r\right)={e}^{c{r}^{2}}$, $c>0$;
(d) Inverse MQs: $\phi \left(r\right)=\frac{1}{\sqrt{{r}^{2}+{c}^{2}}}$, $c>0$,
where $r=\Vert x{x}_{k}\Vert $ represents the Euclidean norm, and $c$ is a free shape parameter. The shape parameter $c$ defined in the MQs and Gaussian basis functions has an important role in the approximation of the meshless methods. Although analysis are available for the shape parameter [24], there is no mathematical theory on choosing its optimal value. In the following section, we will to determine the value by using use TPS RBF which does not require a shape parameter and its solution is not dependent on this parameter. Navier–Stokes equations are of second order, and we want to improve and extend this method to 3D problems, so we must use a higher order TPS [25, 26] and the integer $n$ will be chosen as 2.
Substituting the set of TPS RBF into Eq. (2), the determination of corresponding coefficients for the firstorder derivative is equivalent to solving the following linear equations:
(for simplicity, the notation ${\phi}_{k}\left(X\right)$ is adopted to replace $\phi (\Vert x{x}_{k}\Vert )$) or in the matrix form:
Therefore, if the collocation matrix $\left[A\right]$ is nonsingular, the coefficient vector [19] can be obtained by:
The nonsingularity of the collocation matrix $\left[A\right]$ depends on the properties of RBFs used. According to the work of [25], the matrix $\left[A\right]$ is conditionally positive definite for TPS RBFs. This fact guarantees the nonsingularity of matrix $\left[A\right]$ for distinct supporting points. The coefficient vector $\left\{w\right\}$ can be used to approximate the firstorder derivative in the $x$direction for any locally smooth function at node ${x}_{I}$. The calculation of weighting coefficients for other derivatives can follow the same procedure.
It is clear that the implementation of the DQ method becomes difficult when an irregular region must be considered. The essential reason behind this is the method lacks the upwind mechanism. Therefore, DQ method using to calculate the flow field characteristics with high Reynolds number will bring error. To overcome the difficulties and drawbacks, a local RBFDQ method with upwind mechanism is utilized in this paper. Unlike the RBFDQ method, derivative of a function with respect to a coordinate direction at an internal grid point is expressed as a weighted liner sum of the function values at partial grid points along that coordinate direction rather than at all grid points.
Fig. 2The local support domain model with upwind scheme (al=3, ar=2)
Local RBFDQ method uses local support domain point function to approximate the partial derivatives of a function:
$n$ is the number of the local support domain points.
${a}_{l}$ and ${a}_{r}$ are the upstream and downstream support domain coefficients. ${a}_{l}>{a}_{r}$, ${d}_{x}$ is the mean interval in the $x$ direction.
The local support domain model with upwind scheme is considered as shown in Fig. 2.
2.4. Expansion of the NavierStokes by upwind local RBFDQ method
By using the upwind RBFDQ method described above, Eq. (1) can be discretized as the following formula:
The traditional pressurecorrection methods are widely used to solve NS equations. An example of this approach is the SIMPLE algorithm [27]. This method integrates the NavierStokes equations in time at each timestep by firstly solving the momentum equations using an approximate pressure field to yield an intermediate velocity field that will not, in general, satisfy continuity. A Poisson equation is then solved with the divergence of the intermediate velocity as a source term to provide a pressure correction, which is then used to correct the intermediate velocity field, providing a divergence free velocity. The pressure is updated and integration then proceeds to examine the convergence. If the solution is converged the results are finalized. If not, the iterative process continues from the momentum equation with the new pressure and velocity field.
Kassab and Divo [28] combined the pressurecorrection method with the radial basis collocation method to solve heat diffusion–advection problems. Their approach, which is similar to the SIMPLE algorithm, is summarized as follows: The solution algorithm starts with an initial velocity field, which satisfies continuity, then the velocity field is introduced in the NS equations, and the time derivatives are approximated using backward differencing, transforming the NS to the Helmholtz Poisson equation set for estimating the velocity field. Once the Helmholtz Poisson problem is solved, the velocity field is updated and forced to satisfy continuity. An equation to update the pressure field is obtained by taking the divergence of the momentum equation and applying the continuity equation. The Poisson equation for the pressure field can be solved by imposing a proper and complete set of boundary conditions.
Our approach is distinct from above, instead, we use a leastsquare method, which is inherently iterative, to find primitive variables ($u$, $v$ and $p$) directly. The NS and the boundary condition equations together constitute a system of nonlinear equations, which are to be solved. These nonlinear equations are then solved by one of the leastsquare methods, LevenbergMarquardt. At the end we obtain the solutions of the velocity and pressure field directly. This procedure is described in detail below.
The NavierStokes and boundary condition equations together constitute a set of nonlinear equations, which are to be solved. The expanded velocity and pressure terms and their first and second derivatives are directly inserted into the equations. The boundary condition equations are collocated at the boundary points $(i=1,\dots ,NB)$ and the continuity and momentum equations are collocated at the interior domain points $\left(i=NB+1,\dots ,NB+NI\right)$. Then a nonlinear system of equations can be solved by the leastsquare optimization techniques such as the LevenbergMarquardt method [29]. The LevenbergMarquardt algorithm is an iterative technique that locates the minimum of a multivariate function that is expressed as the sum of squares of nonlinear realvalued functions. It has became a standard technique for nonlinear leastsquares problems. The reader is referred to the literature [29] for further details of this method.
3. Numerical results and discussions
Our objective in the paper is to demonstrate the feasibility of the direct solution approach using the upwind local RBFDQ method. Once the feasibility established in both 2D and 3D cases, we will set about optimization of the technique and code. The parameters are chosen as $Re=100$, $200$ and $300$, the three cases were chosen after systematic analysis of a number of examples over a wide range of $Re$. The 2D channel flow predictions are compared with those obtained using the finite volume method (FVM) which is embodied in the CFDFluent computational package. In order to perform error analysis, the relative percent error ($\epsilon $) is calculated for each node, and then the arithmetic mean ($\mu $) and the standard deviation ($\sigma $) of the relative error are calculated for each case, using the following expressions:
in which $X$ is any generic variable (velocity, pressure, stream function).
3.1. Girdindependence study
In the present study, a number of different grid sizes are considered in order to check the sufficiency of the mesh resolution. Table 2 depicts the comparisons of error and CPU times between RBFDQ method and CFDFluent method in the different mesh density, three different mesh sizes are used for discretization in the three flow types. It can be clearly seen that the error $\mu $ and $\sigma $ decrease as the node density increases.The upwind codes have a tendency to loose stability at higher distance between nodes.
As expected, the CPU time for the meshless method depends on the number of nodes employed, as well as the LevenbergMarquardt convergence criteria. However, a comparison of the results (following section) indicates that the node numbers in Table 2 provide sufficiently accurate results. Considering that there is still some improvement for optimization of the RBF code in the future, the computational cost as indicated in Table could be considered comparable to the commercial codes.
Consider the error and CPU time, the distribution of nodes in the domain for the flow situations considered is shown in Fig. 3. Fig. 3a depicts 2D channels with one obstacle on the lower wall, Fig. 3b represents 2D channels with two obstacles on the lower wall and Fig. 3c involves flow through a 2D channel with three obstructions in staggered pattern. The obstruction is located midway between the inlet and the outlet. The nondimensional length of the obstruction is 1.4, and the height is 0.5 in all cases. Fig. 3 shows that for 2D channels flow with one obstacle, a total of 86 collocation points are used on the boundary and 186 collocation points are used in the interior for RBFDQ expansion. For channel flow with two obstructions, we used a total of 118 collocation points on the boundary and 246 collocation points in the interior. A total of 150 collocation points are used on the boundary and 306 collocation points in the interior channel with three obstructions.
Table 2The comparisons of error and CPU times between RBFDQ method and CFDFluent method in the different mesh density (Re=100)
Flow type  Boundary nodes  Interior nodes  Total nodes  CPU time (s)  CFDFluent (s)  Error  
$\mu $ (%)  $\sigma $ (%)  
Channel flow with one obstruction  43  93  136  380  260  8.90  12.6 
86  186  272  405  280  1.61  3.25  
172  372  544  456  405  1.24  2.51  
Channel flow with two obstructions  59  123  184  412  275  9.65  15.61 
118  246  364  435  310  3.57  4.23  
236  492  728  512  396  3.18  3.96  
Channel flow with three obstructions  75  153  228  434  328  10.51  18.44 
150  306  456  510  350  4.78  5.06  
300  612  912  582  471  4.25  4.81 
Fig. 3Node distribution for 2D channel flow considered
3.2. Influence of upwind support domain on the numerical results
The upwind mechanism is important to the numerical computation of fluid flow. Table 3 shows the comparison of error in the different upwind support domain coefficients, it is shown that the upwind support domain can improve the accuracy of the numerical result. As $Re$ increases, the upstream support domain coefficient is bigger than the downstream support domain coefficient.
Table 3The comparison of error in the different upwind support domain coefficients (channel flow with one obstruction)
Support domain coefficients  $Re=100$  $Re=200$  $Re=300$  
$\mu $ (%)  $\sigma $ (%)  $\mu $ (%)  $\sigma $ (%)  $\mu $ (%)  $\sigma $ (%)  
${a}_{l}=2$, ${a}_{r}=2$  4.56  6.98  6.28  7.94  7.27  9.25 
${a}_{l}=3$, ${a}_{r}=2$  1.61  3.25  2.38  3.77  3.45  4.89 
${a}_{l}=4$, ${a}_{r}=2$  4.28  6.54  4.83  5.22  3.10  4.34 
Fig. 4Predicted velocity vector for 2D channel flow
(a) 2D channel flow with one obstruction
(b) 2D channel flow with two obstructions
(c) 2D channel flow with three obstructions
Fig. 5Predicted pressure distribution along the lower wall for 2D channel flow
(a) 2D channel flow with one obstruction
(b) 2D channel flow with two obstructions
(c) 2D channel flow with three obstructions
3.3. Result on flow characteristics
Table 4 shows the comparison of error between upwind local RBFDQ method and CFDFluent method, it is seen that the numerical solutions converge for varying Reynolds numbers. The predicted velocity vectors for 2D channel flow with different obstructions are shown in Fig. 4a, b and c, respectively. The velocity vectors become crowded above the obstruction, indicating flow acceleration in the narrow section. The minimum velocity occurs just downstream of the obstruction, creating a dead zone. Subsequently the flow starts to develop and recover to its original profile along the channel. As $Re$ increases, the recirculation in this dead zone increases.
Table 4A comparison of error between RBFDQ method and CFDFluent method
Flow type  $Re=100$ $({a}_{l}=3,{a}_{r}=2)$  $Re=200$ $({a}_{l}=3,{a}_{r}=2)$  $Re=300$ $({a}_{l}=4,{a}_{r}=2)$  
$\mu $ (%)  $\sigma $ (%)  $\mu $ (%)  $\sigma $ (%)  $\mu $ (%)  $\sigma $ (%)  
Channel flow with one obstruction  1.61  3.25  2.38  3.77  3.10  4.34 
Channel flow with two obstructions  3.57  4.23  4.16  4.59  4.46  5.64 
Channel flow with three obstructions  4.78  5.06  4.96  5.68  5.46  5.89 
The predicted pressure distributions along the lower wall for channel flow with different obstructions are presented in Fig. 5 for both RBFDQ and CFDFluent models. The two sets of results are in good agreement. The wall pressure decreases slowly from the channel entrance to the obstruction as the flow resistance. Then it drops rapidly on the obstruction due to the bigger velocity gradient on the obstruction. Subsequently the wall shear stress increases slowly as the flow recovers towards the outlet.
4. Conclusions
In this study, we have developed a novel method based on the upwind local RBFDQ method for direction solution of the Navier–Stokes equations. We validated the proposed method by studying the classical flow situations, namely, the 2D channel flow with different obstructions. The predicted results are then compared with the well validated commercial code for the channel flow.
Numerical results show that the upwind local RBFDQ method is well suited to solve irregular regions. Although the method has been developed in the context of twodimensional problem, it can be extended to threedimensional problem is fairly straightforward. It is hoped that the insights presented in this paper will help to spur more interest in and launch more investigations into the use of RBFDQ method for the numerical flow simulation.
References

J. Monaghan A turbulence model for smoothed particle hydrodynamics. European Journal of MechanicsB/Fluids, Vol. 30, 2011, p. 360370.

K. Selim, A. Logg, M. G. Larson An adaptive finite element splitting method for the incompressible Navier–Stokes equations. Computer Methods in Applied Mechanics and Engineering, 2011.

N. Nguyen, J. Peraire, B. Cockburn An implicit highorder hybridizable discontinuous Galerkin method for the incompressible Navier–Stokes equations. Journal of Computational Physics, Vol. 230, 2011, p. 11471170.

P. Guan, S. Chi, J. Chen, T. Slawson, M. Roth SemiLagrangian reproducing kernel particle method for fragmentimpact problems. International Journal of Impact Engineering, Vol. 38, 2011, p. 10331047.

J. Xu, S. Rajendran A partitionofunity based ‘FEMeshfree’QUAD4 element with radialpolynomial basis functions for static analyses. Computer Methods in Applied Mechanics and Engineering, Vol. 200, 2011, p. 33093323.

C. A. Duarte, J. T. Oden Hp cloudsan hp meshless method. Numerical methods for partial differential equations, Vol. 12, 1996, p. 673706.

F. Perazzo, R. Löhner, L. PerezPozo Adaptive methodology for meshless finite point method. Advances in Engineering Software, Vol. 39, 2008, p. 156166.

W. L. Nicomedes, R. C. Mesquita, F. J. S. Moreira The meshless local Petrov–Galerkin method in twodimensional electromagnetic wave analysis. IEEE Transactions on Antennas and Propagation, Vol. 60, 2012, p. 19571968.

B. AtaieAshtiani, L. Farhadi A stable movingparticle semiimplicit method for free surface flows. Fluid Dynamics Research, Vol. 38, 2006, p. 241256.

T. Liszka An interpolation method for an irregular net of nodes. International Journal for Numerical Methods in Engineering, Vol. 20, 2005, p. 15991612.

E. J. Kansa Multiquadrics–A scattered data approximation scheme with applications to computational fluiddynamics–I surface approximations and partial derivative estimates. Computers & Mathematics with Applications, Vol. 19, 1990, p. 127145.

E. J. Kansa Multiquadrics–A scattered data approximation scheme with applications to computational fluiddynamics–II solutions to parabolic, hyperbolic and elliptic partial differential equations. Computers & Mathematics with Applications, Vol. 19, 1990, p. 147161.

G. E. Fasshauer Solving partial differential equations by collocation with radial basis functions. Proceedings of Chamonix, Vanderbilt University Press, Nashville, TN, 1996, p. 18.

G. Demirkaya, C. Wafo Soh, O. Ilegbusi Direct solution of Navier–Stokes equations by radial basis functions. Applied Mathematical Modelling, Vol. 32, 2008, p. 18481858.

T. A. Driscoll, B. Fornberg Interpolation in the limit of increasingly flat radial basis functions. Computers & Mathematics with Applications, Vol. 43, 2002, p. 413422.

Y. Hon, Z. Wu A quasi‐interpolation method for solving stiff ordinary differential equations. International Journal for Numerical Methods in Engineering, Vol. 48, 2000, p. 11871197.

W. Chen, M. Tanaka A meshless, integrationfree, and boundaryonly RBF technique. Computers & Mathematics with applications, Vol. 43, 2002, p. 379391.

C. Chen, C. Brebbia The dual reciprocity method for Helmholtztype operators. Boundary elements, Vol. 10, 1998, p. 495504.

S. Xiang, G. Li, W. Zhang, M. Yang A meshless local radial point collocation method for free vibration analysis of laminated composite plates. Composite Structures, Vol. 93, 2011, p. 280286.

W. Zongmin HermiteBirkhoff interpolation of scattered data by radial basis functions. Approximation Theory and its Applications, Vol. 8, 1992, p. 110.

C. W. Bert, M. Malik Differential quadrature method in computational mechanics: a review. Applied Mechanics Reviews, Vol. 49, 1996.

C. Shu, H. Ding, H. Chen, T. Wang An upwind local RBFDQ method for simulation of inviscid compressible flows. Computer Methods in Applied Mechanics and Engineering, Vol. 194, 2005, p. 20012017.

J. Sun, Z. Zhu Upwind local differential quadrature method for solving incompressible viscous flow. Computer Methods in Applied Mechanics and Engineering, Vol. 188, 2000, p. 495504.

A. H. D. Cheng, M. Golberg, E. Kansa, G. Zammito Exponential convergence and H‐c multiquadric collocation method for partial differential equations. Numerical Methods for Partial Differential Equations, Vol. 19, 2003, p. 571594.

M. D. Buhmann Radial basis functions: theory and implementations. Cambridge University Press, 2003.

M. Griebel, M. A. Schweitzer Meshfree methods for partial differential equations. Springer Verlag, 2003.

S. Patankar Numerical Heat Transfer and Fluid Flow. Hemisphere, Washington, DC, 1980, The FLUENT (trademark of FLUENT Inc.) fluid dynamics analysis package has been used in this study, 1980.

E. Divo, A. J. Kassab A meshless method for conjugate heat transfer problems. Engineering Analysis with Boundary Elements, Vol. 29, 2005, p. 136149.

K. Madsen, H. B. Nielsen, J. Søndergaard Robust Subroutines for NonLinear Optimization. IMM, Informatics and Mathematical Modelling, Technical University of Denmark, 2002.
About this article
This research was supported by the Natural Science Foundation of China (No. 11302133) and Education Fund Item of Liaoning Province (No. L2013071). I am grateful to Dr. Tianzhi Yang for reading the manuscript and for helpful suggestions.