1 Introduction

It is well-known that the finite element approximation of the Naiver-Stokes equations governing the flow of an incompressible fluid may face two problems: spurious oscillations due to the advection-dominated flows and violation of the discrete inf-sup condition. When facing either of the two problems,the standard Galerkin finite element method cannot work well. Several possible remedies for the undesirable behaviors are available,such as the streamline upwind Petrov Galarkin (SUPG) method,the least-square method,and the residual-free bubble method^{[1, 2, 3, 4]}. However,all of these methods are developed with the residuals of the momentum equation with added stabilization terms,requiring approximation of the residual of the second-order derivations for Naiver-Stokes equations,in which pressure and velocity derivatives vanish or are poorly approximated. Recently,stabilized mixed methods involving non-residual stabilization such as the subgrid eddy viscosity model have attracted considerable attention. The eddy viscosity models are firstly proposed by Frisch and Orszag^[5],developed by Iliescu and Layton^[6],and introduced a dissipative mechanism by Smagorinsky^[7]. At present,these models have been further improved by various numerical methods^{[8, 9]}. Some recent models are established by introducing the scale separation subgrid terms based on two-hierarchy mesh structure^{[10, 11]}. Bochev et al.^[12] proposed a variational multiscale method (VMM) in which the diffusion acted only on the finest resolved scales. This VMM is very effective in the convection dominated case. In order to circumvent the so-called inf-sup condition,the pressure projection stabilizations^[13] based on the idea of the subgrid and large eddy simulation (LES) become attractive. Up to now,there are numerous related researches. However,the VMM generally adopts conforming finite element approximations and appears rarely in nonconforming case.

Nonconforming finite element methods are attractive for computational fluid dynamics applications since their edge-oriented degrees of freedom can result in cheap local communications when the methods are parallelized on MIMD-machines^{[14, 15]}. Zhu et al.^[16] proposed a new nonconforming finite element method for stationary Navier-Stokes equations using the P^NC₁ /P1 triangular finite element spaces. Park and Sheen^[17] presented the P1-nonconforming quadrilateral element (denoted by P^NQ₁ ) having only three degrees of freedom to solve the second-order elliptic problem. Recently,Feng et al.^[18] presented the finite element method for Stokes equations adopting the finite elements P^NQ₁ /Q₀ and P^NQ₁ /P^NQ₁ ,respectively. In this paper,we aim to propose a semi-discrete scheme and a full-discrete scheme for a transient Naiver-Stokes problem by using two pairs nonconforming inf-sup unstable finite element spaces,i.e.,P^NC₁ /P^NC₁ triangular and P^NQ₁ /P^NQ₁ quadrilateral finite element spaces. It is attractive due to the simplicity in solving solid and fluid mechanics problems. The full-discrete one has second-order estimations of time,and it is unconditionally stable. Moreover,it has some important features compared with traditional stabilized mixed finite element methods: simple,efficient,and independent of the Reynolds number. Numerical investigations are given to prove the theoretical results.

The rest of this paper is organized as follows. In the next section,we state the abstract functional setting and semi-discrete scheme for Naiver-Stokes equations. In Section 3,we derive the optimal error estimate for the scheme. In Section 4,we analyze the full-discrete one. In Section 5,we give some numerical results to support the theoretical results. 2 Notation and semi-discrete scheme

The flow of an incompressible fluid is governed by the incompressible transient Naiver-Stokes equations

Here,Ω ⊂ R² is a bounded domain with the Lipschitz-continuous boundary Γ,[0,T] is the finite time interval,u(t,x) is the fluid velocity,and p(t,x) is the fluid pressure. ν > 0 is the viscosity,which is inversely proportional to the Reynolds number Re = O(ν⁻¹),the body force f(t,x),and the initial field u₀(x).

Standard notations for Lebesgue and Sobolev spaces are used throughout this paper. We denote by (·,·) the inner product on L²(Ω ) or (L²(Ω ))². If Y denotes a functional space with the norm || · ||Y ,for v = v(x,t) and r > 0,

To introduce the variational formulation,set The continuous bilinear forms a(·,·) and c(·,·) are,respectively,defined by and the trilinear term b(·,·,·) is defined by which is the skew-symmetric form of the convective term.

The mixed variational form of (1)-(2) is to seek (u,p) ∈ X × Q (t > 0) such that for all (v,q) ∈ X × Q,

Let now T_h be a regular Ω into triangular or quadrilateral elements such that h denotes the maximum diameter of the elements in T_h. We define the decomposition T_h in the same way such that H denotes the maximum diameter of the elements in T_h and h 6 H. We denote the boundary of the cell K_i ∈ T_h on ∂Ω by Γ_i = ∂Ω ∩ ∂K_i,the interface between elements K_i ∈ T_h and K_j ∈ T_h by Γ_ij = Γ_ij = ∂K_i ∩ ∂K_j,and the centers of Γ_i and Γ_ij by M_i and M_ij ,respectively.

To describe our finite element space X_h,we need some further notations and definitions. In what follows,εh denotes the set of all edges of T_h. We choose for any face E ∈ εh a unit normal nE with an arbitrary but fixed orientation where nE on the boundary faces is the outer unit normal of Ω . For each element K ∈ T_h,we denote by n_K the outwardpointing unit normal vector on ∂K. For a scalar piecewise continuous function or a vector function ψ,the jump [ψ]E and the average {ψ}E on a face E ∈ T_h are defined by

The construction of the two pairs nonconforming finite element spaces for the velocity and the pressure is described as follows.

(i) When T_h is triangulation,

(ii) When T_h is quadrilateral,we introduce the P1-nonconforming quadrilateral element which was defined in Ref. ^[17] as follows:

In the following paper,we denote X_h × Q_h be the finite element space instead of Xⁱ _h × Qⁱ _h (i = 1,2) for simplicity when there is no confusion.

Now,we introduce the discrete linear terms and the trilinear term as follows:

where the discrete gradient and divergence operators h and h· are defined as follows:

For simplicity,we denote ah(·,·),bh(·,·,·),and ch(·,·) by a(·,·),b(·,·,·),and c(·,·) in the following part. Generic constants which do not depend on the Reynolds number Re and the mesh width h or H are denoted by C,which may take different values at different places.

It is known that the space X_h is not a subspace of X. For any vh in X_h,the following compatibility conditions hold for all i,j^[18]:

It is well-known that the pairs of the finite elements X_h × Q_h do not satisfy the inf-sup condition. Consequently,we introduce the stabilized term G(p_h,q_h) to ensure the stability and convergence of the solutions. Here,G(p_h,q_h) can be defined by

Assume that the operator π_h : Q_h → R_h satisfies

where R_h = {q ∈ L²₀ (Ω) : q|K is constant for all K ∈ T_h}.

We know that for high Reynolds number fluid flows,when the fluid convection dominates fluid flow fields,under the finite-resolution of meshes,the flow becomes very instable. Here,We follow the concept of the variational multi-scale method and choose the following stabilization term.

We assume that LH can be defined as follows:

(i) LH = {LH ∈ (L²(Ω))^2×2 : LH|K ∈ (R_0H(K))^2×2},R0H = {q ∈ L²(Ω ) : q|K ∈ P0, K ∈ T_h};

(ii) LH = {LH ∈ (L²(Ω))^2×2 : LH|K ∈ (P₁(K))^2×2,K ∈ T_h, R E[LH]ds = 0,E ∈ εh}.

Let PH : (L²(Ω))^2×2 → LH be the standard L²-projection with the following properties:

The stabilization term is defined by

where α > 0 is a stabilization parameter. Since PH is an L²-projection,it follows for vh ∈ Xh and ||v_h||₀ > 0 that From 0 6 ||P_hv_h||₀ 6 ||v_h||₀,it follows that

Using the stabilization terms M(u_h,v_h) and G(p_h,q_h),we give the following stabilized semidiscrete scheme for (5)-(7):

Find (u_h,p_h) ∈ X_h × Q_h (t ∈ [0,T]) such that

For the error analysis,we define the following functional space and norms:

Lemma 1 For any q_h ∈ Q_h,there exists a positive constant C₁ independent of h,H,and ν such that

Proof By the subjectivity of the divergence operator,there exists vq ∈ [H¹ ₀ (Ω )]² and a positive constant CF independent of h,H,and ν such that

It is well-known that the pair of spaces X1 h × Rh0 satisfy the uniform inf-sup condition. Then,there exists a positive constant Cπ independent of h,H,ν and a Fortin interpolant π : X → X¹ _h ^[17] such that

Taking the inequality (22) to (21),we have

For the first term in the right-hand side of (23),we have Dividing by ||q_h||₀ in both sides of (23) gives Using |πv_q|_1,h 6 C_FC_π||q_h||₀,there holds

Proof The continuity of A is easy to get by using the older inequality.

In the analysis below,we will often use Young’s inequality: for any real numbers x and y, and δ > 0,

Using Young’s inequality and the stability properties of πvp and vp,we have

Collecting (27)-(29),we have Using the triangular inequality,there holds

Collecting (30) and (31),we deduce that Theorem 1 holds. 3 Error analysis of semi-discrete formulation

We assume that there exists an operator Π₁ : (H¹(Ω ))² → X_h such that

Lemma 2^{[18, 19]} For any s,w ∈ X ∪X_h, For all (v_h,q_h) ∈ X_h × Q_h,the exact solution (u,p) satisfies Subtracting (18) from (34) yields

Theorem 2 Let (u,p) ∈ X × Q be the solution of (1)-(4) and (u_h,p_h) ∈ X_h × Q_h be the solution of (18) and (19). If u_t ∈ L²(0,T ;H¹(Ω )),u ∈ L^∞([0,T] × ) ∩ L²(0,T ;H2(Ω)),and p ∈ L²(0,T ;H1(Ω)),the following inequality holds:

Proof Let vh = Φh in (35). Then,the following equality holds:

We now bound the terms of the right-hand side of (36).

With the Cauchy-Schwarz inequality,Young’s inequality,and the approximation result,the first two terms are bounded as follows:

The first term using Young’s inequality and the inverse inequality is bounded by

The second term is similarly bounded by For the last term,we have Due to Lemma 2,we have With the Cauchy-Schwarz inequality and Young’s inequality,the last term leads to From (36) to (37),we have Integrating over 0 and t and using Gronwall’s lemma yield Then,the theorem can be obtained by the triangle inequality.

Proof From Lemma 1,we can get

Let qh = 0. Subtracting (18) from (34) yields Substituting (39) into (38),we can obtain that all the above terms except the last term in (38) can be handled as those in Theorem 2,and Taking vh = uh and qh = ph in (18) and integrating from 0 to t,we can get Using the triangle inequality for the last term in (38) yields Then,integrating between 0 and T and using the triangle inequality yield Theorem 3. 4 Full discretized formulation

To discretize in time,we divide the interval [0,T] into m equal subintervals,i.e.,Δt = T/m, with grid-points tⁿ = nΔt (0 6 n 6 m). Let uⁿ and uⁿ _h be the values of u and u_h at tⁿ, respectively.

The full discretized formulation reads:

At the first time level,(u¹ _h,p¹ _h) ∈ X_h × Q_h is sought to satisfy

The following identity holds: Firstly,we prove the existence and the uniqueness of the solution of (40). Then,we give the error estimate for the solution.

Lemma 3 (Stability) We have

Proof Let v_h = uⁿ⁺¹ _h in (40) and multiply 4Δt on each side of (40). Then,

Lemma 4 (Uniqueness) When Δt is sufficiently small,(40) has a unique solution.

Proof Due to the theorem of the saddle point of Brouwer and Lemma 1,the problem (40) has at least a solution (uⁿ _h,pⁿ _h). In the following,we give the proof of the uniqueness of the solution. We consider two solutions (u_i,p_i) (i = 1,2). Then,we have

Summing the above result ,we have

where the constant C^* is independent of ν,t,h,and H. The right-hand side of the above inequality holds due to the fact that all the norms are equivalent in finite dimension.

In order to induce the error estimate of (40),we establish the estimate for the solution of one iteration of Euler’s scheme (41).

where C is a positive constant independent of h,H,Δt,and ν.

Proof Due to the regularity assumption of u,we have

In the same way,we denote

In the similar way to Theorem 2,we can deduce the right-hand of the above inequality except the last two terms. For these terms,we have

The theorem is obtained.

Theorem 5 Let (u,p) ∈ X × Q be the solution of (1)-(4) and (uⁿ⁺¹ _h ,pⁿ⁺¹ _h ) ∈ X_h × Q_h be the solution of (40). If u_t ∈ L²(0,T ;H1(Ω)),u ∈ L_∞([0,T] × ) ∩ L²(0,T ;H²(Ω)),and p ∈ L²(0,T ;H¹(Ω)),the following inequality holds:

where C is a positive constant independent of h,H,Δt,and ν.

Proof The following second-order backward finite difference scheme holds:

Subtracting (40) from (34) taken at t = tⁿ⁺¹ and multiplying 4Δt on each side of the reduction, we have Combining with (42),we obtain

Let us study the terms of the right-hand side of (48). The first term e1 is bounded by Theorem 4.

Using (47),we have

The third term is treated as follows:

The fourth term is bounded as follows:

For the non-linear terms,we have

For the stabilized term,we obtain Using Lemma 2,we have and From (48) to (49),we have Applying Gronwall’s lemma yields the following inequality:

Theorem 5 holds. Theorem 6 Let (u,p) ∈ X × Q be the solution of (1)-(4) and (uⁿ⁺¹ _h ,pⁿ⁺¹ _h ) ∈ X_h × Q_h be the solution of (40). If u_t ∈ L²(0,T ;H¹(Ω)),u ∈ L_∞([0,T] × ) ∩ L²(0,T ;H²(Ω)),and p ∈ L²(0,T ;H¹(Ω)),the following inequality holds:

Proof The proof is similar to that of Theorem 3. From Lemma 1,we can get

Let q_h = 0. Subtracting (40) from (34) taken at t = tⁿ⁺¹ yields 5 Numerical computation

We do two experiments to confirm the results. The first one is defined on the unit square with the velocity u = (u₁,u₂) and the pressure p expressed by

Let T = 1,Re = 1 /ν = 10⁶,α = 3h,Δt = h,and H = 2h. The numerical calculation is performed by the nonconforming element pair of v¹_h ×Q¹ _h. The results are presented in Table 1.

The second numerical calculation is defined on the unit square with the velocity u = (u1,u2) and the pressure p expressed by

Let T = 1,Re = 1 /ν = 10⁶,α = h,Δt = 0.01,and H = 2h. The numerical calculation has been performed by the nonconforming element pair of v²_h × Q² _h. The results are presented in Table 2.

As shown in Tables 1 and 2,we notice that for Re = 10⁶,the numerical calculation results of the stabilized method agree well with those of the theoretical analysis.