On a two-level element-free Galerkin method for incompressible fluid flow

Applied Numerical Mathematics 59 (2009) 1894–1904

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

On a two-level element-free Galerkin method for incompressible fluid flow Lin Zhang, Jie Ouyang ∗ , Xiao-Hua Zhang Department of Applied Mathematics, Northwestern Polytechnical University, Xi’an 710072, China

a r t i c l e

i n f o

Article history: Received 29 April 2008 Received in revised form 30 December 2008 Accepted 16 February 2009 Available online 20 February 2009 Keywords: Stabilization Incompressible flows EFG HVM Two-level method

a b s t r a c t In this paper, a new element-free Galerkin method called the two-level element-free Galerkin method is presented for incompressible fluid flow. This method consists of two parts: One at the global level and the other at the local level. The new method is based on the Hughes’ variational multiscale formulation, and arises from a decomposition of the unknown variables into coarse/resolved scale and fine/unresolved scale. In order to find the unresolved part of the solution, the momentum equation is used to define the unresolved governing equation on the local level. Finally the application of the proposed method to incompressible fluid flow is displayed. Numerical solutions obtained from the two-level element-free Galerkin method are present for three benchmark problems, and the results indicate that the proposed method is a robust and accurate method. This is due to the addition of the unresolved scale which adds stability to the discrete problem when equal order basis is used for pressure and velocity. © 2009 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction Accurate and efficient modeling of incompressible flows is an important issue in finite elements. Applications of the Galerkin finite element method (FEM) to incompressible flows were introduced in the early 1970s [29,32]. It was soon recognized that the use of equal-order interpolations for both velocity and pressure variables, which was the most desirable choice from the implementational point of view, generated spurious pressure modes. Solvability of the problem depends on a proper choice of finite element spaces for velocity and pressure. They must satisfy a compatibility condition, namely the so-called LBB (or inf-sup) condition. If this is not the case, alternative formulations (usually depending on a numerical parameter) are devised to circumvent the LBB condition and enable the use of velocity–pressure pairs that are unstable in the standard Galerkin formulation. In general, it is not trivial to verify analytically the LBB condition for a given interpolation of pressure and velocity, and this has spurred the use of numerical inf-sup testing [2,6,7,23]. In the mid-90s Hughes revisited the origins of the stabilization schemes from a variational multiscale view point and presented the variational multiscale method [16,17]. In this method the different stabilization techniques come together as special cases of the underlying subgrid scale modeling concept. Therefore, since this method termed as the Hughes variational multiscale (HVM) method was proposed, the stabilized methods have made great progress. Masud and co-workers developed multiscale/stabilized formulations for the linearized incompressible Navier–Stokes equations [24], the advection– diffusion equation [27], the convective-diffusive heat transfer [1], the Darcy flow equation [26], and the Fokker–Planck equation [25]. Franca et al. proposed a two-level finite element for the convection-diffusion problem [14] as well as the incompressible Navier–Stokes equations [13] and so on.

*

Corresponding author. Tel.: +86 29 88495234; fax: +86 29 88491000. E-mail address: [email protected] (J. Ouyang).

0168-9274/$30.00 © 2009 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2009.02.003

L. Zhang et al. / Applied Numerical Mathematics 59 (2009) 1894–1904

1895

Among the methods mentioned above, the spatial domain where the partial differential governing equations are defined is often discretized into meshes. Generally speaking, the creation of suitable meshes is very essential for acquiring accurate results. However, it is commonly very difficult and needs to consume a lot of time and labor to generate this kind of highquality meshes. Recently, a kind of so-called meshfree or meshless methods has developed quickly. Meshfree methods only use a set of nodes scattered within the problem domain as well as a set of nodes scattered on the boundary, thus can avoid the process of mesh generation. So far, there have existed many meshfree methods such as element-free Galerkin (EFG) method [3,4,20], meshless local Petrov Galerkin (MLPG) method [20] and reproducing kernel particle method (RKPM) [3] and so on. Among these methods, EFG which is proposed by Belytschko et al. [4] possesses a preferable foundation of mathematical theory. The more details of these meshfree methods can refer to [19–21]. It is well-known that the standard Galerkin finite element method is not ideally suited to deal with the incompressible flows when equal order basis is used for pressure and velocity. Unfortunately, the same thing also happens to meshfree methods. In order to simulate the incompressible flows by meshfree methods successfully, many researchers have done a lot of work and given three approaches. Firstly, the stream–vorticity formulation is adopted to describe the incompressible Navier–Stokes equations, e.g. Kim et al. used meshfree point collocation method to simulate the stream-vorticity of 2D incompressible Navier–Stokes and emphasized the novel formulation of effective vorticity condition on no-slip boundaries [18]. Secondly, the meshfree methods which can satisfy LBB condition have also been proposed. For instance, a new stabilized meshfree method for the approximation of the Stokes problem using weighted extended B-splines as shape functions has been developed, in which the web-spline based bilinear velocity-constant pressure element satisfies LBB condition [31]; Choe et al. proposed the moving least-square reproducing kernel (MLSRK) method and obtained the existence of discrete solution and its error estimate, finally they applied the method to solve the driven cavity flow numerically [8]; and in [11] a new EFG formulation is proposed using selective reduced integration. Thirdly, one of the very popular procedures of dealing with the pressure instability in the incompressible flow context is fractional step or projection method. Recent works of the authors have been on combining the convection stabilization developed via higher order time stepping with a fractional step method. The resulted method is widely known as the characteristic based split (CBS) scheme [28,35,36]. However, the above approaches have the shortcomings of their own. Specifically, for the first one, the boundary condition is not easy to designate; for the second one, we must use two set of nodes for velocity and pressure; for the third one, time step is essential to pressure stability and its accuracy is not high. In order to avoid these shortcomings and make use of the advantages of EFG, based on HVM a new kind of numerical method called the two-level element-free Galerkin (TLEFG) method is developed for the incompressible flows in the paper. So far, a lot of work has been done for the combination of FEM and HVM, but few papers have concerned for the combination of EFG and HVM. Thus the original contribution of the paper is to combine EFG with HVM to deal with the flow fluid problems. Consequently, the Stokes problem, the driven cavity flow problem and flow past a circular cylinder at low Reynolds number are solved numerically by the proposed method using the linear basis. The aim is to confirm that the proposed method can avoid the restriction of the LBB condition when equal order basis is used for pressure and velocity, especially the linear basis that is easy to implement. An outline of the paper is as follows: Section 2 presents the fundamental principle of the EFG method. Emphasis in the paper is the description of the TLEFG method, which is presented in Section 3. Section 4 presents the numerical results, and conclusions are drawn in Section 5. 2. The fundamental principle of the EFG method According to the moving least square (MLS) interpolant [3,4], a local approximation u h (x) to the function u (x) in the domain Ωx of influence of node x can be defined by u h (x) =

m 

p i (x)ai (x) = p T (x)a(x)

(1)

i =1

where p (x) is a complete polynomial basis of arbitrary order and a(x) is a coefficient needed to be determined, which as indicated, is the function of the space coordinate x. For the sake of simplicity, linear basis is chosen. Assume that we have known the nodal value u I = u (x I ) for the function u (x) at n nodes x I ( I = 1, 2, . . . , n) in the domain Ω . Then the unknown coefficient a(x) in (1) is obtained at any point x by minimizing the following weighted, discrete error norm J=

n 



w (x − x I ) p T (x I )a(x) − u I

2

(2)

I =1

where w (x − x I ) is a weight function of compact support (often called the domain of influence of node I ). The choice of the weight function is more or less arbitrary, and the spline function is often used in practice together with the exponential function. In the paper the following cubic spline function is chosen

⎧2 2 3 0  z  0.5, ⎨ 3 − 4z + 4z , w ( z) = 4 − 4z + 4z2 − 4 z3 , 0.5 < z  1, ⎩3 3 0 otherwise.

1896

L. Zhang et al. / Applied Numerical Mathematics 59 (2009) 1894–1904

Minimization of (2) with respect to a(x) then yields the following system of linear equations for the coefficient a(x): A (x)a(x) = B (x)u where A (x) and B (x) can easily be obtained from (2), and u is the vector of nodal unknowns. If A is invertible, the coefficient a(x) can be expressed as a(x) = A −1 (x) B (x)u . Substituting the above equation back into (1) leads to u h (x) = p T (x) A −1 (x) B (x)u = φ T (x)u where the shape function of the EFG method is given by

φ T (x) = p T (x) A −1 (x) B (x). The MLS approximation is obtained by a special least squares method, thus the function obtained by the MLS approximation is smooth curve and it does not pass through the nodal values. Therefore, the MLS shape function does not, in general, satisfy the Kronecker delta condition at each node, i.e.

φi (x j ) = δi j which brings great difficulties in exerting Dirichlet boundary condition. Today there have existed many techniques to resolve this problem, such as Lagrange multiplier approaches [3], penalty methods [3], modified variational principles [22], perturbed Lagrangian [9], coupling to finite elements [5] and so on. Among them Lagrange multiplier approaches is adopted in the paper. 3. The two-level element-free Galerkin method Based on the HVM method, we decompose the unknown variables into coarse/resolved scale and fine/unresolved scale, e.g. u = u¯ + u  where u¯ and u  are called the ‘resolved scale’ and ‘unresolved scale’, respectively. Substituting such decomposition into the classical Galerkin weak formulation of the problem, a new Galerkin form containing the unresolved scale is then obtained. Here we name this new Galerkin form as the global problem or the problem on the global level, because it is defined on the whole domain of the problem. However, the unresolved part of the solution in this new form is still unknown. In order to compute the unresolved part, before solving the global problem we define a local problem in each background integral cell of the standard EFG method used in the global problem. Here the solvability of the unresolved scale will also be accomplished by using the standard EFG method. Finally an appropriate solver is invoked to solve the global problem after assembling the unresolved part of the solution into the discretization on the global level. 3.1. Computational algorithm on the global level The Stokes problem is the approximation of the incompressible Navier–Stokes equations in the case of low Reynolds numbers. For the sake of simplicity, here we consider the Stokes problem with a pure Dirichlet boundary. The strong form of the problem reads

⎧ 2 ⎨ −η∇ u + ∇ p = f ∇·u=0 ⎩ u=0

in Ω,

(3) (4) (5)

in Ω, on Γ

where u = (u , v ) is the velocity, p is the kinematic pressure, f is the body force per unit mass and η is the kinematic viscosity. We use the momentum equation (3) to define the governing equation of the unresolved scale [13]. Thus for a given background integral cell Ω cell of the standard EFG method, the unresolved scale u  is zero outside Ω cell and satisfies the following system



−η∇ 2 (u¯ + u  ) + ∇ p¯ = f u = 0

in Ω cell , on ∂Ω

cell

.

(6) (7)

It should be noted that the system (6)–(7) is used to find the unresolved scale of the solution at the pre-processing stage which will be extensively described in the next section. For the system (3)–(5), we can infer that the Galerkin method containing the unresolved scale reads [13]:



¯ , p¯ ) + q¯ , div(u¯ + u  ) = ( w ¯ , f) ¯ , η∇(u¯ + u  ) − (div w ∇w

(8)

where u  is given in terms of u¯ through the relation (6)–(7). From the next section, we know that the approximations for velocity and pressure can be written explicitly in terms of the basis functions as follows:

L. Zhang et al. / Applied Numerical Mathematics 59 (2009) 1894–1904



u¯ =



uai ψai ei ,

(9)

a=1,...,nen i =1,...,nsd



u =







uai ϕai ei +

a=1,...,nen i =1,...,nsd

p¯ =

1897

pa˜

a˜ =1,...,nenp

 i =1,...,nsd

ϕ˜ a˜i ei +



ϕ if ei ,

(10)

i =1,...,nsd

pa˜ ψa˜

(11)

a˜ =1,...,nenp

where uai and pa˜ represent the nodal approximation values, nen is the number of velocity nodes, nenp is the number of pressure nodes (in this paper nenp = nen), nsd is the dimension of the problem (in this paper nsd = 2), ψ ’s are the shape functions of the standard EFG method on the global level, ϕ ’s are the basis functions of the unresolved scale and we will define these ϕ ’s precisely later. Now we insert u¯ , u  and p¯ into (8) so that the discrete variational problem (8) can be written in the following matrix form. Denoting the number of equations by neq, then we get A neq×neq × U neq×1 = B neq×1 where

 A neq×neq =

(12)

j

j

(∇ψb e j , η∇ψai ei ) + (∇ ϕb e j , η∇ψai ei ) j



j

j  j, j =1,...,nsd (∇ ˜ b˜ e

(∇ · ψb e j , ψa˜ ) + (∇ · ϕb e j , ψa˜ ) 

j ( f , ψai ei ) − j =1,...,nsd (∇ ϕ f e j , η∇ψai ei ) B neq×1 = , j − j =1,...,nsd (∇ · ϕ f e j , ψa˜ )



ϕ

η∇ψai ei ) − (ψb˜ , ∇ · ψai ei )

j =1,...,nsd (∇



U neq×1 =

uj pj



j

In the formulation (12), there exist many terms, namely (∇ ϕb e j , η∇ψai ei ), (∇ · ϕb e j , ψa˜ ), j =1,...,nsd (∇

j

· ϕ˜ ˜ e j , ψa˜ ), b



j =1,...,nsd (∇

ϕ

j e , f j

η

∇ψai ei )

and



j =1,...,nsd (∇

·ϕ

,

b

.

j



j

· ϕ˜ ˜ e j , ψa˜ )

j e , ψa˜ ), f j



j  j, j =1,...,nsd (∇ ˜ b˜ e

ϕ

η∇ψai ei ),

which are absent in the standard EFG

method. These terms are actually contributed to the formulation by the unresolved scale and responsible for the stability of the scheme. 3.2. Computational algorithm on the local level To compute the contribution of the unresolved scale to the problem on the global level, we need to calculate the unresolved part in each background integral cell of the standard EFG method. This can be done by solving (6), (7). The system (6), (7) can be written as



−η∇ 2 u  = f + η∇ 2 u¯ − ∇ p¯ , in Ω cell , u = 0 on ∂Ω cell

which can be viewed as a scalar Poisson equation in each component. To find the solution of the problem (13)–(14), we specify the standard EFG method in the integral cell to approximate the unresolved scale, and then use these approximations in the exact solutions of the unresolved scale. Additionally, let L = −η∇ 2 . Since the operator L is decomposing (13)–(14) into its basis components [13]. Thus, to find u  , we solve the basis functions ϕai ’s, ϕ˜ a˜i ’s and ϕ if ’s for the unresolved scale,

(13) (14) interior of each background problem (8) in place of the linear, this can be done by with a varying from one to

the number of the velocity nodes in the background integral cell (nen), a˜ varying from one to the number of the pressure nodes in the background integral cell (nenp) and i varying from one to nsd. Specifically, these basis functions are obtained via the following governing equations, respectively. For ϕai ’s:



L ϕai ei = − L ψai ei

in Ω cell , on ∂Ω cell .

(15) (16)

L ϕ˜ a˜i ei = − ∂ x a˜ ei i

in Ω cell ,

(17)

˜ a˜i ei

on ∂Ω cell

(18)

i i ae

ϕ

=0

For ϕ˜ a˜i ’s:



∂ψ

ϕ

=0

where ψai ’s and ψa˜ ’s are the shape functions of the standard EFG method for u¯ and p¯ on the global level, respectively. For ϕ if ’s:



L ϕ if ei = f i ei

ϕ ei = 0 i f

in Ω cell , on ∂Ω

cell

(19) .

(20)

1898

L. Zhang et al. / Applied Numerical Mathematics 59 (2009) 1894–1904

Thus, if



u¯ =



uai ψai ei ,

a=1,...,nen i =1,...,nsd



p¯ =

pa˜ ψa˜

a˜ =1,...,nenp

then



u =





uai ϕai ei +

a=1,...,nen i =1,...,nsd

pa˜

a˜ =1,...,nenp

 i =1,...,nsd

ϕ˜ a˜i ei +



ϕ if ei .

i =1,...,nsd

So we have demonstrated why the unresolved scale u  is approximated by (10). Substituting these approximations into (8), we get the matrix formulation (12). Now we consider how to solve (15)–(16), (17)–(18) and (19)–(20) by the standard EFG method. Let us denote by Ω ∗,cell an arbitrary background integral cell in the global EFG method and by ψl∗ the basis function for each local EFG’s node in the background integral cell Ω ∗,cell , so our unknown basis functions

⎧ ϕ i = c (a,i) ψl∗ , ⎪ ⎪ ⎨ a l l ϕ˜ a˜i = l cl(˜a,i) ψl∗ , ⎪ ⎪ ⎩ i ( f ,i ) ∗ ϕ f = l cl ψl .

ϕai ’s, ϕ˜ a˜i ’s and ϕ if ’s can be approximated by

(21)

refer to the specific Here l runs over all unknown interior EFG’s nodes in Ω ∗,cell , say N ∗ . The index (a, i ), (˜a, i ) and cl shape function that we are trying to compute. Substituting (21) in (15), (16) and (17), (18), respectively, the matrix forms of these problems read: For each a (from 1 to ( f ,i )

(a,i )

nen), for each a˜ (from 1 to nenp) and for i = 1, . . . , nsd, find cl



(a,i )



(˜a,i )

(˜a,i )

∗

∗

η∇ψl , ∇ψm

cl



l



 ∂ψa˜ ∗ , ψm . = − ∂ xi 

Similarly, we can also account for the basis functions associated with f and



, l = 1, 2, . . . , N ∗ , such that

η∇ψl∗ , ∇ψm∗ = η ψai , ψm∗ ,

cl

l





, l = 1, 2, . . . , N ∗ and cl

( f ,i )





η∇ψl∗ , ∇ψm∗ = f i , ψm∗

cl

ϕ f : Find cl( f ,i) , l = 1, 2, . . . , N ∗ , such that



l

for i ∈ {i , . . . , nsd}.

(a,i )

(˜a,i )

( f ,i )

Here the standard EFG method is used to obtain the coefficients cl , cl and cl in (21). Additionally, non-standard EFG methods [30] can also be used in this pre-processing stage instead of the standard EFG method. (a,i )

i a,

Once the constants cl

˜ a˜i

i f

(˜a,i )

, cl

( f ,i )

and cl

are found, we substitute them in (21) to get the approximate basis functions

ϕ ϕ and ϕ for the unresolved scale. These approximations are then plugged into the global problem with the help of the representing identities in (9)–(11) in place of the exact basis functions. Finally the global problem (8) can be solved numerically by invoking an appropriate solver. 4. Numerical examples 4.1. Analytical test In order to illustrate that the TLEFG method can violate the celebrated LBB condition when equal order basis is used for velocity and pressure, especially the linear basis that is easy to implement, we consider a two-dimensional Stokes problem [15] in the square domain Ω = [0, 1] × [0, 1], which possesses a closed-form analytical solution. The problem consists of determining the velocity field u = (u , v ) and the pressure p such that (3)–(5), where the fluid viscosity is took as η = 1. The components of the body force f are prescribed as





f x = (12 − 24 y )x4 + (−24 + 48 y )x3 + −48 y + 72 y 2 − 48 y 3 + 12 x2



+ −2 + 24 y − 72 y 2 + 48 y 3 x + 1 − 4 y + 12 y 2 − 8 y 3 ,



f y = 8 − 48 y + 48 y 2 x3 + −12 + 72 y − 72 y 2 x2 + 4 − 24 y + 48 y 2 − 48 y 3 + 24 y 4 x − 12 y 2 + 24 y 3 − 12 y 4 .

L. Zhang et al. / Applied Numerical Mathematics 59 (2009) 1894–1904

Fig. 1. The solutions of the Stokes problem: (a) analytical solutions, (b) EFG’s solutions, (c) TLEFG’s solutions, respectively.

1899

1900

L. Zhang et al. / Applied Numerical Mathematics 59 (2009) 1894–1904

Fig. 2. Rate of convergence for (a) u velocity, (b) v velocity and (c) pressure, respectively.

Fig. 3. The node distribution of the driven cavity flow problem: (a) regular node distribution for EFG, (b) regular node distribution for TLEFG, (c) irregular node distribution for TLEFG.

With this prescribed body force, the exact solution is









u (x, y ) = x2 (1 − x)2 2 y − 6 y 2 + 4 y 3 v (x, y ) = − y 2 (1 − y )2 2x − 6x2 + 4x3 p (x, y ) = x(1 − x). Additionally, in order to analyze the convergence rate of the proposed method, the following L 2 error norm is used as error measures:

  u − u h 



L2

=



u − uh

2

 1/ 2

dΩ

,

Ω

and the corresponding relative error norm is u − u h rel = u − uh  L 2 /u  L 2 . L2 We solve this problem with the standard EFG method and the TLEFG method, respectively. These two methods both employ linear basis p T = (1, x, y ) for velocity and pressure. In these two methods, 21 × 21 EFG’s nodes and 20 × 20 background integral cells are distributed uniformly in the domain Ω = [0, 1] × [0, 1]. Additionally, for the TLEFG method 3 × 3 EFG’s nodes and 2 × 2 background integral cells are distributed uniformly in the background integral cell of the standard EFG method on the global level. Fig. 1 represents the analytical solutions (Fig. 1(a)), the numerical solutions via EFG (Fig. 1(b)) and the numerical solutions via TLEFG (Fig. 1(c)), respectively. They all contain four figures, namely (a) u contours, (b) v contours, (c) pressure contours and (d) 3-D pressure. When p T = (1, x, y ) is used for velocity and pressure, from Fig. 1(b) we find that the accuracy of the velocity u = (u , v ) via EFG is very high, however, for the pressure there exist a lot of sharp prickles which imply that this result is not acceptable. Under the same condition, we can obtain simultaneously the satisfactory results for velocity and pressure via TLEFG (see Fig. 1(c)) instead of EFG. In other words, when equal order basis especially the linear basis is used for velocity and pressure, the proposed method can avoid the restriction of the LBB condition for this problem. In the present study, the error norm mentioned above for the velocity and the pressure is plotted with respect to the number of nodes. In order to study the convergence of the method, four regular node distributions with 121, 256, 441 and 676 nodes are used for the MLS approximation of the quantities. Fig. 2 presents the convergence rates for all variables in L 2 and it is observed that the results are consistent with the error estimate theory. 4.2. Driven cavity flow problem This example has become a standard benchmark test for incompressible flows. It models a plane flow of an isothermal fluid in a square lid-driven cavity. The upper side of the cavity moves in its own plane at unit speed, while the other sides

L. Zhang et al. / Applied Numerical Mathematics 59 (2009) 1894–1904

1901

Fig. 4. The solutions of the driven cavity flow problem via: (a) EFG with regular node distribution, (b) TLEFG with regular node distribution, (c) TLEFG with irregular node distribution.

1902

L. Zhang et al. / Applied Numerical Mathematics 59 (2009) 1894–1904

Fig. 5. Problem description for stationary flow around cylinder: (a) problem description, (b) mesh.

Fig. 6. The solutions of flow past a circular cylinder via the TLEFG method (Re = 40).

are fixed. The domain is Ω = [0, 1] × [0, 1], the body force is neglected, namely f = (0, 0) T , and the boundary conditions are formally defined as: g = (0, 0) T on ∂Ω\{ y = 1} and g = (1, 0) T on ∂Ω ∩ { y = 1}. Similarly, this well-known benchmark problem is solved by the standard EFG method and the TLEFG method, respectively. Fig. 3 presents the node distribution adopted in the above mentioned methods. For the TLEFG method, we adopt two kinds of node distribution, namely regular node distribution and irregular node distribution; while for the standard EFG method, we just only use regular node distribution. When regular node distribution is used, the EFG method and the TLEFG both adopt the same parameters as those in the previous problem and employ linear basis p T = (1, x, y ) to approximate velocity and pressure. Fig. 4 depicts the corresponding numerical results, namely (a) the results obtained by the standard EFG method with regular node distribution, (b) the results obtained by the TLEFG method with regular node distribution, (c) the results obtained by the TLEFG method with irregular node distribution. From Fig. 4(a) it is observed clearly that the standard EFG method does not only present spurious pressure modes which vibrate severely and are not acceptable at all, but also presents oscillations on the velocity fields. Nevertheless, we can obtain simultaneously the satisfactory results for velocity and pressure via the TLEFG method (see Fig. 4(b) and Fig. 4(c)) instead of the standard EFG method, in particular no spurious pressure modes are observed. Meanwhile, comparing Fig. 4(b) with Fig. 4(c), we find that the results obtained with regular node distribution are smoother than those obtained with irregular node distribution. Finally from the comparison between the numerical results of these two methods, we can conclude that for this problem the proposed method can violate the LBB condition and has a good stability in spite of equal order basis for velocity and pressure. 4.3. Flow past a circular cylinder at low Reynolds number We finally consider the two-dimensional flow of an incompressible fluid past a circular cylinder. The Reynolds number considered here is 40, for which a steady-state solution exists. The problem description is given in Fig. 5(a). We divide the domain into 40 × 20 four node elements, where along the circular cylinder we divide the boundary into 40 equal parts and along the radius direction we set 21 points in terms of geometric series. A zoomed view of the mesh that consist of four node elements wherein the nodes are concentrated around the cylinder, is presented in Fig. 5(b). Fig. 6 shows

L. Zhang et al. / Applied Numerical Mathematics 59 (2009) 1894–1904

1903

the contours of u velocity, v velocity and pressure, respectively. It is easy to see that we can obtain simultaneously the satisfactory results for velocity and pressure via TLEFG. This also confirms that the TLEFG method can indeed avoid LBB condition when velocity and pressure both adopt linear approximation. Finally, the total drag coefficient (dimensionless) 2 C D = F D /(1/2ρ U ∞ d) is estimated numerically where F D is the total drag force exerted by the fluid on the cylinder. The value of C D is equal to 1.50891 by the proposed method and agrees considerably well with 1.4980 ± 0.0005 in [12], 1.52 in [10,34], 1.53 in [33] and so on. 5. Conclusions In this paper, we have presented the two-level element-free Galerkin method for incompressible fluid flow. Three wellknown benchmark examples are solved by the standard EFG method and the TLEFG method, respectively. The numerical results indicate that the standard EFG method fails the numerical inf-sup test for equal order basis for pressure and velocity. However, the proposed method TLEFG allows the equal order basis, especially the linear basis for pressure and velocity, and stability is observed as refinement progresses. So the TLEFG method is a robust and accurate method due to the addition of the unresolved scale which adds stability to the discrete problem. Finally this method differs from FEM-like approaches by avoiding the need of meshing, a very demanding task for complicated geometry problems. Acknowledgements The support from the National Natural Science Foundation of China (NSFC) (No: 10590353, 10871159) and the National Basic Research Program of China (No: 2005CB321704). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31]

M. Ayub, A. Masud, A new stabilized formulation for convective-diffusive heat transfer, Numer. Heat Transfer 43 (2003) 601–625. K.J. Bathe, D. Hendriana, F. Brezzi, G. Sangalli, Inf-sup testing of upwind methods, Int. J. Numer. Methods Engrg. 48 (2000) 745–760. T. Belytschko, Y. Krongauz, D. Organ, Meshless method: An overview and recent developments, Comput. Methods Appl. Mech. Engrg. 139 (1996) 3–47. T. Belytschko, Y.Y. Lu, L. Gu, Element-free Galerkin methods, Int. J. Numer. Methods Engrg. 37 (1994) 229–256. T. Belytschko, D. Organ, Y. Krongauz, A coupled finite element–element-free Galerkin method, Comput. Mech. 17 (1995) 186–195. F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer Series in Computational Mathematics, vol. 15, Springer-Verlag, New York, 1991. D. Chapelle, K.J. Bathe, The inf-sup test, Comput. Struct. 47 (1993) 537–545. H.J. Choe, D.W. Kim, H.H. Kim, Y. Kim, Meshless method for the stationary incompressible Navier–Stokes equations, Discrete Contin. Dyn. Syst. Ser. B 1 (4) (2001) 495–526. Y.A. Chu, B. Moran, A computational model for nucleation of the solid–solid phase transformations, Modelling Simul. Mater. Sci. Engrg. 3 (1995) 455–471. S.C.R. Dennis, G.Z. Chang, Numerical solution for steady flow past a circular cylinder at Reynolds number up to 100, J. Fluid Mech. 42 (1970) 471–489. J. Dolbow, T. Belytschko, Volumetric locking in the element free Galerkin method, Int. J. Numer. Methods Engrg. 46 (6) (1999) 925–942. B. Fornberg, A numerical study of steady viscous flow past a circular cylinder, J. Fluid Mech. 98 (4) (1980) 819–855. L.P. Franca, A. Nesliturk, On a two-level finite element method for the incompressible Navier–Stokes equations, Int. J. Numer. Meth. Engrg. 52 (2001) 433–453. L.P. Franca, A. Nesliturk, M. Stynes, On the stability of residual-free bubbles for convection–diffusion problems and their approximation by a two-level finite element method, Comput. Methods Appl. Mech. Engrg. 166 (1998) 35–49. A. Huerta, Y. Vidal, P. Villon, Pseudo-divergence-free element free Galerkin method for incompressible fluid flow, Comput. Methods Appl. Mech. Engrg. 193 (2004) 1119–1136. T.J.R. Hughes, Multiscale phenomena: Green’s function, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods, Comput. Methods Appl. Mech. Engrg. 127 (1995) 387–401. T.J.R. Hughes, G.R. Feijoo, M. Luca, Q. Jean-Baptiste, The variational multiscale method- a paradigm for computational mechanics, Comput. Mehtods Appl. Mech. Engrg. 166 (1998) 3–24. Y. Kim, D.W. Kim, S. Jun, J.H. Lee, Meshfree point collocation method for the stream–vorticity formulation of 2D incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg. 196 (2007) 3095–3109. G.R. Liu, Mesh Free Methods: Moving Beyond the Finite Element Method, CRC Press, New York, 2002. G.R. Liu, Y.T. Gu, An Introduction to Meshfree Methods and Their Programming, Springer, Berlin and New York, 2005. G.R. Liu, M.B. Liu, Smoothed Particle Hydrodynamics: A Meshfree Particle Method, World Scientific, Singapore, 2003. Y.Y. Lu, T. Belytschko, L. Gu, A new implementation of the element-free Galerkin method, Comput. Methods Appl. Mech. Engrg. 113 (1994) 397–414. D.S. Malkus, Eigenproblems associated with the discrete LBB condition for incompressible finite elements, Int. J. Engrg. Sci. 19 (1981) 1299–1310. A. Masud, On a stabilized finite element formulation for incompressible Navier–Stokes equations, in: Proceedings of the Fourth US–Japan Conference on Computational Fluid Dynamics, Tokyo, Japan, May 2002, pp. 28–30. A. Masud, L.A. Bergman, Application of multiscale finite element methods to the solution of the Fokker–Planck equation, Comput. Methods Appl. Mech. Engrg. 194 (2005) 1513–1526. A. Masud, T.J.R. Hughes, A stabilized mixed finite element method for Darcy flow, Comput. Methods Appl. Mech. Engrg. 1991 (2002) 4341–4370. A. Masud, R.A. Khurram, A multiscale/stabilized finite element method for the advection-diffusion equation, Comput. Methods Appl. Mech. Engrg. 193 (2004) 1997–2018. P. Nithiarasu, R. Codina, O.C. Zienkiewicz, Characteristic based split (CBS) scheme: a unified approach to fluid dynamics, Int. J. Numer. Methods Engrg. (2006) (special issue). J.T. Oden, L.C. Welloford, Analysis of flow of viscous fluids by the finite element method, AIAA J. 10 (1972) 590–1599. J. Ouyang, L. Zhang, X.H. Zhang, Nonstandard element free Galerkin method for solving unsteady convection dominated problems, Acta Aerodynam. Sinica 25 (2007) 287–293. V.V.K. Srinivas Kumar, B.V. Rathish Kumar, P.C. Das, Weighted extended B-spline method for the approximation of the stationary Stokes problem, J. Comput. Appl. Math. 186 (2006) 335–348.

1904

L. Zhang et al. / Applied Numerical Mathematics 59 (2009) 1894–1904

[32] C. Taylor, P. Hood, A numerical solution of the Navier–Stokes equations using the finite element technique, Comput. Fluids 1 (1973) 73–100. [33] Y.H. Tseng, J.H. Ferziger, A ghost-cell immersed boundary method for flow in complex geometry, J. Comput. Phys. 192 (2003) 593–623. [34] T. Ye, R. Mittal, H.S. Udaykumar, W. Shyy, An accurate Cartesian grid method for viscous incompressible flows with complex immersed boundary, J. Comput. Phys. 156 (1999) 209–240. [35] O.C. Zienkiewicz, R. Codina, A general algorithm for compressible and incompressible flow—Part I: The split, characteristic-based scheme, Int. J. Numer. Methods Fluids 20 (1995) 869–885. [36] O.C. Zienkiewicz, P. Nithiarasu, R. Codina, M. Vazquez, P. Ortiz, The characteristic based split procedure: an efficient and accurate algorithm for fluid problems, Int. J. Numer. Methods Fluids 31 (1999) 359–392.