A finite element variational multiscale method based on two-grid discretization for the steady incompressible Navier–Stokes equations

Accepted Manuscript A finite element variational multiscale method based on two-grid discretization for the steady incompressible Navier-Stokes equations Yueqiang Shang, Jin Qin PII: DOI: Reference:

S0045-7825(15)00366-7 http://dx.doi.org/10.1016/j.cma.2015.11.013 CMA 10757

To appear in:

Comput. Methods Appl. Mech. Engrg.

Received date: 10 September 2015 Revised date: 3 November 2015 Accepted date: 10 November 2015 Please cite this article as: Y. Shang, J. Qin, A finite element variational multiscale method based on two-grid discretization for the steady incompressible Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. (2015), http://dx.doi.org/10.1016/j.cma.2015.11.013 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.

*Manuscript Click here to download Manuscript: 2level-mutiscale2015-10).pdf

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Click here to view linked References

A finite element variational multiscale method based on two-grid discretization for the steady incompressible Navier-Stokes equations I Yueqiang Shanga,∗, Jin Qinb a School

b School

of Mathematics and Statistics, Southwest University, Chongqing 400715, P. R. China of Mathematics and Computer Science, Zunyi Normal college, Zunyi 563002, P. R. China

Abstract By combining the best algorithmic features of two-grid discretization method and a recent variational multiscale method, a two-level finite element variational multiscale method based on two local Gauss integrations for convection dominated incompressible Navier-Stokes equations is proposed and analyzed. In this method, a fully nonlinear NavierStokes problem is first solved on a coarse grid, and then a linear problem is solved on a fine grid to correct the coarse grid solution, where the numerical forms of the Navier-Stokes equations both on coarse and fine grids are stabilized by a stabilization term defined by the difference between a consistent and an under-integrated matrix of the velocity gradient. Error bounds of the approximate solution are analyzed. Algorithmic parameter scalings of the method are derived. Numerical tests are also given to verify the theoretical predictions and demonstrate the efficiency and promise of the method. Keywords: Navier-Stokes equations, finite element, variational multiscale method, two-grid method

1. Introduction For an incompressible Newtonian viscous fluid, the governing equations of its flow are the Navier-Stokes equations. It is generally accepted that incompressible fluid flows, even for high Reynolds numbers, are faithfully modeled by the Navier-Stokes system. Therefore, a computational fluid dynamics model is in general based on the solution of the Navier-Stokes equations and its discretization scheme such as finite element methods, finite volume methods and finite difference methods. As a result, many attentions were put on the development of efficient numerical methods for the Navier-Stokes equations in the last decades; see, for example, the monographs of Temam [1], Girault and Raviart [2, 3], and Elman, Silvester and Wathen[4], among others. Among the successful discretization methods, two-level or multilevel finite element methods get popular since the pioneered work of Xu [5, 6]. For examples, Layton et al. [7, 8, 9, 10, 11], Girault and Lions [12], He et al. [13, 14, 15, 16], Liu and Hou [17, 18] have studied two-level or multilevel methods for the stationary Navier-Stokes I This work was supported by the Natural Science Foundation of China (No. 11361016), the Project Sponsored by the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry, Fundamental Research Funds for the Central Universities (No. XDJK2014C160, SWU113095), and the Science and Technology Foundation of Guizhou Province, China (No. [2013]2212). ∗ Corresponding author. Tel. +86-23-68252350. E-mail: [email protected]

Preprint submitted to Elsevier

November 3, 2015

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

equations, among others. While for the nonstationary Navier-Stokes equations, two-level or multilevel methods have been studied by Girault and Lions [19], Olshanskii [20], He et al. [21, 22, 23, 24, 25, 26], Hou and Mei [27], and Abboud et al. [28]. Based on domain decomposition methods, some parallel two-level discretization methods were also proposed for the Navier-Stokes equations by He et al. [29, 30], Ma et al. [31, 32] and Shang et al. [33, 34, 35, 36]. The key idea of two-level or multilevel methods for the Navier-Stokes equations is to first solve the fully nonlinear Navier-Stokes equations on a coarse grid, and then solve a linear problem on finer grid(s) to correct the solution. Theoretical analysis and numerical tests showed that, for an appropriate choice of the coarse and fine grids sizes, these methods can yield a convergence rate of the same order as the classical one-level finite element method with a substantial reduction in computational time and, consequently, are significantly more efficient than the standard one-level finite element method. However, there is a shortcoming for the two-level or multilevel discretization methods that a fully nonlinear Navier-Stokes problem needs to be solved on the coarsest grid, which limits their applicability in simulation of high Reynolds number flows. It is well known that high Reynolds number flows are characterized by a wide spectrum of sizes of the flow structures (scales) ranging from large ones to very small ones. If the grids used for the discretization of the Navier-Stokes equations are not fine enough to represent all scales of the flow being simulated, spurious oscillations may occur. What is worse, the iterative methods used to solve the nonlinear system may fail to converge, and hence cannot yield an approximate solution at all (cf. [15, 37, 38]). Therefore, to apply two-level or multilevel discretization methods to the simulation of high Reynolds number flows, stabilization techniques are essential. Recently, a finite element variational multiscale (VMS) method based on two local Gauss integrations for the simulation of Navier-Stokes equations was proposed by Zheng et al. in [39]. In this method, the stabilization term is based on two local Gauss integrations at element level. It is defined by the difference between a consistent and an under-integrated matrix only involving the gradient of velocity interpolation, and hence avoids extra storage compared to the usual projection-based variational multiscale methods. The idea of two local Gauss integrations was first proposed by Li and He [40, 41] to offset the discrete pressure space in order to circumvent the inf-sup condition using P1 − P1 finite elements pair to dicretize the Stokes and Navier-Stokes equations. Subsequently, Zheng et al. used

this two local Gauss integrations idea to improve the finite element VMS method [39] and to develop an adaptive technique [42, 43]. This finite element VMS method based on two local Gauss integrations was then applied to the time-dependent Naver-Stokes equations [44] and the stationary conduction-convection problems [45]. It was also combinized with the two-grid discretization method [46], the domain decompositiom method [47] and the partition of unity [48] (for the partition of unity method, we refer the readers to, for example, [49, 50, 51], among others) . In this paper, we consider a synthesis method of the above mentioned finite element VMS method based on two local Gauss integrations with two-level discretization method. Specially, we first solve a nonlinear Navier-Stokes problem on a coarse grid, and then solve a linear problem on a fine grid where the nonlinear convective term is fixed by the coarse grid solutiom. Both on the coarse and fine grids, the numerical forms of the Navier-Stokes equations are stabilized by a stabilization term based on two local Gauss integrations. The significant feature of the 2

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

synthesis method is that it is able to simulate high Reynolds number flows without spurious oscillations compared to the standard two-grid discretization methods; whileas compared to the one-level VMS method, it can save a large amount of computational time. It is also easy to implement based on existing codes. It is noted that some two-level VMS methods have been proposed in, for example, [52, 53, 54, 55] among others. However, there are essential differences between our two-level method and those studied in [52, 53, 54, 55]. Firstly, the purpose and function of the two grids used are different. In our method, the two grids are used to increase the efficiency of our numerical scheme, while in [52, 53, 54, 55], the two grids were used for numerical stabilization. Secondly, our VMS method does not need the construction of any projections, while those in [52, 53, 54, 55] are based on an appropriate projection which defines the large scales. Thirdly, in our method, the coarse grid problem and the fine grid correction are uncoupled in the computational process, which is different from the methods in [52, 53, 54, 55] where the coarse and fine grids variables are coupled by the projection. There are also other projection-based VMS methods such as those in [56, 57, 58, 59]. It is also worthy to mention that although our idea of combination of two-grid discretization method with the finite element VMS method based on two Gauss integrations is similar to that of [46], there is significant novelty for two reasons. First, our method is different from that of [46]. In [46], the linear problem on the fine grid is a Stokes problem, while here is one step of Newton iterations, and consequently, is able to yield a better approximate solution with a higher accuracy. Moreover, in our method, the stabilization term in the fine grid linear problem is explicitly treated by inserting the coarse grid solution, which differs from that of [46]. Second, the numerical analyses are completely different. In [46], all the theoretical results are limited to the so-called uniqueness condition

N k f k−1,Ω ν2

< 1 (see (2.5)

of [46], where ν is the kinematic viscosity of the fluid, f the prescribed body force, Ω the solution domain, and N a postive constant related to the convection term of the Navier-Stokes equations), which is usually invalid for high Reynolds number flows that the VMS methods are supposed to simulate; what is more, their theoretical asymptotical error estimates for the discrete approximate solution are not optimal. While our theoretical analysis does not need the above so-called uniqueness condition, and optimal asymptotical error estimates are obtained. The paper is organized as follows. In the next section, we give mathematical preliminaries for the the NavierStokes equations. Section 3 gives the finite element VMS method based on two local Gauss and two-level finite element VMS method as well as corresponding numerical analysis. Numerical results presented in Section 4, support the analytical results and demonstrate the efficiency of the proposed method. Finally, Section 5 summarizes the results and concludes the article.

3

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

2. Mathematical preliminaries Let Ω be a bounded domain with Lipschitz-continuous boundary ∂Ω in Rd (d=2 or 3). We consider the following incompressible Navier-Stokes problem: −ν∆u + (u · ∇)u + ∇p = f,

in Ω,

(2.1)

∇ · u = 0,

in Ω,

(2.2)

on ∂Ω.

(2.3)

u = 0,

Here, u : Ω → Rd is the velocity, p : Ω → R the pressure, f : Ω → Rd the prescribed body force, and ν the kinematic viscosity. Given a characteristic length L and a characteristic velocity U, the Reynolds number is defined as Re = UL/ν. For the mathematical setting of problem (2.1)-(2.3), we introduce the following Hilbert spaces: Z X = H01 (Ω)d , Y = L2 (Ω)d , M = L02 (Ω) = {q ∈ L2 (Ω) : qdx = 0}, Ω

and denote by k · kk the usual norm of Sobolev space H k (Ω)l for a nonnegative integer k, and by (·, ·) the standard

inner-product of L2 (Ω)l (l = 1, 2, 3). The space H −1 (Ω)l , the dual of H01 (Ω)l , and the corresponding norm k · k−1 will

also be used. Throughout this paper, we shall use the letter c or C to denote a generic positive constant (with or without subscripts) which is independent of mesh parameter and may take on different values on different occurrences. We define trilinear term b(·, ·, ·) as 1 b(u, v, w) = ((u · ∇)v, w) + ((∇ · u)v, w) 2 1 1 = ((u · ∇)v, w) − ((u · ∇)w, v), 2 2

∀ u, v, w ∈ X,

which has the following properties (cf. [1, 3, 60]): b(u, v, w) = −b(u, w, v),

(2.4)

∀ u, v, w ∈ X,

|b(u, v, w)| ≤ Nk∇uk0 k∇vk0 k∇wk0 ,

∀ u, v, w ∈ X,

(2.5)

where N is a positive constant depending only on Ω. The variational formulation of (2.1)-(2.3) reads: find a pair (u, p) ∈ X × M such that ν(∇u, ∇v) + b(u, u, v) − (∇ · v, p) + (∇ · u, q) = ( f, v),

∀(v, q) ∈ X × M.

(2.6)

To introduce the finite element discretization of (2.6), we assume T µ (Ω) = {K} (here µ = H, h with H > h) to

be a shape-regular triangulation (see, e.g., [3, 4]) of Ω into triangles or quadrilaterals (if d = 2), or tetrahedrons or hexahedrons (if d = 3) with mesh size 0 < µ < 1. The fine mesh T h (Ω) can be thought of as generated from the coarse mesh T H (Ω) by a mesh refinement process (see, e.g., [61]), and, therefore, nested. It is not necessary for our 4

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

method, nor needed for our convergence theory to hold. However, we assume them nested since it will simplify our analysis substantially. Throughout this paper, we confine our attention to the context of Taylor-Hood elements pair spaces [62], which, associated with a mesh T µ (Ω), are defined as Xµ = {vµ ∈ H01 (Ω)d : vµ |K ∈ (P2 )d ,

∀K ∈ T µ (Ω)},

Mµ = {qµ ∈ L02 (Ω) ∩ C 0 (Ω) : qµ |K ∈ P1 ,

∀K ∈ T µ (Ω)},

and satisfy the so-called inf-sup inequality (LBB condition): there exists a constant β > 0 such that βkqµ k0 ≤

(∇ · vµ , qµ ) , vµ ∈Xµ ,vµ ,0 k∇vµ k0 sup

(2.7)

∀qµ ∈ Mµ ,

where Pl (l ≥ 1) is the space of polynomials of degree less than or equal to l on the element K. It is well known that the finite element solution (ˆuµ , pˆ µ ) of problem (2.6) with the Taylor-Hood elements satisfies [3, 14, 63]   ku − uˆ µ k0 + µ k∇(u − uˆ µ )k0 + kp − pˆ µ k0 ≤ cµ s+1 (kuk s+1 + kpk s ),

s = 1, 2,

(2.8)

provided the solution (u, p) is regular enough so that the norms at the right-hand side make sense. We define V = {v ∈ X : (∇ · v, q) = 0, ∀q ∈ M},

Vµ = {vµ ∈ Xµ : (∇ · vµ , qµ ) = 0, ∀qµ ∈ Mµ },

and have the following estimate [2, 3] 1 inf k∇(v − vµ )k0 ≤ C(1 + ) inf k∇(v − vµ )k0 , β vµ ∈Xµ

vµ ∈Vµ

∀v ∈ V.

(2.9)

The error analysis we shall perform for our method will be for nonsingular solutions to the Navier-Stokes equations. u is called a nonsingular solution to the Navier-Stokes equations (2.1)-(2.3) if there is a σ = σ(u, ν) > 0 such that [64] inf sup

v∈V w∈V

ν(∇v, ∇w) + b(u, v, w) + b(v, u, w) ≥ σ > 0. k∇vk0 k∇wk0

It is easy to show that a nonsingular solution u to the Navier-Stokes equations is isolated and satisfies (cf. [2]) k∇uk0 ≤ ν−1 k f k−1 ,

|( f, v)| . v∈X,v,0 k∇vk0

k f k−1 = sup

(2.10)

We need the following lemma which was borrowed from [64]: Lemma 2.1. Let u be a nonsingular solution to the Navier-Stokes equations. Then there is a δ > 0 such that for any 0 ≤ α < δ, u0 and u00 ∈ V or Vµ with k∇(u − u0 )k0 < δ, k∇(u − u00 )k0 < δ satisfying inf sup

vµ ∈Vµ wµ ∈Vµ

(ν + α)(∇vµ , ∇wµ ) + b(u0 , vµ , wµ ) + b(vµ , u00 , wµ ) 1 ≥ σ. k∇vµ k0 k∇wµ k0 2 5

(2.11)

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

3. Finite element variational multiscale method 3.1. One-level finite element variational multiscale method The variational multiscale stabilization term is defined by the difference between a consistent and an underintegrated matrices involving the velocity gradient based on two local Gauss integrations at element level as follows:

G(uµ , vµ ) = α

X

K∈T µ (Ω)

here

R

K,s

Z

K,k

∇uµ · ∇vµ dx −

Z

K,1

!

∇uµ · ∇vµ dx ,

∀uµ , vµ ∈ Xµ ,

(3.1)

(·)dx denotes an appropriate Gauss integral over the element K which is exact for polynomials of degree s (

here s = k, 1 with k ≥ 2), and α > 0 is the user-defined stabilization parameter.

The stabilized finite element method based on stabilization term (3.1) for the Navier-Stokes equations is as follows. Method 0. One-level finite element VMS method [39]. Find (uµ , pµ ) ∈ Xµ × Mµ such that ν(∇uµ , ∇vµ ) + b(uµ , uµ , vµ ) − (∇ · vµ , pµ ) + (∇ · uµ , qµ ) + G(uµ , vµ ) = ( f, vµ ),

∀(vµ , qµ ) ∈ Xµ × Mµ .

(3.2)

We define R0 = {v ∈ L2 (Ω) : v|K ∈ P0 ,

∀K ∈ T µ (Ω)},

Lµ = Rd×d 0 ,

L = L2 (Ω)d×d ,

where P0 is the space of constants on the element K, and introduce the standard L2 -orthogonal projection Πµ : L → Lµ with the following properties:

(Πµ ∇u, v) = (∇u, v),

∀u ∈ X, v ∈ Lµ , ∀v ∈ X.

kΠµ ∇vk0 ≤ k∇vk0 ,

(3.3) (3.4)

Noting that k ≥ 2 in (3.1), the stabilization term (3.1) can be rewritten as G(uµ , vµ ) = α((I − Πµ )∇uµ , (I − Πµ )∇vµ ) = α(∇uµ , ∇vµ ) − α(Πµ ∇uµ , ∇vµ ).

(3.5)

For the equivalence of (3.1) and (3.5), we refer to Zheng et al. [39]. It is noted that (3.5) is just used for our theoretical analysis. Remark 3.1. As mentioned before, this paper is confined to the Taylor-Hood elements pair. However, from the definition of the stabilization term (3.1) and its equivalence with (3.5), we can see that the proposed variational multiscale method is applicable to all those finite element pairs where second order polynomials are used for the interpolation of velocity. For the stabilized numerical form (3.2), we have the following lemma [47]. 6

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Lemma 3.1. Assume that (u, p) is a nonsingular solution pair to the Navier-Stokes equations, and α tends to zero as µ tends to zero. Then there exists a µ0 > 0 such that for µ ≤ µ0 the solution (uµ , pµ ) defined by Method 0 satisfies (3.6)

νk∇uµ k0 ≤ k f k−1 , k∇(u − uµ )k0 + kp − pµ k0 ≤ cµ s + Cα,

s = 1, 2.

(3.7)

3.2. Two-level finite element variational multiscale method Now we give our two-level finite element VMS method as follows. Method I. Two-level finite element VMS method. 1. Find a coarse grid solution (uH , pH ) ∈ XH × MH such that ν(∇uH , ∇vH ) + b(uH , uH , v) − (∇ · vH , pH ) + (∇ · uH , qH ) + G(uH , vH ) = ( f, vH ),

(3.8)

∀(vH , qH ) ∈ XH × MH ,

2. Find a fine grid solution (uh , ph ) ∈ Xh × Mh such that ν(∇uh , ∇vh ) + b(uH , uh , vh ) + b(uh , uH , vh ) − (∇ · vh , ph ) + (∇ · uh , qh ) +G∗ (uh , vh ) = ( f, vh ) + b(uH , uH , vh ), Here, the stabilized term G∗ (uh , vh ) is defined as Z X Z ∗ h h G (u , vh ) = α ∇u · ∇vh dx − K∈T h (Ω)

K,k

K,1

(3.9)

∀(vh , qh ) ∈ Xh × Mh . !

∇uH · ∇vh dx ,

∀vh ∈ Xh ,

Remark 3.2. In the definition (3.10) of G∗ (uh , vh ) , that we use uH rather than uh in

R

K,1

(3.10)

∇uH · ∇vh dx is just for

simplicity of theoretical analysis and computations. One can definitely define a scheme using uh as [46] has done.

Theorem 3.1. Under the conditions of Lemma 3.1, the solution (uh , ph ) ∈ Xh × Mh defined by Method I satisfies ! 2 −1 h k∇(u − u )k0 ≤ C(β) 1 + (ν + α + 2ν Nk f k−1 ) inf k∇(u − wh )k0 wh ∈Xh σ ! √ 2 2 −1 + d inf kp − λh k0 + Nk∇(u − uH )k0 + 2αν k f k−1 , (3.11) λh ∈Mh σ √ d 1 kp − ph k0 ≤ (1 + ) inf kp − λh k0 + (ν + α + 2ν−1 Nk f k−1 )k∇(u − uh )k0 β λh ∈Mh β N 2α + k∇(u − uH )k20 + k f k−1 . (3.12) β νβ Proof. Setting (eh , ηh ) = (u − uh , p − ph ) and subtracting (3.9) from (2.6), we have (ν + α)(∇eh , ∇vh ) + b(eh , uH , vh ) + b(uH , eh , vh ) + b(u − uH , u − uH , vh ) −(∇ · vh , ηh ) + (∇ · eh , qh ) = α(∇u, ∇vh ) − α(Πh ∇uH , ∇vh ), ∀(vh , qh ) ∈ Xh × Mh . 7

(3.13)

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Let (wh , λh ) ∈ Vh × Mh is an approximation of (u, p). By setting qh = 0 in (3.13), limiting vh ∈ Vh , and splitting eh = χ − ϕh with χ = u − wh , ϕh = uh − wh , we get

(ν + α)(∇ϕh , ∇vh ) + b(ϕh , uH , vh ) + b(uH , ϕh , vh ) = (ν + α)(∇χ, ∇vh ) + b(χ, uH , vh ) + b(uH , χ, vh ) + b(u − uH , u − uH , vh ) − (∇ · vh , p − λh ) − α(∇u, ∇vh ) + α(Πh ∇uH , ∇vh ),

∀vh ∈ Vh ,

which, together with (2.5), (2.10), (3.4) and (3.6), leads to 1 [(ν + α)(∇ϕh , ∇vh ) + b(ϕh , uH , vh ) + b(uH , ϕh , vh )] k∇vh k0

≤ (ν + α)k∇χk0 + 2Nk∇uH k0 k∇χk0 + Nk∇(u − uH )k20 √ + dkp − λh k0 + α(k∇uk0 + k∇uH k0 )

≤ (ν + α + 2ν−1 Nk f k−1 )k∇χk0 + Nk∇(u − uH )k20 √ + dkp − λh k0 + 2αν−1 k f k−1 . The application of Lemma 2.1 yields 1 σk∇ϕh k0 ≤(ν + α + 2ν−1 Nk f k−1 )k∇χk0 2 √ + Nk∇(u − uH )k20 + dkp − λh k0 + 2αν−1 k f k−1 . By using the triangle inequality, taking the infimum over wh ∈ Vh and λh ∈ Mh , and using (2.9), we get the required

result (3.11).

On the other hand, taking qh = 0 in (3.13), applying (2.5), (2.10), (3.4) and (3.6), we get (∇ · vh , ph − λh ) = (∇ · vh , p − λh ) − (ν + α)(∇eh , ∇vh ) − b(eh , uH , vh ) − b(uH , eh , vh ) − b(u − uH , u − uH , vh ) + α(∇u, ∇vh ) − α(Πh ∇uH , ∇vh ) √  ≤ dkp − λh k0 + (ν + α + 2ν−1 Nk f k−1 )k∇eh k0 k∇vh k0   + Nk∇(u − uH )k20 + 2αν−1 k f k−1 k∇vh k0 ,

which, combing with the inf-sup condition and the triangle inequality, yields √ d 1 h kp − p k0 ≤ (1 + )kp − λµ k0 + (ν + α + 2ν−1 Nk f k−1 )k∇eh k0 β β N 2α 2 + k∇(u − uH )k0 + k f k−1 . β νβ Taking the infimum over λµ ∈ Mµ , we obtain the required result (3.12). 

Corollary 3.1. Under the conditions of Theorem 3.1, the solution (uh , ph ) ∈ Xh × Mh defined by Method I satisfies k∇(u − uh )k0 + kp − ph k0 ≤ ch s + C1 H 2s + C2 α, 8

s = 1, 2.

(3.14)

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

This follows by inserting the estimates (2.8) and (3.7) for the coarse grid solution into the right-hand side of estimators (3.11)-(3.12) and using the fact that α2 ≤ α. Remark 3.3. Corollary 3.1 indicates the typical choices of the algorithmic parameters: H = O(h1/2 )

α = O(h s ),

(3.15)

which ensure optimal accuracy of the solution (uh , ph ). 4. Numerical results In this section, we shall perform some numerical tests to verify the theory predictions and demonstrate the efficiency of the proposed method. The public domain software FreeFem++ [65] is used. In all experiments, the iterative tolerance for the nonlinear iterations on the coarse grid is set as 10−6 . When the number of nonlinear iterations exceeds 1000, we declare the failure of the iterative method. Within the nonlinear iterations, the stabilization term in the one and two-level VMS methods is treated explicitly as X Z G(uµj , vµ ) = α K∈T µ (Ω)

K,k

∇uµj

· ∇vµ dx −

Z

K,1

∇uµj−1

!

· ∇vµ dx ,

∀uµj , uµj−1 , vµ ∈ Xµ ,

where j denotes the number of nonlinear iterations. 4.1. Analytical solutions 4.1.1. A 2D example In this example, we set the boundary conditions and the body force f such that the exact solution to the NavierStokes equations is given by u1 = 10x2 (x − 1)2 y(y − 1)(2y − 1),

u2 = −10y2 (y − 1)2 x(x − 1)(2x − 1),

(4.1)

p = 10(3x2 + 3y2 − 2),

where the solution domain Ω = [0, 1] × [0, 1] ⊂ R2 .

Firstly, we compute the finite element solutions by our two-level VMS method (Method I) on uniform meshes

with size h = n−2 (n = 5, 6, · · · , 10) and corresponding H satisfying H = h1/2 , where ν = 0.01. For the solution given

by (4.1), the estimate (2.8) holds for s = 2. Therefore, according to Remark 3.3, the stabilization parameter is set as α = 0.1h2 . The Newton iterative method (cf. [66]) is used to solve the nonlinear Navier-Stokes problem on the coarse grid. The numerical results are listed in Table 1, where it is the nonlinear iterations count satisfying the stopping condition. From Table 1, we can see that optimal convergence rate both for the velocity and pressure is gotten by our method. The numerical results support theoretical predictions very well. 9

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Table 1: Errors of the approximate solutions by present two-level VMS method (Method I) for 2D example.

uH 1 rate

pL2 rate

it

CPU(s)

0.00499998

kp−ph k0 kpk0

0.0004

-

-

4

0.393797

1/6

0.002344

0.000192901

2.07758

2

3

0.765838

1/49

1/7

0.00124887

0.000104123

2.04223

2

3

1.41608

1/64

1/8

0.000727318

6.10352e-005

2.02435

2

3

2.41467

1/81

1/9

0.000452474

3.8104e-005

2.01486

2

3

4.03517

1/100

1/10

0.000296273

2.5e-005

2.00953

2

3

6.2094

h

H

1/25

1/5

1/36

k∇(u−uh )k0 k∇uk0

Table 2: Errors of the approximate solutions by one-level VMS method [39] for 2D example. k∇(u−uh )k0 k∇uk0

h

kp−ph k0 kpk0

uH 1 rate

pL2 rate

it

CPU(s)

1/25

0.00481599

0.0004

-

-

4

1.44564

1/36

0.00229973

0.000192901

2.02705

2

3

2.33488

1/49

0.00123586

0.000104123

2.01435

2

3

4.39916

1/64

0.000722851

6.10352e-005

2.00821

2

3

7.66604

1/81

0.00045074

3.8104e-005

2.00501

2

3

12.9069

1/100

0.00029553

2.5e-005

2.00322

2

3

20.2422

Table 3: Errors of the approximate solutions by the standard two-level finite element method [7] for 2D example.

uH 1 rate

pL2 rate

it

CPU(s)

0.00470373

kp−ph k0 kpk0

0.0004

-

-

3

0.386772

1/6

0.0022729

0.000192901

1.99454

2

3

0.707862

1/49

1/7

0.00122794

0.000104123

1.99714

2

3

1.39101

1/64

1/8

0.000720111

6.10352e-005

1.99836

2

3

2.36935

1/81

1/9

0.000449668

3.8104e-005

1.999

2

3

3.86676

1/100

1/10

0.000295068

2.5e-005

1.99935

2

3

6.02701

h

H

1/25

1/5

1/36

k∇(u−uh )k0 k∇uk0

Table 4: Errors of the approximate solutions by the two-level VMS method of [46] for 2D example.

uH 1 rate

pL2 rate

it

CPU(s)

0.00639076

kp−ph k0 kpk0

0.000400006

-

-

4

0.303886

1/6

0.00338126

0.000192904

1.74583

2

3

0.582545

1/49

1/7

0.00200278

0.000104125

1.6987

2

3

1.06057

1/64

1/8

0.00128418

6.10361e-05

1.66408

1.99999

3

1.81408

1/81

1/9

0.000873253

3.81048e-05

1.63714

1.99998

3

2.9921

1/100

1/10

0.000621236

2.50005e-05

1.61595

2.00001

3

4.66977

h

H

1/25

1/5

1/36

k∇(u−uh )k0 k∇uk0

10

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Table 5: Comparison of the two-level VMS methods for 2D example.

Present two-level VMS method

ν k∇(u−uh )k0 k∇uk0

kp−ph k0 kpk0

The two-level VMS method of [46] CPU(s)

k∇(u−uh )k0 k∇uk0

kp−ph k0 kpk0

CPU(s)

1

0.00227274

0.000192902

0.764227

0.00227287

0.000192904

0.586508

0.1

0.00227346

0.000192901

0.739596

0.00228663

0.000192904

0.571381

0.01

0.002344

0.000192901

0.765838

0.00338126

0.000192904

0.582545

0.001

0.00607635

0.000192901

0.810667

0.0260391

0.000192909

0.629731

0.0001

0.0561135

0.000192916

0.994974

0.320494

0.000193027

0.796023

Secondly, to evaluate the performance of our method, we also compute the finite element solutions by the one-level VMS method [39], the standard two-level finite element discretization method [7], and the two-level VMS method in [46] with the same stabilization parameter and on the same uniform grid meshes. The numerical results are listed in Tables 2 - 4. Comparing Tables 1 - 3, we see that the accuracy of the computed solutions by our two-level VMS method, the one-level VMS method [39] and the standard two-level finite element discretization method [7] is very comparable; in particular, for the pressure, the computed solutions are the same. However, compared to the one-level VMS method [39], our two-level VMS method takes much less computational time; it saves 72.8%, 67.2%, 67.8%, 68.5%, 68.7% and 69.3% CPU time when h = 1/25, 1/36, 1/49, 1/64, 1/81 and 1/100, respectively. While compared to the standard two-level discretization method [7], the stabilization used in our method does not degrade the accuracy of the approximate solutions. However, from Table 4, we can see that the two-level VMS method of [46] was not able to yield an optimal convergence rate for H 1 -error of the velocity, although it spent the least CPU time among the methods being compared. What is more, accuracy of the computed velocities by the method of [46] is about half lower than our method when the fine grid size h is smaller than 1/49. To further compare our method with that of [46], we compute the finite element solutions with different values of the viscosity ν = 1, 0.1, 0.01, 0.001 and 0.0001, where the meshes with h = 1/36, H = 1/6 are used. The numerical results are shown in Table 5, from which we can see that for the computed pressures, there is no obvious difference for all values of viscosity being tested between the two methods; however, when the viscosity is small (smaller than 0.01 in this numerical test), the computed velocities by our method are much better than those computed by the method of [46]. Thirdly, to test the coarse and fine grid scalings, we compute the finite element solutions by our method with √ √ different coarse grid sizes satisfying H = 21 h and H = 14 h. The numerical results are listed in Table 6 which show that optimal convergence rates are also obtained. Comparing Table 1 with Table 6, we see that there is no obvious increase in accuracy of the computed velocities, while accuracy of the computed pressures is the same as the coarse grid size gets smaller.

11

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Table 6: Errors of the approximate solutions by present two-level VMS method (Method I) with different grids scalings for 2D example.

uH 1 rate

pL2 rate

it

CPU(s)

0.00482324

kp−ph k0 kpk0

0.0004

-

-

4

0.566191

1/12

0.0023027

0.000192901

2.02763

2

3

0.970667

1/49

1/14

0.00123648

0.000104123

2.01691

2

3

1.66415

1/64

1/16

0.000723156

6.10352e-005

2.00851

2

3

2.71215

1/81

1/18

0.00045083

3.8104e-005

2.00595

2

3

4.29813

1/100

1/20

0.000295582

2.5e-005

2.00334

2

3

6.74215

1/25

1/20

0.00481605

0.0004

-

-

4

1.25366

1/36

1/24

0.00229977

0.000192901

2.02704

2

3

1.73613

1/49

1/28

0.00123588

0.000104123

2.01436

2

3

2.74356

1/64

1/32

0.000722866

6.10352e-005

2.00818

2

3

4.13202

1/81

1/36

0.000450745

3.8104e-005

2.00505

2

3

6.13894

1/100

1/40

0.000295532

2.5e-005

2.00323

2

3

8.73467

h

H

1/25

1/10

1/36

k∇(u−uh )k0 k∇uk0

4.1.2. A 3D example In this example, the exact solution of the Navier-Stokes equations is given by u1 = x2 (x − 1)2 [2y(y − 1)(2y − 1)z2 (z − 1)2 − 2y2 (y − 1)2 z(z − 1)(2z − 1)],

u2 = y2 (y − 1)2 [−2x(x − 1)(2x − 1)z2 (z − 1)2 + 2x2 (x − 1)2 z(z − 1)(2z − 1)],

u3 = z2 (z − 1)2 [2x(x − 1)(2x − 1)y2 (y − 1)2 − 2x2 (x − 1)2 y(y − 1)(2y − 1)],

(4.2)

p = (2x − 1)(2y − 1)(2z − 1),

where Ω = [0, 1] × [0, 1] × [0, 1].

We compute the finite element solutions by our two-level VMS method and the standard two-level finite element

discretization method [7]. In this test, the stabilization parameter is set as α = 0.1h2 , and the Newton iterative method is used for the coarse grid nonlinear problem, where the viscosity is chosen as ν = 0.1. The numerical results are shown in Table 7 and Table 8, respectively. From Table 7 we can see that the numerical results are in good agreement with the theoretical predictions. A comparison of Table 7 with Table 8 shows that there is no big difference in accuracy of the computed solutions by our two-level VMS method and the standard two-level finite element discretization method [7]. 4.2. Lid-driven cavity flow 4.2.1. 2D case Here we consider the 2D lid-driven cavity flow problem which is defined in the unit square. With zero source external force, velocities are zero on all boundaries except the top one (the lid), which has a driving horizontal velocity 12

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Table 7: Errors of the approximate solutions by present two-level VMS method (Method I) for 3D example.

h

H

k∇(u − uh )k0

kp − ph k0

uH 1 rate

pL2 rate

it

CPU(s)

1/4

1/2

0.0416349

0.016257

-

-

5

0.627475

1/6

1/3

0.0141966

0.0072019

2.65359

2.00801

4

2.10665

1/8

1/4

0.00625251

0.00404505

2.85043

2.00517

4

6.68291

1/10

1/5

0.00328606

0.00258691

2.88286

2.00333

3

17.5246

1/12

1/6

0.00194906

0.00179576

2.86495

2.00213

3

47.1655

Table 8: Errors of the approximate solutions by the standard two-level finite element method [7] for 3D example.

h

H

k∇(u − uh )k0

kp − ph k0

uH 1 rate

pL2 rate

it

CPU(s)

1/4

1/2

0.0410773

0.0161897

-

-

2

0.482868

1/6

1/3

0.0136797

0.00719583

2.71181

1.99986

2

1.80637

1/8

1/4

0.00608648

0.00404418

2.81506

2.00517

2

6.05043

1/10

1/5

0.00322319

0.00258672

2.84884

2.00269

2

17.9948

1/12

1/6

0.00192166

0.00179571

2.83665

2.0019

2

46.057

set to unity; see Figure 1. The Reynolds number for this problem is defined as Re = UL/ν, where U is the velocity of the top lid and L is the length of the side wall. We adopt the Oseen iterative method (cf. [66, 36]) to solve the nonlinear problem on the coarse grid. Uniform coarse and fine meshes with size H = 1/64 and h = 1/128 are used. For this problem, the estimate (2.8) hold for s = 1. Therefore, the stabilization parameter is chosen as α = 0.1h. We compute an approximate solution at Re = 5000, 7500 and 10000 by our two-level VMS method (Method I), the one-level VMS method in [39], and the standard two-level discretization method [7], and compare our numerical results with those of Ghia et al. [67]. The numerical test shows that the standard two-level finite element discretization method [7] is not able to yield a solution due to the divergence of the nonlinear iterations on the coarse grid, while the one and two-level VMS stabilization methods work well. Figures 2-4 plot the computed u1 -velocity along the vertical centerline and u2 -velocity along the horizontal centerline, compared with those of Ghia et al. [67] where a much finer 257 × 257 grid mesh was used. From Figures 2-4, we can see that the accuracy of the computed solutions is comparable to those of Ghia et al. [67]. There is no obvious

difference between the numerical results computed by the one-level finite element VMS method and our two-level VMS method. However, our two-level VMS method saves a large amount of computational time (see Table 9). It is noted that we also have performed numerical experiment with different values of stabilization parameter α = h, 0.5h, 0.1h, 0.05h and 0.01h. The numerical tests show that at Re = 5000, our two-level VMS method works well for all values of the stabilization parameter being tested, while at Re = 7500 and Re = 10000, the nonlinear 13

u1=1, u2=0

L=1

u1=0, u2=0

u1=0, u2=0

u1=0, u2=0

L=1

Figure 1: Schematic diagram of the lid-driven cavity flow.

Re=5000

Re=5000

1

0.6

0.9 0.4 0.8 0.2

0.7 u −velocity

0.6 0.5

2

y

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

0.4 0.3

−0.2

−0.4 One−level VMS Two−level VMS Ghia et al.

0.2 0.1 0 −0.5

0

0

0.5

One−level VMS Two−level VMS Ghia et al.

−0.6

−0.8

1

0

0.2

0.4

u1−velocity

0.6

0.8

1

x

Figure 2: Comparison of u1 -velocity profiles along the vertical centerline (left) and u2 -velocity profiles along the horizontal centerline (right) for 2D lid-driven cavity flow at Re = 5000.

Table 9: Computaional time of the VMS methods for 2D lid-driven cavity flow problem.

Re

5000

7500

One-level VMS [39]

638.798

940.683

1282.07

Present two-level VMS

151.195

213.281

276.99

Saved time

76.3%

77.3%

78.4%

14

10000

Re=7500

Re=7500

1

0.6

0.9 0.4 0.8 0.2

0.7 u −velocity

0.5

2

y

0.6

0.4 0.3

0.1 0 −0.5

0

−0.2

−0.4 One−level VMS Two−level VMS Ghia et al.

0.2

0

0.5

One−level VMS Two−level VMS Ghia et al.

−0.6

−0.8

1

0

0.2

0.4

0.6

u −velocity

0.8

1

x

1

Figure 3: Comparison of u1 -velocity profiles along the vertical centerline (left) and u2 -velocity profiles along the horizontal centerline (right) for 2D lid-driven cavity flow at Re = 7500.

Re=10000

Re=10000

1

0.6

0.9 0.4 0.8 0.2

0.7 u −velocity

0.6 0.5

2

y

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

0.4 0.3

−0.2

−0.4

One−level VMS Two−level VMS Ghia et al.

0.2

0

One−level VMS Two−level VMS Ghia et al.

−0.6

0.1 0 −0.5

0

0.5

−0.8

1

u1−velocity

0

0.2

0.4

0.6

0.8

1

x

Figure 4: Comparison of u1 -velocity profiles along the vertical centerline (left) and u2 -velocity profiles along the horizontal centerline (right) for 2D lid-driven cavity flow at Re = 10000.

15

Table 10: Numbers of nonlinear iterations of present two-level VMS method for 2D lid-driven cavity flow with different values of stabilization parameter.

h

0.5h

0.1h

0.05h

Re = 5000

453

253

74

49

456

Re = 7500

679

380

108

64

-

Re = 10000

918

513

142

80

-

α

1

0.8

0.8

0.6

0.6

0.01h

Y

1

Y

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

0.4

0.4

0.2

0.2

0

0

0.2

0.4

0.6

0.8

0

1

0

X

0.2

0.4

0.6

0.8

1

X

Figure 5: Computed streamlines (left) and pressure contours (right) for 2D lid-driven cavity flow at Re = 5000.

iterations on the coarse grid fail to converge when α = 0.01h. There is no big difference between the successfully computed results with different values of stabilization parameter although the number of nonlinear iterations are largely different (see Table 10 for detailed information). Figures 5-7 depict the computed streamlines and pressure contours at different Reynolds numbers with α = 0.05h, which demonstrated the effectiveness of our two-level VMS method. 4.2.2. 3D case The fourth numerical example is the 3D lid-driven cavity flow problem defined in the unit cube [0, 1]×[0, 1]×[0, 1]. In this problem, the lid face is imposed an unit driving horizontal velocity as boundary conditions, while on all other faces, homogeneous Dirichlet boundary conditions are set. We set h = 1/10 and H = 1/5, and compute the finite element solutions by our two-level VMS method, where the stabilization parameter is chosen as α = 0.1h. Figures 8-10 plot the cut 2D pictures at y = 0.2, 0.5 and 0.8 of the computed velocity fields and pressure contours at Re = 100, 500 and 1000, respectively, which illustrate the effectiveness of our proposed method. 16

0.8

0.8

0.6

0.6

Y

1

Y

1

0.4

0.4

0.2

0.2

0

0

0.2

0.4

0.6

0.8

0

1

0

0.2

0.4

0.6

X

0.8

1

X

Figure 6: Computed streamlines (left) and pressure contours (right) for 2D lid-driven cavity flow at Re = 7500.

1

0.8

0.8

0.6

0.6

Y

1

Y

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

0.4

0.4

0.2

0.2

0

0

0.2

0.4

0.6

0.8

0

1

0

0.2

0.4

0.6

X

0.8

1

X

Figure 7: Computed streamlines (left) and pressure contours (right) for 2D lid-driven cavity flow at Re = 10000. cut y = 0.2

cut y = 0.5

cut y = 0.8

Figure 8: Cut 2D pictures of the computed velocity fields and pressure contours for 3D lid-driven cavity flow at Re = 100: y = 0.2, 0.5 and 0.8 (from left to right).

17

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

cut y = 0.2

cut y = 0.5

cut y = 0.8

Figure 9: Cut 2D pictures of the computed velocity fields and pressure contours for 3D lid-driven cavity flow at Re = 500: y = 0.2, 0.5 and 0.8 (from left to right). cut y = 0.2

cut y = 0.5

cut y = 0.8

Figure 10: Cut 2D pictures of the computed velocity fields and pressure contours for 3D lid-driven cavity flow at Re = 1000: y = 0.2, 0.5 and 0.8 (from left to right).

4.3. Backward-facing step flow In this test case, we consider the backward-facing step flow problem defined on a long channel [0, 30]×[0.5, 0.5]. The boundary conditions are shown in Figure 11. The Reynolds number for this problem is defined as Re = Uave L/ν, where Uave (= 1) is the average velocity at the inlet boundary and L(= 1) is the channel height. The uniform meshes sizes are set as H = 1/12, h = 1/24, and the stabilization parameter is chosen as α = 0.004h, where Oseen iterative method (cf. [66]) is used for the nonlinear problem on the coarse grid. We compute an approximate solution at Re = 800 by our two-level VMS and compare the numerical results with those of Gartling [68]. In Figure 12, the computed velocity and pressure across the channel at x = 7 and 15

(0, 0.5)

u1=24y (0.5-y) u2=0

u1=u2=0 L=1

u1=u2=0 (0, -0.5)

(30, 0.5) -p+® w u1/ w x=0 u2=0

u1=u2=0

Figure 11: Schematic diagram of the backward-facing step flow.

18

(30, -0.5)

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.2

0.1

0.1

0.1

−0.1

0.5 Two−level VMS (x=7) Two−level VMS (x=15) Gartling (x=7) Gartling (x=15)

0.4

0

0

−0.1

−0.1

−0.2

−0.2

−0.2

−0.3

−0.3

−0.3

−0.4

−0.4

−0.4

−0.5 −0.2

−0.5 −20

0

0.2

0.4 0.6 u1−velocity

0.8

1

1.2

−15

−10

Two−level VMS (x=7) Two−level VMS (x=15) Gartling (x=7) Gartling (x=15)

0.3

y

Two−level VMS (x=7) Two−level VMS (x=15) Gartling (x=7) Gartling (x=15)

0

y

y

−5 u2−velocity

0

−0.5

5

0.16

0.17

0.18

−3

x 10

0.19

0.2 0.21 pressure

0.22

0.23

0.24

0.25

Figure 12: Comparison of u1 -velocity, u2 -velocity and pressure (from left to right) profiles at various locations for backward-facing step flow.

0.01

0.2

0.005

0.15

0 y

0.25

y

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

−0.005

0.1 Lower wall Upper wall

Lower wall Upper wall −0.01

0.05

0

0

5

10

15 x

20

25

−0.015

30

0

5

10

15 x

20

25

30

Figure 13: Pressure (left) and shear stress (right) profiles along upper and lower channel walls for backward-facing step flow.

compared with those of Gartling [68] are shown. From Figure 12, we can see that our numerical results coincide well with those of Gartling [68]. It is worth mentioning that due to the different solutions to uniquely determining the approximate pressure, there is a constant difference between our computed pressure with that of Gartling’s. For the sake of comparison, our pressure data presented in Figure 12 were adjusted by making the computed pressure equal to that of Gartling at the lower channel wall point (x, y) = (7.0, −0.5).

Figure 13 describes the computed pressure and shear stress along the upper and lower channel walls, while Figure

14 plots the computed streamlines and pressure contours, which are also in good agreement with those of Gartling [68].

19

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Figure 14: Computed streamlines (top) and pressure contours (bottom) for backward-facing step flow.

5. Conclusions Combining the best algorithmic features of two-grid discretization approch and a recent VMS method based on two local Gauss integrations, we proposed a two-level finite element VMS method for the convection dominated incompressible Navier-Stokes equations. Convergence theory and algorithmic parameter scalings of the method were developed. Numerical results verified the theoretical results and demonstrated the efficiency of the proposed method. Compared to the one-level VMS method, our two-level VMS method can save a large amount of CPU time, while compared to the standard two-level finite element discretization method, our method can simulate high Reynolds number flows for which the the standard two-level discretization method fails to work. References [1] R. Temam, Navier-Stokes Equations: Theory and Numerical Analysis, North-Holland, Amsterdam, 1984. [2] V. Girault, P.A. Raviart, Finite Element Approximation of the Navier-Stokes Equations, Springer-Verlag, Berlin, 1979. [3] V. Girault, P.A. Raviart, Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms, Springer-Verlag, Berlin Heidelberg, 1986. [4] H.C. Elman, D.J. Silvester, A.J. Wathen, Finite Elements and Fast Iterative Solvers: with Applications in Incompressible Fluid Dynamics, Oxford University Press, Oxford, 2005. [5] J.C. Xu, A novel two-grid method for semilinear elliptic equations, SIAM J. Sci. Comput. 15(1994) 231-237. [6] J.C. Xu, Two-grid discretization techniques for linear and nonlinear PDEs, SIAM J. Numer. Anal. 33(5)(1996) 1759-1777. [7] W. Layton, A two level discretization method for the Navier-Stokes equations, Comput. Math. Appl. 5(26)(1993) 33-38. [8] W. Layton, H.W.J. Leferink, Two-level Picard and modified Picard methods for the Navier-Stokes equations, Appl. Math. Comput. 69(1995) 263-274. [9] W. Layton, H.W.J. Lenferink, A multilevel mesh independence principle for the Navier-Stokes equations, SIAM J. Numer. Anal. 33(1)(1996) 17-30. [10] W. Layton, H.K. Lee, J. Peterson, Numerical solution of the stationary Navier-Stokes equations using a multilevel finite element method, SIAM J. Sci. Comput. 20(1)(1998) 1-12. [11] W. Layton, L. Tobiska, A two-level method with backtraking for the Navier-Stokes equations, SIAM J. Numer. Anal. 35(1998) 2035-2054. [12] V. Girault, J.L. Lions, Two-grid finite element scheme for the steady Navier-Stokes equations in polyhedra, Portugal. Math. 58(2001) 25-57. [13] Y.N. He, K.T. Li, Two-level stabilized finite element methods for the steady Navier-Stokes problem, Computing 74(2005) 337-351.

20

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[14] Y.N. He, A.W. Wang, A simplified two-level method for the steady Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 197(2008) 1568-1576. [15] Y.N. He, J. Li, Two-level methods based on three correction for the 2D/3D steady Navier-Stokes equations, Int. J. Numer. Anal. Model. Ser. B 2(1)(2011) 42-56. [16] Y.N. He, Y. Zhang, Y.Q. Shang, H. Xu, Two-level Newton iterative method for the 2D/3D steady Navier-Stokes equations, Numer. Methods Partial Differential Equations 28(5)(2012) 1620-1642. [17] Q.F. Liu, Y.R. Hou, A two-level defect-correction method for Navier-Stokes equations, Bull. Aust. Math. Soc. 81(2010) 442-454. [18] Q.F. Liu, Y.R. Hou, A two-level finite element method for the Navier-Stokes equations based on a new projection, Appl. Math. Modelling 34(2)(2010) 383-399. [19] V. Girault, J.L. Lions, Two-grid finite element scheme for the transient Navier-Stokes problem, Math. Model. Numer. Anal. 35(2001) 945-980. [20] M.A. Olshanskii, Two-level method and some a priori estimates in unsteady Navier-Stokes calculations, J. Comp. Appl. Math. 104(1999) 173-191. [21] Y.N. He, Two-level method based on finite element and Crank-Nicolson extrapolation for the time-dependent Navier-Stokes equations, SIAM J. Numer. Anal. 41(2003) 1263-1285. [22] Y.N. He, A two-level finite element Galerkin method for the nonstationary Navier-Stokes equations, I: Spatial discretization, J. Comput. Math. 22(2004) 21-32. [23] Y.N. He, H.L. Miao, C.F. Ren, A two-level finite element Galerkin method for the nonstationary Navier-Stokes equations, II: Time discretization, J. Comput. Math. 22(2004) 33-54. [24] Y.N. He, K.M. Liu, A multi-level finite element method in space-time for the Navier-Stokes equations, Numer. Methods Partial Differential Equations 21(2005)1052-1078. [25] Y.N. He, K.M. Liu, W.W. Sun, Multi-level spectral Galerkin method for the Navier- Stokes equations I: spatial discretization, Numer. Math. 101(2005) 501-522. [26] Y.N. He, K.M. Liu, Multi-level spectral Galerkin method for the Navier-Stokes equations II: time discretization, Adv. Comput. Math. 25(2006) 403-433. [27] Y.R. Hou, L.Q. Mei, Full discrete two-level correction scheme for Navier-Stokes equations, J. Comput. Math. 26(2008) 209-226. [28] H. Abboud, V. Girault, T. Sayah, A second order accuracy for a full discretized time-dependent Navier-Stokes equations by a two-grid scheme, Numer. Math. 114(2009) 189-231. [29] Y.N. He, J.C. Xu, A.H. Zhou, Local and parallel finite element algorithms for the Navier-Stokes problem, J. Comput. Math. 24(3)(2006) 227-238. [30] Y.N. He, L.Q. Mei, Y.Q. Shang, J. Cui, Newton iterative parallel finite element algorithm for the steady Navier-Stokes equations, J. Sci. Comput. 44(1)(2010) 92-106. [31] Y.C. Ma, Z.M. Zhang, C.F. Ren, Local and parallel finite element algorithms based on two-grid discretization for the stream function form of Navier-Stokes equations, Appl. Math. Comput. 175(2006) 786-813. [32] F.Y. Ma, Y.C. Ma, W.F. Wo, Local and parallel finite element algorithms based on two-grid discretization for steady Navier-Stokes equations, Appl. Math. Mech. -Engl. Ed. 28(1)(2007) 27-35. [33] Y.Q. Shang, Z.D. Luo, A parallel two-level finite element method for the Navier-Stokes equations, Appl. Math. Mech. -Engl. Ed. 31(11)(2010) 1429-1438. [34] Y.Q. Shang, A parallel two-level linearization method for incompressible flow problems, Appl. Math. Lett. 24(2011) 364-369. [35] Y.Q. Shang, Y.N. He, Z.D. Luo, A comparison of three kinds of local and parallel finite element algorithms based on two-grid discretizations for the stationary Navier-Stokes equations, Comput. Fluids 40(2011) 249-257. [36] Y.Q. Shang, Y.N. He, A parallel Oseen-linearized algorithm for the stationary Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 209-212(2012) 172-183. [37] Y.Q. Shang, A two-level subgrid stabilized Oseen iterative method for the steady Navier-Stokes equations, J. Comput. Phys. 233 (2013)

21

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

210-226. [38] E. Erturk, T. Corke, C. Gokcol, Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers, Int. J. Numer. Methods Fluids 48(2005) 747-774. [39] H.B. Zheng, Y.R. Hou, F. Shi, L.N. Song, A finite element variational multiscale method for incompressible flows based on two local gauss integrations, J. Comput. Phys. 228(2009) 5961-5977. [40] J. Li, Y.N. He, A stabilized finite element method based on two local Gauss integrations for the Stokes equations, J. Comp. Appl. Math. 214 (2008) 58-65. [41] Y.N. He, J. Li, A stabilized finite element method based on local polynomial pressure projection for the stationary Navier-Stokes equations, Appl. Numer. Math. 58 (2008) 1503-1514. [42] H. B. Zheng, Y. R. Hou, F. Shi, Adaptive finite element variational multiscale method for incompressible flows based on two local Gauss integrations, J. Comput. Phys. 229 (2010) 7030-7041. [43] F. Shi, H.B. Zheng, J. Yu, Y. Li, On the convergence of Variational multiscale methods based on Newton’s iteration for the incompressible flows, Appl. Math. Model. 38(2014) 5726-5742. [44] Y.Q. Shang, Error analysis of a fully discrete finite element variational multiscale method for time-dependent incompressible Navier-Stokes equations, Numer. Methods Partial Differential Equations 29(6)(2013) 2025-2046. [45] Y. Jiang, L. Q. Mei, H. M. Wei, W. J. Tian, J. T. Ge, A finite element variational multiscale method based on two local Gauss integrations for stationary conduction-convection problems, Math. Problems Eng. (2012), article ID 747391. [46] Y. Li, L.Q. Mei, Y. Li, K. Zhao, A two-level variational multiscale method for incompressible flows based on two local Gauss integrations, Numer. Methods Partial Differential Equations, 29 (2013) 1986-2003. [47] Y.Q. Shang, A parallel two-level finite element variational multiscale method for the Navier-Stokes equations, Nonlin. Anal. 84 (2013) 103-116. [48] C. Xie, H. B. Zheng, A parallel variational multiscle method for incompressible flows based on the partition of unity, Int. J. Numer. Anal. Model. 11(4)(2014) 854-865. [49] I. Babu sˇka, J. Melenk, The partition of unity method, Int. J. Numer. Methods Eng. 40(1997) 727-758. [50] H.B. Zheng, J.P. Yu, F. Shi, Local and parallel finite element method based on the partition of unity for incompressible flow, J. Sci. Comput. 65(2)(2015) 512-532. [51] J.P. Yu, F. Shi, H.B. Zheng, Local and parallel finite element method based on the partition of unity for the Stokes problem, SIAM J. Sci. Comput. 36(5)(2014) C547-C567. [52] W. Layton, A connection between subgrid-scale eddy viscosity and mixed methods, Appl. Math. Comput. 133(1)(2002) 147-157. [53] S. Kaya, B. Rivi`ere, A two-grid stabilization method for solving the steady-state Navier-Stokes equations, Numer. Methods Partial Differential Equations 22(2005) 728-743. [54] V. John, S. Kaya, A finite element variational multiscale method for the Navier-Stokes equations, SIAM J. Sci. Comput. 26 (2005) 1485-1503. [55] K.J. Galvin, New subgrid artificial viscosity Galerkin methods for the Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 200(2011) 242-250. [56] J. Blasco, R. Codina, Stabilized finite element method for the transient Navier-Stokes equations based on a pressure gradient projection, Comput. Methods Appl. Mech. Engrg. 182 (2000) 277-300. [57] M. Braack, E. Burman, Local projection stabilization for the Oseen problem and its interpretation as a variational multiscale method, SIAM J. Numer. Anal. 43 (2000) 2544-2566. [58] T. Chacon ´ Rebollo, R. Lewandowski, Mathematical and Numerical Foundations of Turbulence Models and Applications, Birkh¨auser Basel, 2014. [59] S. Ganesan, G. Matthies, L. Tobiska, Local projection stabilization with equal order interpolation applied to the Stokes problem, Math. Comp. 77 (2008) 2039-2060. [60] J.G. Heywood, R. Rannacher, Finite element approximation of the nonstationary Navier-Stokes problem I: Regularity of solutions and second-

22

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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

order error estimates for spatial discretization, SIAM J. Numer. Anal. 19(2)(1982) 275-311. [61] J. Maubach, Local bisection refinement for n-simplicial grids generated by reflection, SIAM J. Sci. Comput. 16(1995) 210-227. [62] P. Hood, C. Taylor, A numerical solution of the Navier-Stokes equations using the finite element technique, Comput. Fluids 1(1973) 73-100. [63] A. Quarteroni, A. Valli, Numerical Approximation of Partial Differential Equations, Spring-Verlag, Berlin, 1997. [64] S. Kaya, W. Layton, B. Rivi`ere, Subgrid stabilized defect correction methods for the Navier-Stokes equations, SIAM J. Numer. Anal. 44 (2006) 1639–1654. [65] F. Hecht, New development in FreeFem++, J. Numer. Math. 20 (3-4)(2012) 251-265. [66] Y. He, J. Li, Convergence of three iterative methods based on finite element discretization for the stationary Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 198(2009) 1351-1359. [67] U. Ghia, K. Ghia, C. Shin, High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method, J. Comput. Phys. 48(1982) 387-411. [68] D.K. Gartling, A test problem for outflow boundary conditions-flow over a backward-facing step, Int. J. Numer. Methods Fluids 11(1990) 953-967.

23