Localized method of approximate particular solutions for solving unsteady Navier–Stokes problem

Accepted Manuscript

Localized method of approximate particular solutions for solving unsteady Navier-Stokes problem Xueying Zhang, Muyuan Chen, C.S. Chen, Zhiyong Li PII: DOI: Reference:

S0307-904X(15)00582-X 10.1016/j.apm.2015.09.048 APM 10747

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

7 August 2013 3 June 2015 23 September 2015

Please cite this article as: Xueying Zhang, Muyuan Chen, C.S. Chen, Zhiyong Li, Localized method of approximate particular solutions for solving unsteady Navier-Stokes problem, Applied Mathematical Modelling (2015), doi: 10.1016/j.apm.2015.09.048

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • The LMAPS is proposed to solve Navier-Stokes equations. • A global and sparse matrix is reformulated.

AC

CE

PT

ED

M

AN US

• Fractional step algorithm is adopted.

CR IP T

• The weighting coefficients are obtained by local MQ-RBFs.

1

ACCEPTED MANUSCRIPT

Localized method of approximate particular solutions for solving unsteady Navier-Stokes problem

CR IP T

Xueying Zhanga, 1 , Muyuan Chenb , C. S. Chena,c , Zhiyong Lib a

College of Science, Hohai University, Nanjing, Jiangsu 210098, China Laboratory of Biomechanics, School of Biological Engineering and Medical Engineering, Southeast University, Nanjing, Jiangsu 210096, China c Department of Mathematics, University of Southern Mississippi, Hattiesburg, MS 39406, U.S.A

AN US

b

Abstract

PT

ED

M

The localized method of approximate particular solutions (LMAPS) is proposed to solve two-dimensional transient incompressible Navier-Stokes systems of equations in primitive variables. The equations contain the Laplacian operator. In avoiding ill-conditioning problem, the weight coefficients of linear combination with respect to the function values and its derivatives can be obtained by solving low-order linear system within local supporting domain in which five nearest neighboring points and multiquadrics are used for interpolation. Then local matrices are reformulated in the global and sparse matrix. The obtained large sparse linear systems can be directly solved instead of using more complicated iterative method. The method is assessed on driven cavity problem and flow around cylinder. The numerical experiments show that the newly developed LMAPS is suitable for solving incompressible Navier-Stokes equations with high accuracy and efficiency.

AC

CE

Keywords: Navier-Stokes equations; Local method; Particular solution to Laplace equation; Radial basis functions(RBFs).

1

Corresponding author: Xueying Zhang, College of Science, Hohai University, Nanjing, Jiangsu 210098, China. E-mail address: [email protected]

Preprint submitted to Applied Mathematical Modelling

October 5, 2015

ACCEPTED MANUSCRIPT

1. Introduction

CR IP T

In the computational fluid dynamics realm, many computational methods have been successfully applied to solve the N-S equations, prominent among which are the finite difference method (FDM), the finite element method (FEM) and the finite volume method (FVM). However, these conventional numerical methods encountered many challenges as their numerical results are greatly dependent on mesh generation, which is a major hurdle for the high dimensional and complexgeometry problems.

PT

ED

M

AN US

In order to reduce the high cost of grid generation, a number of meshless methods have been developed in the past two decades. These methods include the smoothed particle hydrodynamics method (SPH)[1], the element-free Galerkin method [2], the finite point method [3] and the differential quadrature (DQ) method [4, 5], MLPG [6], the method of approximate particular solutions (MAPS) [7, 8, 9], etc. Among them radial basis functions (RBFs) have been actively investigated in various areas of applications. The RBF methods have the advantage of not requiring any particular grid system consequently yielding more flexibility when approximating the function in an irregular domain and yielding a fast convergence if the function considered is smooth enough. Initially, the RBFs can be dated from the fact that Hardy first successfully applied multiquadrics RBFs to function interpolation based on multivariate data [10]. However, their truly mesh-free property first motivated Kansa to deal with partial differential equations [11, 12]. Since then, an increasing number of researchers cast their sights on this class of numerical methods and applied them to solve linear and non-linear PDEs in various engineering and academic cases [13, 14, 15, 16].

AC

CE

It’s worth noting that RBFs are used in the collocation technique in global manner, so the resultant coefficient matrix may be subjected to singular, dense and even ill-conditioned problems which restricts our abilities to solve large-scale problems. Consequently a new class of effective constructions were proposed to overcome these shortcomings such as domain decomposition [23], the improved truncated singular valued decomposition [24], the localized formulations [25, 26], fast multipole expansion techniques [27], etc. Among these new techniques, local scheme based on local support interpolations appears to more efficient in handling a large number of collocation points and less sensitive to changes of shape parameter [9, 17, 18, 19, 20, 21, 22, 28, 29]. Motivated by the idea of local RBFs based meshless methods, we will present 3

ACCEPTED MANUSCRIPT

CR IP T

readers a novel method, the local method of approximate particular solutions (LMAPS), and apply it to solve two-dimensional, time-dependent, incompressible Navier-Stokes equations in the primitive variables form by employing the fractional step procedure. In this method the momentum equations and the Poisson equation for pressure are solved separately at each time step.

AN US

The organization of the paper is as follows. In section 2, the LMAPS is proposed to directly approximate the derivatives, in which RBFs are approximated by using local supporting nodes. Then the formulations of the LMAPS are listed for solving two-dimensional, time-dependent Navier-Stokes equations in detail. In section 3, numerical examples are given to demonstrate the effectiveness of the proposed procedure. Conclusions are drawn in Section 4. 2. Formulation of the LMAPS for unsteady N-S equations 2.1. Governing equations

∇ · u = 0,

(1)

1 2 ∂u + u · ∇u = −∇p + ∇ u, ∂t Re

(2)

ED

Continuity equation:

M

In this paper, we consider viscous, unsteady, incompressible flows which is governed by the following Navier-Stokes equations in primitive variable formulation.

PT

Momentum equation:

AC

CE

where p is the pressure, u is the velocity, u = (u, v), Re = v`/ν is Reynolds number which denotes the ratio of the inertial terms to the viscous terms. The computational domain is an unified space-time domain Ω with the initial and boundary conditions. v indicates a characteristic velocity, ` is a characteristic length and ν represents the kinematic viscosity of the fluid. When the Re increases, the inertial forces will predominate in the flow field. As a result, a thin boundary layer is developed with large gradient. It becomes increasingly difficult to compute nonlinear and instable flow in the region of large gradient.

4

ACCEPTED MANUSCRIPT

2.2. Local approximation of derivatives

CR IP T

The ill-conditioning of the interpolation matrix and the dense matrix have limited the application of global RBFs based meshless methods for solving largescale problems. To mitigate these difficulties, localized RBF schemes have been developed with less sensitivity to the selection of the shape parameter of RBFs and avoiding the ill-conditioning problem.

f (x p ) ' fˆ(x p ) = where

ns X j=1

AN US

Let x = (x1 , x2 ) and {x j }Nj=1 be the interpolation points inside a bounded and connected domain Ω. For each point x p ∈ Ω, we create an influence domain Ω p containing a region form by the n s closest neighboring interpolation points ns {x[p] j } j=1 to x p . Assume that the value of a function f (x p ) can be approximated by a linear combination of n s radial basis functions (RBFs) in the following form   [p] α[p] Φ kx − x k , p j j

xp ∈ Ωp,

∇2 Φ(r) = ϕ(r),

(3)

√

(4)

ED

M

2 2 r = kx p − x[p] j k, and ϕ is a RBF. In this paper, we choose MQ-RBF ( r + c ) as the basis function. MQ has been proved as an effective RBF and widely used for solving various types of problems. By direct integration from Eq. (4), we have [7]

Φ(r) =

 √ √ c3  1 2 4c + r2 r2 + c2 − ln c + r2 + c2 . 9 3

(5)

AC

CE

PT

It is known that the accuracy of the MQ-RBF based method greatly depends on the choice of the shape parameter of MQ. How to select the optimal shape parameter of MQ is still an outstanding research topic for global formulation. For local RBF methods, Vertnik and Sarler [17] proposed to adopt the dimensionless shape parameter for local RBF collocation method and produced excellent results. In their approach, the MQ-RBF is changed to the following form p cr0 )2 (6) ϕ(r) = r2 + (e

where r0 is the maximum distance from x p to all nodal points in the influence domain Ω p . In other words, r0 = max kx p − x[p] j k. 1≤ j ≤ n s

5

(7)

ACCEPTED MANUSCRIPT

Note that when the density of nodal points is high, r0 becomes small. Thus, we should choose large shape parameter and vice versa. By adopting such dimensionless shape parameter, the selection of shape parameter is much more predictable and easier to manage.

fˆ(x p ) = Ψ[p] α[p] where

CR IP T

We note that Eq. (3) can be rewritten in the following matrix form

(8)

h    i [p] Ψ[p] = Φ kx p − x[p] k , · · · , Φ kx − x k p ns 1

where fˆ[p]

AN US

[p] [p] T [p] n s and α [p] = [α[p] 1 , α2 , . . . , αn s ] . Since {x j } j=1 ⊂ Ω p , Eq. (3) holds for every x[p] j , 1 ≤ j ≤ n s . By collocation method, we have the following matrix form

Φ[p]α [p] = fˆ[p] , h  i ˆ [p] T [p] = Φ kx[p] − x[p] k = [ fˆ(x[p] j i 1 ), . . . , f (xn s )] , Φ

(9)

16i, j6n s

. It follows that

α [p] = (Φ[p] )−1 fˆ[p] .

M

(10)

PT

ED

The above equation shows that α[p] in each influence domain Ω p can be expressed as the product of inverse of interpolation matrix and function values at nodal points. In this way, we don’t have to store these coefficients α[p] , which is tedious, in each influence domain. Other than simplifying the solution process, Eq. (10) is necessary for the formulation of localized RBF methods. As we shall see, all the α[p] will be replaced by the right hand side of Eq. (10) in the whole solution process. Note that ∂ f (x p )/∂xk , k = 1, 2, and ∇2 f (x p ) can be approximated by

CE

ns  ∂ f (x p ) ∂ fˆ(x p ) X ∂Φ  [p] [p] = Ξ[p] ' = α[p] kx − x k p xk α , j j ∂xk ∂xk ∂x k j=1

xp ∈ Ωp,

(11)

xp ∈ Ωp,

(12)

AC

and

∇ f (x p ) ' ∇ fˆ(x p ) = 2

2

where

Ξ[p] xk

ns X j=1

  [p] α[p] ϕ kx − x k = Θ[p]α [p] , p j j

"

 # ∂Φ  ∂Φ  [p] [p] = kx p − x1 k , · · · , kx p − xns k , k = 1, 2, ∂xk ∂xk 6

ACCEPTED MANUSCRIPT

and

  i h  [p] k , · · · , ϕ kx − x k . Θ[p] = ϕ kx p − x[p] p ns 1

Substituting Eq. (10) into Eqs. (11)∼(12), we have

CR IP T

To be more specific, let x j = (x1j , x2j ). Then, ∂Φ/∂xk in Eq. (11) can be written as follows: !  1 xk − xkj ∂Φ  c3 2 2 r + 2c − , k = 1, 2. (13) kx − x j k = √ √ ∂xk 3 c + r2 + c2 r 2 + c2

AN US

∂ fˆ(x p )  [p] [p] −1  [p] = Ξ xk (Φ ) fˆ , k = 1, 2, ∂xk   ∇2 fˆ(x p ) = Θ[p] (Φ[p] )−1 fˆ[p] .

(14)

(15)

M

In the next subsection, we will channel the local formulation in Eqs. (14)∼(15) into the governing equation of the N-S equation and its boundary conditions to formulate a global sparse linear system of equations. 2.3. The discretization of N-S equations

PT

ED

For a time increment ∆t = t(n+1) −t(n) , discretization of the momentum equation Eq. (2) by Euler method gives # " 1 2 (n) (n+1) (n) (n) (n) (n) ∇u . (16) u = u + ∆t −∇p − u · ∇u + Re

AC

CE

Due to the lack of an independent equation for the pressure and non-existence of a dominant variable in the continuity equation, it is a challenge to directly solve above Navier-Stokes equations. In order to circumvent these difficulties, the socalled projection method or fractional step algorithm is the most popular which is originally proposed by Chorin [30]. The projection procedures for each interior node in the domain are described in the following steps: Firstly, the pressure gradient terms are dropped from the momentum equations. Meanwhile, a provisional velocity u∗ can be predicted by solving the resultant advection-diffusion equations " # 1 2 (n) ∗ (n) (n) (n) (17) u = u + ∆t −u · ∇u + ∇u . Re 7

ACCEPTED MANUSCRIPT

CR IP T

where (n) is the time step index. Secondly, the velocity components are corrected by the pressure gradient term so as to enforce the satisfaction of the continuity equation by solving   u(n+1) = u∗ + ∆t −∇p(n+1) . (18)

At each time step, the solution u(n+1) is expected to satisfy the continuity equation Eq. (1), which implies ∇ · u(n+1) = 0. (19) Taking the divergence of Eq. (18), we obtain   ∇ · u(n+1) = ∇ · u∗ − ∆t∇p(n+1) .

(20)

AN US

Substituting into the resultant Eq. (19), the Eq. (20) can be simplified to pressure Poisson equation

1 ∇ · u∗ . (21) ∆t can be obtained by solving Eq. (21) with Neumann boundary ∇2 p(n+1) =

∂p(n+1) −n = 0 ∂→

(22)

p(n+1) = p0

(23)

M

The pressure p(n+1) condition

ED

or Dirichlet boundary condition

CE

PT

−n is the normal direction to the boundary surface. Since the boundary where → values of velocity components are available, boundary condition Eq. (21) for pressure can be simplified to [16]  1  ∗ ∂p(n+1) (n+1) = u − u . (24) −n ∆t ∂→

AC

When we obtain intermediate velocity u∗ and the pressure p(n+1) by Eqs. (17) and (21), respectively. The new velocity u(n+1) at the next time step can be determined by Eq. (18). For two dimensional N-S equations, let u = (u, v). Linear coefficients obtained by Eqs. (14)∼(15) are used to discrete the corresponding spatial derivatives in Eqs. (17)∼(21), it yields   ns ns ns X X X   1 (n)  (n) (n) (n)  , −u(n) l u (25) u∗i = u(n) a u − v b u + + ∆t i j i j i j j  i j i j   i Re j=1 j=1 j=1 8

ACCEPTED MANUSCRIPT

u(n+1) = u∗i − ∆t i

n  ns s X  1 X ∗ ∗  ai j u j + = bi j v j  . ∆t j=1 j=1

ns X j=1

ai j p(n+1) , v(n+1) = v∗i − ∆t j i

ns X

(26) bi j p(n+1) , j

(27)

CR IP T

∇2 p(n+1) i

j=1

where n s is the number of the interpolation points inside the influence domain Ωi for any node xi ∈ Ω, i = 1, 2, . . . , N.

AN US

We note that the coefficients li j , ai j and bi j in Eqs. (25)∼(27) can be calculated by the coefficient matrices in Eqs. (12)∼(13). Once the weighting coefficients are obtained, we reformulate local coefficient matrices in global and sparse form by inserting zeros in proper positions. For further details of the reformulation, we refer readers to Reference [9]. Thus, reformulating local coefficient matrices, we have the following discretization for advection-diffusion equation (25) 1 dU ∗ = LU (n) − (UA + VB) U (n) , dt Re

(28)

M

where L, A and B are N × N sparse matrices, U and V are N × N diagonal matrices in which all of the diagonal elements come from the elements of the vectors U (n) and V (n) . Further, we obtain P(n+1) from Eq. (26) as follows

ED

P(n+1) =

1 −1 L (A + B)U ∗ . ∆t

(29)

PT

Ultimately, Eqs. (27) can be transformed into matrix formulation U (n+1) = U ∗ − ∆tAP(n+1) , V (n+1) = V ∗ − ∆tBP(n+1) ,

(30)

CE

where U and V are N × 1 vectors. The resultant equations can be solved with implementation of accurate boundary conditions. The general solution procedures of the LMAPS for N-S equations are shown below

AC

1. Discretize the governing equations (17)∼(21) into Eqs. (25)∼(27) with the computed weighting coefficients for the related derivatives, then reformulate them in global form, see Eqs. (28)∼(30). 2. Suppose the initial-boundary conditions of velocity u and pressure p are known at t(n) , solve the time-dependent advection-diffusion equation (28) to obtain the intermediate velocity u∗ . Determine the updated boundary values for intermediate velocity u∗ . 9

ACCEPTED MANUSCRIPT

CR IP T

3. Solve the discretized pressure Poisson equation (29) to obtain the pressure p(n+1) . 4. Solve equation (30) to obtain new velocity u(n+1) , also update the velocity boundary conditions; 5. Go to step 2 until a required time step is reached. 3. Numerical results 3.1. Lid-driven cavity flow

AN US

In this section, the proposed LMAPS is applied to solve a unsteady incompressible Navier-Stokes problem for Lid-driven cavity flow and flow around a circular cylinder to examine the efficiency and accuracy of the method. As has been pointed out in [28], local five-noded subdomains are used for each interior node of uniform nodal distributions, which shows high accuracy and practicability, especially near the boundary. So overlapping five-noded subdomains and uniform nodal points are used in this example and next. 1

M

1 0.9 0.8

0.9 0.8 0.7

ED

0.7

y

0.5

0.3 0.2 0.1

0.2

CE

0 0

PT

0.4

0.4

0.6

y

0.6

0.5 0.4 0.3 0.2 0.1

0.6

0.8

0 0

1

x

0.4

x

0.6

0.8

1

(b) Streamlines, Re = 1000

AC

(a) Sreamlines, Re = 100

0.2

Figure 1: The streamlines of Lid-driven flow in square cavity for Re=100 and Re=1000.

The computational domain Ω ∪ ∂Ω for the problem is 0 ≤ x ≤ 1, 0 ≤ y ≤ 1 and 0 < t < +∞. Ω ∪ ∂Ω is a regular domain where the interior and boundary nodes are uniformly distributed. The boundary conditions of the Eqs. (28)∼(30) are subjected to no slip boundary conditions, i.e. on top wall u = 1, v = 0, and 10

ACCEPTED MANUSCRIPT

. +

*

0.8

LMAPS, Re=100, 101x101 LMAPS, Re=100, 51x51 LMAPS, Re=100, 17x17

. +

0.15

*

0.1

LMAPS, Re=100, 101x101 LMAPS, Re=100, 51x51 LMAPS, Re=100, 17x17

AN US

1

CR IP T

on all other walls u = 0, v = 0. The pressure boundary condition for the Eq. (21) is taken as Neumann boundary condition ∂∂p → −n = 0 on all walls. Similar to the discretization of the governing equations, each boundary node is also discretized using more nodes within a 3r0 radius of local supporting subdomain as defined in Eq. (7) for r0 . Since the boundary values of velocity components are available, the updating of boundary condition for pressure can be derived from Eq. (29). As the Reynolds number increases, the higher velocity gradients near the boundary walls result in vorticity generations at more number of boundary points.

0.05

0.6

V

U

0

0.4

−0.05

0.2

−0.1

−0.15

0

−0.2

0

0.2

0.4

0.6 Y

0.8

1

0

M

−0.2

0.8

0.3

LMAPS Re=1000, 101x101

0.2

CE

0

AC

0

0.2

+ LMAPS,Re=1000,101x101 LMAPS,Re=1000,51x51

0

V

−0.1

−0.2

−0.3

−0.2

−0.4

1

0.1

PT

U

0.4

0.8

.

0.2

LMAPS Re=1000, 51x51

0.6

0.6

(b) v-velocity, Re = 100

ED

. +

0.4

X

(a) u-velocity, Re = 100 1

0.2

−0.4

0.4

0.6

0.8

−0.5

1

Y

0

0.2

0.4

0.6

0.8

1

X

(c) u-velocity, Re = 1000

(d) v-velocity, Re = 1000

Figure 2: u-velocity component along vertical line x = 0.5 and v-velocity component along horizontal line y = 0.5 through the geometric center.

We have solved the incompressible flow in a square driven cavity numerically. 11

ACCEPTED MANUSCRIPT

CR IP T

Time step ∆t is takes as 0.00005 for 51 × 51 uniform nodal distributions and Re = 100. The shape parameter c is taken as 20r0 . The streamlines are shown in the Figure 1 for Re = 100 and Re = 1000. Figure 2 shows the velocity components for u along the vertical centerline and v along horizontal centerline for 17 × 17, 51 × 51 and 101 × 101 uniform nodal distributions. As shown in Figure 2 (a) and (b) for Re = 100 and Figure 2 (c) and (d) for Re = 1000, the accuracy seems to be improved as nodal distributions become more denser. The minimum and maximum velocities through the center of the cavity, calculated with the LMAPS are compared to the results obtained Ghia [31] for Re = 1000.

ED

M

AN US

The minimum u component, the minimum and maximum v component of velocities with the LMAPS are basically consistent with those of Mramor [21], Ghia et al [31] as shown in the Figure 3. In comparison with the data from above literatures, the accuracy of the minimum and the maximum values of velocities along horizontal and vertical lines that pass through the geometric center of the cavity is slightly poorer as shown in Figure 3. This may result from fewer nodes and different boundary treatment. We note that the increase in the Reynolds number from 100 to 1000 results in sharp u and v velocity gradients near the thin boundary layers. It can be observed from Figure 1 that the axis of the primary vortex starts in the upper right half cavity in a stable state, then gradually moves towards the square center as the Reynolds number increases. 1

0.4

0.3

Ghia 129x129 Re=1000 LMAPS 51x51 Re=1000

0

CE

U

0.4

0.1

V

0.6

Ghia 129x129 Re=1000 LMAPS 51x51 Re=1000

0.2

PT

0.8

0.2

−0.1

−0.2

−0.3

0

−0.4

AC

−0.2

−0.4 0

0.1

0.2

0.3

−0.5

0.4

0.5 Y

0.6

0.7

0.8

0.9

−0.6 0

1

(a) u-velocity component along x = 0.5

0.1

0.2

0.3

0.4

0.5 X

0.6

0.7

0.8

0.9

1

(b) v-velocity component along y = 0.5

Figure 3: Comparison of velocities for Re = 1000, calculated with LMAPS and results by Ghia [31]

12

ACCEPTED MANUSCRIPT

2

1.8

1.8

1.6

1.6

1.4

1.4

AN US

2

1.2

1.2

y

1

y

CR IP T

When the aspect ratio of the driven cavity becomes 1 : 2, two obvious vorticities are generated at more number of points along all the three sides apart from the top side as shown in Figure 4. It is noted that the upper eddy gradually moves towards the cavity center, and the lower eddy extends progressively to most of the lower cavity for the case of Re = 100 and Re = 1000.

1

0.8

0.8

0.6

0.6

0.4

0.2 0 0

1

ED

0.5 x

M

0.4

(a) Streamlines, Re = 100

0.2 0 0

0.5 x

1

(b) Streamlines, Re = 1000

PT

Figure 4: The streamlines of Lid-driven flow in rectangular cavity with Re=100 and Re=1000.

CE

3.2. Flow around a circular cylinder

AC

In this subsection, the proposed LMAPS is applied to compute the unsteady flow around a circular cylinder which has served as a test case for various numerical approaches. The configuration of the problem is shown in Figure 5. No slip boundary conditions are set at all the solid surfaces. A parallel flow with unit speed is specified on inflow boundary, while right boundary subjects to outflow boundary conditions. The no-slip boundary condition is set at the rest of the boundaries of the problem domain. The pressure boundary condition is taken as Neumann 13

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 5: Configuration of the flow around a circular cylinder.

ED

M

boundary condition, expect that the outlet is given as Dirichlet boundary condition p = 0. It is complex flow phenomena including flow separation, wake flow and vortex-shedding whose properties and shapes change with the Reynolds number. The value of the shape parameter c is set to 20r0 . The streamlines are plotted in Figure 6 for Re = 30 and Re = 100. It is observed that two symmetrical vortexes are generated behind the circular cylinder which are not separated till steady state for the case of Re = 30 as shown in Figure 6(a). As the Reynolds number is increased to 100, typical Karmann Streets are well organized, and coherent vortexshedding is exhibited in Figure 6(b). The boundary layer is extremely thin and laminar in the attached flow region. The separation point is delayed, and the wake is narrower with the increase of the Reynolds number.

PT

4. Concluding remarks

AC

CE

In this paper, the LMAPS has been proposed to solve two-dimensional, timedependent, incompressible Navier-Stokes equations. In order to overcome the instability of the ill-posed problem, the weight coefficients can be obtained by solving the low order linear systems within local supporting domain. A two-step fractional step algorithm is required to circumvent the difficulties arising from lack of an independent equation for the pressure. The pressure Poisson equation can be directly solved by the LMAPS throughout the whole domain instead of conventional Successive Over-Relaxation (SOR) method. Two classic examples, namely Lid-driven cavity flow and flow around a circular cylinder are implemented to validate the effectiveness of the proposed method. In general, the presented LMAPS appears to be attractive for being used to solve Navier-Stokes system containing 14

ACCEPTED MANUSCRIPT

8

4 2 2

4

6

8

10

12

14

16

8

Y

6 4 2 4

6

8

20

10

12 X

14

16

18

20

22

22

M

2

18

AN US

(a) Streamlines, Re=30

CR IP T

6

(b) Streamlines, Re=100

ED

Figure 6: The streamlines of flow around a circular cylinder for Re=30 and Re=100.

AC

CE

PT

second-order Laplacian operator. However, some issues still remain pending. For example, the choice of optimal value of shape parameter relies on the tests or experiences. The robust implementation for the boundary conditions of pressure Poisson equation and intermediate velocity still needs further study. Moreover, the accuracy by the LMAPS may not be as accurate as that by high precision finite difference method. In future work, it is hoped that the LMAPS can be extended to complex unsteady flows based on the upwind strategy, specially using appropriate nonuniform nodal distribution and large Reynolds numbers. Acknowledgement The first author thanks the support of the Scientific Research Foundation for Returned Overseas Personnel (No.20155003412). The third author thanks the support of Fellowship for Distinguished Overseas Scholar by the Ministry of Ed15

ACCEPTED MANUSCRIPT

CR IP T

ucation of China. The fourth author thanks the partially supported of the National 973 Basic Research Program of China [No. 2013CB733800] and the National Natural Science Foundation of China (NSFC) (No. 11272091). References

[1] L.B. Lucy, A numerical approach to the testing of the fission hypothesis, Astro. J., 8 (1977) 1013–1024. [2] T. Belytschko, Y.Y. Lu, L. Gu, Element-free Galerkin methods, Int. J. Numer. Meth. Engrg., 37 (1994) 229–256.

AN US

[3] E. Onate, S. Idelsohn, O.C. Zienkiewicz, R.L. Taylor, A finite point method in computational mechanics, Application to convective transport and fluid flow, Int. J. Numer. Meth. Engrg., 39 (1996) 3839–3866.

M

[4] C. Shu, H. Ding, K.S. Yeo, Local radial basis function-based differential quadrature method and its application to solve two-dimensional incompressible Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg., 192 (2003)941-954.

ED

[5] C. Shu, Q. Yao, K.S. Yeo, Block-marching in time with DQ discretization: an efficient method for time-dependent problems, Comput. Methods Appl. Mech. Engrg., 191 (2002) 4587–4597.

PT

[6] S.N. Atluri, T. Zhu, A new meshless Local Petrov-Galerkin (MLPG) approach in computational mechanics, Computational Mechanics, 22 (1998) 117–127.

CE

[7] C.S. Chen, C.M. Fan, P.H. Wen, The method of particular solutions for solving elliptic problems with variable coefficients, International Journal of Computational Methods, 8 (2011) 545–559.

AC

[8] P.H. Wen, C.S. Chen, The method of particular solutions for solving scalar wave equations, International Journal for Numerical Methods in Biomedical Engineering, 26 (2010) 1878–1889. [9] Guangming Yao, C.S. Chen, Joseph Kolibal, A localized approach for the method of approximate particular solutions, Computers and Mathematics with Applications, 61 (2011) 2376–2387. 16

ACCEPTED MANUSCRIPT

[10] R.L. Hardy, Multiquadric equations of topography and other irregular surfaces, J. Geophys. Res., 76 (1971) 1905–1915.

CR IP T

[11] E.J. Kansa, Multiquadrics, A scattered data approximation scheme with applications to computational fluiddynamics. I. Surface approximations and partial derivative estimates, Comput. Math. Appl., 19 (1990) 127–145.

[12] E.J. Kansa, Multiquadrics, a scattered data approximation scheme with applications to computational fuid dynamics. II. Solutions to parabolic, hyperbolic and elliptic partial differential equations, Comput. Math. Appl., 19 (1990) 147–161.

AN US

[13] B. Sˇarler, G. Kuhn, Primitive variable dual reciprocity boundary element method solution of incompressible Navier-Stokes euqaitons, Engineering Analysis with Boundary Elements 23 (1999) 443–455.

[14] I. Kovaˇcevi´c, A. Poredo sˇ, B. Sˇarler, Solving the Stefan problem with the radial basis function collocation method, Numerical Heat Transfer, Part B, 44 (2003) 575–599.

M

[15] B. Sˇarler, J. Perko, C.S. Chen, Radial basis function collocation method solution of natural convection in porous media, International Journal of Numerical Methods for Heat & Fluid Flow, 14 (2004), 187–212.

ED

[16] Y.V.S.S. Sanyasiraju, G. Chandhini, Local radial basis function based gridfree scheme for unsteady incompressible viscous flows, Journal of Computational Physics, 227 (2008) 8922–8948.

CE

PT

[17] R. Vertnik, B. Sˇarler, Meshless local radial basis function collocation method for convective-diffusive solid-liquid phase change problems, International Journal of Numerical Methods for Heat and Fluid Flow, 16 (2006) 617–640. [18] G. Kosec, B. Sˇarler, Local RBF collocation method for Darcy flow, Computer modeling in engineering & sciences, 3 (2008) 197–208.

AC

[19] R. Vertnik, B, Sˇarler, Solution of incompressible turbulent flow by a meshfree method, Computer Modeling in Engineering & Sciences, 44 (2009) 65– 96. [20] R. Vertnik, B. Sˇarler, Local collocation approach for solving turbulent combined forced and natural convection problems, Advances in Applied Mathematics and Mechanics, 3 (2011) 259–279. 17

ACCEPTED MANUSCRIPT

[21] K. Mramor, R. Vertnik, B. Sˇarler, Low and intermediate re solution of lid driven cavity problem by local radial basis function collocation method, Computers, materials & continua, 36 (2013), 1–21.

CR IP T

[22] K. Mramor, R. Vertnik, B. Sˇarler, Simulation of natural convection influenced by magnetic field with explicit local radial basis function collocation method. Computer modeling in engineering & sciences, 92 (2013), 327–352. [23] E.J. Kansa, Y.C. Hon, Circumventing the ill-conditioning problem with multiquadric radial basis functions: applications to elliptic partial differential equations, Computcomputer and Mathematics with Applications, 39 (2000) 123–137.

AN US

[24] N.A. Libre, A. Emdadi, E.J. Kansa, M. Rahimian, M. Shekarchi, A stabilized RBF collocation scheme for neumann type boundary value problems, Computer Modeling in Engineering and Science, 24 (2008) 61–80. [25] C.K. Lee, X. Liu, S.C. Fan, Local multiquadric approximation for solving boundary value problems, Computational Mechanics, 30 (2003) 396–409.

M

[26] B. Sˇarler, R. Vertnik, Meshfree explicit local radial basis function collocation method for diffusion problems, Computers and Mathematics with Applications, 21 (2006) 1269–1282.

ED

[27] R.K. Beatson, L. Greengard, A short course on fast multipole methods, in: M. Ainsworth, J. Levesley, W. Light, M. Marletta (Eds.), Wavelets, Multilevel Methods and Elliptic PDEs, Oxford University Press, 1997, pp.1–37.

CE

PT

[28] X. Zhang, H. Zhu, L.H. Kuo, A comparison study of the LMAPS method and the LDQ method for time-dependent problems, Engineering Analysis with Boundary Elements 37 (2013) 1408–1415.

AC

[29] X. Zhang, H. Tian, W. Chen, Local method of approximate particular solutions for two-dimensional unsteady Burgers’ equations, Computers and Mathematics with Applications 66 (2014) 2425–2432. [30] A.J. Chorin, Numerical solution of the Navier-Stokes equations. Math Comput, 22 (1968) 745–762. [31] U. Ghia, K. N. Ghia, C. T. Shin, High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method. J Comput Phys, 48 (1982) 387–411. 18