1 Introduction

We consider the following generalized Oseen equations with a high Reynolds number Re:

(1)

where is a bounded domain with the polygonal boundary ∂Ω, u=(u^[1], u^[2]) denotes the velocity, p is the pressure, and f=(f^[1], f^[2]) ∈ L²(Ω) is the body force. ε=Re^-1 is the viscosity coefficient, and 0 < ε≪1. Equation (1) can be seen as the linearization of the stationary (σ=0) and the nonstationary (σ>0) time discretization Navier-Stokes equations. Here, we assume that b=(b^[1], b^[2])∈ (W^{1, ∞}(Ω))² is divergence-free, i.e., ∇·b=0.

The Oseen equations are important for the development of robust and efficient numerical methods in solving the incompressible Navier-Stokes equations. The finite element method (FEM) of the Oseen equations need stability for advective-dominated flows and compatibility between the approximating spaces of velocity and pressure. The streamline diffusion FEM (SDFEM) is an efficient method, which was first proposed by Hughes and Brooks^[1], and has been successfully used to solve compressible and incompressible flow problems, advection-diffusion problems, etc.^[2-9]. The analysis of the SDFEM, including the case of equal order interpolation, was performed in Ref. [8], in which various terms were added to the weak formulation. The residual-based stabilisation FEMs were studied in Refs. [10] and [11]. The residual-free bubble method of the nonconforming P₁^nc/P₀ element pair was investigated in Ref. [2]. An important feature of the method is that the stabilization does not generate an additional coupling between the mass equation and the momentum equation, which will appear when the SDFEM is applied to equal-order interpolation. Different local projection type stabilization methods have also been applied to the Oseen equations^[12-15]. However, the superconvergent results about the FEMs for this problem have not been found up to now.

Recently, Chen et al.^[16] used the integral identities of triangular linear conforming elements to study the two-dimensional time-dependent advection-diffusion equations, and obtained a uniform optimal-order error estimate. Some superconvergent estimates under the ε weighted energy norm were also derived through introducing an interpolation postprocessing operator. Inspired by this idea, we will use low order rectangular Bernardi-Raugel (B-R) elements to approximate the generalized Oseen equations under the inf-sup condition. Then, with the help of the Bramble-Hilbert lemma, we will derive the corresponding optimal error estimates about the velocity u, the pressure p, and the superconvergent results about the pressure p in L²-norm when the exact solutions are reasonably smooth.

The rest of this paper is organized as follows. In Section 2, we briefly introduce some notations and preliminaries. In Section 3, we present the corresponding streamline diffusion (SD) scheme with B-R elements, and provide the stability analysis and error estimates of this scheme. In Section 4, an interpolation postprocessing operator of pressure is constructed to improve the convergent behavior. In Section 5, a slightly modified SD scheme is proposed to improve the range of the parameter ε. In Section 6, a numerical experiment is carried out to confirm our theoretical results.

We will use the notations for the Sobolev spaces W^{m, p}(Ω) with norm ‖·‖_{m, p} and semi-norm |·|_{m, p}, and W^{m, p}(K) with norm ‖·‖_{m, p, K} and semi-norm |·|_{m, p, K}, where m and p are nonnegative integer numbers. Especially, when p=2, p will be omitted in the above notations.

2 Notations and preliminaries

The variational formulation for Eq. (1) is written as follows: find (u, p)∈V × M such that

(2)

where

The existence and uniqueness of the solution for Eq. (2) can be verified from the following two conditions^[17]:

(3)

(4)

where Z = {v∈V:∇·v = 0}, and α and β are positive constants independent of h.

Let T_h be a rectangular partition of the convex polygon domain Ω. For a given element K ∈ T_h with the center point (x_K, y_K), its four vertices are denoted by a₁(x₁, y₁), a₂(x₂, y₂), a₃(x₃, y₃), and a₄(x₄, y₄), and its four edges are denoted by Let h_K denote the diameter of K, and . Let be the reference element in the ξη-plane with the nodes and and the edges .

We define the affine mapping by

(5)

where h_{x, K} and h_{y, K} are the half lengths of the element K along the x-and y-directions, respectively.

When

we define the associated interpolation operators as follows:

which satisfy

(6)

The B-R finite element space pair is defined by

Now, we introduce a subspace Z_h of V_h as follows:

The global interpolation operators Π_h:(H¹(Ω))²↦ V_h and I_h: L²(Ω)↦ M_h are defined by

When u ∈ (H²(Ω)∩H₀¹(Ω))², q∈ H¹(Ω)∩L₀²(Ω), and Π_hu ∈ Z_h, we have

(7)

(8)

Throughout the paper, C (with or without subscripts) denotes a positive constant independent of the mesh parameter h, and may be different at each occurrence.

Now, we give the following useful lemmas.

Lemma 1^[18]  When u ∈ (H³(Ω))² and p ∈ H²(Ω), for all v_h ∈ V_h, we have

(9)

(10)

Lemma 2  When u ∈ (H³(Ω)∩ H₀¹(Ω))² and b ∈ (W^{1, ∞}(Ω))², we have

(11)

Proof  For all v_h=(v_h^[1], v_h^[2]) ∈ V_h, we have

(12)

Now, we start to estimate each term of A_i (i=1, 2, 3, 4).

For the term A₁, we consider the functional

(13)

according to the Sobolev embedding theorem and the inverse inequality as follows:

(14)

A direct computation shows that

(15)

which gives

(16)

Hence,

(17)

The line integrals in Eq. (17) will be canceled due to the continuousness of u_xx^[1] v_h^[1] across l_i (i=2, 4) of each element K ∈ T_h. As a result, we obtain

(18)

Let

Then, we have

(19)

For the adjacent elements K_i and K_j, there holds |b^[1]|_Ki-b^[1]|_Kj|≤Ch|b^[1]|_{1, ∞}.

Therefore, we can derive

(20)

(21)

From Eqs. (19)-(21), we have

(22)

Similarly, for the term A₂, we consider the following functional:

(23)

It is easy to check that

(24)

which gives

(25)

With the Bramble-Hilbert lemma, we can derive

Applying the same techniques as for Eqs. (19)-(21) yields

(26)

Similar to the estimations of A₁ and A₂, we have

which leads to the desired result.

Lemma 3   For p ∈ H³(Ω) and v_h ∈ V_h, there holds

(27)

Proof  From Ref. [18], we have

Using the inverse inequality, we can derive the desired result (27) directly.

3 SDFEM and error estimates

The SDFE scheme of Eq. (1) is as follows: find (u_h, p_h) ∈ V_h × M_h such that for all (v_h, q_h)∈ V_h × M_h,

(28)

where

(29)

(30)

We set the parameter δ_K=δ₀h² with a positive constant δ₀. We define the norm in V_h as follows:

(31)

Lemma 4  There exist positive constants and independent of h such that

(32)

(33)

Proof  Since v_h∈V_h is continuous across the inner edges of K and ∇·b=0, we have

(34)

which implies that

(35)

Using the inverse inequality and noting 0 < ε≪1, we have

(36)

which leads to the desired result (32).

Since ‖v_h‖₀≤ C₁|v_h|₁, we have

(37)

where

Then, from Eqs. (4), (6), and (37), we have

(38)

which is the desired result (33) with . The proof is complete.

Now, we state the following error estimate of the velocity u.

Theorem 1   Assume that (u, p)∈V × M and (u_h, p_h) ∈ V_h × M_h are the solutions of Eq.(1) and Eq. (28), respectively. If

we have

(39)

Proof  Let

Then, for all θ ∈ Z_h, from Eqs. (1) and (28), we have

(40)

Now, we start to estimate I_i (i=1, 2, …, 11).

First, we can see that

(41)

(42)

(43)

An application of Lemmas 2 and 3 gives

(44)

(45)

where is an arbitrary small constant.

Second, from Eq. (34), we have

(46)

Furthermore, we have

(47)

(48)

(49)

Then, using the inverse inequality and noting δ_K=δ₀h², we can derive

(50)

where , which means that

Last, we use Green's formula and Young's inequality to obtain

(51)

Since the terms

in the above estimates can be adsorbed into the left-hand side of Eq. (40), from Eqs. (40)-(51), we have

(52)

Let

Then, the above inequality becomes

(53)

which implies the desired result (39). The proof is complete.

Consequently, we present the superclose estimate for the pressure p.

Theorem 2  Under the assumptions of Theorem 1, we have

(54)

Proof  Let e=u-u_h. For all v_h ∈ V_h, from Eqs. (1) and (28), we can derive

(55)

Obviously, we obtain

(56)

(57)

(58)

Using the inverse inequality, we can derive

(59)

Integrating the term ((b·∇)θ, v_h) by parts, we have

(60)

Combining the above estimates yields

(61)

Moreover,

(62)

Since the spaces V_h and M_h satisfy the discrete inf-sup condition^[19], we have

(63)

Thus,

(64)

which is the desired result.

4 Global superconvergent result about pressure

Now, we will introduce a proper interpolation postprocessing operator to get the global superconvergent result about pressure. For this purpose, we furthermore assume that T_2h has been obtained from T_h by dividing each element into four congruent rectangles. As in Ref. [18], we let

and define the interpolation operator I_2h on the partition T_2h by

(65)

where Q₁=span{1, x, y, xy}.

It has been shown in Ref. [18] that the interpolation operator I_2h defined by Eq. (65) is well posed and satisfies the following properties:

(66)

(67)

Theorem 3  Under the assumptions of Theorem 2, we can get the following global superconvergent result:

(68)

Proof  Since

(69)

from Eq. (66), we have

(70)

Consequently, it follows from Eqs. (67) and (54) that

(71)

The desired result can be derived from Eqs. (69)-(71).

5 A modified SDFEM

When the SDFE scheme (28) is used to solve the problem (1), we find that the numerical results are not convergent when ε is not small enough (see Table 6). This phenomenon may be caused by the asymmetric term (εΔu_h, δ_K(b·∇) v_h) in Eq. (29). Since Δu_h≠ 0 for the B-R element pair, the condition has to be added in the error estimates, which restricts the range of the parameter ε.

Now, we remove the term (εΔu_h, δ_K(b·∇) v_h) in the classical SDFE scheme (28), and consider the following modified SDFE scheme: to find (u_h, p_h) ∈ V_h × M_h such that

(72)

where

(73)

For the modified SDFE scheme (72), it is easy to check that the results of Theorems 1 and 2 are also valid if their proofs are modified slightly.

6 Numerical experiment

To observe the behavior of the numerical method, we consider the following example. In the computation, we let Ω= [0, 1]² and δ_K=h². The exact solutions u and p are given, and the function f can be calculated by Eq. (1).

In Tables 1-5, we present the errors of the velocity u in |||·|||-norm and the superclose and superconvergent results of the pressure p in L²-norm for the SDFE scheme (72) with different mesh generations, where N denotes the number of the elements. We can see that for a fixed σ, e.g., σ=1, the numerical results are uniform convergent even ε is very small.

The numerical results for the SDFE scheme (28) when σ=1 and ε=10^-6 are listed in Table 6. We can see from Table 6 that they are not convergent.

In order to examine the effects of the parameters σ and ε, we also give the corresponding results with h=1/32 when these two parameters take different values in Fig. 1. We can see that the numerical results are stabilized. This shows the properties of the scheme (72) are good.

In order to show the effect when the parameter σ becomes large, we set ε=10^-6 and σ=100. The choice of σ corresponds to a length of the time step of 0.01 in the nonstationary Navier-Stokes equations. The numerical results are showed in Fig. 2. We can see that, though the accuracies of ||p_h-I_hp||₀ and ||p-I_2hp_h||₀ are a little smaller in coarse meshes, they are soon restored to O(h²) when the meshes are refined.