Stabilized Multiscale Nonconforming Finite Element Method for the Stationary Navier-Stokes Equations

Abstract

We consider a stabilized multiscale nonconforming finite element method for the two-dimensional stationary incompressible Navier-Stokes problem. This method is based on the enrichment of the standard polynomial space for the velocity component with multiscale function and the nonconforming lowest equal-order finite element pair. Stability and existence uniqueness of the numerical solution are established, optimal-order error estimates are also presented. Finally, some numerical results are presented to validate the performance of the proposed method.

1. Introduction

As the development of science and technology, finite element method has became an important and powerful tool for the complex fluid problems, such as for the Navier-Stokes equations. It is well known that the pressure and velocity pairs satisfy the discrete Inf-Sup condition [1] that plays the key role for simulating the Navier-Stokes equations. However, some unstable mixed finite element pairs which violate the so-called Inf-Sup condition are also popular, see [2–4]. In order to overcome this restriction, various of stabilized methods have been proposed, including the bubble condensation-based methods [5], pressure projection method (PPM) [6–8], the local Gauss integration method (LGIM) [9–11], multiscale method [12, 13], macroelement stabilized method [3, 14], and so on. Most of these stabilized methods necessarily need to introduce the stabilization parameters either explicitly or implicitly. In addition, some of these techniques are conditionally stable or are of suboptimal accuracy. Therefore, the development of mixed finite element methods free from stabilization parameters has become increasingly important.

In 2005, Franca et al. gave a new multiscale method for the reaction-diffusion equation in [15]. The chief characteristic of their method is to use the Petrov-Galerkin approach to split the solution into two parts, and the trial function space is enriched with an unstable bubble-like function, which is the solution to a local problem. Later, Barrenechea and Valentin [13] considered the relationship between the enriched multiscale method and stabilized techniques for generalized Stokes problem based on the 𝑃₁-𝑃₁ pair. By enriching the velocity space with an unusual bubble function, Araya et al. established the convergence for the Stokes problem in [16], their method is different from usual residual free bubble method in [5], in which one should choose local basis functions to enrich the standard finite element spaces by solving some local problem analytically. Furthermore, the method proposed in [15] can also be used to treat the unsteady reaction-diffusion problem (see [17]).

Compared with conforming finite element method, the nonconforming finite element methods are more popular due to their simplicity and small support sets of basis functions. Crouzeix and Ravizrt in [18] used the nonconforming piecewise linear velocity and a piecewise constant pressure to solve the Stokes equations. In this paper, motivated by the ideas of [13, 15, 16], we will use the Petrov-Galerkin approach based on the nonconforming velocity space to handle with the steady Navier-Stokes equations. The main differences between [13, 15, 16, 19] and this work lie in the following: (i) the finite element spaces of velocity are different; it is nonconforming element in this paper; (ii) the treated problems are different; we consider the nonlinear problem; (iii) the finite element pairs are different; the NCP₁-𝑃₁ pair is used in this paper.

The outline of this work is arranged as follows. In the following section, the abstract functional setting for steady Navier-Stokes equations is recalled. Section 3 is devoted to derive the general form of enriched multiscale method based on the NCP₁-𝑃₁ pair. After providing the stability and existence uniqueness for the approximation solution, the optimal error estimates are established in Section 4. In Section 5, Some numerical results are presented to verify the established theoretical analysis. Finally, Some conclusions are made in Section 6.

2. Preliminaries

Let Ω be an open bounded domain of ℝ² with Lipschitz continuous boundary 𝜕Ω and satisfy a further condition stated in (A1) below. The incompressible stationary Navier-Stokes equations with the homogeneous Dirichlet boundary condition are ()𝑢|−𝜈Δ𝑢+𝑢⋅∇𝑢+∇𝑝=𝑓inΩ,div𝑢=0,inΩ,_𝜕Ω=0on𝜕Ω, (2.1) where 𝑢=(𝑢₁(𝑥),𝑢₂(𝑥))^𝑇 represents the velocity, 𝑝=𝑝(𝑥) the pressure, 𝑓=𝑓(𝑥)∈𝐿²(Ω)² the prescribed body force, 𝜈>0 the viscosity coefficient.

In order to introduce the variational formulation for problem (2.1), we set 𝑋=𝐻¹₀(Ω)²,𝑌=𝐿²(Ω)²,𝐷(𝐴)=𝐻²(Ω)²∩𝑋,𝑀=𝐿²₀(Ω)=𝑞∈𝐿²(Ω)∶_Ω𝑞𝑑𝑥=0. (2.2)

The standard notations of Sobolev space 𝑊^𝑚,𝑟(Ω) are used. To simplify, we use 𝐻^𝑚(Ω) instead of 𝑊^𝑚,𝑟(Ω) as 𝑟=2 and ‖⋅‖_𝑚 for ‖⋅‖_𝑚,2. The spaces 𝐿²(Ω)^𝑚(𝑚=1,2) are endowed with the usual 𝐿²-scalar product (⋅,⋅) and 𝐿²-norm ‖⋅‖₀. The spaces 𝐻¹₀(Ω) and 𝑋 are equipped with the scalar product (∇𝑢,∇𝑣) and the norm |𝑢|²_1,Ω, 𝑢, 𝑣∈𝐻¹₀(Ω)(or𝑋).

Define 𝐴𝑢=−Δ𝑢 is the operator associated with the Navier-Stokes problem, it is positive self-adjoint operator from 𝐷(𝐴) onto 𝑌.

Introducing the bilinear operator 𝐵()=()1𝑢,𝑣𝑢⋅∇𝑣+2(div𝑢)𝑣∀𝑢,𝑣∈𝑋, (2.3) and defining a trilinear form on 𝑋×𝑋×𝑋 as follows: 𝑏(𝑢,𝑣,𝑤)=⟨𝐵(𝑢,𝑣),𝑤⟩_{𝑋′×𝑋}=))+1((𝑢⋅∇𝑣,𝑤2))=1((div𝑢𝑣,𝑤2))−1((𝑢⋅∇𝑣,𝑤2((𝑢⋅∇)𝑤,𝑣). (2.4)

The variational formulation of problem (2.1) reads as: find (𝑢,𝑝)∈(𝑋,𝑀) such that for all (𝑣,𝑞)∈(𝑋,𝑀)𝑎(𝑢,𝑣)−𝑑𝑣,𝑝+𝑑𝑢,𝑞+𝑏(𝑢,𝑢,𝑣)=𝑓,𝑣, (2.5) where 𝑎()()=,𝐵𝑢,𝑣=𝜈∇𝑢,∇𝑣,𝑑𝑣,𝑞=−𝑣,∇𝑞𝑞,div𝑣₀𝑢,𝑝;𝑣,𝑞=𝑎(𝑢,𝑣)−𝑑𝑣,𝑝+𝑑𝑢,𝑞. (2.6)

Clearly, the bilinear forms 𝑎(⋅,⋅) and 𝑑(⋅,⋅) are continuous on 𝑋×𝑋 and 𝑋×𝑀, respectively. Moreover, 𝑑(⋅,⋅) also satisfies (see [20]): sup_{0≠𝑣∈𝑋}𝑑𝑣,𝑞|𝑣|_1,Ω≥𝛽𝑞_0,Ω, (2.7) where 𝛽 is a positive constant depending only on Ω.

It is easy to verify that 𝐵₀ satisfies the following important properties for all (𝑢,𝑝), (𝑣,𝑞)∈(𝑋,𝑀) (see [1]): 𝐵₀𝑢,𝑝;𝑢,𝑝=𝜈‖𝑢‖²_1,Ω,𝐵₀𝑢,𝑝;𝑣,𝑞≤𝐶‖𝑢‖_1,Ω+𝑝_0,Ω‖𝑣‖_1,Ω+𝑞_0,Ω𝛽₀,‖𝑢‖_1,Ω+𝑝_0,Ω≤sup_{(𝑣,𝑞)∈(𝑋,𝑀)}𝐵₀𝑢,𝑝;𝑣,𝑞‖𝑣‖_1,Ω+𝑞_0,Ω, (2.8)(2.9)(2.10) where 𝛽₀>0 is a constant. Here and below, the letter 𝐶 (with or without subscript) denotes a generic positive constant, depending at most on the data 𝜈, Ω and 𝑓. Furthermore, the following estimates about 𝑏(⋅,⋅,⋅) are hold [1, 20]: 𝑏()(),|𝑏(≤1𝑢,𝑣,𝑤=−𝑏𝑢,𝑤,𝑣𝑢,𝑣,𝑤)|2𝑐₀‖𝑢‖_0,Ω^1/2|𝑢|_1,Ω^1/2|𝑣|_1,Ω‖𝑤‖^1/2_0,Ω|𝑤|^1/2_1,Ω+‖𝑣‖^1/2_0,Ω|𝑣|^1/2_1,Ω|𝑤|_1,Ω, (2.11)(2.12) for all 𝑢, 𝑣, 𝑤∈𝑋 and |𝑏(𝑢,𝑣,𝑤)|+|𝑏(𝑣,𝑢,𝑤)|+|𝑏(𝑤,𝑢,𝑣)|≤𝐶|𝑢|_1,Ω‖𝐴𝑣‖_0,Ω‖𝑤‖_0,Ω, (2.13) for all 𝑢∈𝑋,𝑣∈𝐷(𝐴), 𝑤∈𝑌.

As mentioned above, a further assumption about Ω is needed (see [1]).

(A1) Assume that Ω is regular so that the unique solution (𝑣,𝑞)∈(𝑋,𝑀) of the steady Stokes equations −𝜈Δ𝑣+∇𝑞=𝑔,div𝑣=0inΩ,𝑣|_𝜕Ω=0, (2.14) for a prescribed 𝑔∈𝑌 exists and satisfies ‖𝐴𝑣‖_0,Ω+𝑞_1,Ω≤𝐶𝑔_0,Ω. (2.15)

Under the assumption of (A1), if 𝜕Ω is of 𝐶² or Ω is a two-dimensional convex polygon, it has been shown that (see [20]) ‖𝑣‖_0,Ω≤𝛾₀|𝑣|_1,Ω,∀𝑣∈𝑋,|𝑣|_1,Ω≤𝛾₀‖𝐴𝑣‖_0,Ω,∀𝑣∈𝐷(𝐴), (2.16) where 𝛾₀ is a positive constant only depending on Ω.

The following existence and uniqueness results about problem (2.5) are classical (see [1, 20]).

Theorem 2.1. Assume that 𝜈 and 𝑓∈𝑌 satisfy the following uniqueness condition: 𝑐1−⁰𝛾²₀𝜈²𝑓_0,Ω>0. (2.17) Then problem (2.5) admits a unique solution (𝑢,𝑝)∈(𝐷(𝐴),𝐻¹(Ω)∩𝑀) with div𝑢=0 such that |𝑢|_1,Ω≤𝛾₀𝜈𝑓_0,Ω,‖𝐴𝑢‖_0,Ω+𝑝_1,Ω≤𝐶𝑓_0,Ω. (2.18)

3. Enriched Nonconforming Finite Element Method

Let 𝒯_ℎ be a regular triangulation of Ω into element {𝐾_𝑗}∶Ω=∪𝐾_𝑗, that is, |𝐾_𝑗|≃𝐶ℎ²_𝐾𝑗, where |𝐾_𝑗| is the area of the element 𝐾_𝑗 and ℎ_𝐾𝑗 is the diameter of 𝐾_𝑗; the mesh parameter ℎ is given by ℎ=max{ℎ_𝐾𝑗∶𝐾_𝑗∈𝒯_ℎ}. Denote the boundary segment and the interior boundary by 𝛾_𝑗=𝜕Ω∩𝜕𝐾_𝑗 and 𝛾_𝑗𝑘=𝛾_𝑘𝑗=𝜕𝐾_𝑗∩𝜕𝐾_𝑘, respectively. Let Γ_ℎ and Γ_𝑗 be the sets of 𝛾_𝑗𝑘 and 𝛾_𝑗. The centers of 𝛾_𝑗 and 𝛾_𝑗𝑘 are indicated by 𝜉_𝑗 and 𝜉_𝑗𝑘, respectively. The finite element spaces investigated in this paper are the following mixed finite element spaces: NCP₁=𝑣∈𝑌∶𝑣|_𝐾𝑗∈𝑃₁𝐾_𝑗²,𝑣𝜉_𝑗𝑘=𝑣𝜉_𝑘𝑗,𝑣𝜉_𝑗=0∀𝑗,𝑘,𝐾_𝑗∈𝒯_ℎ,𝑃₁=𝑞∈𝐻¹(Ω)∶𝑞_𝐾𝑗∈𝑃₁𝐾_𝑗,∀𝐾_𝑗∈𝒯_ℎ, (3.1) where 𝑃₁(𝐾_𝑗) is the set of line polynomials on 𝐾_𝑗, and noting that the nonconforming finite element space NCP₁ is not a subspace of 𝑋. Defining the energy norm ‖𝑣‖_1,ℎ=_𝐾𝑗|𝑣|²_1,𝐾𝑗^1/2,∀𝑣∈NCP₁. (3.2) The finite element spaces NCP₁ and 𝑃₁ satisfy the following approximation property (see [4, 21]): for (𝑣,𝑞)∈𝐻²(Ω)×𝐻¹(Ω), there are two approximations 𝑣_𝐼∈NCP₁ and 𝑞_𝐼∈𝑃₁ such that 𝑣−𝑣_𝐼0,Ω+ℎ𝑣−𝑣_𝐼1,ℎ+𝑞−𝑞_𝐼0,Ω≤𝐶ℎ²‖𝐴𝑣‖_0,Ω+𝑞_1,Ω, (3.3) and the compatibility conditions hold for all 𝑗 and 𝑘: _𝛾𝑗𝑘[𝑣]𝑑𝑠=0,_Γ𝑗𝑣𝑑𝑠=0∀𝑣∈NCP₁, (3.4) where [𝑣]=𝑣_𝛾𝑗𝑘−𝑣_𝛾𝑘𝑗 denotes the jump of the function 𝑣 across the boundary 𝛾_𝑗𝑘.

Set ⟨⋅,⋅⟩_𝑗=(⋅,⋅)_𝜕𝐾𝑗 and |⋅|_𝑚,𝑗=|⋅|_{𝑚,𝐾𝑗}. Then for all 𝑢, 𝑣∈𝐻¹(𝐾_𝑗)², 𝑞∈𝐿²(Ω), the discrete bilinear forms are 𝑎_ℎ(𝑢,𝑣)=_𝐾𝑗𝜈(∇𝑢,∇𝑣)_𝐾𝑗,𝑑_ℎ𝑣,𝑝=_𝐾𝑗div𝑣,𝑝_𝐾𝑗. (3.5) For the nonconforming space NCP₁, we define a local operator Π_𝑗∶𝐻¹𝐾_𝑗²→NCP₁𝐾_𝑗, (3.6) satisfying _𝜕𝐾𝑗𝑣−Π_𝑗𝑣𝑑𝑠=0. (3.7) Then the local operator Π_𝑗 satisfies (see [21]) 𝑣−Π_𝑗𝑣_1,𝐾𝑗≤𝐶ℎ^𝑖|𝑣|_{𝑖+1,𝐾𝑗},𝑣∈𝐻^𝑖+1𝐾_𝑗,𝑖=0,1,Π_𝑗𝑣_1,𝐾𝑗≤𝐶‖𝑣‖_1,𝐾𝑗. (3.8) The global operator Π_ℎ∶𝑋→NCP₁ is defined as Π_ℎ𝑣|_𝑗=Π_𝑗𝑣,𝑣∈𝑋.

As noted, the choice NCP₁-𝑃₁ is an unstable pair that does not satisfy the discrete Inf-Sup condition. Therefore, we need to introduce the enrichment multiscale method to overcome this restriction.

Let 𝐸_ℎ be a finite dimensional space, called multiscale space, such that 𝐸_ℎ∈𝐻¹𝒯_ℎ²,𝐸_ℎ∩NCP₁={0},where𝐻¹𝒯_ℎ²=𝑣∈𝑌∶𝑣|_𝐾𝑗∈𝐻¹𝐾_𝑗². (3.9) The discrete weak formulation of the Stokes equations is to find 𝑢_ℎ+𝑢_𝑒∈NCP₁⊕𝐸_ℎ and 𝑝_ℎ∈𝑃₁, such that 𝑎_ℎ𝑢_ℎ+𝑢_𝑒,𝑣_ℎ−𝑑_ℎ𝑣_ℎ,𝑝_ℎ+𝑑_ℎ𝑢_ℎ+𝑢_𝑒,𝑞_ℎ=𝑓,𝑣_Ω, (3.10) for all 𝑣_ℎ∈NCP₁⊕𝐻¹₀(𝒯_ℎ)² and 𝑞_ℎ∈𝑃₁. Let 𝑢_𝑒|_𝐾𝑗=𝑢^𝐾𝑗_𝑒+𝑢^𝜕𝐾𝑗_𝑒, we can solve it through the following local problem: −𝜈Δ𝑢^𝐾𝑗_𝑒=𝑓+𝜈Δ𝑢_ℎ−∇𝑝_ℎin𝐾_𝑗,𝑢^𝐾𝑗_{𝑒𝜕𝐾𝑗}=0,−𝜈Δ𝑢^𝜕𝐾𝑗_𝑒=0in𝐾_𝑗,𝑢^𝜕𝐾𝑗_𝑒=𝑔_𝑒on𝜕𝐾_𝑗,−𝜈𝜕_𝑠𝑠𝑔_𝑒=1ℎ_𝑒𝜈𝜕_𝑛𝑢_ℎ+𝑝_ℎ𝐼⋅𝑛_𝐸,𝑔_𝑒=0atthenodes, (3.11) where ℎ_𝑒 denotes the length of the edge 𝐸∈𝜕𝐾_𝑗; 𝑛 the normal outward vector on 𝜕𝐾_𝑗;𝜕_𝑠,𝜕_𝑛 are the tangential and normal derivative operators, respectively; 𝐼 is the ℝ^2×2 identity matrix. Equation (3.11) is well posed, that is, 𝑢_𝑒 can be expressed by 𝑢_ℎ, 𝑝_ℎ, and 𝑓 on each element 𝐾_𝑗. For convenience, we define two local operators ℳ_𝐾𝑗∶𝐿²(𝐾_𝑗)²→𝐻¹₀(𝐾_𝑗)² and ℋ_𝐾𝑗∶𝐿²(𝜕𝐾_𝑗)²→𝐻¹(𝐾_𝑗)² by 𝑢^𝐾𝑗_𝑒=1𝜈ℳ_𝐾𝑗𝑓+𝜈Δ𝑢_ℎ−∇𝑝_ℎ,∀𝐾_𝑗∈𝒯_ℎ,𝑢_𝑒^𝜕𝐾=1𝜈ℋ_𝐾𝑗𝜈𝜕_𝑛𝑢_ℎ+𝑝_ℎ𝐼⋅𝑛_𝐸,∀𝐾_𝑗∈𝒯_ℎ,𝐸∈Γ_ℎ. (3.12) With Green formulation and (3.12), for all (𝑢_ℎ,𝑝_ℎ), (𝑣_ℎ,𝑞_ℎ)∈NCP₁×𝑃₁, (3.10) can be rewritten as _𝐾𝑗𝜈∇𝑢_ℎ,∇𝑣_ℎ𝐾𝑗−𝑝_ℎ,∇⋅𝑣_ℎ𝐾𝑗+𝑞_ℎ,∇⋅𝑢_ℎ𝐾𝑗+_𝐾𝑗1𝜈ℳ_𝐾𝑗−𝜈Δ𝑢_ℎ+∇𝑝_ℎ−ℋ_𝐾𝑗𝜈𝜕_𝑛𝑢_ℎ+𝑝_ℎ𝐼⋅𝑛_𝐸,𝜈Δ𝑣_ℎ+∇𝑞_ℎ𝐾𝑗+_𝐸∈Γℎ1𝜈ℋ_𝐾𝑗𝜈𝜕_𝑛𝑢_ℎ+𝑝_ℎ𝐼⋅𝑛_𝐸,𝜈𝜕_𝑛𝑣_ℎ+𝑞_ℎ𝐼⋅𝑛_𝐸=_𝐾𝑗𝑓,𝑣_ℎ𝐾𝑗+1𝜈ℳ_𝐾𝑗𝑓,𝜈Δ𝑣_ℎ+∇𝑞_ℎ𝐾𝑗. (3.13) With the help of (3.13), the enriched nonconforming finite element method for the stationary Navier-Stokes equations (2.1) is rewritten as follows: find (𝑢_ℎ,𝑝_ℎ)∈NCP₁×𝑃₁ such that 𝐵𝑢_ℎ,𝑝_ℎ;𝑣_ℎ,𝑞_ℎ+𝑏𝑢_ℎ,𝑢_ℎ,𝑣_ℎ=𝐹𝑣_ℎ,𝑞_ℎ (3.14) for all (𝑣_ℎ,𝑞_ℎ)∈NCP₁×𝑃₁, where 𝐵𝑢_ℎ,𝑝_ℎ;𝑣_ℎ,𝑞_ℎ=𝐵_ℎ𝑢_ℎ,𝑝_ℎ;𝑣_ℎ,𝑞_ℎ+_𝐾𝑗1𝜈ℳ_𝐾𝑗∇𝑝_ℎ,∇𝑞_ℎ𝐾𝑗+_𝐸∈Γℎ1𝜈ℋ_𝐾𝑗𝜈𝜕_𝑛𝑢_ℎ,𝜈𝜕_𝑛𝑣_ℎ𝐸≜𝐵_ℎ𝑢_ℎ,𝑝_ℎ;𝑣_ℎ,𝑞_ℎΩ+_𝐾𝑗𝜏_𝐾𝑗∇𝑝_ℎ,∇𝑞_ℎ𝐾𝑗+_𝐸∈Γℎ𝜏_𝐸𝜈𝜕_𝑛𝑢_ℎ,𝜈𝜕_𝑛𝑣_ℎ𝐸,𝐹𝑣_ℎ,𝑞_ℎ=_𝐾𝑗𝑓,𝑣_ℎ𝐾𝑗+1𝜈ℳ_𝐾𝑗𝑓,∇𝑞_ℎ𝐾𝑗,𝐵_ℎ𝑢_ℎ,𝑝_ℎ;𝑣_ℎ,𝑞_ℎ=𝑎_ℎ𝑢_ℎ,𝑣_ℎ−𝑑_ℎ𝑣_ℎ,𝑝_ℎ+𝑑_ℎ𝑢_ℎ,𝑞_ℎ. (3.15)

By applying the technique to one used in [16], we can obtain that (𝑏_𝐾𝑗,1)_𝐾𝑗/|𝐾_𝑗|⋍𝐶ℎ²_𝐾𝑗, (𝑎_𝐾𝑗,1)_𝐸/ℎ_𝑒⋍ℎ_𝑒/12, 𝜏_𝐾𝑗⋍𝐶ℎ²_𝐾𝑗 and 𝜏_𝐸⋍ℎ_𝑒/(12𝜈). Moreover, if 𝑓 is a piecewise constant, then we have ℳ_𝐾𝑗(𝑓)=𝑏_𝐾𝑗𝑓, ℳ_𝐾𝑗𝑓,∇𝑞_ℎ𝐾𝑗=𝑏_𝐾𝑗,1_𝐾𝑗𝐾_𝑗𝑓,∇𝑞_ℎ𝐾𝑗≃𝐶ℎ²_𝐾𝑗𝑓,∇𝑞_ℎ𝐾𝑗. (3.16) Define the mesh-dependent norms as follows: ‖|𝑢|‖²_ℎ=𝜈‖𝑢‖²_1,ℎ+_𝐸∈Γℎ𝜏_𝐸𝜈𝜕_𝑛𝑢²_0,𝐸,𝑞²_ℎ=_𝐾𝑗𝜏_𝐾𝑗𝑞²_1,𝐾𝑗. (3.17)

Remark 3.1. The assumption of piecewise constant 𝑓 is made simply to analyze the problem (3.14), but this assumption does not affect the precision of this method, and (3.14) may be implemented as it is presented for a general function 𝑓∈𝐿²(Ω)². Here, we do not give the detail proof about this fact; readers can visit Appendix  B of the paper [16] for 𝑓∈𝐻¹(Ω)².

Remark 3.2. Generally speaking, the following linear algebra equations can be obtained from the discrete system of original problem: 𝐷𝐴−𝐷_𝑇0𝑈𝑃=𝐹0, (3.18) where the matrices 𝐴 and 𝐷 are deduced from the diffusion, convection, and incompressible terms; 𝐹 is the variation of the source term. The norm of matrix 𝐴 gets smaller as the convection increases; therefore, some unnecessary oscillations will be created. In order to eliminate these oscillations, we introduce the stabilized term, in this case, the coefficient matrix of discrete formulation transforms into 𝐷𝐴−𝐷_𝑇𝐺, (3.19) where 𝐺 is derived from the stabilized term, that is, the term of (∇𝑝_ℎ,∇𝑞_ℎ). As the considered problem has strong convection, in order to obtain a good behavior of matrix 𝐷𝐴−𝐷_𝑇𝐺, we should choose a proper 𝐺. In this way, the singularly perturbed problem can be treated effectively. The reason that we treat the convection term not use enriched function technique is to simply the theoretical analysis and computation, and the discrete convection term has no influence about the stabilized term 𝐺(⋅,⋅).

Lemma 3.3. Let (𝑣_ℎ,𝑞_ℎ)∈NCP₁×𝑃₁, then, 𝐵𝑣_ℎ,𝑞_ℎ;𝑣_ℎ,𝑞_ℎ=𝑣_ℎ²_ℎ+𝑞_ℎ²_ℎ. (3.20)

Proof. The results follow from the definition of (3.15) and the mesh-dependent norms in each 𝐾∈𝒯_ℎ.
Before establishing the stability of scheme (3.14), we introduce the local trace theorem (see [1]). There exists 𝐶>0, independent of ℎ, such that ‖𝑢‖²_{0,𝜕𝐾𝑗}≤𝐶ℎ_𝐾⁻¹_𝑗‖𝑢‖²_0,𝐾𝑗+ℎ_𝐾𝑗|𝑢|²_1,𝐾𝑗,∀𝑢∈𝐻¹𝐾_𝑗. (3.21)

Theorem 3.4. There exist two positive constants 𝐶, 𝛽 depending on 𝜈, for all (𝑢_ℎ,𝑝_ℎ), (𝑣_ℎ,𝑞_ℎ)∈NCP₁×𝑃₁ such that 𝐵𝑢_ℎ,𝑝_ℎ;𝑣_ℎ,𝑞_ℎ≤𝐶𝑢_ℎ1,ℎ+𝑝_ℎ0,Ω𝑣_ℎ1,ℎ+𝑞_ℎ0,Ω,sup_{0≠(𝑣ℎ,𝑞ℎ)∈NCP1,𝑃1}𝐵𝑢_ℎ,𝑝_ℎ;𝑣_ℎ,𝑞_ℎ𝑣_ℎ1,ℎ+𝑞_ℎ0,Ω≥𝛽𝑢_ℎ1,ℎ+𝑝_ℎ0,Ω. (3.22)(3.23)

Proof. It follows from (𝑢_ℎ,𝑝_ℎ), (𝑣_ℎ,𝑞_ℎ)∈NCP₁×𝑃₁, inverse inequality, (3.15), and (3.21) that 𝐵𝑢_ℎ,𝑝_ℎ;𝑣_ℎ,𝑞_ℎ≤𝜈𝑢_ℎ1,ℎ𝑣_ℎ1,ℎ+𝑢_ℎ1,ℎ𝑞_ℎ0,Ω+𝑝_ℎ0,Ω𝑣_ℎ1,ℎ+_𝐾𝑗𝜏_𝐾𝑗𝑝_{ℎ1,𝐾𝑗}𝑞_{ℎ1,𝐾𝑗}+_𝐸∈Γℎ𝜏_𝐸𝜈𝜕_𝑛𝑢_ℎ0,𝐸𝜈𝜕_𝑛𝑣_ℎ0,𝐸≤𝜈𝑢_ℎ1,ℎ𝑣_ℎ1,ℎ+𝑢_ℎ1,ℎ𝑞_ℎ0,Ω+𝑝_ℎ0,Ω𝑣_ℎ1,ℎ+𝐶₁𝑝_ℎ0,Ω𝑞_ℎ0,Ω+𝐶_{2𝐸∈Γℎ}𝜏_𝐸ℎ_𝐾^−1/2_𝑗∇𝑢_{ℎ0,𝐾𝑗}+ℎ_𝐾^1/2_𝑗∇𝑢_{ℎ1,𝐾𝑗}ℎ_𝐾^−1/2_𝑗∇𝑣_{ℎ0,𝐾𝑗}+ℎ_𝐾^1/2_𝑗∇𝑣_{ℎ1,𝐾𝑗}≤𝐶(𝜈)𝑢_ℎ1,ℎ𝑣_ℎ1,ℎ+𝑢_ℎ1,ℎ𝑞_ℎ0,Ω+𝑝_ℎ0,Ω𝑣_ℎ1,ℎ+𝐶₁𝑝_ℎ0,Ω𝑞_ℎ0,Ω≤𝐶𝑢_ℎ1,ℎ+𝑝_ℎ0,Ω𝑣_ℎ1,ℎ+𝑞_ℎ0,Ω, (3.24) that is, the continuity result (3.22) holds.
From the properties of the nonconforming finite element given in [18], for all 𝑝_ℎ∈𝐿²(𝐾_𝑗), there exists a function 𝑤∈𝐻¹(𝐾_𝑗)², such that ‖𝑤‖_1,ℎ=‖𝑝_ℎ‖₀ and ∇⋅𝑤,𝑝_ℎ𝐾𝑗=𝑝_ℎ²_0,𝐾𝑗,‖𝑤‖_1,ℎ≤𝐶𝑝_ℎ0,Ω. (3.25) Using the Cauchy-Schwartz inequality and (3.25), we have 𝐵𝑢_ℎ,𝑝_ℎ;(−𝑤,0)=−𝜈∇𝑢_ℎ,∇𝑤+𝑝_ℎ,∇⋅𝑤−_𝐸∈Γℎ𝜏_𝐸𝜈𝜕_𝑛𝑢_ℎ,𝜈𝜕_𝑛𝑤_𝐸≥−𝜈𝑢_ℎ1,ℎ‖𝑤‖_1,ℎ+𝐶₀‖𝑤‖_1,ℎ𝑝_ℎ0,Ω−_𝐸∈Γℎ𝜏_𝐸[𝜈𝜕_𝑛𝑢_ℎ]_0,𝐸[𝜈𝜕_𝑛𝑤]_0,𝐸≥−𝜈𝑢_ℎ²_1,ℎ+_𝐸∈Γℎ𝜏_𝐸𝜈𝜕_𝑛𝑢_ℎ²_0,𝐸^1/2𝜈‖𝑤‖²_1,ℎ+_𝐸∈Γℎ𝜏_𝐸𝜈𝜕_𝑛𝑤²_0,𝐸^1/2+𝐶₀‖𝑤‖_1,ℎ𝑝_ℎ0,ℎ. (3.26) Using (3.21) and inverse inequality, we obtain that 𝜏_𝐸𝜈𝜕_𝑛𝑤²_0,𝐸≤ℎ_𝑒ℎ12𝜈^𝐾−1‖𝜈∇𝑤⋅𝑛‖²_0,𝐾+ℎ_𝐾|𝜈∇𝑤⋅𝑛|²_1,𝐾≤ℎ_𝑒𝜈12ℎ_𝐾|𝑤|²_1,𝐾+𝐶_𝐾𝜈ℎ_𝑒12ℎ_𝐾|𝑤|²_1,𝐾≤𝜈1+𝐶_𝐾|𝑤|12^21,𝐾. (3.27) Combining (3.26) with (3.27) yields 𝐵𝑢_ℎ,𝑝_ℎ;()√−𝑤,0≥−𝐶𝜈‖𝑤‖_1,ℎ𝜈𝑢_ℎ²_1,ℎ+_𝐸∈Γℎ𝜏_𝐸𝜈𝜕_𝑛𝑢_ℎ²_0,𝐸^1/2+𝐶₀‖𝑤‖_1,ℎ𝑝_ℎ0,Ω=−√𝐶𝜈‖𝑤‖_1,ℎ𝑢_ℎℎ+𝐶₀𝑝_ℎ²_0,Ω≥−𝐶𝜈𝛾₁⁻¹𝑢_ℎ²_ℎ+𝐶₀−𝛾₁𝑝_ℎ²_0,Ω, (3.28) where 𝐶=(1+𝐶₀)/12 with 𝐶₀=max_{𝑘∈𝒯ℎ}𝐶_𝐾, and 𝛾₁ is chosen small enough. Let 𝑣_ℎ,𝑞_ℎ=𝑢_ℎ−𝛿𝑤,𝑝_ℎ,𝛿>0. (3.29) Using (3.26) and Lemma 3.3 we have 𝐵𝑢_ℎ,𝑝_ℎ;𝑣_ℎ,𝑞_ℎ=𝐵𝑢_ℎ,𝑝_ℎ;𝑢_ℎ,𝑝_ℎ+𝛿𝐵𝑢_ℎ,𝑝_ℎ;()≥𝑢−𝑤,0_ℎ2ℎ+𝑝_ℎ²_ℎ+𝛿−𝐶𝜈𝛾₁⁻¹𝑢_ℎ²_ℎ+𝐶₀−𝛾₁𝑝_ℎ²_0,Ω≥1−𝐶𝛿𝜈𝛾₁⁻¹𝑢_ℎ²_ℎ+𝑝_ℎ²_ℎ+𝛿𝐶₀−𝛾₁𝑝_ℎ²_0,Ω≥𝜈1−𝐶𝛿𝜈𝛾₁⁻¹𝑢_ℎ²_1,ℎ+𝛿𝐶₀−𝛾₁𝑝_ℎ²_0,Ω, (3.30) provided that 0<𝛿<𝛾₁/(𝐶𝜈) and 0<𝛾₁<𝐶₀. Denote 𝐶(𝜈)≜min𝜈1−𝐶𝛿𝜈𝛾₁⁻¹,𝛿𝐶₀−𝛾₁,𝐶(𝛿)≜max2,1+2𝛿². (3.31) Then we have 𝑣_ℎ²_1,ℎ+𝑞_ℎ²_0,Ω=𝑢_ℎ−𝛿𝑤²_1,ℎ+𝑝_ℎ²_0,Ω≤2𝑢_ℎ²_1,ℎ+1+2𝛿²𝑝_ℎ²_0,Ω≤𝐶(𝛿)𝑢_ℎ²_1,ℎ+𝑝_ℎ²_0,Ω. (3.32) Taking 𝛽=𝐶(𝜈)/(𝐶(𝛿)), we obtain the desired result (3.23).

Theorem 3.5. Under the assumptions of Theorem 2.1 and the following condition: 𝑡ℎ𝑒𝑠𝑡𝑟𝑜𝑛𝑔𝑢𝑛𝑖𝑞𝑢𝑒𝑛𝑒𝑠𝑠𝑐𝑜𝑛𝑑𝑖𝑡𝑖on:1−𝑐₀𝛾₀𝛾₀+𝜈^1/2𝜈²𝑓_0,Ω>0. (3.33) Problem (3.14) admits a unique solution (𝑢_ℎ,𝑝_ℎ)∈(NCP₁,𝑃₁), and satisfying 𝑢_ℎ1,ℎ≤𝛾₀+𝜈^1/2𝜈𝑓_0,Ω,𝑝_ℎ0,Ω≤𝛽⁻¹𝛾₀+𝐶𝜈𝑓_0,Ω+𝑐₀𝛾₀𝜈⁻²𝛾₀+𝜈^1/22𝑓²_0,Ω. (3.34)

Proof. Let Hilbert space 𝐻_ℎ=(NCP₁,𝑃₁) be with the scalar product and norm 𝑣,𝑞;(𝑤,𝑟)_𝐻ℎ=_𝐾𝑗(∇𝑣,∇𝑤)_𝐾𝑗+𝑞,𝑟, (3.35) and 𝐾_ℎ be a nonvoid, convex, and compact subset of 𝐻_ℎ defined by 𝐾_ℎ=𝑣,𝑞∈𝐻_ℎ∶‖𝑣‖_1,ℎ≤𝛾₀+𝜈^1/2𝜈𝑓₀,𝑞_0,Ω≤𝜈𝛾₀+𝐶𝑓𝛽𝜈^0,Ω+𝑐₀𝛾₀𝛾₀+𝜈^1/22𝛽𝜈²𝑓²_0,Ω. (3.36)
Defining a continuous mapping from 𝐾_ℎ into 𝐻_ℎ as follows: given (𝑣,𝑞)∈𝐾_ℎ,for all (𝑤,𝑟)∈𝐻_ℎ, find (𝑣,𝑞)=𝐹(𝑣,𝑞) such that 𝐵𝑣,𝑞;(𝑤,𝑟)+𝑏(𝑣,𝑣,𝑤)=𝑓,𝑤+_𝐾𝑗𝜏_𝐾𝑗𝑓,∇𝑟_𝐾𝑗. (3.37) Taking (𝑤,𝑟)=(𝑣,𝑞), using (2.8)–(2.13) and inverse inequality yields 𝜈‖𝑣‖²_1,ℎ+𝑞²_ℎ≤‖|𝑣|‖²_ℎ+𝑞²_ℎ≤𝛾₀𝑓_0,Ω‖𝑣‖_1,ℎ+_𝐾𝑗𝜏_𝐾𝑗𝑓_0,𝐾𝑗∇𝑞_0,𝐾𝑗≤𝛾²₀+ℎ2𝜈2²𝑓²_0,Ω+𝜈2‖𝑣‖²_1,ℎ+12_𝐾𝑗𝜏_𝐾𝑗∇𝑞²_0,𝐾𝑗=𝛾²₀+ℎ2𝜈2²𝑓²_0,Ω+𝜈2‖𝑣‖²_1,ℎ+12𝑞²_ℎ. (3.38) As a consequence, we have ‖𝑣‖_1,ℎ≤𝛾₀+𝜈^1/2𝜈𝑓_0,Ω. (3.39) Using again (2.17), (3.23), (3.37), and inverse inequality, we arrive at 𝛽‖𝑣‖_1,ℎ+𝑞_0,Ω≤𝑓,𝑤+∑_𝐾𝑗𝜏_𝐾𝑗𝑓,∇𝑟_𝐾𝑗‖𝑤‖_1,ℎ+‖𝑟‖₀+𝑐₀𝛾₀‖𝑣‖_1,ℎ‖𝑣‖_1,ℎ≤𝛾₀𝑓_0,Ω+𝐶ℎ𝜈𝑓_0,Ω+𝑐₀𝛾₀𝜈⁻²𝛾₀+𝜈^1/22𝑓²_0,Ω≤𝛾₀𝜈+𝐶ℎ𝜈𝑓_0,Ω+𝑐₀𝛾₀𝜈⁻²𝛾₀+𝜈^1/22𝑓²_0,Ω. (3.40) Hence, the two estimates imply (𝑣,𝑞)=𝐹(𝑣,𝑞)∈𝐾_ℎ, thanks to the fixed point theorem, the mapping (𝑣,𝑞)=𝐹(𝑣,𝑞)∈𝐾_ℎ has at least one fixed point (𝑢_ℎ,𝑝_ℎ)∈𝐾_ℎ; namely, (𝑢_ℎ,𝑝_ℎ)∈𝐾_ℎ is a numerical solution of problem (3.14).
Next, we shall prove that the problem (3.14) has a unique solution (𝑢_ℎ,𝑝_ℎ). In fact, if (𝑣_ℎ,𝑞_ℎ) also satisfies (3.14), then for all (𝑤,𝑟)∈(NCP₁,𝑃₁) we have 𝐵𝑢_ℎ−𝑣_ℎ,𝑝_ℎ−𝑞_ℎ;(𝑤,𝑟)=𝑏𝑣_ℎ−𝑢_ℎ,𝑢_ℎ,𝑤+𝑏𝑣_ℎ,𝑣_ℎ−𝑢_ℎ,𝑤. (3.41) Taking (𝑤,𝑟)=(𝑢_ℎ−𝑣_ℎ,𝑝_ℎ−𝑞_ℎ) in (3.41) and using again (2.8)–(2.13), Lemma 3.3, it follows that 𝜈𝑢_ℎ−𝑣_ℎ²_1,ℎ≤𝑐₀𝛾₀𝑢_ℎ1𝑢_ℎ−𝑣_ℎ²_1,ℎ≤𝑐₀𝛾₀𝛾₀+𝜈^1/2𝜈𝑓_0,Ω𝑢_ℎ−𝑣_ℎ²_1,ℎ, (3.42) Which, together with the strong uniqueness condition 𝜈−𝑐₀𝛾₀𝛾₀+𝜈^1/2𝜈𝑓_0,Ω=𝜈1−𝑐₀𝛾₀𝛾₀+𝜈^1/2𝜈²𝑓_0,Ω>0, (3.43) gives 𝑢_ℎ=𝑣_ℎ. Using again (2.13), (3.23), and (3.41), we obtain 𝛽‖𝑝_ℎ−𝑞_ℎ‖²_0,Ω≤0 which implies 𝑝_ℎ=𝑞_ℎ.

4. Error Estimates

In order to derive the error estimates of the numerical solution (𝑢_ℎ,𝑝_ℎ), we introduce the Galerkin projection (𝑅_ℎ,𝑄_ℎ):(𝑋,𝑀)→(NCP₁,𝑃₁) defined as follows: for all (𝑣_ℎ,𝑞_ℎ)∈(NCP₁,𝑃₁)𝐵𝑅_ℎ𝑣,𝑞,𝑄_ℎ𝑣,𝑞;𝑣_ℎ,𝑞_ℎ=𝐵₀𝑣,𝑞;𝑣_ℎ,𝑞_ℎ. (4.1) Noting the Theorem 3.4, (𝑅_ℎ(𝑣,𝑞),𝑄_ℎ(𝑣,𝑞)) is well defined.

By using a similar argument to the one used in [14, 22], we have the following lemma.

Lemma 4.1. Let (𝑢,𝑝)∈𝐷(𝐴)×(𝐻¹(Ω)∩𝑀); under the assumptions of Theorems 3.4 and 3.5, the projection operator (𝑅_ℎ,𝑄_ℎ) satisfies 𝑢−𝑅_ℎ(𝑢,𝑝)_0,Ω+ℎ𝑢−𝑅_ℎ(𝑢,𝑝)_1,ℎ+𝑝−𝑄_ℎ(𝑢,𝑝)_0,Ω≤𝐶ℎ²‖𝐴𝑢‖_0,Ω+𝑝_1,Ω. (4.2)

Proof. From (𝑢,𝑝)∈[𝐻²(Ω)²∩𝑋]×[𝐻¹(Ω)∩𝑀], we have [𝜈𝜕_𝑛𝑢]_𝐸=0. For all (𝑣_ℎ,𝑞_ℎ)∈NCP₁×𝑃₁, using (4.1) yields 𝐵𝑅_ℎ𝑢,𝑝,𝑄_ℎ𝑢,𝑝𝑣_ℎ,𝑞_ℎ=𝐵₀𝑢,𝑝;𝑣_ℎ,𝑞_ℎ=𝐵𝑢,𝑝;𝑣_ℎ,𝑞_ℎ−_𝐾𝑗𝜏_𝐾𝑗∇𝑝,∇𝑞_ℎ𝐾𝑗. (4.3) From the definition of (𝑅_ℎ(𝑢,𝑝),𝑄_ℎ(𝑢,𝑝)), (3.3), combining Theorem 3.4, (4.3), the triangular with inverse inequalities, we arrive at 𝑢−𝑅_ℎ(𝑢,𝑝)_1,ℎ+𝑝−𝑄_ℎ(𝑢,𝑝)_0,Ω≤𝑢−𝑢_𝐼1,ℎ+𝑝−𝑝_𝐼0,Ω+𝑢_𝐼−𝑅_ℎ(𝑢,𝑝)_1,ℎ+𝑝_𝐼−𝑄_ℎ(𝑢,𝑝)_0,Ω≤𝑢−𝑢_𝐼1,ℎ+𝑝−𝑝_𝐼0,Ω+𝛽⁻¹sup_{(𝑣ℎ,𝑞ℎ)∈NCP1×𝑃1}𝐵𝑢_𝐼−𝑅_ℎ𝑢,𝑝,𝑝_𝐼−𝑄_ℎ𝑢,𝑝;𝑣_ℎ,𝑞_ℎ𝑣_ℎ1,ℎ+𝑞_ℎ0,Ω≤𝛽⁻¹sup_{(𝑣ℎ,𝑞ℎ)∈NCP1×𝑃1}𝐵𝑢_𝐼−𝑢,𝑝_𝐼−𝑝;𝑣_ℎ,𝑞_ℎ+∑_𝐾𝑗𝜏_𝐾𝑗∇𝑝,∇𝑞_ℎ𝐾𝑗𝑣_ℎ1,ℎ+𝑞_ℎ0,Ω+𝑢−𝑢_𝐼1,ℎ+𝑝−𝑝_𝐼0,Ω. (4.4) It is easy to check that 𝐵𝑢_𝐼−𝑢,𝑝_𝐼−𝑝;𝑣_ℎ,𝑞_ℎ≤𝐶𝑢−𝑢_𝐼1,ℎ+𝑝−𝑝_𝐼0,Ω+ℎ𝑝₁𝑣_ℎ1,ℎ+𝑞_ℎ0,Ω. (4.5) Combining (4.4), (4.5), and inverse inequality yields 𝑢−𝑅_ℎ(𝑢,𝑝)_1,ℎ+𝑝−𝑄_ℎ(𝑢,𝑝)_0,Ω≤𝐶𝑢−𝑢_𝐼1,ℎ+𝑝−𝑝_𝐼0,Ω+ℎ𝑝₁+𝛽⁻¹sup_{(𝑣ℎ,𝑞ℎ)∈NCP1×𝑃1}∑_𝐾𝑗𝜏_𝐾𝑗∇𝑝,∇𝑞_ℎ𝐾𝑗𝑣_ℎ1,ℎ+𝑞_ℎ0,Ω≤𝐶₁ℎ‖𝐴𝑢‖_0,Ω+𝑝_1,Ω+𝐶₂ℎsup_{(𝑣ℎ,𝑞ℎ)∈NCP1×𝑃1}∑_𝐾𝑗ℎ_𝐾𝑗∇𝑝_0,𝐾𝑗∇𝑞_{ℎ0,𝐾𝑗}𝑣_ℎ1,ℎ+𝑞_ℎ0,Ω≤𝐶ℎ‖𝐴𝑢‖_0,Ω+𝑝_1,Ω. (4.6) In order to derive the estimate in the 𝐿²-norm, we consider the following dual problem with (𝑒,𝜂)=(𝑢−𝑅_ℎ(𝑢,𝑝),𝑝−𝑄_ℎ(𝑢,𝑝)): Φ|−ΔΦ+∇Ψ=𝑒inΩ,divΦ=0inΩ,_𝜕Ω=0on𝜕Ω. (4.7)(4.8)(4.9) Based on the assumption of (A1), (4.7)–(4.9) have a unique solution and satisfy ‖𝐴Φ‖_0,Ω+|Ψ|_1,Ω≤𝐶𝑢−𝑅_ℎ(𝑢,𝑝)_0,Ω. (4.10) Multiplying (4.7) and (4.8) by 𝑒 and 𝜂, respectively, integrating over Ω, and using (4.3) with (𝑣_ℎ,𝑞_ℎ)=(Φ_𝐼,Ψ_𝐼), we see that 𝑢−𝑅_ℎ(𝑢,𝑝)²₀=𝐵₀𝑢−𝑅_ℎ𝑢,𝑝,𝑝−𝑄_ℎ𝑢,𝑝;(Φ,Ψ)−_𝐾𝑗𝜕Φ𝜕𝑛,𝑢−𝑅_ℎ𝑢,𝑝_𝑗+_𝐾𝑗𝑢−𝑅_ℎ𝑢,𝑝⋅𝑛,Ψ_𝑗=𝐵𝑢−𝑅_ℎ𝑢,𝑝,𝑝−𝑄_ℎ𝑢,𝑝;(Φ,Ψ)−_𝐾𝑗𝜏_𝐾𝑗∇𝑝−𝑄_ℎ𝑢,𝑝,∇Ψ_𝐾𝑗−_𝐾𝑗𝜕Φ𝜕𝑛,𝑢−𝑅_ℎ𝑢,𝑝_𝑗+_𝐾𝑗𝑢−𝑅_ℎ𝑢,𝑝⋅𝑛,Ψ_𝑗=𝐵𝑢−𝑅_ℎ𝑢,𝑝,𝑝−𝑄_ℎ𝑢,𝑝;Φ−Φ_𝐼,Ψ−Ψ_𝐼−_𝐾𝑗𝜏_𝐾𝑗∇𝑝−𝑄_ℎ𝑢,𝑝,∇Ψ_𝐾𝑗+_𝐾𝑗𝜏_𝐾𝑗∇𝑝,∇Ψ_𝐼𝐾𝑗−_𝐾𝑗𝜕Φ𝜕𝑛,𝑢−𝑅_ℎ𝑢,𝑝_𝑗+_𝐾𝑗𝑢−𝑅_ℎ𝑢,𝑝⋅𝑛,Ψ_𝑗, (4.11) where (Φ_𝐼,Ψ_𝐼) is the finite element interpolation of (Φ,Ψ) in (NCP₁,𝑃₁) and satisfies (3.3). For each 𝐸∈𝜕𝐾_𝑗, we define the mean value of 𝑢−𝑅_ℎ(𝑢,𝑝) and Ψ on 𝐸𝑢−𝑅_ℎ1(𝑢,𝑝)=ℎ_𝑒𝐸𝑢−𝑅_ℎ𝑢,𝑝|_𝐾𝑗𝑑𝑠;1Ψ=ℎ_𝑒𝐸Ψ|_𝐾𝑗𝑑𝑠. (4.12) Note that each interior edge appears twice in the sum of (4.11); 𝑢−𝑅_ℎ(𝑢,𝑝) and Ψ are constants. Then it follows from (4.11) that 𝑢−𝑅_ℎ(𝑢,𝑝)²_0,Ω=𝐵𝑢−𝑅_ℎ𝑢,𝑝,𝑝−𝑄_ℎ𝑢,𝑝;Φ−Φ_𝐼,Ψ−Ψ_𝐼+_𝐾𝑗𝜏_𝐾𝑗∇𝑝,∇Ψ_𝐼𝐾𝑗−_𝐾𝑗𝜏_𝐾𝑗∇𝑝−𝑄_ℎ𝑢,𝑝,∇Ψ_𝐾𝑗+_{𝐾𝑗𝐸∈𝜕𝐾𝑗}𝑢−𝑅_ℎ𝑢,𝑝⋅𝑛,Ψ−Ψ_𝐸−_{𝐾𝑗𝐸∈𝜕𝐾𝑗}𝜕Φ𝜕𝑛,𝑢−𝑅_ℎ𝑢,𝑝−𝑢−𝑅_ℎ𝑢,𝑝_𝐸. (4.13) Combining (3.3) with Lemma 4.1, we deduce that 𝐵𝑢−𝑅_ℎ𝑢,𝑝,𝑝−𝑄_ℎ𝑢,𝑝;Φ−Φ_𝐼,Ψ−Ψ_𝐼≤𝐶ℎ²‖𝐴𝑢‖_0,Ω+𝑝_1,Ω‖𝐴Φ‖_0,Ω+|Ψ|_1,Ω,_𝐾𝑗𝜏_𝐾𝑗∇𝑝,∇Ψ_𝐼𝐾𝑗−_𝐾𝑗𝜏_𝐾𝑗∇𝑝−𝑄_ℎ𝑢,𝑝,∇Ψ_𝐾𝑗≤_𝐾𝑗𝜏_𝐾𝑗𝑝_1,𝐾𝑗Ψ_{𝐼1,𝐾𝑗}+_𝐾𝑗𝜏_𝐾𝑗𝑝−𝑄_ℎ𝑢,𝑝_1,𝐾𝑗|Ψ|_1,𝐾𝑗. (4.14)(4.15) With the help of (4.12), we have _𝐸𝑢−𝑅_ℎ𝑢,𝑝−𝑢−𝑅_ℎ𝑢,𝑝𝑑𝑠=0. (4.16) Combining the definition of Π_𝑗, (4.16), and local trace theorem (3.21) with the standard argument for the nonconforming element (see [21]), we see that _{𝐾𝑗𝐸∈𝜕𝐾𝑗}𝜕Φ𝜕𝑛,𝑢−𝑅_ℎ𝑢,𝑝−𝑢−𝑅_ℎ𝑢,𝑝_𝐸=_{𝐾𝑗𝐸∈𝜕𝐾𝑗}𝜕Φ−𝜕Π𝜕𝑛^𝑗Φ𝜕𝑛,𝑢−𝑅_ℎ𝑢,𝑝−𝑢−𝑅_ℎ𝑢,𝑝_𝐸≤_{𝐾𝑗𝐸∈𝜕𝐾𝑗}∇(Φ−Π_𝑗Φ)_𝐿2(𝐸)𝑢−𝑅_ℎ(𝑢,𝑝)−𝑢−𝑅_ℎ(𝑢,𝑝)_𝐿2(𝐸)≤𝐶ℎ²‖𝐴𝑢‖_0,Ω+𝑝_1,Ω‖𝐴Φ‖_0,Ω. (4.17) In a similar way, we have _{𝐾𝑗𝐸∈𝜕𝐾𝑗}𝑢−𝑅_ℎ𝑢,𝑝⋅𝑛,Ψ−Ψ_𝐸≤𝐶ℎ²‖𝐴𝑢‖_0,Ω+𝑝₁|Ψ|_1,Ω. (4.18) By combining (4.13)–(4.15) with (4.17)-(4.18), we deduce that 𝑢−𝑅_ℎ(𝑢,𝑝)_0,Ω≤𝐶ℎ²‖𝐴𝑢‖_0,Ω+𝑝_1,Ω, (4.19) which, together with (4.6). We finish the proof.

Theorem 4.2. Assume that the conditions of Theorems 3.4 and 3.5 are valid; let (𝑢,𝑝), (𝑢_ℎ,𝑝_ℎ) be the solutions of (2.1) and (3.14), respectively, then 𝑢−𝑢_ℎ1,ℎ+𝑝−𝑝_ℎ0,Ω≤𝐶ℎ. (4.20)

Proof. We get the following error equation by combining (2.1) with (3.14), for all (𝑣_ℎ,𝑞_ℎ)∈(NCP₁,𝑃₁)𝐵𝑢−𝑢_ℎ,𝑝−𝑝_ℎ;𝑣_ℎ,𝑞_ℎ+𝑏𝑢−𝑢_ℎ,𝑢,𝑣_ℎ+𝑏𝑢_ℎ,𝑢−𝑢_ℎ,𝑣_ℎ−_𝐾𝑗𝜕𝑢𝜕𝑛,𝑣_ℎ𝑗+_𝐾𝑗𝑣_ℎ⋅𝑛,𝑝_𝑗=_𝐾𝑗𝜏_𝐾𝑗∇𝑝,∇𝑞_ℎ𝐾𝑗−_𝐾𝑗𝜏_𝐾𝑗𝑓,∇𝑞_ℎ𝐾𝑗. (4.21) With (4.3), (4.21) can be rewritten as 𝐵𝑒_ℎ,𝜂_ℎ;𝑣_ℎ,𝑞_ℎ+𝑏𝑢−𝑅_ℎ𝑢,𝑝+𝑒_ℎ,𝑢,𝑣_ℎ+𝑏𝑢_ℎ,𝑢−𝑅_ℎ𝑢,𝑝+𝑒_ℎ,𝑣_ℎ−_𝐾𝑗𝜕𝑢𝜕𝑛,𝑣_ℎ𝑗+_𝐾𝑗𝑣_ℎ⋅𝑛,𝑝_𝑗=−_𝐾𝑗𝜏_𝐾𝑗𝑓,∇𝑞_ℎ𝐾𝑗, (4.22) where 𝑒_ℎ=𝑅_ℎ(𝑢,𝑝)−𝑢_ℎ and 𝜂_ℎ=𝑄_ℎ(𝑢,𝑝)−𝑝_ℎ.
From Theorem 3.5 and (4.22), we get that 𝜂_ℎ0,Ω≤𝛽⁻¹sup_{0≠(𝑣ℎ,𝑞ℎ)∈(NCP1,𝑃1)}𝐵𝑒_ℎ,𝜂_ℎ;𝑣_ℎ,𝑞_ℎ𝑣_ℎ1,ℎ+𝑞_ℎ0,Ω≤𝛽⁻¹sup_{0≠(𝑣ℎ,𝑞ℎ)∈NCP1,𝑃1}1𝑣_ℎ1,ℎ+𝑞_ℎ0,Ω×𝑏𝑢−𝑅_ℎ𝑢,𝑝+𝑒_ℎ,𝑢,𝑣_ℎ+_𝐾𝑗𝜏_𝐾𝑗𝑓,∇𝑞_ℎ𝐾𝑗+𝑏𝑢_ℎ,𝑢−𝑅_ℎ𝑢,𝑝+𝑒_ℎ,𝑣_ℎ+_𝐾𝑗𝜕𝑢𝜕𝑛,𝑣_ℎ𝑗+_𝐾𝑗𝑣_ℎ⋅𝑛,𝑝_𝑗. (4.23) Again, with (2.13), Theorem 2.1, inverse inequality, and Lemma 4.1, we have 𝑏𝑢−𝑅_ℎ𝑢,𝑝,𝑢,𝑣_ℎ+𝑏𝑢_ℎ,𝑢−𝑅_ℎ𝑢,𝑝,𝑣_ℎ≤𝑐₀𝛾₀|𝑢|_1,Ω+𝑢_ℎ1,ℎ𝑢−𝑅_ℎ𝑢,𝑝_1,ℎ𝑣_ℎ1,ℎ≤𝐶ℎ𝑣_ℎ1,ℎ,𝑏𝑒_ℎ,𝑢,𝑣_ℎ+𝑏𝑢_ℎ,𝑒_ℎ,𝑣_ℎ≤𝑐₀𝛾₀|𝑢|_1,Ω+𝑢_ℎ1,ℎ𝑒_ℎ1,ℎ𝑣_ℎ1,ℎ,_𝐾𝑗𝜏_𝐾𝑗𝑓,∇𝑞_ℎ𝐾𝑗≤_𝐾𝑗𝜏_𝐾𝑗𝑓_0,𝐾𝑗∇𝑞_{ℎ0,𝐾𝑗}≤_𝐾𝑗𝐶ℎ²_𝐾𝑗𝑓_0,𝐾𝑗∇𝑞_{ℎ0,𝐾𝑗}≤𝐶(𝜈)ℎ𝑓_0,Ω𝑞_ℎ0,Ω, (4.24)(4.25)(4.26) and using the similar arguments as for (4.17)-(4.18) yields _𝐾𝑗𝜕𝑢𝜕𝑛,𝑣_ℎ𝑗+_𝐾𝑗𝑣_ℎ⋅𝑛,𝑝_𝑗≤𝐶ℎ²‖𝐴𝑢‖_0,Ω+𝑝_1,Ω𝑣_ℎ1,ℎ. (4.27) Combining (4.23)–(4.27) with Theorem 3.5, we arrive at 𝜂_ℎ0,Ω≤𝐶ℎ+2𝑐₀𝛾₀𝛾₀+𝜈^1/2𝜈𝑓_0,Ω𝑒_ℎ1,ℎ. (4.28) Choosing (𝑣_ℎ,𝑞_ℎ)=(𝑒_ℎ,𝜂_ℎ) in (4.22), we obtain that 𝐵𝑒_ℎ,𝜂_ℎ;𝑒_ℎ,𝜂_ℎ+𝑏𝑒_ℎ,𝑢,𝑒_ℎ+_𝐾𝑗𝜏_𝐾𝑗𝑓,∇𝜂_ℎ𝐾𝑗=−𝑏𝑢−𝑅_ℎ𝑢,𝑝,𝑢,𝑒_ℎ−𝑏𝑢_ℎ,𝑢−𝑅_ℎ𝑢,𝑝,𝑒_ℎ+_𝐾𝑗𝜕𝑢𝜕𝑛,𝑣_ℎ𝑗−_𝐾𝑗𝑒_ℎ⋅𝑛,𝑝_𝑗. (4.29) Using (2.13), (2.17), Theorem 2.1, and Lemma 3.3, we get 𝐵𝑒_ℎ,𝜂_ℎ;𝑒_ℎ,𝜂_ℎ−𝑏𝑒_ℎ,𝑢,𝑒_ℎ=𝑒_ℎ²_ℎ+𝜂_ℎ²_ℎ−𝑐₀𝛾₀|𝑢|₁𝑒_ℎ²_1,ℎ≥𝜈𝑒_ℎ²_1,ℎ−𝑐₀𝛾₀|𝑢|₁𝑒_ℎ²_1,ℎ≥𝜈1−𝑐₀𝛾²₀𝜈⁻²𝑓₀𝑒_ℎ²_1,ℎ>0. (4.30) Combining (4.24)–(4.28) with (4.29) yields: 𝑒_ℎ1,ℎ≤𝐶ℎ. (4.31) From (4.28) and (4.31), we obtain that ‖𝜂_ℎ‖_0,Ω≤𝐶ℎ. Furthermore, we finish the proof by combining triangles inequality with Lemma 4.1, (4.28), and (4.31).

Theorem 4.3. Let (𝑢,𝑝) and (𝑢_ℎ,𝑝_ℎ) be the solutions of (2.1) and (3.14), respectively, then we have 𝑢−𝑢_ℎ0,Ω≤𝐶ℎ²‖𝐴𝑢‖_0,Ω+𝑝_1,Ω. (4.32)

Proof. Using the duality argument for a linearized stationary Navier-Stokes problem; for some given 𝑔∈𝑌 and the solution (𝑢,𝑝) of (2.1), defining (Φ,Ψ)∈(𝑋,𝑀) by 𝐵()()𝑢|−𝜈ΔΦ+∇Ψ+𝑢,Φ−𝐵𝑢,Φ=𝑔,inΩ,divΦ=0inΩ,_𝜕Ω=0on𝜕Ω, (4.33)(4.34)(4.35) where 𝐵(𝑢,Φ) is defined as ⟨𝑣,𝐵(𝑢,Φ)⟩_{𝑋×𝑋′}=𝑏(𝑣,𝑢,Φ), for all 𝑣∈𝑋; multiplying (4.33) and (4.34) by 𝑣∈𝑋 and 𝑞∈𝑀, respectively; integrating over Ω, from (2.8)–(2.11), it is easily to see that the bilinear form 𝑎(Φ,𝑣)−𝑑(𝑣,Ψ)+𝑑(Φ,𝑞) is continuity and 𝑋×𝑀 coercive, by using the Lax-Milgram’s Lemma, (4.33)–(4.35) have a unique solution (Φ,Ψ).
Multiplying (4.33) and (4.34) by Φ and Ψ, respectively, using (2.13) and Theorem 2.1, we have 𝜈|Φ|²_1,Ω−𝑏(Φ,𝑢,Φ)≥𝜈|Φ|²_1,Ω−𝑐₀𝛾₀|𝑢|_1,Ω|Φ|²_1,Ω≥𝜈1−𝑐₀𝛾²₀𝜈⁻²𝑓_0,Ω|Φ|²_1,Ω>0. (4.36) On the other hand, estimating the right term yields Φ,𝑔≤‖Φ‖_0,Ω𝑔_0,Ω≤𝛾₀|Φ|_1,Ω𝑔_0,Ω. (4.37) By using (4.36) and (4.37), we arrive at |Φ|_1,Ω≤𝐶𝑔_0,Ω. (4.38) Setting Ψ=0 and taking the scalar product of (4.33) with 𝐴Φ in 𝑌 yields 𝜈‖𝐴Φ‖²₀+𝑏(𝐴Φ,𝑢,Φ)−𝑏(𝑢,Φ,𝐴Φ)=𝑔,𝐴Φ. (4.39) Using the Gagliardo-Nirenberg inequality yields ‖𝑣‖²_𝐿4≤𝐶‖𝑣‖_0,Ω|𝑣|_1,Ω‖‖,∀𝑣∈𝑋,∇𝑣²_𝐿4≤𝐶‖𝐴𝑣‖_0,Ω|𝑣|_1,Ω,∀𝑣∈𝐻²(Ω)²∩𝑋. (4.40) With the help of the Agmon’s inequality, we have ‖𝑣‖²_𝐿^∞≤𝐶‖𝑣‖_0,Ω‖𝐴𝑣‖_0,Ω,∀𝑣∈𝐻²(Ω)²∩𝑋. (4.41) Furthermore, the following estimates are hold: 𝑏(𝐴Φ,𝑢,Φ)≤𝐶‖𝐴Φ‖_0,Ω‖∇𝑢‖_𝐿⁴‖Φ‖_𝐿⁴+‖𝑢‖_𝐿^∞‖∇Φ‖_0,Ω≤𝜈4‖𝐴Φ‖²_0,Ω+𝐶𝜈⁻¹‖𝐴𝑢‖²_0,Ω‖∇Φ‖²_0,Ω𝑏()≤𝜈),𝑢,Φ,𝐴Φ4‖𝐴Φ‖²_0,Ω+𝐶𝜈⁻¹‖𝐴𝑢‖²_0,Ω‖∇Φ‖²_0,Ω≤𝜈),𝑔,𝐴Φ4‖𝐴Φ‖²_0,Ω+𝐶𝜈⁻¹𝑔_0,Ω. (4.42) Combining above inequalities with (4.38) and (4.39), we arrive at ‖𝐴Φ‖²_0,Ω≤𝑐𝜈⁻²1+‖𝐴𝑢‖²_0,Ω𝑔²_0,Ω. (4.43) Applying the continuous Inf-Sup condition (2.8) yields |Ψ|_1,Ω≤𝑐𝜈‖𝐴Φ‖_0,Ω+𝑐‖𝐴𝑢‖_0,Ω|Φ|_1,Ω+𝑐𝑔_0,Ω. (4.44) Combining (4.43)-(4.44) with (2.18), (4.38), we arrive at ‖𝐴Φ‖_0,Ω+|Ψ|_1,Ω≤𝐶𝑔_0,Ω. (4.45) Taking 𝑔=𝑅_ℎ(𝑢,𝑝)−𝑢_ℎ, multiplying (4.33) and (4.34) by 𝑅_ℎ(𝑢,𝑝)−𝑢_ℎ and 𝑄_ℎ(𝑢,𝑝)−𝑝_ℎ, respectively, using (4.1) yields 𝑅_ℎ(𝑢,𝑝)−𝑢_ℎ²_0,Ω=𝐵𝑅_ℎ𝑢,𝑝−𝑢_ℎ,𝑄_ℎ𝑢,𝑝−𝑝_ℎ;(Φ,Ψ)+𝑏𝑢,𝑅_ℎ𝑢,𝑝−𝑢_ℎ,Φ−_𝐾𝑗𝜏_𝐾𝑗∇𝑄_ℎ𝑢,𝑝−𝑝_ℎ,∇Ψ_𝐾𝑗+_𝐾𝑗Ψ,𝑅_ℎ𝑢,𝑝−𝑢_ℎ⋅𝑛,_𝑗+𝑏𝑅_ℎ𝑢,𝑝−𝑢_ℎ,𝑢,Φ−_𝐾𝑗𝜕Φ𝜕𝑛,𝑅_ℎ𝑢,𝑝−𝑢_ℎ𝑗. (4.46) Setting (𝑣_ℎ,𝑞_ℎ)=(Φ_𝐼,Ψ_𝐼) in (4.22), and using (4.46), we obtain that 𝑅_ℎ(𝑢,𝑝)−𝑢_ℎ²_0,Ω=𝐵𝑅_ℎ𝑢,𝑝−𝑢_ℎ,𝑄_ℎ𝑢,𝑝−𝑝_ℎ;Φ−Φ_𝐼,Ψ−Ψ_𝐼−_𝐾𝑗𝜏_𝐾𝑗∇𝑄_ℎ𝑢,𝑝−𝑝_ℎ,∇Ψ_𝐾𝑗+𝑏𝑢−𝑢_ℎ,𝑅_ℎ𝑢,𝑝−𝑢_ℎ,Φ+_𝐾𝑗𝜕𝑢𝜕𝑛,𝑣_ℎ𝑗+𝑏𝑢−𝑅_ℎ𝑢,𝑝,𝑢,Φ_𝐼+𝑏𝑅_ℎ𝑢,𝑝−𝑢_ℎ,𝑢,Φ−Φ_𝐼+𝑏𝑢_ℎ,𝑢−𝑅_ℎ𝑢,𝑝,Φ−Φ_𝐼+𝑏𝑢_ℎ,𝑅_ℎ𝑢,𝑝−𝑢_ℎ,Φ−Φ_𝐼−𝑏𝑢_ℎ,𝑢−𝑅_ℎ𝑢,𝑝,Φ+_𝐾𝑗Ψ,𝑅_ℎ𝑢,𝑝−𝑢_ℎ⋅𝑛,_𝑗−_𝐾𝑗𝜏_𝐾𝑗𝑓,∇Ψ_𝐼𝐾𝑗−_𝐾𝑗𝜕Φ𝜕𝑛,𝑅_ℎ𝑢,𝑝−𝑢_ℎ𝑗−_𝐾𝑗𝑣_ℎ⋅𝑛,𝑝_𝑗. (4.47) We now estimate the right terms of (4.47). Thanks to (3.3); Theorems 2.1, 3.4-3.5, and 4.2; inverse inequality; Lemma 4.1; we know that 𝐵𝑅_ℎ𝑢,𝑝−𝑢_ℎ,𝑄_ℎ𝑢,𝑝−𝑝_ℎ;Φ−Φ_𝐼,Ψ−Ψ_𝐼≤𝐶𝑅_ℎ𝑢,𝑝−𝑢_ℎ1,ℎ+𝑄_ℎ𝑢,𝑝−𝑝_ℎ0,ΩΦ−Φ_𝐼1,ℎ+Ψ−Ψ_𝐼0,Ω+ℎ|Ψ|_1,Ω≤𝐶ℎ²‖𝐴𝑢‖_0,Ω+𝑝_1,Ω‖𝐴Φ‖_0,Ω+|Ψ|_1,Ω,𝑏𝑅_ℎ𝑢,𝑝−𝑢_ℎ,𝑢,Φ−Φ_𝐼+𝑏𝑢_ℎ,𝑅_ℎ𝑢,𝑝−𝑢_ℎ,Φ−Φ_𝐼≤2𝑐₀|𝑢|_1,Ω+𝑢_ℎ1,ℎ𝑅_ℎ(𝑢,𝑝)−𝑢_ℎ1,ℎΦ−Φ_𝐼1,ℎ≤𝐶ℎ²‖𝐴𝑢‖_0,Ω+𝑝_1,Ω‖𝐴Φ‖_0,Ω,_𝐾𝑗𝜏_𝐾𝑗∇𝑄_ℎ𝑢,𝑝−𝑝_ℎ,∇Ψ_𝐾𝑗≤_𝐾𝑗𝐶ℎ²_𝐾∇(𝑄_ℎ(𝑢,𝑝)−𝑝_ℎ)_0,𝐾𝑗‖∇Ψ‖_0,𝐾𝑗≤𝐶ℎ²‖𝐴𝑢‖_0,Ω+𝑝₁|Ψ|_1,Ω,𝑏𝑢−𝑅_ℎ𝑢,𝑝,𝑢,Φ_𝐼≤𝐶𝑢−𝑅_ℎ(𝑢,𝑝)_0,Ω‖𝐴𝑢‖_0,ΩΦ_𝐼1,ℎ≤𝐶ℎ²‖𝐴𝑢‖_0,Ω+𝑝_1,Ω|Φ|_1,Ω,𝑏𝑢−𝑢_ℎ,𝑅_ℎ𝑢,𝑝−𝑢_ℎ,Φ≤𝑐₀𝑢−𝑢_ℎ1,ℎ𝑅_ℎ(𝑢,𝑝)−𝑢_ℎ1,ℎ|Φ|_1,Ω≤𝐶ℎ²‖𝐴𝑢‖_0,Ω+𝑝₁|Φ|_1,Ω,𝑏𝑢_ℎ,𝑢−𝑅_ℎ𝑢,𝑝,Φ≤𝐶𝑢_ℎ1,ℎ𝑢−𝑅_ℎ𝑢,𝑝_0,Ω‖𝐴Φ‖_0,Ω,_𝐾𝑗𝜏_𝐾𝑗𝑓,∇Ψ_𝐼𝐾𝑗≤_𝐾𝑗𝐶ℎ²_𝐾𝑗𝑓_0,𝐾𝑗∇Ψ_{𝐼0,𝐾𝑗}≤𝐶ℎ²𝑓_0,Ω|Ψ|_1,Ω,𝑏𝑢_ℎ,𝑢−𝑅_ℎ𝑢,𝑝,Φ−Φ_𝐼≤𝑐₀𝑢_ℎ1,ℎ𝑢−𝑅_ℎ(𝑢,𝑝)_1,ℎΦ−Φ_𝐼1,ℎ≤𝐶ℎ²‖𝐴𝑢‖_0,Ω+𝑝_1,Ω‖𝐴Φ‖_0,Ω. (4.48) Applying the argument used for (4.18)-(4.21), (4.28) gives _𝐾𝑗𝜕Φ𝜕𝑛,𝑅_ℎ𝑢,𝑝−𝑢_ℎ𝑗+_𝐾𝑗Ψ,𝑅_ℎ𝑢,𝑝−𝑢_ℎ⋅𝑛,_𝑗≤𝐶ℎ²‖𝐴𝑢‖_0,Ω+𝑝_1,Ω‖𝐴Φ‖_0,Ω+|Ψ|_1,Ω,_𝐾𝑗𝜕𝑢𝜕𝑛,Φ_𝐼+_𝐾𝑗Φ_𝐼⋅𝑛,𝑝_𝑗≤𝐶ℎ²‖𝐴𝑢‖_0,Ω+𝑝_1,ΩΦ_𝐼1,ℎ≤𝐶ℎ²‖𝐴𝑢‖_0,Ω+𝑝_1,ΩΦ_𝐼1,Ω. (4.49) Combining above inequalities, Theorems 2.1, 3.5, (4.46), and (4.51), we get that 𝑅_ℎ(𝑢,𝑝)−𝑢_ℎ²_0,Ω≤𝐶ℎ²‖𝐴𝑢‖_0,Ω+𝑝_1,Ω‖𝐴Φ‖_0,Ω+|Ψ|_1,Ω+𝐶𝑢_ℎ1,ℎ𝑢−𝑅_ℎ(𝑢,𝑝)_0,Ω‖𝐴Φ‖_0,Ω≤𝐶ℎ²‖𝐴𝑢‖_0,Ω+𝑝_1,Ω𝑅_ℎ(𝑢,𝑝)−𝑢_ℎ1,Ω𝛾+𝐶⁰+𝜈^1/2𝜈𝑓_0,Ω𝑢−𝑅_ℎ(𝑢,𝑝)²_0,Ω. (4.50) Choosing the appropriate 𝜈,Ω and 𝑓 such that 1−𝐶((𝛾₀+𝜈^1/2)/𝜈)‖𝑓‖₀>0, then we have 𝑅_ℎ(𝑢,𝑝)−𝑢_ℎ0,Ω≤𝐶ℎ²‖𝐴𝑢‖_0,Ω+𝑝_1,Ω. (4.51) By applying the triangles inequality, Lemma 4.1 and (4.51), we finish the proof.

Lemma 4.4. Under the assumptions of Theorem 4.2, the following estimate about 𝑢−𝑢_ℎ in mesh-dependent norm 𝑢−𝑢_ℎℎ≤𝐶ℎ, (4.52) holds, where 𝑢 and 𝑢_ℎ are the solution of problem (2.1) and (3.14), respectively.

Proof. According to the definition of ‖|⋅|‖_ℎ, 𝑢_ℎ∈𝑋_ℎ, with Theorems 2.1 and 4.2, inverse and local trace inequalities (3.21), we have 𝑢−𝑢_ℎ²_ℎ=𝜈𝑢−𝑢_ℎ²_1,ℎ+_𝐸∈Γℎ𝜏_𝐸𝜈𝜕_𝑛(𝑢−𝑢_ℎ)²_0,𝐸=𝜈𝑢−𝑢_ℎ²_1,ℎ+_𝐸∈Γℎℎ_𝑒12𝜈𝜈𝜕_𝑛(𝑢−𝑢_ℎ)²_0,𝐸≤𝐶ℎ²+𝐶_𝐸∈Γℎℎ_𝑒ℎ_𝐾⁻¹_𝑗∇𝑢−𝑢_ℎ⋅𝑛²_0,𝐾𝑗+ℎ_𝐾𝑗∇𝑢−𝑢_ℎ⋅𝑛²_1,𝐾𝑗≤𝐶ℎ². (4.53)

5. Numerical Validations

In this section, we provide two numerical examples to illustrate the theoretical analysis of the method (3.14). In all experiments, we consider the domain Ω to be the square [0,1]×[0,1]. The mesh consists of triangular elements that are obtained by dividing Ω into subsquares of equal size and then drawing the diagonal in each sub-square, see Figure 1. The software Freefem++, developed by Hecht et al. [23], is used in our experiments.

Figure 1 

Uniform triangulation of Ω into triangulars.

5.1. An Analytical Solution: Convergence Validation

For this test, our purpose is to verify the theoretical analysis which has been established in the previous section by setting the viscosity coefficient 𝜈=1 and 𝑓 is given by the exact solution 𝑢₁=10𝑥²(𝑥−1)²𝑦,𝑢𝑦−12𝑦−1₂=−10𝑥(𝑥−1)(2𝑥−1)𝑦²𝑦−1²,𝑝=10(2𝑥−1)2𝑦−1. (5.1)

In our numerical validation, The experimental rates of convergence with respect to the mesh size ℎ are calculated by the formula (log(𝐸_𝑖/𝐸_𝑖+1))/(log(ℎ_𝑖/ℎ_𝑖+1)), where 𝐸_𝑖 and 𝐸_𝑖+1 are the relative errors corresponding to the meshes of sizes ℎ_𝑖 and ℎ_𝑖+1.

In order to show the efficiency of the enriched multiscale method, we compare the numerical results obtained by using different methods, which are shown in Tables 1, 2, 3, and 4. The compared methods include the pressure projection method (PPM) in [6], the local Gauss integration method (LGIM) in [7, 10, 19], and the standard Galerkin method (SGM) with MINI element (see [1]), respectively. From these tables, we can see that the stabilized multiscale method has good precision for pressure, and the precision of velocity, worse than other methods. Table 5 explains the CPU times that needed for solving the steady Navier-Stokes equations in different mesh sizes. From these data, we know that our method takes less time than other methods. Furthermore, from Tables 1–4, we can see that the numerical results reproduce the established theoretical analysis and show an 𝒪(ℎ) order of convergence for |𝑢−𝑢_ℎ|_1,Ω and ‖𝑝−𝑝_ℎ‖_0,Ω, and an 𝒪(ℎ²) convergence for ‖𝑢−𝑢_ℎ‖_0,Ω.

Table 1 

Numerical results for enriched multiscale method with NCP₁-𝑃₁ element.

Table 2 

Numerical results for PPM in [6] with NCP₁-𝑃₁ element.

Table 3 

Numerical results for SGM with 𝑃₁𝑏-𝑃₁ element.

Table 4 

Numerical results for LGIM in [4, 7, 19] with NCP₁-𝑃₁ element.

Table 5 

CPU time for solving the steady Navier-Stokes equations with different mesh sizes.

5.2. Lid-Driven Cavity Problem

In this test, we consider the incompressible lid-driven cavity flow problem defined on the unit square. Setting 𝑓=0 and the boundary condition 𝑢=0 on [{0}×(0,1)]∪[(0,1)×{0}]∪[{1}×(0,1)] and 𝑢=(1,0)^𝑇 on (0,1)×{1}, see Figure 2. The mesh consists of triangular element and the mesh size ℎ=1/50.

Figure 2 

Lid-driven cavity flow.

Figure 3 shows the pressure contours at different Reylond numbers, where the stopping criterion ‖𝑢_ℎ^𝑛+1−𝑢^𝑛_ℎ‖_0,Ω/‖𝑢_ℎ^𝑛+1‖_0,Ω≤10⁻⁶ is employed, where 𝑢_ℎ^𝑛+1 is the approximation of 𝑢_ℎ at the 𝑛+1 Newton iterative. From Figure 3, we can see that the oscillations are absented for the pressure isovalues by the NCP₁-𝑃₁ approximations, and compared with the results given in [24], we can see that our method has the effect to stabilized flow field (see Figure 4). In this sense, we say that the stabilized multiscale nonconforming finite element method is effective for the stationary Navier-Stokes problem.

(a)

(b)

Figure 3 

Pressure contours of the driven cavity by using the stabilized multiscale method with NCP1-𝑃1 at different Reynolds numbers. (a) Re = 1, (b) Re = 100.

(a)

(b)

Figure 4 

The computed velocity profiles through the geometric center at Re = 100. (a) Horizontal velocity, (b) vertical velocity.

6. Conclusion

In this paper we have derived a theoretical analysis of enriched multiscale nonconforming finite element method for the steady incompressible Navier-Stokes equations. The analysis has extended the work in [16] from the linear problem to the nonlinear problem. The discretization uses nonconforming and conforming piecewise linear finite elements for velocity and pressure over triangles elements, respectively. Numerical tests show that this stabilized method is computationally efficient, and it can be performed locally at the element level with minimal additional cost; at the same time, our numerical results obtained are in good agreement with the established theoretical results.

Acknowledgments

This work was supported by the Natural Science Foundation of China (no. 11126117) and the Doctor Fund of Henan Polytechnic University (no. 648689).