The Postprocessed Mixed Finite-Element Method for the Navier--Stokes Equations

SIAM J. NUMER. ANAL.c 2005 Society for Industrial and Applied Mathematics Vol. 43, No. 3, pp. 1091–1111 THE POSTPROCESSED MIXED FINITE-ELEMENT METHOD FOR THE NAVIER–STOKES EQUATIONS∗ BLANCA AYUSO†, BOSCO GARC´ IA-ARCHILLA‡,AND JULIA NOVO† Abstract. A postprocessing technique for mixed finite-element methods for the incompressible Navier–Stokes equations is studied. The technique was earlier developed for spectral and standard finite-element methods for dissipative partial differential equations. The postprocessing amounts to solving a Stokes problem on a finer grid (or higher-order space) once the time integration on the coarser mesh is completed. The analysis presented here shows that this technique increases the convergence rate of both the velocity and the pressure approximations. Numerical experiments are presented that confirm both this increase in the convergence rate and the corresponding improvement in computational efficiency. Key words. Navier–Stokes equations, mixed finite-element methods AMS subject classifications. 65M60, 65M20, 65M15, 65M12 DOI. 10.1137/040602821 1. Introduction. This paper in a sense culminates the development of a postprocessing technique to increase the accuracy and computational efficiency of Galerkin methods for dissipative partial differential equations introduced in [18]. We turn to the equations which gave rise to this postprocessing technique, the incompressible Navier–Stokes equations, and we address those Galerkin methods for these equations which, when complex-shaped bodies are present, are acknowledged to be of wider applicability, mixed finite-element (MFE) methods. The postprocessing technique we study here was originally developed for spectral methods [18], [19]. At that moment, either its analysis and understanding or its development seemed to depend heavily on the properties of the Fourier modes, although this was not a shortcoming to prove its usefulness in the study of nonlinear shell vibrations [27]. In later works [13], [14], the dependence on the Fourier modes was overcome. Of particular importance to the present work, besides [14], has been the development of the postprocessing technique for finite-element methods in [20], [15]. In [20], it was devised how to carry out the postprocessing without the help of an approximate inertial manifold [11], [12], a concept more suited to spectral methods and eigenfunction expansion. In [15], it is shown what gains can be expected when postprocessing low-order elements. As is usually the case with MFE methods, it is the experience and understanding gained in previous works (see [14], [15], [16], [17], [18], [19], [20], and the references cited therein) with simpler equations and methods which has allowed the present one to be written. Furthermore, although for simplicity we focus on Hood–Taylor [26] elements, the postprocessing technique can be easily adapted to other kinds of mixed ∗Received by the editors January 2, 2004; accepted for publication (in revised form) March 18, 2005; published electronically September 23, 2005. http://www.siam.org/journals/sinum/43-3/60282.html †Departamento de Matem´aticas, Universidad Aut´onoma de Madrid, Madrid, Spain (blanca. a[email protected], julia.no[email protected]). The research of the first author was supported by project HPRNCT-2002-00284. The research of the third author was supported by DGI-MCYT under project MTM2004-02847 (cofinanced by FEDER funds) and by JCYL under project VA044/03. ‡Departamento de Matem´atica Aplicada II, Universidad de Sevilla, Sevilla, Spain (bosco@ matina.us.es). This author’s research was supported by DGICYT project BFM2003-00336. 1091 Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 1092 BLANCA AYUSO, BOSCO GARC´ IA-ARCHILLA, AND JULIA NOVO elements. In fact, in [3] (see also [5]) the so-called mini-element is shown to render similar gains as Hood–Taylor elements when postprocessed if the provisions in [15] are taken into account. Let us describe what this postprocessing technique is. We consider the incompressible Navier–Stokes equations, which, in appropriate dimensionless variables, can be written as ut−νΔu+(u·∇)u+∇p=f,(1.1) div(u)=0 in a bounded domain Ω ⊂Rd(d=2,3) with smooth boundary subject to homogeneous Dirichlet boundary conditions u=0on∂Ω. In (1.1), uis the velocity field, pthe pressure, and fa given force field. Suppose that for the solution uand p corresponding to a given initial condition u(·,0) = u0;(1.2) we are interested in its value at a certain time T>0. We first compute MFE approximations uhand phto the velocity and pressure, respectively, by integrating in time the corresponding discretization of (1.1)–(1.2) from t=0tot=T. Then, in the postprocessing step, we obtain an approximation to the solution ˜u,˜pof the Stokes problem −νΔ˜u+∇˜p=f−d dt uh(T)−(uh(T)·∇)uh(T) div(˜u)=0 in Ω, ˜u=0 on∂Ω. (1.3) The MFE of this last step is either the same-order Hood–Taylor element over a finer grid or a higher-order Hood–Taylor element over the same grid. The rate of convergence of the discrete velocity and pressure in the resulting method is proved to be the same as the rate of convergence of the MFE used in the postprocessed step. The overcost of the postprocessed procedure is nearly negligible since the Stokes problem using the enhanced MFE is solved only once, when the time integration has been completed. In this respect, it radically differs from some other research [2], [32], with low-order MFEs for the Navier–Stokes equations that also developed from the ideas in [11] and [12], since in [2] and [32] computations with the enhanced element or on the finer grid are carried out all the way through the interval (0,T]. Some superconvergence results are obtained in the paper and are used as a tool to get the rate of convergence of the postprocessed method. In particular, we derive a superconvergence result for the error between the MFE approximation to the velocity and the discrete Stokes projection introduced in [24]. For simplicity of analysis, we derive these results under the strong regularity hypotheses in (2.2), which, as pointed out in [24], are unrealistic in practical situations. In a more practical setting, assumptions (2.2) should be assumed from some positive time t0>0 onwards, and, as we comment in section 2, computations (and their analysis) up to this time should take into account the lower regularity at t=0. Finally, we remark that recent research [16], [17] has shown the usefulness of the postprocessing technique in obtaining efficient a posteriori error estimators in partial differential equations of evolution, a field much less developed than in the case of steady problems. The application of the postprocessing technique to get a posteriori error estimates for Navier–Stokes equations using the results obtained in this paper will be the subject of future work. Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php POSTPROCESSED MFE METHOD FOR NAVIER–STOKES EQUATIONS 1093 The rest of the paper is as follows. In section 2 we recall some properties of MFE methods and collect some inequalities to be used later. In section 3 we first specify the postprocessing technique and then carry out the convergence analysis. Finally, in section 4 numerical experiments are presented to assess the capabilities of the new technique. 2. Preliminaries and notations. Let Ω be a bounded domain in Rd,d=2,3, not necessarily convex, but of class Cm,m≥3, and let Hand Vbe the Hilbert spaces H={u∈L2(Ω))d,|div(u)=0,u·n|∂Ω=0},V={u∈H1 0(Ω))d,|div(u)=0}, endowed with the inner product of L2(Ω)dand H1 0(Ω)d,respectively. For 1 ≤q≤∞ and l≥0, we consider the standard Sobolev spaces, Wl,q(Ω)d, of functions with derivatives up to order lin Lq(Ω), and Hl(Ω)d=Wl,2(Ω)d. The norm in Hl(Ω)dwill be denoted by ·lwhile ·−lwill represent the norm of its dual space. We consider also the quotient spaces Hl(Ω)/Rwith norm pHl/R= inf{p+cl|c∈R}. We shall frequently use the following Sobolev’s imbeddings [1]. There exists a constant C=C(Ω,q) such that for q∈[1,∞),q <∞, it holds that vLq(Ω)d≤CvWs,q(Ω)d,1 q≥1 q≥1 q−s d>0,v∈Ws,q(Ω)d.(2.1) For q=∞, (2.1) holds with 1 q<s d. Let Π : L2(Ω)d−→ Hbe the Leray projector that maps each function in L2(Ω)d onto its divergence-free part. We denote by Athe Stokes operator in Ω: A:D(A)⊂H−→ H, A=−ΠΔ,D(A)=H2(Ω)d∩V. Applying Leray’s projector to (1.1), the equations can be written in the form ut+νAu+B(u, u)=Πfin Ω, where B(u, u) = Π((u·∇)u). In what follows we will assume that the solution (u, p) of (1.1)–(1.2) satisfies max 0≤t≤Tu(t)r+p(t)Hr−1/R<∞,max 0≤t≤Tut(t)r+pt(t)Hr−1/R<∞.(2.2) We refer the reader to [30] for a study about the regularity of the solutions of the Navier–Stokes equations. Notice, however, that, as pointed out in [24], it is unrealistic to assume such a strong regularity up to time t= 0. The assumption in (2.2) is for simplicity in the analysis. In a more realistic setting, t= 0 should be replaced by some positive time t0, and error bounds requiring less regularity such as those in [24] and [25] should be considered from t=0tot=t0. In order to maintain the accuracy levels that a higher regularity would allow from t0onwards, computations up to t=t0 should be carried out on an adequate finer grid. Notice also that among the conditions to ensure (2.2) (see, e.g., Theorem 4 in [23]) is that Ω is of class Cr. Let Th=(τh i,φ h i)i∈Ih,h>0, be a family of partitions of suitable domains Ωh, where the parameter his the maximum diameter of the elements τh i∈T hand φh i are the mappings of the reference simplex τ0onto τh i. We restrict ourselves to quasiuniform and regular meshes Th. Let r≥2, we consider the finite-element spaces  Sh,r =χh∈C0(Ωh)|χh|τh i◦φh i∈Pr−1(τ0)⊂H1(Ωh), ◦ Sh,r =χh∈C0(Ωh)|χh|τh i◦φh i∈Pr−1(τ0),χ h(x)=0∀x∈∂Ωh⊂H1 0(Ωh), Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 1094 BLANCA AYUSO, BOSCO GARC´ IA-ARCHILLA, AND JULIA NOVO where Pr−1(τ0) denotes the space of polynomials of degree at most r−1onτ0.Asa consequence of restricting our study to quasi-uniform partitions, the following inverse inequality holds (see, e.g., [9, Theorem 3.2.6]) ∀τ=τh i∈T h, with diam(τ)=hτ≤h, vh∈(◦ Sh,r)d: (2.3) vhWm,q(τ)d≤Chl−m−d(1 q−1 q)vhWl,q(τ)d,0≤l≤m≤2,1≤q≤q≤∞. In order to guarantee convergence of the MFE approximation, we choose a stable combination of two finite-element spaces (see [7]). We introduce the finite-element spaces in which our MFE approximation to (u, p) will be carried out. We shall denote by (Xh,r,Q h,r−1) the so-called Hood–Taylor element, where Xh,r =◦ Sh,rd,Q h,r−1= Sh,r−1∩L2(Ωh)/R,r≥3. For this mixed element a uniform inf-sup condition is satisfied (see [26], [6]), that is, there exists a constant β>0 independent of the mesh grid size hsuch that inf qh∈Qh,r−1 sup vh∈Xh,r (qh,∇·vh) vh1qhL2/R≥β.(2.4) The approximate velocity solution belongs to the discretely divergence-free space Vh,r =Xh,r ∩χh∈H1 0(Ωh): Ωh qhdiv(χh)=0 ∀qh∈Qh,r−1. We observe that for the Hood–Taylor element, Vh,r is not a subspace of V. For any v∈C0(Ω)d, we consider the standard interpolant operator Ih:C0(Ω)d−→ Xh,r. Let v∈Hr(Ω)d∩H1 0(Ω)d; it is well known that Ihsatisfies v−Ih(v)L2(Ω∩Ωh)d+hv−Ih(v)H1(Ω∩Ωh)d≤ChrvHr(Ω)d.(2.5) We briefly discuss next under what circumstances (2.5) can be extended to a global estimate (i.e., to an estimate in Ω and not just in Ω∩Ωh). The interpolation operator Ih(v) is extended by zero in Ω \Ωh, and defining δ(h) = maxx∈∂Ωhdist(x, ∂Ω), one obtains v−Ih(v)L2(Ω)d+hv−Ih(v)H1(Ω∩Ωh)d≤C(hr+δ(h))vHr(Ω)d.(2.6) For x∈Ω∩Ωh, (2.5) (and so (2.6)) follows from standard theory of interpolation and the Bramble–Hilbert lemma (see, e.g., [9, p. 192]). For x∈Ω\Ωh,v(x) can be bounded by means of the mean-value theorem, v−Ih(v)L2(Ω\Ωh)d=vL2(Ω\Ωh)d≤δ(h)∇vL2(Ω)d. We observe that using isoparametric elements δ(h)≤Chr, and so in (2.6) the righthand side is further bounded by ChrvWr,q (Ω)d(see [9, section 4.4]). As regards the global estimate for the gradient, isoparametric modification is not enough to preserve the optimal approximability properties of the finite-element space. Following [3], we shall assume in what follows the use of superparametric elements at the boundary. By this type of approximation we mean that δ(h)≤Ch2r−2so that the outside effects will Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php POSTPROCESSED MFE METHOD FOR NAVIER–STOKES EQUATIONS 1095 not pollute the optimal estimate. Under these assumptions [3], [4], the interpolant Ih satisfies v−Ih(v)L2(Ω)d+hv−Ih(v)H1(Ω)d≤ChrvHr(Ω)d.(2.7) Notice then that the condition δ(h)≤Ch2r−2allows us to forget about the discrepancies between Ω and Ωhin most of the arguments that follow. Observe, however, that one must then assume that Ω is piecewise of class C2r−2. For each fixed time t∈[0,T] the solution (u, p) of (1.1)–(1.2) is also the solution of a Stokes problem with right-hand side f−ut−(u·∇)u. We will denote by (sh,q h)∈(Xh,r,Q h,r−1),its MFE approximation satisfying ν(∇sh,∇φh)−(qh,∇·φh)=ν(∇u, ∇φh)−(p, ∇·φh) =(f−ut−(u·∇u),φ h)∀φh∈Xh,r,(2.8) (∇·sh,ψ h)=0 ∀ψh∈Qh,r−1. We observe that sh=Sh(u):V−→ Vh,r is the so-called discrete Stokes projection of the solution (u, p) of (1.1)–(1.2) (see [24]) and satisfies (∇Sh(u),∇χh)=(∇u, ∇χh)−(p, ∇·χh)=(f−ut−(u·∇)u, χh)∀χh∈Vh,r. The following bound holds for 2 ≤l≤r: u−sh0+hu−sh1≤Chlul+pHl−1/R.(2.9) The proof of (2.9) for Ω = Ωhcan be found in [25]. For the general case superparametric approximation at the boundary is assumed; see [3], [4]. Under the same conditions, the bound for the pressure is [21] p−qhL2/R≤Cβhl−1ul+pHl−1/R,(2.10) where the constant Cβdepends on the constant βin the inf-sup condition (2.4). Since we are assuming that Ω is of class Cmwith m≥3 (and that δ(h)≤Ch2r−2) using standard duality arguments and (2.9), one obtains [3], [4] u−sh−s≤Chr+s(ur+pHr−1/R),0≤s≤min(r−2,1).(2.11) Let Πh,r :L2(Ω)d−→ Vh,r be the discrete Leray’s projection defined by demanding that (Πh,r(u),χ h)=(u, χh)∀χh∈Vh,r. By definition, the projection is stable in the L2norm. For divergence-free functions, by writing Πh,ru=(Π h,ru−Sh(u))+Sh(u) and using the quasi-uniformity of the meshes, one easily shows that Πh,ru1≤Cu1∀u∈V.(2.12) We will denote by Ahthe discrete Stokes operator defined by (∇vh,∇φh)=(Ahvh,φ h)=A1/2 hvh,A1/2 hφh∀vh,φ h∈Vh,r. Since Ahis a discrete self-adjoint operator, it is easy to show that, for each 0 ≤α<1, there exists a positive constant Cα, which is independent of h, such that Aα he−tAh0≤Cαt−α∀0≤α<1.(2.13) Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 1096 BLANCA AYUSO, BOSCO GARC´ IA-ARCHILLA, AND JULIA NOVO In our analysis we shall frequently use the following relations for f∈L2(Ω)d: A−s/2 hΠh,rf0≤Chsf0+A−s/2Πf0,s=1,2,(2.14) A−s/2Πf0≤Chsf0+A−s/2 hΠh,rf0,s=1,2.(2.15) These inequalities are readily deduced from the estimates A−s/2−A −s/2 hΠh,r0≤ Chsfor s=1,2 [29]. Similarly, since ∀vh∈Vh,r,(A−1/2 hΠh,rf,vh)=(f,A−1/2 hvh), it follows that A−1/2 hΠh,rf0≤Cf−1,(2.16) and since ∀v∈V, we have (A−1/2Π(Πh,rf),v)=(Π h,rf,A−1/2v)=(f,Πh,rA−1/2v), from (2.12) it follows that A−1/2Π(Πh,rf)0≤Cf−1,f∈L2(Ω)2.(2.17) 2.1. The suggested method. Let us suppose that we want to approximate the solution of (1.1)–(1.2) at time T.Ford= 3, the final time Tis assumed to satisfy 0<T<T ∗, where T∗is the critical time until which the existence and uniqueness of a strong solution of (1.1)–(1.2) has been proven. The postprocessing technique can be seen as a two-level method. We first compute the MFE approximation to (1.1)–(1.2) at time T. Given uh(0) an initial approximation to u(0), we find that uh:[0,T]−→ Xh,r and ph:[0,T]−→ Qh,r−1satisfy (˙uh,φ h)+ν(∇uh,∇φh)+bh(uh,u h,φ h)+(∇ph,φ h)=(f,φh)∀φh∈Xh,r,(2.18) (∇·uh,ψ h)=0 ∀ψh∈Qh,r−1,(2.19) where bh(·,·,·) is a suitable discrete approximation to its continuous counterpart. As an initial condition we will take uh(0) = Sh(u0), although other choices are possible. In the second step, the discrete velocity and pressure (uh(T),p h(T)) are postprocessed. Basically, we enhance this approximation by solving a single discrete Stokes problem, via MFE. The MFE in this step, denoted by (  X,  Q), is either •the same-order Hood–Taylor element over a finer grid (  X,  Q)=(X˜ h,r,Q ˜ h,r−1), r≥3,˜ h<h,or •a higher-order Hood–Taylor element over the same grid (  X,  Q)=(X˜ h,r+1,Q ˜ h,r), r≥3, ˜ h=h. That is, we shall search for (˜uh,˜ph)∈( X,  Q) satisfying ν∇˜u˜ h,∇˜ φ+∇˜p˜ h,˜ φ=(f, ˜ φ)−b˜ h(uh(T),u h(T),˜ φ)−(˙uh(T),˜ φ)∀˜ φ∈ X,(2.20) ∇·˜u˜ h,˜ ψ=0 ∀˜ ψ∈ Q.(2.21) We will denote by  Vthe corresponding discretely divergence-free space that can be either  V=V˜ h,r or  V=Vh,r+1 depending on the selection of the postprocessed space. The discrete Leray’s projection into  Vwill be denoted by  Π˜ h, and we will represent by  A˜ hthe discrete Stokes operator acting on functions in  V. The postprocessed Hood–Taylor approximation to the velocity, ˜u˜ h, is the solution of the pressure-free formulation ν∇˜u˜ h,∇˜χh=f, ˜χh−b˜ huh(T),u h(T),˜χh−˙uh(T),˜χh∀˜χh∈ V.(2.22) Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php POSTPROCESSED MFE METHOD FOR NAVIER–STOKES EQUATIONS 1097 In the next section, we show that the solution (˜uh,˜ph) of (2.20)–(2.21) is a more accurate approximation to the solution of (1.1)–(1.2) than the Galerkin MFE approximation (uh,p h) that solves (2.18)–(2.19). For the discrete approximation to the nonlinear term, following [24], we define bh in the following way: bh(uh,v h,φ h)=((uh·∇)vh,φ h)+1 2(div(uh)vh,φ h)∀uh,v h,φ h∈Xh,r ⊂H1 0(Ω)d. For all u, v ∈H1 0(Ω)d, the corresponding continuous operator will be denoted by F(u, v)=(u·∇)v+(1/2) div(u)v. Extending the definition of bhto functions in H1 0(Ω)d(not necessarily in Xh,r), we observe that ∀u, v, w ∈H1 0(Ω)d,b h(u, v, w)= (F(u, v),w).It is straightforward to verify that bhenjoys the skew-symmetry property bh(u, v, w)=−bh(u, w, v)∀u, v, w ∈H1 0(Ω)d.(2.23) Let us observe that B(u, v)=ΠF(u, v)ifu∈V. Finally, we shall denote by Bh(u, v)=Π h,rF(u, v)∀u, v ∈H1 0(Ω)d. 3. Analysis of the postprocessed method. This section is devoted to the analysis of convergence of the postprocessed MFE method. Our first aim will be to show a superconvergence result for the error between the MFE approximation to the velocity uhand the Stokes projection of the velocity field u,sh. This superconvergence behavior occurs for both the L2and H1norms, as will be shown in Theorem 3.7 and Corollary 3.8, respectively. In the first part of the section, we shall concentrate our efforts in Theorem 3.7. It will be achieved by a stability plus consistency argument (Propositions 3.2 and 3.6, respectively). For the purpose of analysis, we shall mainly be concerned with the pressure-free formulation associated with (2.18)–(2.19). If (uh,p h) is the MFE approximation to the solution (u, p) of (1.1)–(1.2), then uh∈Vh,r is the solution of (˙uh,χ h)+ν(∇uh,∇χh)+bh(uh,u h,χ h)=(f,χh)∀χh∈Vh,r,(3.1) which can also be expressed in abstract operator form as ˙uh+νAhuh+Bh(uh,u h)=Π h,rf.(3.2) The Stokes projection shsatisfies the abstract equation ˙sh+νAhsh+Bh(sh,s h)=Π h,rf+Th,(3.3) where Th(t) is the truncation error, defined as Th(t)= ˙sh−Πh,r(ut)+Bh(sh,s h)−Bh(u, u).(3.4) Let us now consider mappings vh:[0,T]−→ Vh,r satisfying the following threshold condition: sh(t)−vh(t)0≤cτh2∀t∈[0,t 1],0<t 1≤T.(3.5) We define their truncation error as  Th=˙vh+νAhvh+Bh(vh,v h)−Πh,rf.(3.6) Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 1098 BLANCA AYUSO, BOSCO GARC´ IA-ARCHILLA, AND JULIA NOVO Prior to establishing the stability restricted to the threshold (3.5) (Proposition 3.2), we prove a lemma which provides some estimates for the convective term. Lemma 3.1. Let (u, p)be the solution of the Navier–Stokes problem (1.1)–(1.2). Let sh=Sh(u)be the discrete Stokes projection of the velocity field uand let vh: [0,T]−→ Vh,r satisfy the threshold condition (3.5). Then, there exists a constant K>0, independent of t1in (3.5), such that ∀t∈[0,t 1], F(sh(t),s h(t)) −F(vh(t),v h(t))0≤Ksh(t)−vh(t)1,(3.7) F(sh(t),s h(t)) −F(vh(t),v h(t))]−1≤Ksh(t)−vh(t)0,(3.8) where the constant K=Kcτ,max0≤t≤T(u(t)2+p(t)H1/R). Proof. In order to simplify the notation, we shall omit the dependence on tin the proof. Denote by eh=vh−sh. We proceed by standard duality arguments, using the splitting F(vh,v h)−F(sh,s h)=F(vh,e h)+F(eh,s h).(3.9) We start by showing (3.7). We first observe that F(eh,s h)0= sup φ0=1 (eh·∇sh,φ)+1 2((∇·eh)sh,φ) ≤CehL2d/(d−1)(Ω)d∇shL2d(Ω)d+Ceh1sh∞. Let us show that both sh∞,∇shL2d(Ω)dare bounded. Since, by virtue of Sobolev’s imbeddings (2.1), we have sh∞≤C∇shL2d(Ω)d, we only need to bound the second term. Application of the inverse inequality (2.3) and the error estimates (2.9) and (2.7) together with (2.1) give ∇shL2d(Ω)d≤Ch −(1+d) 2(sh−u0+u−Ihu0)+∇IhuL2d(Ω)d (3.10) ≤Ch(3−d)/2(u2+pH1/R)+CuW1,2d(Ω)d≤K. Using again (2.1) we obtain ehL2d/(d−1)(Ω)d≤Ceh1/2≤Ceh1, and so F(eh,s h)0≤Keh1.As regards the other term in (3.9), the same arguments lead to F(vh,e h)0= sup φ0=1 (vh·∇eh,φ)+1 2((∇·vh)eh,φ) ≤Cvh∞eh1+C∇vhL2d(Ω)dehL2d/(d−1)(Ω)d. As before, to conclude we must show that the above norms of vhare bounded. We only need to handle ∇vhL2d(Ω)d. Using the inverse inequality (2.3) and the threshold conditions (3.5) and (3.10), we find ∇vhL2d(Ω)d≤h −(1+d) 2vh−sh0+∇shL2d(Ω)d≤cτh(3−d)/2+K≤K. Therefore, (3.7) follows. We now show (3.8). Applying (3.9), we find F(vh,v h)−F(sh,s h)−1≤F(vh,e h)−1+F(eh,s h)−1,(3.11) Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php POSTPROCESSED MFE METHOD FOR NAVIER–STOKES EQUATIONS 1099 so that the proof is reduced to estimate each of the above negative norms on the right-hand side. Using the skew-symmetry property (2.23), one gets for the first term: F(vh,e h)−1= sup φ1=1 −((vh·∇)φ, eh)−1 2((∇·vh)φ, eh) ≤sup φ1=1 eh0vh∞φ1+eh0∇·vhL2d/(d−1) φL2d(Ω)d)≤Keh0. Regarding the other term in (3.11), integrating by parts, we obtain F(eh,s h)−1= sup φ1=1  1 2(eh·∇)sh,φ −1 2(eh·∇)φ, sh ≤sup φ1=1 (eh0∇shL2d/(d−1)(Ω)dφL2d(Ω)d+eh0φ1sh∞)≤Keh0. This finishes the proof of (3.8). Proposition 3.2 (stability). Let T>0be fixed; let sh=Sh(u)be the discrete Stokes projection of the velocity field usolution of (1.1)–(1.2) and let vh:[0,T]−→ Vh,r satisfy the threshold condition (3.5). Then, there exists a positive constant Ks>0 such that ∀t1≤T, the following estimate holds: max 0≤t≤t1sh(t)−vh(t)0≤eKst1sh(0) −vh(0)0 (3.12) + max 0≤t≤t1   t 0 e−ν(t−s)Ah[Th(s)− Th(s)]ds   0, where Th(s)and  Th(s)are the truncation errors given in (3.4) and (3.6), respectively. Proof. We denote by eh=sh−vh. Subtracting (3.6) from (3.3), it follows that ehsatisfies the error equation ˙eh(t)+νAheh(t)=Bh(vh(t),v h(t)) −Bh(sh(t),s h(t)) + Th(t)− Th(t). Then, by integrating the above error equation from time 0 up to time t, we find that eh(t)=e−νtAhΠh,reh(0) + t 0 e−ν(t−s)AhΠh,r[Bh(vh,v h)−Bh(sh,s h)]ds +t 0 e−ν(t−s)AhΠh,r[Th(s)− Th(s)]ds. Since {e−νtAhΠh,r}t>0is a contraction e−νtAhΠh,reh(0)0≤eh(0)0.As regards the second term, estimates (2.13), (2.16), and (3.8) from Lemma 3.1 lead to    t 0 e−ν(t−s)Ah[Bh(sh,s h)−Bh(vh,v h)]ds   0 ≤C1/2 √νt 0 A−1/2 hΠh,rF(sh,s h)−Πh,rF(vh,v h) 0 √t−sds ≤KC1/2 √νt 0 eh(s)0 √t−sds. Then, eh(t)0≤eh(0)0+KC1/2 √νt 0 eh(s)0 √t−sds +t 0 e−ν(t−s)Ah[Th(s)− Th(s)]ds. Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 1106 BLANCA AYUSO, BOSCO GARC´ IA-ARCHILLA, AND JULIA NOVO Then, applying  A−1 ˜ hto both sides of the above equation, we obtain ˜u˜ h− S˜ h(u)0≤1 ν  A−1 ˜ h Π˜ h[F(u(T),u(T)) −F(uh(T),u h(T))] 0 +  A−1 ˜ h Π˜ h[ut(T)−˙uh(T)] 0. Thus, our aim is reduced to estimate each of the above norms. As regards the nonlinear term, taking into account (2.14), with s=2,wefind   A−1 ˜ h Π˜ h[F(u, u)−F(uh,u h)] 0≤C˜ h2 F(u, u)−F(uh,u h) 0 + A−1Π[F(u, u)−F(uh,u h)] 0. Now, using estimates (3.7) from Lemma 3.1 and (3.13) from Lemma 3.4, we get  A−1 ˜ h Π˜ h[F(u, u)−F(uh,u h)] 0≤C˜ h2u−uh1+C(u−uh−1+u−uh0u−uh1). To conclude, we shall estimate each term in both sums. The required estimates in the L2and H1norms are granted by Corollary 3.9. As regards the estimate in the H−1 norm, note that by means of (2.11) and (3.17), one readily finds u−uh−1≤u−sh−1+sh−uh−1≤u−sh−1+sh−uh0 ≤Chr+1(ur+pHr−1/R)+Khr+1|log(h)|. Then, we finally get   A−1 ˜ h Π˜ h[F(u, u)−F(uh,u h)] 0≤Khr+1|log(h)|. We next deal with the estimate for the time derivative. Applying again (2.14) with s= 2 together with estimates (3.20) and (3.21) from Lemma 3.10, we reach   A−1 ˜ h Π˜ h[ut(T)−˙uh(T)] 0≤˜ h2ut(T)−˙uh(T)0+A−1Π[ut(T)−˙uh(T)]0 ≤K˜ h2hr−1|log(h)|+Khr+1|log(h)|≤Khr+1|log(h)|. Hence the proof for the L2norm is also finished. Theorem 3.15. Let T>0be fixed. Let (uh,p h)be the MFE approximation to the solution (u, p)of (1.1)–(1.2).Let(˜u˜ h,˜p˜ h)be the postprocessed MFE approximation at time T. Then, there exists a constant K(u, p, ν)such that (i) if the postprocessing element is ( X,  Q)=(X˜ h,r,Q ˜ h,r−1), then p(T)−˜p˜ hL2/R≤Cβ(˜ h)r−1u(T)r+p(T)Hr−1/R (3.34) +K(u, p, ν, β)hr|log (h)|; (ii) if at time Tthe solution (u(T),p(T)) belongs to (Hr+1(Ω)d∩V)×Hr(Ω)/R, and the postprocessing element is ( X,  Q)=(Xh,r+1,Q h,r), then p(T)−˜p˜ hL2/R≤Cβhru(T)r+1 +p(T)Hr/R+K(u, p, ν, β)hr|log (h)|.(3.35) Proof. Let us denote by ˜q˜ hthe MFE approximation to the pressure p(T) obtained by solving the Stokes problem (2.8) at time Tin the postprocessed space (  X,  Q). Adding and subtracting ˜q˜ h,weget p(T)−˜p˜ hL2/R≤p(T)−˜q˜ hL2/R+˜q˜ h−˜p˜ hL2/R. Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php POSTPROCESSED MFE METHOD FOR NAVIER–STOKES EQUATIONS 1107 The first term can easily be estimated applying (2.10): p(T)−˜q˜ hL2/R≤⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ Cβ(˜ h)r−1u(T)r +p(T)Hr−1/R, X,  Q=(X˜ h,r,Q ˜ h,r−1), Cβhru(T)r+1 +p(T)Hr/R, X,  Q=(Xh,r+1,Q ˜ h,r). Let us now bound the second term. Using the equations that satisfy ˜p˜ hand ˜q˜ h((2.20), (2.8), respectively), we deduce ˜p˜ h−˜q˜ h,∇·˜ φ=ν∇˜u˜ h− S˜ h(u),∇˜ φ+F(uh,u h) −F(u, u),˜ φ+˙uh−ut,˜ φ∀˜ φ∈ X. Using the inf-sup condition (2.4), we obtain β˜p˜ h−˜q˜ hL2/R≤ν˜u˜ h− S˜ h1+F(uh,u h)−F(u, u)−1+uh−ut−1. Taking into account (3.33), (3.8) from Lemma 3.1, and (3.22) from Lemma 3.10, we reach ˜p˜ h−˜q˜ hL2/R≤1 β(Khr|log(h)|+u−uh0+Khr|log(h)|), so that, applying Corollary 3.8, we have completed the proof. Remark 3.1. Observe that for the velocity we used piecewise polynomials of degree at least 2. In general, the postprocessed method does not increase the rate of convergence in the L2norm in the linear case although an improvement in the energy norm is obtained. The application of the postprocessing technique to the mini-element approximation to Navier–Stokes equations is studied in [3], [5]. 4. Numerical experiments. In this section, we present some numerical experiments in order to support the analysis developed in the paper and to assess the merit of the postprocessed method when compared with the standard MFE method. We consider the Navier–Stokes equations (1.1) over the domain Ω = [0,1] ×[0,1] subject to homogeneous Dirichlet boundary conditions. The value of the viscosity in the experiments is ν= 1, and the final time is T=1.2. We set to zero the initial velocity field u0(1.2) and choose the external force fso that the exact solution is u1(x, y, t)=−6·[1 −cos(πt)]sin3(πx) sin2(πy) cos(πy),(x, y, t)∈Ω×[0,T], u2(x, y, t)=6·1−cos(πt)sin2(πx) sin3(πy) cos(πx),(x, y, t)∈Ω×[0,T], p(x, y, t) = (sin(2πt)/2)sin4(πx) + sin3(πy)−p0,(x, y, t)∈Ω×[0,T], where p0denotes the mean of the pressure. In spite of the simplicity of this solution and its lack of physical meaning, we remark that our main interest has been to check the improvement in the rate of convergence achieved with the postprocessing technique and whether this also increases the efficiency of the standard MFE approximation. In our calculations we take the so-called regular pattern triangulations of Ω, which are induced by the set of nodes (i/N, j/N), 0 ≤i, j ≤N, where N=|Ω|/h is an integer. The MFE approximation to (1.1)–(1.2) is carried out using the Hood–Taylor element (Xh,3,Q h,2) that we will denote by P2P1. That is, we use Lagrange quadratic Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 1108 BLANCA AYUSO, BOSCO GARC´ IA-ARCHILLA, AND JULIA NOVO 10 1 10 −7 10 −6 10 −5 10 −4 10 −3 10 −2 N=|Ω|/h L 2 −norm velocity error slope=−3.7521 slope=−4.0632 slope=−3.0694 10 1 10 −5 10 −4 10 −3 10 −2 10 −1 N=|Ω|/h H 1 −norm velocity error slope=−2.9489 slope=−2.9671 slope=−1.9714 Fig. 4.1.Convergence diagrams for the first component of the velocity with P2P1(continuous line), P3P2(dashed-dotted line), and the postprocessed method with P3P2(dashed line). On the left the errors are measured in the L2norm (circles ◦) and on the right in the H1norm (diamonds ♦). elements for the approximation to the velocity and linear elements to approximate the pressure. For the postprocessing step, due to the smoothness of the solution (u, p), we perform the experiments not only with the same MFE over a finer grid, (Xh,3,Q h,2), h<h, but also with the higher-order Hood–Taylor element over the same grid, (Xh,4,Q h,3); i.e., Lagrange cubic for the velocity and Lagrange quadratic for the pressure. This element will be denoted by P3P2. For the time integration we use the well-known semi-implicit method where linear terms are approximated by the implicit midpoint rule (i.e., the Crank–Nicolson method) and nonlinear terms by the two-step explicit Adams formula (see, e.g., [8, p. 105]). The modified Stokes problems that arise at each step are solved by means of a standard projection method [31, pp. 27–28] (see also [3, section 4.6]). For each hused in the triangulations of Ω, every experiment was carried out with different values of the time step dt. There is always a point, depending on h, at which further reduction of the time step dt does not reduce the errors anymore. This means that the error arising from the time discretization is smaller than the error arising from the MFE discretization. To avoid wrong conclusions from our numerical experiments, we have been careful to ensure that the dominant error in all the computations presented here is the spatial discretization error. For the computational cost in the efficiency diagrams shown here, we use the largest time step among those in which the spatial discretization error is dominant. In what follows, we use the same symbols in all the plots to represent the relative errors. For the velocity we plot the errors in the first component. Similar errors are obtained for the second. The different methods are distinguished by the line used to join the symbols. For the MFE-P2P1 approximation, we use continuous line, and for the MFE-P3P2 dashed-dotted line. The MFE-P2P1 has been postprocessed in two different ways: using P3P2 (dashed line) and refining the mesh (dotted line). In Figure 4.1 we present two convergence diagrams showing the errors committed by the methods when used with h=|Ω|/N, N =8,16,32,64, both in the L2norm (left) and the H1norm (right). We have plotted the errors of the MFE-P2P1 and P3P2 methods and the postprocessed errors with P3P2. One can observe that the postprocessing technique with P3P2 provides an approximate velocity with about the Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php POSTPROCESSED MFE METHOD FOR NAVIER–STOKES EQUATIONS 1109 10 1 10 −4 10 −3 10 −2 10 −1 N=|Ω|/h L 2 /ℜ−norm pressure error slope=−2.9076 slope=−2.9864 slope=−1.9714 101 10−4 10−3 10−2 10−1 N=|Ω|/h H 1 0 −norm velocity error h´=h/2 h´ ≈ h 3/2 Fig. 4.2.Left: convergence diagram for the pressure approximation with P2P1(continuous line), P3P2(dashed-dotted lines), and the postprocessed method with P3P2(dashed lines). Right: convergence diagram for the first component of the velocity approximation with P2P1(continuous line) and the postprocessed P2P1over finer grids. same accuracy as that corresponding to the MFE-P3P2 method. This is especially true for the H1norm, in which the two methods produce virtually the same errors. Measures of the slopes of the plots confirm the rates predicted by the theory (i.e., the errors in the plots decrease like Nslope = const.h−slope). Similar conclusions can be reached from the errors of the approximations to the pressure in Figure 4.2 (left). Except for the first point, which correspond to h=1/8, the postprocessed errors lies on a line (almost) parallel to the one joining the MFE- P3P2 errors. The rate of convergence of these two methods is one unit larger than that of the MFE-P2P1 in agreement with what the theory predicts. In Figure 4.2 (right), we plot the errors obtained postprocessing the MFE-P2P1 refining the grid. We have represented the errors measured in the H1norm; similar results have been obtained for the L2norm. In view of Theorem 3.14, in order to get a gain of one order of convergence in the H1norm, we should use a mesh of size h≈h3/2. The improvement in the rate of convergence of the postprocessed method can be observed in the figure. We can also observe in the plot that using a refined mesh of size h=h/2 (only one regular refinement), the errors are considerably reduced. In fact, observe that the postprocessed error with h=h/2 is almost the same as that of the standard MFE-P2P1 carried out using a mesh of size h/2 over the full interval [0,T]. This fact can be of interest when the cost of the postprocessing step with a refined mesh of size a power of his not affordable for computational reasons. The relevant question now is whether the improvement in the rate of convergence also implies improved efficiency. In Figure 4.3, we have represented the same errors as in Figure 4.1 (right) and Figure 4.2 (left) against the smallest amount of time needed to achieve them. We have also plotted the errors of the postprocessed method refining the mesh (Figure 4.2 (right)). In the plot we observe that the efficiency of the two postprocessing procedures is very similar. We can conclude that the postprocessed method really improves the efficiency of the standard MFE method for both approximations to the velocity and to the pressure. For any error that we may demand, the postprocessed method achieves that error in less computing time than the standard P2P1 and P3P2-MFE methods. The reason for this improvement is that the error of Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 1110 BLANCA AYUSO, BOSCO GARC´ IA-ARCHILLA, AND JULIA NOVO 10 1 10 2 10 3 10 4 10 5 10 6 10 7 10 −4 10 −3 10 −2 10 −1 cputime H 1 0 −norm velocity error 10 1 10 2 10 3 10 4 10 5 10 6 10 7 10 −3 10 −2 10 −1 cputime L 2 /ℜ−norm pressure error Fig. 4.3.Efficiency diagrams for the first component of the velocity in the H1norm (left) and the pressure (right) with P2P1(continuous line), P3P2(dashed-dotted line) and the postprocessed method with P3P2(dashed line) and refining the grid (dotted line). the MFE-P2P1 method is reduced when the postprocessing is done, but this is done at very little cost: that of solving a single discrete Stokes problem at the final time. All numerical experiments were carried out on a Pentium IV, with 1 GB of Rimm memory, under the Solaris8 (Intel) operating system, with SUN Workshop 5 compilers. The programs were written in Fortran 77. REFERENCES [1] R. Adams,Sobolev Spaces, Academic Press, New York, 1975. [2] A. Ait Ou Ammi and M. Marion,Nonlinear Galerkin methods and mixed finite elements: Two-grid algorithms for the Navier–Stokes equations, Numer. Math., 68 (1994), pp. 189– 213. [3] B. Ayuso,The Postprocessed Mixed Finite Element Method for the Navier–Stokes Equations, Ph.D. thesis, Universidad Aut´onoma de Madrid, Madrid, Spain, 2003. [4] B. Ayuso and B. Garc´ ıa-Archilla,Regularity constants of the Stokes problem. Application to finite-element methods on curved domains, Math. Models Methods Appl. Sci., 15 (2005), pp. 437–470. [5] B. Ayuso, J. de Frutos, and J. Novo,Improving the accuracy of the mini-element approximation to Navier–Stokes equations, submitted for publication. [6] F. Brezzi and R. S. Falk,Stability of higher-order Hood–Taylor methods, SIAM J. Numer. Anal., 28 (1991), pp. 581–590. [7] F. Brezzi and M. Fortin,Mixed and Hybrid Finite Element Methods, Springer-Verlag, New York, 1991. [8] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang,Spectral Methods in Fluid Dynamics, Springer Ser. Comput. Math., Springer-Verlag, Berlin, 1988. [9] P. G. Ciarlet,The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. [10] P. Constantin and C. Foias,Navier–Stokes Equations, Chicago Lectures in Math., The University of Chicago, Chicago, 1988. [11] C. Foias, O. Manley, and R. Temam,Modelling of the interaction of small and large eddies in two dimensional turbulent flows, RAIRO Mod´el. Math. Anal. Num´er., 22 (1988), pp. 93– 118. [12] C. Foias, G. R. Sell, and E. S. Titi,Exponential tracking and approximation of inertial manifolds for dissipative nonlinear equations, J. Dynam. Differential Equations, 1 (1989), pp. 199–243. [13] J. de Frutos, B. Garc´ ıa-Archilla, and J. Novo,A postprocessed Galerkin method with Chebyshev or Legendre polynomials, Numer. Math., 86 (2000), pp. 419–442. Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php POSTPROCESSED MFE METHOD FOR NAVIER–STOKES EQUATIONS 1111 [14] J. de Frutos and J. Novo,A spectral element method for the Navier–Stokes equations with improved accuracy, SIAM J. Numer. Anal., 38 (2000), pp. 799–819. [15] J. de Frutos and J. Novo,Postprocessing the linear-finite-element method, SIAM J. Numer. Anal., 40 (2002), pp. 805–819. [16] J. de Frutos and J. Novo,A posteriori error estimation with the p-version of the finite element method for nonlinear parabolic differential equations, Comput. Methods Appl. Mech. Engrg., 191 (2002), pp. 4893–4904. [17] J. de Frutos and J. Novo,Element-wise a posteriori estimates based on hierarchical bases for non-linear parabolic problems, Internat. J. Numer. Methods Engrg., 63 (2005), pp. 1146–1173. [18] B. Garc´ ıa-Archilla, J. Novo, and E. S. Titi,Postprocessing the Galerkin method: A novel approach to approximate inertial manifolds, SIAM J. Numer. Anal., 35 (1998), pp. 941–972. [19] B. Garc´ ıa-Archilla, J. Novo, and E. S. Titi,An approximate inertial manifold approach to postprocessing Galerkin methods for the Navier–Stokes equations, Math. Comp., 68 (1999), pp. 893–911. [20] B. Garc´ ıa-Archilla and E. S. Titi,Postprocessing the Galerkin method: The finite-element case, SIAM J. Numer. Anal., 37 (2000), pp. 470-499. [21] V. Girault and P. A. Raviart,Finite Element Methods for Navier–Stokes Equations, Springer-Verlag, Berlin, 1986. [22] D. Henry,Geometric Theory of Semilinear Parabolic Equations, Springer-Verlag, Berlin, 1991. [23] J. G. Heywood,The Navier–Stokes equations: On the existence, regularity and decay of solutions, Indiana Univ. Math. J., 29 (1980), pp. 639–681. [24] J. G. Heywood and R. Rannacher,Finite element approximation of the nonstationary Navier–Stokes problem. Part I: Regularity of solutions and second-order error estimates for spatial discretization, SIAM J. Numer. Anal., 19 (1982), pp. 275–311. [25] J. G. Heywood and R. Rannacher,Finite element approximation of the nonstationary Navier–Stokes problem. Part III: Smoothing property and higher order error estimates for spatial discretization, SIAM J. Numer. Anal., 25 (1988), pp. 489–512. [26] P. Hood and C. Taylor,A numerical solution of the Navier–Stokes equations using the finite element technique, Comput. Fluids, 1 (1973), pp. 73–100. [27] C. R. Laing, A. McRobie, and J. M. T. Thompson, The post-processed Galerkin method applied to non-linear shell vibrations, Dynam. Stability Systems, 14 (1999), pp. 163–181. [28] J. C. L´ opez Marcos and J. M. Sanz-Serna,Stability and convergence in numerical analysis, III. Linear investigation of nonlinear stability, IMA J. Numer. Anal., 8 (1988), pp. 71–84. [29] H. Okamoto,On the semidiscrete finite element approximation for the nonstationary Stokes equation, J. Fac. Sci. Univ. Tokyo Sect. IA Math., 29 (1982), pp. 241–260. [30] R. Temam,Navier Stokes Equations and Nonlinear Functional Analysis, CBMS-NSF Regional Conf. Ser. in Appl. Math. 41, SIAM, Philadelphia, PA, 1983. [31] S. Turek,Efficient Solvers for Incompressible Flow Problems. An Algorithmic and Computational Approach, Springer-Verlag, Berlin, 1999. [32] O. Walsh,On Approximate Inertial Manifolds for the Navier–Stokes Equations Using Finite Elements, Ph.D. thesis, University of British Columbia, Vancouver, BC, Canada, 1994. Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php