Static two-grid mixed finite-element approximations to the Navier-Stokes equations

# Static two-grid mixed finite-element approximations to the Navier-Stokes equations

Javier de Frutos Departamento de Matemática Aplicada, Universidad de Valladolid. Spain. Research supported by Spanish MEC under grant MTM2010-14919 and by JCyL under grant VA001A10-1 (frutos@mac.uva.es)    Bosco García-Archilla Departamento de Matemática Aplicada II, Universidad de Sevilla, Sevilla, Spain. Research supported by Spanish MEC under grant MTM2009-07849 (bosco@esi.us.es)    Julia Novo Departamento de Matemáticas, Universidad Autónoma de Madrid, Instituto de Ciencias Matemáticas CSIC-UAM-UC3M-UCM, Spain. Research supported by Spanish MEC under grant MTM2010-14919 (julia.novo@uam.es)
###### Abstract

A two-grid scheme based on mixed finite-element approximations to the incompressible Navier-Stokes equations is introduced and analyzed. In the first level the standard mixed finite-element approximation over a coarse mesh is computed. In the second level the approximation is postprocessed by solving a discrete Oseen-type problem on a finer mesh. The two-level method is optimal in the sense that, when a suitable value of the coarse mesh diameter is chosen, it has the rate of convergence of the standard mixed finite-element method over the fine mesh. Alternatively, it can be seen as a postprocessed method in which the rate of convergence is increased by one unit with respect to the coarse mesh. The analysis takes into account the loss of regularity at initial time of the solution of the Navier-Stokes equations in absence of nonlocal compatibility conditions. Some numerical experiments are shown.

## 1 Introduction

We consider the incompressible Navier–Stokes equations

 ut−νΔu+(u⋅∇)u+∇p = f, (1) div(u) = 0,

in a bounded domain () with a smooth boundary subject to homogeneous Dirichlet boundary conditions on . In (1), is the velocity field, the pressure, the diffusion coefficient and  a given force field.

In this paper we study the following two-grid mixed finite-element method for the spatial discretization of the above equations. First, for the solution of the fully nonlinear Navier-Stokes equations (1) corresponding to a given initial condition

 u(⋅,0)=u0, (2)

the mixed finite-element approximation  over a coarse mesh of diameter is computed. Then, for any time , the postprocessed approximation is obtained as the mixed finite-element approximation over a finer mesh () to the following steady Oseen-type problem:

 −νΔ~u+(uH(t)⋅∇)~u+∇~p=f−ddtuH(t)\@@LTX@noalign\vskip3.0ptplus3.0ptminus1.0pt\omitdiv(~u)=0⎫⎪ ⎪⎬⎪ ⎪⎭\rm in~{}Ω,\omit\span\omit\span\@@LTX@noalign\vskip6.0ptplus3.0ptminus1.0pt\omit~u=0,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\,\rm on~{}% ∂Ω.

In this paper we prove that, in terms of the fine mesh diameter , this two-grid technique is of optimal order in the sense that, for appropriate choices of the coarse mesh diameter , the method has the same rate of convergence of standard mixed finite element approximations in the fine mesh. On the other hand, for a suitable value of the discretization parameter , the rate of convergence of the postprocessed approximation in terms of increases by one unit the rate of convergence of the coarse standard approximation. The improvement in precision is achieved in both the norm for the velocity and the norm for the pressure in the case of linear, quadratic and cubic elements. For other than linear elements the rate of convergence in the norm of the velocity is also increased by one unit. We remark that time evolution is performed only at the coarse mesh whereas at the fine grid the time appears only as a parameter (see equation (1)), thus the name of static two-grid method.

Two-grid or two-level methods are a well established technique for nonlinear steady problems, see [34]. In [25], [26] several two-level methods are considered to approximate the steady Navier-Stokes equations. They require solving a nonlinear system over a coarse mesh and, depending on the algorithm chosen, one Stokes problem, one linear Oseen problem or one Newton step over the fine mesh. The corresponding algorithms obtain the optimal rate of convergence in the fine mesh for appropriate choices of the coarse mesh diameter .

In the case of nonlinear evolutionary equations, two-grid techniques have been proposed and studied in [1], [22], [24], [14]. In these methods, as opposed to the method studied in the present paper, time evolution is also performed over the fine mesh. The advantage of the method studied in the present paper is that since the time integration is only carried out on the coarse mesh, computations on the fine grid can be done at selected target time levels where an improved approximation is desired, with the corresponding reduction of computing time, specially if these target time levels are sufficiently spaced in time. For this reason, although some of the two-grid methods that incorporate the evolution in time of the fine mesh approximation are more accurate, the method we present can still be more efficient in terms of computational effort for a given error level.

Two-grid techniques that integrate in time only on the coarse level have previously been developed in [16], [17] (see also [27]) for spectral methods, and later extended to mixed finite-element formulations in [3], [4], [10]. In all these works the two grid method is referred to as postprocessed Galerkin method, and, instead of (1), the approximation  is found as an approximation to the following Stokes problem

 −νΔ~u+∇~p=f−ddtuH(t)−(uH(t)⋅∇)uH(t)\@@LTX@noalign\vskip3.0ptplus3.0ptminus1.0pt\omitdiv(~u)=0⎫⎪ ⎪⎬⎪ ⎪⎭\rm in~{}Ω,\omit\span\omit\span\@@LTX@noalign\vskip6.0ptplus3.0ptminus1.0pt\omit~u=0,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\,\rm on~{}% ∂Ω.

This two-grid method will be termed standard postprocessed method, to differentiate it to that studied in the present paper, which will be termed new postprocessed method. Both, the standard and the new postprocessed methods, have the same rate of convergence. However, as already noted in [12] for nonlinear convection-diffusion problems, the new postprocessing technique produces more accurate approximations than the standard postprocessed method, for moderate to small values of the diffusion parameter . This will also be the case in the numerical experiments in the present paper for moderate values of the Reynolds number.

In the present paper we take into account the loss of regularity suffered by the solutions of the Navier-Stokes equations at the initial time in the absence of nonlocal compatibility conditions. Thus, for the analysis, we do not assume the solution  to have more than second-order spatial derivatives bounded in up to initial time , since demanding further regularity requires the data to satisfy nonlocal compatibility conditions unlikely to be fulfilled in practical situations [19], [20]. Due to the loss of regularity at , the best error bound that we can obtain is . For this reason we do not analyze higher than cubic finite elements. The same limit in the rate of convergence was found in [20] for standard mixed finite-element approximations and in [10], [14] for two-grid schemes.

In practice, any method to numerically solve evolutionary equations needs of some time discretization procedure. For brevity reasons, we have preferred to present the method in a semidiscrete manner without reference to any particular time discretization. However, we emphasize that being static, the method we present can be applied exactly in the same form, with any time discretization. The analysis of fully discrete procedures can be developed along the same lines that appear in [11], [13].

The rest of the paper is as follows. In Section 2 we introduce some preliminaries and notation. In Section 3 we carry out the error analysis of the new method. Finally, some numerical experiments are shown in the last section.

## 2 Preliminaries and notations

We will assume that is a bounded domain in , of class , for . When dealing with linear elements ( below) may also be a convex polygonal or polyhedral domain. We consider the Hilbert spaces

 H ={u∈L2(Ω)d∣div(u)=0,u⋅n|∂Ω=0}, V ={u∈H10(Ω)d∣div(u)=0},

endowed with the inner product of and , respectively. For integer and , we consider the standard spaces, , of functions with derivatives up to order in , and . We will denote by the norm in , and  will represent the norm of its dual space. We consider also the quotient spaces with norm .

We recall the following Sobolev’s imbeddings [2]: For , there exists a constant such that

 ∥v∥Lq′≤C∥v∥Ws,q,1q′≥1q−sd>0,q<∞,v∈Ws,q(Ω)d. (3)

For , (3) holds with .

The following inf-sup condition is satisfied (see [18]), there exists a constant such that

 infq∈L2(Ω)/Rsupv∈H10(Ω)d(q,∇⋅v)∥v∥1∥q∥L2/R≥β. (4)

Let be the projection onto . We denote by the Stokes operator on :

 A:D(A)⊂H⟶H,A=−ΠΔ,D(A)=H2(Ω)d∩V.

We shall assume that is a strong solution up to time , so that

 ∥u(t)∥1≤M1,∥u(t)∥2≤M2,0≤t≤T, (5)

for some constants and . We shall also assume that there exists a constant such that

 ∥f∥1+∥ft∥1+∥ftt∥1≤~M2,0≤t≤T. (6)

Finally, we shall assume that for some

 sup0≤t≤T∥∥∂⌊k/2⌋tf∥∥k−1−2⌊k/2⌋+⌊(k−2)/2⌋∑j=0sup0≤t≤T∥∥∂jtf∥∥k−2j−2<+∞,

so that, according to Theorems 2.4 and 2.5 in [19], there exist positive constants and such that the following bounds hold:

 ∥u(t)∥k+∥ut(t)∥k−2+∥p(t)∥Hk−1/R≤Mkτ(t)1−k/2, (7) ∫t0σk−3(s)(∥u(s)∥2k+∥us(s)∥2k−2+∥p(s)∥2Hk−1/R+∥ps(s)∥2Hk−3/R)ds≤K2k, (8)

where and  for some . Observe that for , we can take  and . For simplicity, we will take these values of and .

Let , be a family of partitions of suitable domains , where is the maximum diameter of the elements , and are the mappings of the reference simplex onto .

Let , we consider the finite-element spaces

 Sh,r={χh∈C(¯¯¯¯Ωh)|χh|τhi∘ϕhi∈Pr−1(τ0)}⊂H1(Ωh), S0h,r=Sh,r∩H10(Ωh),

where denotes the space of polynomials of degree at most on .

We shall denote by the so-called Hood–Taylor element [7, 21], when , where

 Xh,r=(S0h,r)d,Qh,r−1=Sh,r−1∩L2(Ωh)/R,r≥3,

and the so-called mini-element [8] when , where , and . Here, is spanned by the bubble functions , , defined by , if  and 0 elsewhere, where denote the barycentric coordinates of . For these elements a uniform inf-sup condition is satisfied (see [7]), that is, there exists a constant independent of the mesh grid size such that

 infqh∈Qh,r−1supvh∈Xh,r(qh,∇⋅vh)∥vh∥1∥qh∥L2/R≥β. (9)

The approximate velocity belongs to the discrete divergence-free space

 Vh,r=Xh,r∩{χh∈H10(Ωh)d∣(qh,∇⋅χh)=0∀qh∈Qh,r−1},

which is not a subspace of .

Let be the solution of a Stokes problem with right-hand side , we will denote by the so-called Stokes projection (see [20]) defined as the velocity component of the solution of the following problem: find such that

 ν(∇sh,∇ϕh)+(∇qh,ϕh) =(g,ϕh) ∀ϕh∈Xh,r, (10) (∇⋅sh,ψh) =0 ∀ψh∈Qh,r−1. (11)

The following bound holds for :

 ∥u−sh∥0+h∥u−sh∥1≤Chl(∥u∥l+∥p∥Hl−1/R). (12)

The proof of (12) for can be found in [20]. The bound for the pressure is [18]

 ∥p−qh∥L2/R≤Cβhl−1(∥u∥l+∥p∥Hl−1/R), (13)

where the constant depends on the constant in the inf-sup condition (9).

We consider the semi-discrete finite-element approximation  to , solution of (1)–(2). That is, given , we compute and , , satisfying

 (˙uH,ϕH)+ν(∇uH,∇ϕH)+b(uH,uH,ϕH)+(∇pH,ϕH) =(f,ϕH) ∀ϕH∈XH,r, (14) (∇⋅uH,ψH) =0 ∀ψH∈QH,r−1, (15)

where for any .

For , provided that (12)–(13) hold for , and (7)–(8) hold for , then we have

 ∥u(t)−uH(t)∥0+H∥u(t)−uH(t)∥1≤CHrt(r−2)/2,0≤t≤T, (16)

(see, e.g., [10, 19, 20]), and also,

 ∥p(t)−pH(t)∥L2/R≤CHr−1t(r′−2)/2,0≤t≤T, (17)

where if and if .

## 3 The new postprocessed method

The postprocessing technique we propose is a two-level or two-grid method. In the first level, we choose a coarse mesh of size and compute the mixed finite-element approximation   to  defined by (14)-(15). In the second level, the discrete velocity and pressure are postprocessed by solving the following linear Oseen problem: find , , satisfying for all and

 ν(∇~uh(t),ϕh)+((uH(t)⋅∇)~uh(t),ϕh)+(∇~ph(t),ϕh) = (f(t)−˙uH(t),ϕh), (18) (∇⋅~uh(t),ψh) = 0. (19)

Equations (18)-(19) can also be solved over a higher order mixed finite-element space over the same grid. For simplicity in the exposition we will only consider the case in which we refine the mesh at the postprocessing step.

Let us observe that projecting equation (18) over the discretely-free space , and avoiding for simplicity the dependence on in the notation, we get that satisfies

 ν(∇~uh,vh)+((uH⋅∇)~uh,vh)=(f−˙uH,ϕh),∀vh∈Vh,r. (20)

We now prove that equation (20) is well-posed, i.e., for small enough there exists a unique function solving (20). Let us denote by the bilinear form defined by

 BH(uh,vh)=ν(∇uh,∇vh)+((uH⋅∇)uh,vh),uh,vh∈Vh,r. (21)

We proceed to show that is coercive which implies that there exists a unique function satisfying (20). Let us also observe that once a unique is found, using the inf-sup condition (9) one easily obtains the existence and uniqueness of the pair satisfying (18)-(19).

###### Lemma 1

Let be the bilinear form defined in (21). Then, there exists a constant such that for the following bound holds:

 ∣∣BH(vh,vh)∣∣≥(ν−CHr−1+γt(r−2)/2)∥vh∥21,∀vh∈Vh,r, (22)

where if the dimension  is , and if .

###### Proof

To prove the coercivity of we follow [25, p. 2042]. Let us first observe that for any

 BH(vh,vh)=ν∥∇vh∥20−12(∇⋅uH,vh⋅vh).

Let be the orthogonal projection of over , so that applying standard finite-element theory [9] and interpolation theory on Hilbert spaces (see e. g. [31, § II.2] we have , for . Taking into account that the velocity satisfies then

 BH(vh,vh)=ν∥∇vh∥20−12(∇⋅(uH−u),vh⋅vh−qH). (23)

And then

 |(∇⋅(uH−u),vh⋅vh−qH)|≤C∥uH−u∥1∥vh⋅vh−qH∥L2(Ω)/R.

Following [25, p. 2042] we get

 ∥vh⋅vh∥γ≤C∥vh∥21, (24)

where if , and if . Using (24) together with (16) we get

 |(∇⋅(uH−u),vh⋅vh−qH)|≤CHr−1t(r−2)/2Hγ∥vh∥21.

Finally, going back to (23) we reach (22).

Let us observe that, for and , as a consequence of Lemma 1, there exists a unique satisfying (20).

We introduce now a linearized problem that will be used in the proof of Theorem 1 where we state the rate of convergence of the new method. Let be the velocity in the solution of (1)-(2). We will denote by the solution of the following linearized problem

 −νΔv+(u⋅∇)v+∇j = d (25) div(v) = 0

in the domain subject to homogeneous Dirichlet boundary conditions. Let us observe that since the divergence of is zero the bilinear form:

 B(v,w)=ν(∇v,∇w)+((u⋅∇)v,w),v,w∈V.

associated to this problem is continuous and coercive. Since the solution of (25) satisfies

 B(v,w)=(d,w),∀w∈V

by the Lax-Milgram theorem there exists a unique solution . Due to (4) there exists also a unique pressure .

We will assume in the sequel that both problem (25) and its dual problem satisfy the regularity assumption

 ∥v∥2+∥j∥H1(Ω)/R≤C∥d∥0. (26)

The regularity assumption (26) can be proved by using the analogous regularity of the Stokes problem and a bootstrap argument, see [25, Remark 2.1].

In the following lemma we state the rate of convergence of the mixed finite-element approximation to the solution of (25) defined as follows: find such that

 ν(∇vh,∇ϕh)+((u⋅∇)vh,ϕh)+(∇jh,ϕh) = (d,ϕh), ∀ϕh∈Xh,r, (27) (∇⋅vh,ψh) = 0,  ∀ψh∈Qh,r−1. (28)
###### Lemma 2

Let be the solution of (25) and let be its mixed finite-element approximation. Then, the following bounds hold for

 ∥v−vh∥0+h∥v−vh∥1 ≤ Chl(∥v∥l+∥j∥Hl−1/R), (29) ∥j−jh∥L2/R ≤ Chl−1(∥v∥l+∥j∥Hl−1/R). (30)
###### Proof

Let us denote by the Stokes projection of . More precisely, will be the solution of (10)-(11) with right-hand-side . Let us denote by . Then, from (27) and (10) we get

 ν(∇eh,∇wh)+((u⋅∇)eh,wh)=((u⋅∇)(sh−v),wh),∀wh∈Vh,r. (31)

Taking in (31) and using (3) we get

 ν∥eh∥21≤C∥u∥L2d/(d−1)∥sh−v∥1∥eh∥L2d≤C∥u∥1/2∥sh−v∥1∥eh∥1,

so that

 ∥eh∥1≤C∥sh−v∥1. (32)

Since applying (12) we conclude is bounded by the righ-hand side of (29). The bound (30) for the pressure is readily obtained by means of the auxiliary value . Subtracting (27) from (29) and applying the inf-sup condition (9) one easily gets

 β∥kh∥L2/R≤ν∥eh∥1+C∥u∥1/2∥sh−v∥1,

so that due to (32) and (12) it follows that is bounded by the right-hand side of (30). Since , applying (13) we finally prove (30).

We are left with the task of proving the bound for the norm of the error in the velocity. We will argue by duality. Let us observe that

 ∥eh∥0=supφ∈L2 φ≠0|(eh,φ)|∥φ∥0. (33)

Let us fix and let us denote by the solution of the linearized dual problem

 −νΔw−(u⋅∇)w+∇k=φ,div(w)=0,}\rm in Ω,u=0,\rm on ∂Ω. (34)

As stated before we assume that this problem satisfies the regularity assumption (26), so that

 ∥w∥2+∥k∥H1(Ω)/R≤C∥φ∥0. (35)

We will denote by the mixed finite-element approximations to . Reasoning exactly as before and applying (35) we obtain

 ∥w−wh∥1 ≤ Ch(∥w∥2+∥k∥H1/R)≤Ch∥φ∥0, (36) ∥k−kh∥L2/R ≤ Ch(∥w∥2+∥k∥H1/R)≤Ch∥φ∥0. (37)

Integrating by parts we reach

 (eh,φ) = ν(∇eh,∇w)+((u⋅∇)eh,w)−((∇⋅eh),k) = ν(∇eh,∇(w−wh))+((u⋅∇)eh,w−wh)−((∇⋅eh),k−kh) +ν(∇eh,∇wh)+((u⋅∇)eh,wh).

And then, applying (36) and (37) we reach

 |(eh,φ)| ≤ Cν∥eh∥1h∥φ∥0+C∥u∥1/2∥eh∥1h∥φ∥0+C∥eh∥1Ch∥φ∥0 (38) +|ν(∇eh,∇wh)+((u⋅∇)eh,wh)|.

Then, to conclude, it only remains to bound which by (31) is equal to . Let us decompose

 |((u⋅∇)(sh−v),wh)|≤|((u⋅∇)(sh−v),wh−w)|+|((u⋅∇)(sh−v),w)|.

Then, integrating by parts in the last term

 |((u⋅∇)(sh−v),wh)|≤C∥u∥1/2∥sh−v∥1∥wh−w∥1+|((u⋅∇)w,sh−v)|,

and the bound for the first term on the right hand side above concludes by applying (12) and (36). Finally, since

 |((u⋅∇)w,sh−v)|≤C∥u∥L2d/(d−1)∥∇w∥L2d∥sh−v∥0.

Applying Sobolev inequality (3) together with (35) and (12) we reach

 |((u⋅∇)w,sh−v)|≤C∥u∥1/2∥φ∥0hl(∥v∥l+∥j∥Hl−1/R),

so that the proof is finished.

We now state some results that will be use to get the rate of convergence of the new postprocessed method. The proof of the following lemma can be found in [15, Lemma 4] for the case and in [10, Lemma 5.1] for .

###### Lemma 3

Let be the solution of (1)–(2) and let be the mixed finite-element approximation to . Then, there exists a positive constant such that

 ∥ut(t)−˙uH(t)∥−1 ≤ Ct(r−1)/2Hr|log(H)|r′, t∈(0,T], r=2,3,4, (39) ∥A−1Π(ut(t)−˙uH(t))∥0 ≤ Ct(r−1)/2Hr+1|log(H)|, t∈(0,T], r=3,4, (40)

where when and otherwise.

The proof of the following lemma can be found in [10, p. 226].

###### Lemma 4

Let be the solution of (1)–(2) and let be the mixed finite-element approximation to . Then, there exists a positive constant such that

 ∥u(t)−uH(t)∥−1 ≤ Ct(r−1)/2Hr+1|log(H)|, t∈(0,T],r=3,4. (41)

We end this section with a theorem that states the rate of convergence of the new postprocessed method.

###### Theorem 1

Let be the solution of (1)–(2) and for let (7)–(8) hold with . Then, there exist a positive constant such that the new postprocessed approximation defined by (18)-(19) satisfies the following bounds for and small enough:

 ∥u(t)−~uh(t)∥1 ≤ Ch+Ct1/2H2|log(H)|2,r=2, (42) ∥u(t)−~uh(t)∥j ≤ Ct(r−2)/2hr−j+Ct(r−1)/2Hr+1−j|log(H)|, j=0,1, r=3,4. (43) ∥p(t)−~ph(t)∥L2/R ≤ Ct(r−2)/2hr−1+Ct(r−1)/2Hr|log(H)|r′,r=2,3,4, (44)

where for and otherwise.

###### Proof

Let us consider the linearized problem (25) with right hand side . Then, the solution of (25) is the solution of (1)–(2). Let us denote by its mixed finite-element approximation, that is the solution of (27)–(28). This approximation satisfy the error bounds (29) and (30) for . Let us decompose and . To bound the first terms in these two decompositions we will apply (29) and (30). In the rest of the proof we deal with the other two terms.

Let us denote by . Subtracting (18) from (27) it is easy to see that satisfies

 ν(∇eh,∇ϕh)+((uH⋅∇)eh,ϕh)=(˙uH−ut,ϕh)+(((uH−u)⋅∇)vh,ϕh),

for all . Taking in the above equation and applying (22) we get that for there exists a constant such that

 ∥eh∥1≤C(∥ut−˙uH∥−1+∥uH−u∥1∥vh−u∥1+∥uH−u∥0∥u∥3/2),

and applying (39) from Lemma 3, (16) and (29) we get

 ∥eh∥1≤Ct(r−1)/2Hr|log(H)|r′+Ct(r−2)/2(hHr−1