A posteriori error analysis of an augmented mixed method for the Navier-Stokes equations with nonlinear viscosity

A posteriori error analysis of an augmented mixed method for the Navier-Stokes equations with nonlinear viscosity ∗ Gabriel N. Gatica† ,Ricardo Ruiz-Baier‡ ,Giordano Tierra§ Abstract In this work we develop the a posteriori error analysis of an augmented mixed finite element method for the 2D and 3D versions of the Navier-Stokes equations when the viscosity depends nonlinearly on the module of the velocity gradient. Two different reliable and efficient residual-based a posteriori error estimators for this problem on arbitrary (convex or non-convex) polygonal and polyhedral regions are derived. Our analysis of reliability of the proposed estimators draws mainly upon the global inf-sup condition satisfied by a suitable linearization of the continuous formulation, an application of Helmholtz decomposition, and the local approximation properties of the Raviart- Thomas and Cl´ement interpolation operators. In addition, differently from previous approaches for augmented mixed formulations, the boundedness of the Cl´ement operator plays now an interesting role in the reliability estimate. On the other hand, inverse and discrete inequalities, and the localization technique based on triangle-bubble and edge-bubble functions are utilized to show their efficiency. Finally, several numerical results are provided to illustrate the good performance of the augmented mixed method, to confirm the aforementioned properties of the a posteriori error estimators, and to show the behaviour of the associated adaptive algorithm. 1 Introduction In the recent work [11], a new dual-mixed finite element method for the Navier-Stokes equations with constant density and variable viscosity – depending nonlinearly on the gradient of velocity – has been introduced and analyzed. More precisely, the approach in [11] employs a technique previously applied to the Navier-Stokes equations with constant viscosity (see [12] and [32]), which is based on the introduction of a modified pseudostress tensor involving the diffusive and convective terms, and the pressure. The latter is then eliminated thanks to the incompressibility condition, and the nonlinear viscosity is handled by incorporating the gradient of velocity as an auxiliary tensor unknown. In addition, the fact that the convective term forces the velocity to live in a smaller space motivates the augmentation of the variational formulation with suitable penalty terms arising mainly from the constitutive and equilibrium equations, and the relation defining the aforementioned gradient. As a consequence, the resulting augmented scheme can be written equivalently as a fixed point equation, ∗This research was partially supported by BASAL project CMM, Universidad de Chile, by Centro de Investigaci´on en Ingenier´ıa Matem´atica (CI2MA), Universidad de Concepci´on, by CONICYT project Anillo ACT1118 (ANANUM), by the Swiss National Science Foundation through the research grant SNSF PP00P2-144922, by the Elsevier Mathematical Sciences Sponsorship Fund (MSSF), and by the Ministry of Education, Youth and Sports of the Czech Republic through the ERC-CZ project LL1202. †CI2MA and Departamento de Ingenier´ıa Matem´atica, Universidad de Concepci´on, Casilla 160-C, Concepci´on, Chile, email: [email protected] ‡Mathematical Institute, University of Oxford, A. Wiles Building, Radcliffe Observatory Quarter, Woodstock Road, Oxford OX2 6GG, United Kingdom, email: [email protected] §Department of Mathematics, Temple University, Philadelphia, PA 19122, USA, e-mail: [email protected] 1 and hence the well-known Schauder and Banach theorems, together with classical results on monotone operators, yield the well-posedness of the continuous and discrete systems. In particular, the usual discrete inf-sup conditions are not required anymore, and therefore the unique solvability of the Galerkin scheme is guaranteed with arbitrary finite element subspaces of the respective continuous spaces. The purpose of the present paper is to develop a residual-based a posteriori error analysis for the problem and method from [11]. Regarding this goal, we begin by remarking that standard Galerkin procedures such as finite element and mixed finite element methods inevitably lose accuracy when applied to nonlinear problems on quasi-uniform discretisations. The main reason for this fact is that, in general, no a priori hints on how to mesh the domains are available in these cases, and hence one normally employs adaptive algorithms that are based on a posteriori error estimates. While the list of references on a posteriori error analysis for linear and nonlinear problems is nowadays quite extensive, and includes some important contributions in recent years, some of which are mentioned in what follows, it is important to remark that most of the main ideas and associated techniques can be found in the early works [3], [38] and the references therein. Now, aiming at as complete as possible bibliographic discussion on the main topics involved in this paper, we refer first to [34] and [35], where reliable a posteriori error estimates for finite element approximations of primal variational formulations of the compressible Euler and Navier-Stokes equations have been developed. The main tools utilized there include the element residual method and the introduction of a special norm to measure the error in the velocity and pressure. Further contributions in this direction for the steady incompressible Navier-Stokes equations are found in [6] and [8], where an abstract estimate yielding Hood-Taylor and modified Hood-Taylor finite element as particular examples, and the Becker-Rannacher “dual-weighted-residual method” for optimal-control problems, are employed, respectively. Other related works dealing with the classical velocity-pressure formulation of the steady incompressible Navier-Stokes equations include [39], [37], and [33]. In particular, residual type a posteriori error estimates for a stabilised finite element method are developed in [33] by applying a general framework established by Verf¨urth for nonlinear equations. Moreover, a simple error estimator in L2−L2norm is also presented in [33] by using a duality argument. Furthermore, a posteriori error analysis for discontinuous Galerkin approximations of the Navier-Stokes equations are developed in [30] and [31]. More precisely, a duality argument of weighted type is applied in [30] to derive a posteriori bounds and an adaptive mesh design for the interior penalty discontinuous Galerkin finite element approximation of the compressible Navier-Stokes equations. In turn, upper and lower bounds for the velocity-pressure error measured in terms of the energy norm of the discretisation for the two-dimensional stationary incompressible Navier-Stokes equations (in the case of small data) are obtained in [31]. On the other hand, a number of dual mixed approaches for the Navier-Stokes equations and the derivation of corresponding a posteriori error estimates, have begun to appear only in the last decade in the literature. For instance, we can refer to [20] where quasi-optimal a priori estimates and a posteriori error estimates for a mixed finite element approximation of this system in a polygonal domain of the plane are provided. Dirichlet boundary conditions for the velocity and the velocity gradient as an additional unknown are considered in [20]. In addition, the main tools for the a posteriori error analysis developed there include an abstract nonlinear theory and the posteriori error estimates for the Stokes equations from [19]. Other contributions in this direction include [26], [27], and the already mentioned work [32], which, as explained above, suggested one of the techniques employed in [11], which, in turn, motivated the present paper. Concerning the aforementioned references, we remark that the key aspects of the analysis in [27], which deals with a velocity-pseudostress approach for a class of quasi-Newtonian Stokes flows, are a global inf-sup condition for a linearised version of the resulting nonlinear twofold saddle point operator equation, and a conveniently constructed Helmholtz 2 decomposition of the space containing the stresses of the fluid, together with its discrete counterpart. The technique from [27] was then extended in [21] to any nonlinear twofold saddle point variational formulation. Nevertheless, irrespective of it, and as we show below in the present paper, the approach from [27] and [21] can also be applied to other nonlinear problems not showing a twofold saddle point structure. Finally, and just for sake of completeness, we mention that a posteriori error estimates for mixed finite element discretisations of the evolutionary Navier-Stokes equations are derived in [10] and [16], and that an a posteriori error analysis for an iterative algorithm solving Navier-Stokes has been proposed very recently in [9]. According to the above discussion, in this paper we adapt the strategy proposed in [27] and [21] to derive reliable and efficient residual-based a posteriori error estimators for the model and augmented mixed method introduced in [11]. This means that our analysis begins with a global infsup condition for the linearisation arising from the use of the Gˆateaux derivatives of the nonlinear terms of the formulation. Then, our remaining tools include a suitable handling of the corresponding convective term of the Navier-Stokes equations, continuous and discrete Helmholtz’s decompositions, local approximation properties of the Raviart-Thomas and Cl´ement interpolation operators, inverse inequalities, and the localization technique based on triangle-bubble and edge-bubble functions. The rest of the work is organised in the following manner. Both strong and weak forms of the velocity gradient-stress-velocity formulation of the Navier-Stokes equations with nonlinear viscosity are stated in Section 2. We recall the main hypotheses on the viscosity nonlinearity and define a Galerkin scheme, detailing a particular choice of the discrete spaces and stating the a priori error bounds established in [11]. Next, in Section 3 we consider the 2D case, introduce two a posteriori error indicators, and derive the corresponding theoretical bounds yielding reliability and efficiency of each estimator. The analysis and results from Section 3 are then extended to the 3D case in Section 4. Finally, in Section 5 we collect several numerical examples illustrating the good performance and good effectivity indexes of both error estimators under diverse scenarios in 2D and 3D, and confirming the satisfactory behaviour of the corresponding adaptive refinement strategies. 2 The Navier–Stokes equations with nonlinear viscosity 2.1 Preliminaries Let us denote by Ω ⊆Rn,n∈ {2,3}, a given open, bounded, and connected domain with polygonal (polyhedral) boundary Γ, and denote by νthe outward unit normal vector on Γ. Standard notation will be adopted for Lebesgue spaces Lp(Ω) and Sobolev spaces Hs(Ω) with norm k·ks,Ωand seminorm | · |s,Ω. In addition, by Mand Mwe will refer to the corresponding vector and tensorial counterparts of the generic scalar functional space M, whereas k·k,with no subscripts, will stand for the natural norm of either an element or an operator in any product functional space. On the other hand, for any vector fields v= (vi)i=1,n and w= (wi)i=1,n, we set the gradient, divergence, and tensor product operators, as ∇v:= ∂vi ∂xji,j=1,n ,div v:= n X j=1 ∂vj ∂xj ,and v⊗w:= (viwj)i,j=1,n . Furthermore, for any tensor fields τ= (τij)i,j=1,n and ζ= (ζij)i,j=1,n, we let div τbe the divergence operator div acting along the rows of τ, and define the transpose, the trace, the tensor inner product, and the deviatoric tensor, respectively, as τt:= (τji)i,j=1,n,tr(τ) := n X i=1 τii,τ:ζ:= n X i,j=1 τijζij ,and τd:= τ−1 ntr(τ)I. 3 Finally, we define the following tensorial functional spaces: H0(div; Ω) := ζ∈H(div; Ω) : ZΩ tr(ζ) = 0 , and L2 tr(Ω) := ns∈L2(Ω) : tr s= 0o. 2.2 The boundary value problem We consider the stationary Navier-Stokes equations with nonlinear viscosity, that is −divµ(|∇u|)∇u+u· ∇u+∇p=fin Ω , div u= 0 in Ω , u=gon Γ , (2.1) where the unknowns are the velocity uand the pressure pof a non-Newtonian fluid occupying the region Ω. The flow is characterised by the nonlinear function µ:R+−→ Rdescribing the viscosity field, a volume force f∈L2(Ω), and the boundary velocity g∈H1/2(Γ) satisfying the compatibility condition RΓg·ν= 0. Uniqueness of a pressure solution of (2.1) is ensured in the space L2 0(Ω) = q∈L2(Ω) : ZΩ q= 0. Furthermore, we assume that µis of class C1, and that there exist constants µ1, µ2>0, such that µ1≤µ(s)≤µ2and µ1≤µ(s) + s µ0(s)≤µ2∀s≥0,(2.2) which imply Lipschitz-continuity and strong monotonicity of the nonlinear operator induced by µ, in the sense of the following result. Lemma 2.1 Let Lµ:= max µ2,2µ2−µ1, where µ1and µ2are the bounds of µgiven in (2.2). Then for each r,s∈L2(Ω) there holds kµ(|r|)r−µ(|s|)sk0,Ω≤Lµkr−sk0,Ω,(2.3) ZΩnµ(|r|)r−µ(|s|)so:r−s≥µ1kr−sk2 0,Ω.(2.4) Proof. See [28, Theorem 3.8] for details. Typical examples of functions satisfying (2.2) are e.g. the classical power law and the Carreau parametrisation: µ(s) := 2 + (1 + s)−1and µ(s) := α0+α1(1 + s2)(β−2)/2, where α0, α1>0 and β∈[1,2]. Now, after introducing the additional tensor unknown, characterised by the constitutive law σ:= µ(|∇u|)∇u−(u⊗u)−pIin Ω , 4 and applying the incompressibility constraint, we obtain µ(|∇u|)∇u−(u⊗u)d=σd, p =−1 ntr( σ+u⊗u) in Ω ,(2.5) which implies that the pressure can be eliminated from (2.1). Moreover, we also introduce the velocity gradient as an auxiliary variable t:= ∇uin Ω. Consequently, we end up with the following system, expressed in terms of the unknowns t,σ, and u: ∇u=tin Ω , µ(|t|)t−(u⊗u)d=σdin Ω , −div σ=fin Ω , u=gon Γ , RΩtr(σ+u⊗u) = 0 . (2.6) Notice from (2.5) that the required belonging of the pressure field to L2 0(Ω) is indeed guaranteed by the last equation of (2.6). 2.3 The augmented mixed formulation Let us now recall that the augmented (through the incorporation of suitable Galerkin redundant terms) mixed formulation for (2.6), introduced in [11] and motivated by velocity regularity requirements, reads: Find ~ t:= (t,σ,u)∈H:= L2 tr(Ω) ×H0(div; Ω) ×H1(Ω) such that A+Bu(~ t),~ s=F,~ s∀~ s:= (s,τ,v)∈H,(2.7) where ·,·stands for the duality pairing between H0and H,A:H−→ H0is the nonlinear operator A(~ t),~ s:= ZΩ µ(|t|)t:s−ZΩ σd:s+ZΩ τd:t+ZΩ u·div τ−ZΩ v·div σ +κ1ZΩnσd−µ(|t|)to:τd+κ2ZΩ div σ·div τ +κ3ZΩn∇u−to:∇v+κ4ZΓ u·v, the bounded linear functional F:H−→ R is defined as F,~ s:= hτ ν,gi+ZΩ f·nv−κ2div τo+κ4ZΓ g·v, where h·,·i denotes the duality pairing of H−1/2(Γ) and H1/2(Γ) with respect to the L2(Γ)-inner product, and for each z∈H1(Ω), Bz:H−→ H0is the bounded linear operator Bz(~ t),~ s:= −ZΩ (z⊗u)d:nκ1τd+so, for all ~ t:= (t,σ,u),~ s:= (s,τ,v)∈H. The coefficients κ1, κ2, κ3, κ4are positive parameters assuming the following values, dictated by the stability analysis of the augmented formulation: κ1∈ 0,2δµ1 Lµand κ3∈0,2e δµ1−κ1Lµ 2δ, with δ∈0,2 Lµand e δ∈(0,2), and κ2,κ4>0. 5 The unique solvability of (2.7) has been established in [11] using fixed point arguments. We omit details and refer the reader to [11, Sections 3.2 and 3.3]. Only for further use throughout the rest of the paper, we now let α(Ω) be the strong monotonicity constant of the nonlinear operator A(cf. [11, eq. (3.23)]), and let c1(Ω) and cTbe the boundedness constants specified in [11, eqs. (3.6), (3.19)]. In addition, we set (cf. [11, eq. (3.25)]) ρ0:= α(Ω) 2c1(Ω) (κ2 1+ 1)1/2.(2.8) 2.4 The Galerkin scheme Let us consider arbitrary finite dimensional subspaces Ht h,Hσ hand Hu hof the continuous spaces L2 tr(Ω), H0(div; Ω), and H1(Ω), respectively. As usual, hdenotes the size of a regular triangulation Thof Ω made up of triangles K(when n= 2) or tetrahedra K(when n= 3) of diameter hK, that is h:= max nhK:K∈ Tho. Then, the Galerkin scheme associated with the nonlinear problem (2.7) reads: Find ~ th:= (th,σh,uh)∈Hh:= Ht h×Hσ h×Hu hsuch that A+Buh(~ th),~ sh=F,~ sh∀~ sh:= (sh,τh,vh)∈Hh.(2.9) Its linearised counterpart is defined via the discrete fixed point operator Th:Hu h−→ Hu h: Th(zh) := uh∀zh∈Hu h, where uhis the third component of the unique solution (to be confirmed by Theorem 2.2) of the discrete problem: Find ~ th:= (th,σh,uh)∈Hhsuch that A+Bzh(~ th),~ sh=F,~ sh∀~ sh:= sh,τh,vh)∈Hh.(2.10) The well-posedness of the discrete problem is given in the following result Theorem 2.2 Given ρ∈(0, ρ0), with ρ0defined by (2.8), we let Wh ρ:= nzh∈Hu h:kzhk1,Ω≤ρo, and assume that the data satisfy (cf. [11, eq. (3.27)]) cTnkfk0,Ω+kgk0,Γ+kgk1/2,Γo≤ρ . Then there exists a unique solution ~ th∈Hhto (2.9) with uh∈Wh ρ, and there holds k~ thk ≤ cTnkfk0,Ω+kgk0,Γ+kgk1/2,Γo.(2.11) Proof. We refer the reader to [11, Lemma 4.1] for a detailed proof. In what follows, given an integer k≥0 and T∈ Th, we let Pk(T) (resp. e Pk(T)) be the space of polynomials on Tof degree ≤k(resp. of degree = k), and, as indicated at the beginning of Section 2.1, we set Pk(T) := [Pk(T)]nand Pk(T) := [Pk(T)]n×n. In turn, the global version of Pk(T) and Pk(T) are defined, respectively, as Pk(Th) := nvh∈L2(Ω) : vhT∈Pk(T)∀T∈ Tho and Pk(Th) := nsh∈L2(Ω) : shT∈Pk(T)∀T∈ Tho. 6 In addition, denoting by xa generic vector in Rn, we let RTk(T) be the local Raviart–Thomas space of order k, that is RTk(T) := Pk(T)⊕e Pk(T)x, and define the tensor version of the corresponding global Raviart–Thomas space of order kas RTk(Th) := nτh∈H(div; Ω) : ctτT∈RTk(T)∀c∈Rn,∀T∈ Tho. Then, examples of specific finite element subspaces Ht h,Hσ hand Hu hfulfilling Theorem 2.2 are approximations of t,σand uby piecewise polynomial tensors of degree ≤k, by tensor Raviart– Thomas elements of order k, and by continuous piecewise polynomial vectors of degree ≤k+ 1, respectively, that is Ht h:= L2 tr(Ω) ∩Pk(Th), Hσ h:= H0(div; Ω) ∩RTk(Th), Hu h:= nvh∈C(Ω) : vhT∈Pk+1(T)∀T∈ Tho. (2.12) In such a framework, the following a priori error bounds are available (see the derivation in [11]): Theorem 2.3 Suppose that there exists s > 0such that t∈Hs(Ω),σ∈Hs(Ω),div σ∈Hs(Ω), and u∈Hs+1(Ω), and that the finite element subspaces are defined by (2.12). Then, there exists C > 0, independent of h, such that for each h > 0there holds k~ t−~ thk+kp−phk0,Ω≤C hmin{s,k+1}nktks,Ω+kσks,Ω+kdiv σks,Ω+kuks+1,Ωo. 3 A posteriori error analysis: The 2D case In this section we derive two reliable and efficient residual based a posteriori error estimators for the two-dimensional version of (2.10). The corresponding a posteriori error analysis for the 3D case, which follows from minor modifications of the one to be presented next, will be addressed in Section 4. 3.1 Preliminaries Let Ehbe the set of all edges of Th, and E(T) denotes the set of edges of a given T∈ Th. Then Eh=Eh(Ω) ∪ Eh(Γ), where Eh(Ω) := {e∈ Eh:e⊆Ω},Eh(Γ) := {e∈ Eh:e⊆Γ}. Moreover, hestands for the length of a given edge e. Also, for each edge e∈ Ehwe fix a unit normal vector νe:= (ν1, ν2)t, and let se:= (−ν2, ν1)tbe the corresponding fixed unit tangential vector along e. However, when no confusion arises, we simple write νand sinstead of νeand se, respectively. Now, let v∈L2(Ω) such that v|T∈C(T) on each T∈ Th. Then, given T∈ Thand e∈ E(T)∩ Eh(Ω), we denote by [v·s] the tangential jump of vacross e, that is [v·s] := (v|T−v|T0)|e·s, where T and T0are the triangles of Thhaving eas a common edge. Similar definitions hold for the tangential jumps of scalar and tensor fields ϕ∈L2(Ω) and τ∈L2(Ω), respectively, such that ϕ|T∈C(T) and τ|T∈C(T) on each T∈ Th. Moreover, given scalar, vector, and tensor valued fields v,ϕ:= (ϕ1, ϕ2)t, and τ=: (τij), respectively, we denote curl v:= ∂v ∂x2 −∂v ∂x1!,curl(ϕ) := curl(ϕ1)t curl(ϕ2)t!and curl(τ) :=     ∂τ12 ∂x1 −∂τ11 ∂x2 ∂τ22 ∂x1 −∂τ21 ∂x2    . 7 Now, let Ih:H1(Ω) →Xhbe the vector version of the usual Cl´ement interpolation operator (cf. [15]), where Xh:= nϕh∈C(¯ Ω) : ϕh|T∈P1(T)∀T∈ Tho, and let Πh:H1(Ω) →RTk(Th) be the Raviart-Thomas interpolator, which, according to its characterisation properties (see e.g. [22, Section 3.4.1]), verifies div(Πh(τ)) = Ph(div τ)∀τ∈H1(Ω) ,(3.1) where Phis the L2(Ω)-orthogonal projector onto Pk(Th). Further approximation properties of Ihand Πhare summarised in the following lemmas (see a proof in e.g. [15] and [22, Lemmas 3.16 and 3.18], respectively). Lemma 3.1 There exist c1, c2>0, independent of h, such that for all ϕ∈H1(Ω) there holds kϕ−Ih(ϕ)k0,T ≤c1hTkϕk1,∆(T)∀T∈ Th, and kϕ−Ih(ϕ)k0,e ≤c2h1/2 ekϕk1,∆(e)∀e∈ Eh(Ω) ∪ Eh(Γ) , where ∆(T) := ∪{T0∈ Th:T0∩T6=∅} and ∆(e) := ∪{T0∈ Th:T0∩e6=∅}. Lemma 3.2 There exist C1, C2>0, independent of h, such that for all τ∈H1(Ω) there holds kτ−Πh(τ)k0,T ≤C hTkτk1,T ∀T∈ Th, and k(τ−Πh(τ)) νk0,e ≤C h1/2 ekτk1,Te∀e∈ Eh(Ω) ∪ Eh(Γ) , where Teis a triangle of Thcontaining the edge eon its boundary. 3.2 Error estimators and statement of our main results Let ~ th:= (th,σh,uh)∈Hhbe the unique solution of (2.10) and let us set the following residuals r1(~ th;0) := σd h−µ(|th|)th+ (uh⊗uh)din Ω , r2(~ th;f) := f+div σhin Ω , r3(~ th;0) := ∇uh−thin Ω , r4(~ th;g) := g−uhon Γ . (3.2) Then, we define for each T∈ Ththe (local) a posteriori error indicators θ2 1,T := kr1(~ th;0)k2 0,T +kr2(~ th;f)k2 0,T +kr3(~ th;0)k2 0,T +X e∈ E(T)∩Eh(Γ) kr4(~ th;g)k2 0,e ,(3.3) and θ2 2,T := θ2 1,T +h2 Tkcurl(th)k2 0,T +kf− Phfk2 0,T +X e∈ E(T)∩Eh(Ω) he ths 2 0,e +X e∈ E(T)∩Eh(Γ) he    dg ds−ths    2 0,e ,(3.4) 8 so that the global a posteriori error estimators are given, respectively, by θ1:= X T∈Th θ2 1,T +kr4(~ th;g)k2 1/2,Γ1/2 and θ2:= X T∈Th θ2 2,T 1/2 .(3.5) The main results of this section are as follows. Theorem 3.3 Assume that f∈L∞(Ω) and g∈H1(Γ), and that there holds (cf. (2.8)) cTnkfk+kgk0,Γ+kgk1/2,Γo≤ρ0 2.(3.6) In addition, let ~ t∈Hand ~ th∈Hhbe the unique solutions of the continuous and discrete formulations (2.7) and (2.10), respectively. Then, there exist constants Crel >0and Ceff >0, independent of h, such that Ceff θ1≤ k~ t−~ thk ≤ Crel θ1.(3.7) Theorem 3.4 Under the assumptions of Theorem 3.3, there exist constants crel >0and ceff >0, independent of h, such that ceff θ2≤ k~ t−~ thk ≤ crel θ2.(3.8) The upper and lower bounds in (3.7), (3.8) are known as the reliability and efficiency estimates, respectively, which are derived below in Sections 3.4 and 3.5 under the assumption that fand g are piecewise polynomials on Thand the induced triangulation on Γ, respectively, for each h > 0. Otherwise, higher order terms arising from polynomial approximations of these functions would appear in (3.7) and (3.8). The analysis of the first a posteriori error estimator is straightforward, taking advantage of the fact that u∈H1(Ω), which allows us to integrate by parts in some terms. For the second estimator, we exploit the properties of the Helmholtz decomposition jointly with the Cl´ement and Raviart-Thomas interpolation operators. In this case, new terms that capture the jumps between triangles will appear. 3.3 A general a posteriori error estimate In order to establish the reliability estimates of the a posteriori error estimators θ1and θ2, that is the upper bounds in (3.7) and (3.8), and as announced in the Introduction, here we follow the strategy recently proposed in [21], which is based on the definition of a linear operator depending on the Gˆateaux derivatives of the nonlinear terms of the formulation. More precisely, we begin with the following main result. Lemma 3.5 Assume that the data satisfy (3.6), and let ~ t∈Hand ~ th∈Hhbe the unique solutions of the continuous and discrete formulations (2.7) and (2.10), respectively. Then, there exists a constant C > 0, independent of h, such that k~ t−~ thk ≤ CnkR1kL2 tr(Ω)0+kR2kH0(div;Ω)0+kR3kH1(Ω)0o,(3.9) where R1∈L2 tr(Ω)0,R2∈H0(div; Ω)0, and R3∈H1(Ω)0, are defined by R1(s) := ZΩ r1(~ th;0) : s, R2(τ) := hτ ν,gi − ZΩ uh·div τ−ZΩ th:τd −κ1ZΩ r1(~ th;0) : τd−κ2ZΩ r2(~ th;f)·div τ, R3(v) := ZΩ r2(~ th;f)·v−κ3ZΩ r3(~ th;0) : ∇v+κ4ZΓ r4(~ th;g)·v, (3.10) 9 As a straightforward consequence of Lemmas 3.9 and 3.10, the identity (3.32), and the stability estimate of the Helmholtz decomposition (cf. (3.27)), we conclude the required upper bound for kR2kH0(div;Ω)0, that is kR2kH0(div;Ω)0≤C  X T∈The θ2 2,T +X T∈Thb θ2 2,T    1/2 ,(3.44) where C > 0 is a constant independent of h. We are now ready to establish the reliability estimate for θ2. Theorem 3.11 Assume that the data satisfy (3.6), and let ~ t∈Hand ~ th∈Hhbe the unique solutions of the continuous and discrete formulations (2.7) and (2.10), respectively. Then there exists crel >0, independent of h, such that k~ t−~ thk ≤ crel θ2.(3.45) Proof. Having in mind the general estimate provided by (3.9) (cf. Lemma 3.5), we begin by recalling, as announced before, that the reliability of θ2arises from the same upper bounds of kR1kL2 tr(Ω)0and kR3kH1(Ω)0given in (3.24) and (3.25), respectively, together with the new estimate of kR2kH0(div;Ω)0 provided by (3.44). In this way, and observing that the terms h2 Tkr1(~ th;0)k2 0,T ,h2 Tkr3(~ th;0)k2 0,T , and hekr4(~ th;g)k2 0,e, which form part of b θ2 2,T (cf. (3.41)), are dominated by kr1(~ th;0)k2 0,T ,kr3(~ th;0)k2 0,T , and kr4(~ th;g)k2 0,e, respectively, which appear in the aforementioned bounds, we conclude (3.45) after replacing (3.24), (3.25), and (3.44) back into (3.9). 3.5 Efficiency of the a posteriori error estimators In this section we prove the lower estimates announced in (3.7) (cf. Theorem 3.3) and (3.8) (cf. Theorem 3.4). Most of the corresponding analysis consists of deriving suitable upper bounds depending on the true errors for each one of the terms defining the local error indicators θ2 1,T and θ2 2,T . For this purpose, we make extensive use of the original system of equations (2.6), which is recovered from the augmented continuous formulation (2.7) by choosing suitable test functions and integrating by parts backwardly the corresponding equations. We begin with the efficiency estimate for θ1. Theorem 3.12 Let ~ t∈Hand ~ th∈Hhbe the unique solutions of the continuous and discrete formulations (2.7) and (2.10), respectively. Then there exists Ceff >0, independent of h, such that Ceff θ1≤ k~ t−~ thk.(3.46) Proof. We first deduce from (3.2) and the second equation of (2.6) that r1(~ th;0) = σd h−σd−µ(|th|)th−µ(|t|)t+(uh⊗uh)d−(u⊗u)d, which, applying (2.3) (cf. Lemma 2.1), and noting that kτdk0,Ω≤ kτk0,Ωfor each τ∈L2(Ω), yields kr1(~ th;0)k2 0,Ω≤3nkσh−σk2 0,Ω+L2 µkt−thk2 0,Ω+k(uh⊗uh)−(u⊗u)k2 0,Ωo.(3.47) 16 In turn, subtracting and adding uin the first term involving the tensor product ⊗, and using Cauchy- Schwarz’s inequality for the resulting expressions, we find that k(uh⊗uh)−(u⊗u)k2 0,Ω≤2nkuhk2 L4(Ω) +kuk2 L4(Ω)oku−uhk2 L4(Ω) , from which, employing that the injection ic:H1(Ω) →L4(Ω) is compact (and hence continuous) (see Rellich-Kondrachov compactness Theorem in [1, Theorem 6.3] or [36, Theorem 1.3.5]), and using from [11, Theorems 3.9 and 4.2] that kuk1,Ωand kuhk1,Ωare bounded above by ρ0(cf. (2.8)), we deduce that k(uh⊗uh)−(u⊗u)k2 0,Ω≤4ρ2 0kick4ku−uhk2 1,Ω.(3.48) On the other hand, it is readily seen from (3.2) and the remaining equations of (2.6) that kr2(~ th;f)k2 0,Ω≤ kdiv σ−σhk2 0,Ω, kr3(~ th;0)k2 0,Ω≤2n|u−uh|2 1,Ω+kt−thk2 0,Ωo, kr4(~ th;g)k2 0,Γ≤Cku−uhk2 1,Ω, and kr4(~ th;g)k2 1/2,Γ≤Cku−uhk2 1,Ω,(3.49) where the last two inequalities make use of the trace inequalities in L2(Γ) and H1/2(Γ), respectively. In this way, the required efficiency estimate (3.46) follows straightforwardly from the definition of θ1 (cf. (3.5)) and the inequalities (3.47), (3.48) – (3.49). We now aim to establish the efficiency of θ2. More precisely, thanks to the previous theorem, it only remains to prove the corresponding upper bounds for the other four terms defining θ2 2,T (cf. (3.4)). For this purpose, we proceed as in [13, 24, 4], and apply the localization technique based on triangle-bubble and edge-bubble functions, together with extension operators, discrete trace and inverse inequalities. According to the above, we now introduce additional notations and preliminary results. Given T∈ Thand e∈ E(T), we let ψTand ψebe the usual triangle-bubble and edge-bubble functions, respectively (see [38, eqs. (1.5) and (1.6)]), for which there hold: i) ψT∈Pn+1(T), ψT= 0 on ∂T, supp(ψT)⊆T, and 0 ≤ψT≤1 in T. ii) ψe|T∈Pn(T), ψe= 0 on ∂T \e, supp(ψe)⊆we:= ∪{T0∈ Th:e∈ E(T0)}, and 0 ≤ψe≤1 in we. It is well-known that, given k∈N∪ {0}, there exists an extension operator L:C(e)→C(T) that satisfies L(p)∈Pk(T) and L(p)|e=pfor all p∈Pk(e). Further properties of ψT,ψeand Lare stated in the following lemma (see e.g. [38]). Lemma 3.13 Given k∈N∪ {0}, there exist positive constants c1,c2, and c3, depending only on k and the shape regularity of the triangulations (minimum angle condition), such that for each T∈ Th and e∈ E(T), there hold , kpk2 0,e ≤c2kψ1/2 epk2 0,e ∀p∈Pk(e), and kψ1/2 eL(p)k2 0,T ≤c3hekpk2 0,e ∀p∈Pk(e). 17 The following inverse and discrete trace inequalities are also employed. Lemma 3.14 Let k,l,m∈N∪ {0}such that l≤m. Then there exists c > 0, depending only on k, l,mand the shape regularity of the triangulations, such that for each T∈ Ththere holds |q|m,T ≤c hl−m T|q|l,T ∀q∈Pk(T). Proof. See [14, Theorem 3.2.6]. Lemma 3.15 There exists C > 0, depending only on the shape regularity of the triangulations, such that for each T∈ Thand e∈ E(T), there holds kvk2 0,e ≤Ch−1 ekvk2 0,T +he|v|2 1,T ∀v∈H1(T). Proof. See [2, Theorem 3.10]. In turn, the following lemma, whose proof makes use of Lemmas 3.13 and 3.14, is applied next to bound the terms involving curl and the tangential jumps across the edges of Th. Lemma 3.16 Let ρh∈L2(Ω) be a piecewise polynomial tensor of degree k≥0on each T∈ Th, and let ρ∈L2(Ω) be such that curl(ρ) = 0 in Ω. Then, there exist c, ec > 0, independent of h, such that kcurl(ρh)k0,T ≤c h−1 Tkρ−ρhk0,T ∀T∈ Th,(3.50) and k[ρhs]k0,e ≤ec h−1/2 ekρ−ρhk0,ωe∀e∈ Eh(Ω) .(3.51) Proof. For the proof of (3.50) we refer to [7, Lemma 4.3], whereas (3.51) is a slight modification of the proof of [7, Lemma 4.4]. We omit further details. The following lemma provides the required upper bounds for the second and fourth terms on the right hand side of (3.4). Lemma 3.17 There exists C1, C2>0, independent of h, such that h2 T curl(th) 2 0,T ≤C1kt−thk2 0,T ∀T∈ Th, and he ths 2 0,e ≤C2kt−thk2 0,we∀e∈ Eh(Ω) . Proof. It suffices to apply Lemma 3.16 with ρh=thand ρ=t=∇u. Now, the third and fifth terms on the right hand side of (3.4) are estimated as follows. Lemma 3.18 There holds kf− Phfk0,T ≤2kdiv σ−div σhk0,T ∀T∈ Th.(3.52) 18 Proof. Using that f=−div σ, and then adding and subtracting Ph(div σh) = div σh, we readily find that kf− Phfk0,T ≤ kPhdiv σ−div σhk0,T +kdiv σh−div σk0,T , which yields (3.52) and ends the proof. Lemma 3.19 There exists C3>0, independent of h, such that he    dg ds−ths    2 0,e ≤C3kt−thk2 0,Te∀e∈ Eh(Γ) , where Teis the triangle of Thhaving eas and edge. Proof. It follows similarly to the proof of [26, Lemma 4.15]. We omit further details. As a consequence of Theorem 3.12 and Lemmas 3.17, 3.18, and 3.19, we are now in position to state the efficiency of θ2. Theorem 3.20 Let ~ t∈Hand ~ th∈Hhbe the unique solutions of the continuous and discrete formulations (2.7) and (2.10), respectively. Then there exists ceff >0, independent of h, such that ceff θ2≤ k~ t−~ thk. 4 A posteriori error analysis: The 3D case In this section we extend the results from Section 3 to the three-dimensional version of (2.10). Similarly as in Section 3, given a tetrahedron T∈ Th, we let E(T) be the set of its faces, and let Ehbe the set of all faces of the triangulation Th. Then, we write Eh=Eh(Ω)∪Eh(Γ), where Eh(Ω) := {e∈ Eh:e⊆Ω} and Eh(Γ) := {e∈ Eh:e⊆Γ}. Also, for each face e∈ Ehwe fix a unit normal νeto e, so that given τ∈L2(Ω) such that τ|T∈C(T) on each T∈ Th, and given e∈ Eh(Ω), we let τ×νebe the corresponding jump of the tangential traces across e, that is τ×νe:= (τ|T−τ|T0)|e×νe, where T and T0are the elements of Thhaving eas a common face. In what follows, when no confusion arises, we simple write νinstead of νe. Now, we recall that the curl of a 3D vector v:= (v1, v2, v3) is the 3D vector curl(v) = ∇ × v:= ∂v3 ∂x2 −∂v2 ∂x3 ,∂v1 ∂x3 −∂v3 ∂x1 ,∂v2 ∂x1 −∂v1 ∂x2, and that, given a tensor function τ:= (τij)3×3, the operator curl denotes curl acting along each row of τ, that is, curl(τ) is the 3 ×3 tensor whose rows are given by curl(τ) :=   curl(τ11, τ12, τ13) curl(τ21, τ22, τ23) curl(τ31, τ32, τ33) . In addition, τ×νstands for the 3 ×3 tensor whose rows are given by the tangential components of each row of τ, that is, τ×ν:=   (τ11, τ12, τ13)×ν (τ21, τ22, τ23)×ν (τ31, τ32, τ33)×ν . 19 Having introduced these notations, and employing the same definitions of the residuals given in (3.2), we now set for each T∈ Ththe local a posteriori error indicators θ2 1,T (exactly as in (3.3)), and θ2 2,T := θ2 1,T +h2 Tkcurl(th)k2 0,T +kf− Phfk2 0,T +X e∈ E(T)∩Eh(Ω) he th×ν 2 0,e +X e∈ E(T)∩Eh(Γ) hek∇g×ν−th×νk2 0,e .(4.1) In this way, the corresponding global a posteriori error estimators are defined as in (3.5), that is θ1:= X T∈Th θ2 1,T +kr4(~ th;g)k2 1/2,Γ1/2 and θ2:= X T∈Th θ2 2,T 1/2 , and the main estimates, which are the analogue of Theorems 3.3 and 3.4, are stated as follows. Theorem 4.1 Assume that f∈L∞(Ω) and g∈H1(Γ), and that there holds (3.6). In addition, let ~ t∈Hand ~ th∈Hhbe the unique solutions of the continuous and discrete formulations (2.7) and (2.10), respectively. Then, there exist constants Crel >0,Ceff >0,crel >0, and ceff >0, independent of h, such that Ceff θ1≤ k~ t−~ thk ≤ Crel θ1. and ceff θ2≤ k~ t−~ thk ≤ crel θ2. The proof of Theorem 4.1 follows very closely the analysis of Section 3, except a few issues to be described throughout the following discussion. Indeed, we first observe that the general a posteriori error estimate given by Lemma 3.5 is also valid in 3D, and that the corresponding upper bounds of kR1kL2 tr(Ω)0,kR3kH1(Ω)0, and kR2kH0(div;Ω)0yielding the reliability of θ1are the same as those given in (3.24), (3.25), and (3.26), respectively. In turn, for the reliability of θ2, we need to use a 3D version of the stable Helmholtz decomposition provided by Lemma 3.7, which is valid only for 2D. This required result and a further extension of it to a particular case of Neumann boundary conditions were established recently for arbitrary polyhedral domains in [23, Theorem 3.1]. Before it, the stability of the Helmholtz decomposition in 3D was known only for convex regions (see, e.g. [40, Proposition 4.52]), whose proof is consequence of some results from [5]. We remark that the analysis in [23] also makes use of several estimates available in [5]. Next, the associated discrete Helmholtz decomposition and the functional R2are set and rewritten exactly as in (3.31) and (3.32), respectively. Furthermore, in order to derive the new upper bound of kR2kH0(div;Ω)0, we now need the 3D analogue of the integration by parts formula on the boundary given by (3.33) (cf. Lemma 3.8). In fact, by applying again the identities from [29, Chapter I, eq. (2.17) and Theorem 2.11], we deduce that in this case there holds hcurl χ ν,φi=− h∇φ×ν,χi ∀ χ∈H1(Ω) ,∀φ∈H1(Ω) . In addition, the integration by parts formula on each tetrahedron T∈ Th, which is employed in the proof of the 3D analogue of Lemma 3.9, becomes (cf. [29, Chapter I, Theorem 2.11]) ZT curl τ:χ−ZT τ:curl χ=hτ×ν,χi∂T ∀τ∈H(curl; Ω) ,∀χ∈H1(Ω) , where h·,·i∂T is the duality pairing between H−1/2(T) and H1/2(T), and, as usual, H(curl; Ω) is the space of tensors in L2(Ω) whose curl lie also in L2(Ω). Note that the foregoing identities explain the 20 D.o.f. h e(t)r(t)e(σ)r(σ)e(u)r(u)eff(θ1)eff(θ2) iter Augmented P0−RT0−P1scheme 74 1.4142 4.7503 – 91.5264 – 11.5238 – 1.0219 0.6711 4 258 0.7071 3.1161 0.6082 53.1640 0.7837 5.2062 1.1462 1.0093 0.6720 5 962 0.3535 1.5991 0.9624 32.6029 0.7054 2.6013 1.0009 1.0073 0.6732 5 3714 0.1767 0.8003 0.9986 17.7405 0.8779 1.2612 1.0443 1.0069 0.6778 5 14594 0.0883 0.3986 1.0055 8.9096 0.9936 0.6210 1.0221 1.0070 0.6828 4 57858 0.0441 0.1990 1.0021 4.4821 0.9911 0.3090 1.0070 1.0069 0.6849 5 230402 0.0220 0.0994 1.0007 2.2497 0.9944 0.1542 1.0019 1.0069 0.6856 5 919554 0.0110 0.0497 1.0002 1.1259 0.9985 0.0771 1.0005 1.0069 0.6858 5 Augmented P1−RT1−P2scheme 173 1.4142 3.9997 – 73.6908 – 5.4413 – 0.9627 0.2804 5 589 0.7071 1.7202 1.2172 30.2571 1.4738 2.0013 1.4429 0.9792 0.2899 5 2165 0.3535 0.3219 2.4176 12.1352 1.3180 0.4759 2.0720 0.9790 0.2939 5 8293 0.1767 0.0734 2.1328 4.7582 1.7127 0.1125 2.0797 0.9807 0.2964 5 32453 0.0883 0.0184 1.9905 1.2888 1.8223 0.0260 2.1110 0.9794 0.3030 5 128389 0.0441 0.0048 1.9320 0.3577 1.8492 0.0064 1.9940 0.9799 0.3016 5 510725 0.0220 0.0012 1.9834 0.0913 1.9691 0.0016 1.9944 0.9788 0.2996 6 2037253 0.0110 0.0003 1.9987 0.0227 1.9786 0.0004 1.9953 0.9800 0.2996 5 Table 5.1: Test 1: convergence history and Newton iteration count for the Pk−RTk−Pk+1 approximations of the Navier–Stokes problem, with k= 0,1. appearing of the expressions th×νand ∇g×ν−th×νin the 3D definition of θ2 2,T (cf. (4.1)). The rest of the proof of the reliability of θ2and the entire analysis yielding the efficiency of both θ1and θ2 proceed as in Sections 3.4 and 3.5, respectively. Just for sake of completeness, we remark that most of the details concerning the 3D version of the efficiency estimates given in Lemmas 3.16, 3.17, and 3.19, can be found in (or derived from) [25, Lemmas 4.9, 4.10, 4.11, and 4.13]. 5 Numerical results We numerically investigate the performance and accuracy of the proposed augmented finite element scheme along with the properties of the a posteriori error estimators θ1and θ2derived in Section 3 and 4, respectively. In this regard, we remark that for purposes of adaptivity, which requires to have locally computable indicators, we now use that kr4(~ th;g)k2 1/2,Γ≤ckr4(~ th;g)k2 1,Γ=cX e∈ Eh(Γ) kr4(~ th;g)k2 1,e , and redefine θ1as θ1:= X T∈Th θ2 1,T 1/2 , where θ2 1,T := kr1(~ th;0)k2 0,T +kr2(~ th;f)k2 0,T +kr3(~ th;0)k2 0,T +X e∈ E(T)∩Eh(Γ) kr4(~ th;g)k2 1,e . 21 Figure 5.1: Test 1: P1−RT1−P2approximation of velocity gradient components (top panels), Cauchy stress first components and postprocessed pressure field (center panels), and viscosity with velocity components and vector directions (bottom row). Under this redefinition θ1is certainly still reliable, but efficient only up to all its terms, except for the last one, associated to the boundary Γ. Nevertheless, the numerical results to be displayed below allow us to conjecture that this modified θ1actually verifies both properties. Furthermore, a Newton-Raphson procedure is employed to linearise (2.9) (where we recall that the convective term is already linearised via a fixed-point strategy), and impose a fixed tolerance of 1E-6 on the norm of the incremental solutions. Linear systems were solved with the iterative GMRES method preconditioned with a multilevel incomplete-LU factorisation. In all subsequent cases, the viscosity is set as µ(s) := α0+α1(1 + s2)−1/2, µ1=α0< µ(s)< α1=µ2,for s≥0, and the stabilisation coefficients are chosen according to their optimal values suggested by the analysis of Section 2, that is, Lµ= max{µ2,2µ2−µ1},δ=L−1 µ,κ1=δµ1/Lµ,κ2=κ1,κ3=µ1−κ1Lµ 2δ, κ4=µ1/4. In addition, the mean value of tr σhover Ω is fixed via a penalisation strategy. 22 D.o.f. e(t)r(t)e(σ)r(σ)e(u)r(u)e(p)r(p)eriter eff(θi) Augmented P0−RT0−P1scheme with uniform refinement 794 6.3200 – 450.6412 – 11.0093 – 17.2742 – 450.8528 – 6 – 1876 5.1547 0.4741 413.9833 0.1953 7.7599 0.8135 14.6693 0.3802 414.0881 0.1946 5 – 4880 4.4945 0.2867 386.1393 0.1456 5.2116 0.8327 9.9066 0.8212 386.2006 0.1458 6 – 15480 2.8735 0.7749 260.9739 0.6787 3.8975 0.7170 6.9631 0.6108 261.0058 0.6788 6 – 55622 1.8025 0.7291 177.1064 0.6061 2.2694 0.7904 3.7081 0.9852 177.1201 0.6062 6 – 206242 1.1616 0.7588 98.3499 0.7977 1.4177 0.6963 2.2597 0.7558 98.3555 0.7977 7 – 799732 0.8884 0.7997 51.9209 0.9114 1.0376 0.7494 1.1785 0.8915 53.9234 0.8714 7 – Augmented P0−RT0−P1scheme with adaptive refinement via θ1 638 4.5173 – 411.7368 – 7.9698 – 16.911 – 411.8387 – 6 1.0013 1308 3.0896 1.0582 339.4973 0.5374 4.2947 1.4880 7.5044 2.2634 339.5373 0.5378 6 1.0007 2418 1.8071 1.7456 168.9206 2.2721 3.2822 1.4399 4.2425 1.8564 168.9419 1.5720 6 1.0009 6199 0.8863 1.5134 80.7837 1.5670 1.5660 1.2590 2.4891 1.1328 80.7932 1.3670 5 1.0012 21914 0.5149 0.9601 34.9716 1.3260 0.9218 1.0924 1.2172 1.1302 34.9793 1.2259 5 1.0018 110799 0.2231 1.0321 14.6106 1.0771 0.3215 1.0573 0.5409 1.0008 14.6140 1.0771 6 1.0020 668283 0.0940 0.9614 5.8537 1.0180 0.1245 1.0072 0.2217 0.9927 5.8552 1.0179 6 1.0022 Augmented P0−RT0−P1scheme with adaptive refinement via θ2 615 4.3007 – 402.5308 – 7.3464 – 16.1987 – 402.6208 – 6 0.6988 1243 2.9774 1.0451 320.5792 0.6597 4.3928 1.4616 8.6691 1.7768 330.6218 0.5599 6 0.7001 2614 1.6849 1.5318 165.1663 1.4669 1.8411 2.3396 4.0207 2.0671 165.1852 1.8669 6 0.7041 6529 0.8727 1.4373 76.9355 1.3692 0.8917 1.5838 2.3679 1.1567 76.9456 1.3691 6 0.7045 24145 0.4898 0.9832 33.0832 1.2906 0.4837 0.9553 1.1645 1.0852 33.0904 1.2904 6 0.7038 126045 0.2138 1.0030 13.3836 1.0949 0.2139 0.9873 0.5099 0.9994 13.3901 1.0949 5 0.7032 795283 0.0861 0.9968 5.3079 1.0043 0.0857 0.9921 0.2019 1.0057 5.3093 1.0043 6 0.7033 Table 5.2: Test 2: convergence history and Newton iteration count for the Pk−RTk−Pk+1 approximations of the Navier-Stokes problem with k= 0, and convergence of the postprocessed pressure field. Values computed under uniform (top rows) and adaptive (bottom) refinement. Test 1. Our first example concentrates on the accuracy of the augmented method when a manufactured solution of (2.7) is given by the smooth functions u=−cos(πx1) sin(πx2) sin(πx1) cos(πx2),t=∇u,˜ σ=µ(|∇u|)∇u−u⊗u−(x2 1−x2 2)I,σ=˜ σ−1 8(ZΩ tr˜ σ)I, defined on the square domain Ω = (−1,1)2. The external load fand the boundary datum gare determined from these solutions. The viscosity parameters are set as α0= 3, α1= 4. A sequence of successively refined uniform triangulations of the domain is used to present the error history displayed in Table 5.1. Discrete norms and the convergence rate between two consecutive meshes of size hand ˆ hare defined as e(t) = kt−thk0,Ω, e(σ) = kσ−σhkdiv,Ω, e(u) = ku−uhk1,Ω, r(·) = log(e(·)/ˆe(·))[log(h/ˆ h)]−1. The results reported in Table 5.1 are in accordance with the theoretical bounds established in Theorem 2.3. In addition, we also compute the global a posteriori error indicators θ1,θ2and measure their reliability and efficiency with the efficiency index. For the two orders tested, these estimators remain always bounded. The approximate solutions computed with k= 1 are displayed in Figure 5.1. 23 D.o.f. Individual error 103104105106 10−3 10−2 10−1 100 101 102 103 Ch e(t) uniform e(t) adaptive after θ1 e(t) adaptive after θ2 e(σ) uniform e(σ) adaptive after θ1 e(σ) adaptive after θ2 e(u) uniform e(u) adaptive after θ1 e(u) adaptive after θ2 D.o.f. Total error 103104105106 100 101 102 103 Ch uniform refinement adaptive refinement after θ1 adaptive refinement after θ2 Figure 5.2: Test 2: Convergence test against analytical solutions employing the lowest order family k= 0. Individual contributions (left) and global errors (right) computed on a sequence of uniformly and adaptively refined meshes according to the a posteriori error estimators θ1and θ2. Test 2. Our next numerical experiment is aimed at testing the features of adaptive mesh refinement after the a posteriori error estimators θ1,θ2. Removing the unit square from the computational domain considered in Test 1 we end up with an L-shaped domain Ω = (−1,1)2\[0,1]2. Let us further consider the same exact velocity and velocity gradient as in Test 1, but now the volumetric part of the Cauchy stress (the pressure field) assumes a different expression, yielding ˜ σ=µ(|∇u|)∇u−u⊗u−(1 −x1)(x1−1/20)2+ (x2−1/20)2−1I. Under uniform refinement, a hampering of the convergence rates is expected due to the stress singularity near the origin (where the domain has a reentrant corner). For Tests 2 and 3 we compute individual convergence rates as r(·) := −2 log(e(·)/ˆe(·))[log(N/ ˆ N)]−1, where Nand ˆ Ndenote the corresponding degrees of freedom at each triangulation, and we also define the total error, its convergence rate, and the effectivity index associated to the scheme where the mesh refinement has been applied according to a given global estimator θi,i∈ {1,2}as e:= [e(t)]2+ [e(σ)]2+ [e(u)]21/2,r:= −2 log(e/ˆe)[log(N/ ˆ N)]−1,eff(θi) := eθ−1 i. We carry out a classical adaptive mesh refinement based on the equi-distribution of the local error indicators: the diameter of each element in the new fine mesh, which is contained in a generic element Ton the initial coarse mesh, is proportional to the diameter of the old element times the ratio ˜ θi θi,T , where ˜ θistands for the average of θiover the old triangulation (cf. [38]). Table 5.2 along with Figure 5.2 report on the error history under uniform and adaptive refinement. Sub-optimal rates are observed in the first case (not only for the total error, but for each field), whereas adaptive refinement according to either a posteriori error indicator yield optimal convergence and stable effectivity indexes. Notice also that, even if the total errors have practically the same rate of convergence and the same orders of magnitude for both error estimators, θ2delivers lower errors for the velocity. Approximate solutions computed with an augmented P0−RT0−P1family are 24 Figure 5.3: Test 2: Approximate solutions computed with a lowest order family k= 0. shown in Figure 5.3, and examples of some adapted meshes generated using θ1and θ2are collected in Figure 5.4. We can observe a clear clustering of elements near the reentrant corner (where pressure and stresses exhibit high gradients), and also in the zones where the pseudo-stress and the viscosity show sharp profiles. Test 3. To conclude, we replicate Test 2 in a three-dimensional setting. This time the manufactured exact solutions adopt the form u=  cos(x1) sin(x2) sin(x3) sin(x1) cos(x2) sin(x3) −2 sin(x1) sin(x2) cos(x3) ,t=∇u, ˜ σ=µ(|∇u|)∇u−u⊗u−(1 −x2 1−x2 2−x2 3)(x1−1/20)2+ (x2−1/20)2+ (x3−1.01)2−1I, and the domain consists on the polygon Ω = (−1,1)2\[0,1]2×(0,1), so σ=˜ σ−1 9(RΩtr˜ σ)I. All remaining parameters and functions are taken as in the previous test. Again, from the computations we observe a disturbed convergence under uniform refinement, and we confirm the recovering of optimal convergence rates when using adaptive refinement guided by the a posteriori error estimator θ1(see a summary in Table 5.3). Resulting approximations after five mesh refinement steps are collected in Figure 5.5, whereas snapshots of intermediate meshes are shown in Figure 5.6. References [1] R.A. Adams and J.J.F. Fournier,Sobolev Spaces. Second edition. Pure and Applied Mathematics (Amsterdam), 140. Elsevier/Academic Press, Amsterdam, 2003. 25