Convergence to Suitable Weak Solutions for a Finite Element Approximation of the Navier–Stokes Equations with Numerical Subgrid Scale Modeling

CONVERGENCE TO SUITABLE WEAK SOLUTIONS FOR A FINITE ELEMENT APPROXIMATION OF THE NAVIER-STOKES EQUATIONS WITH NUMERICAL SUBGRID SCALE MODELING SANTIAGO BADIA†AND JUAN VICENTE GUTI´ ERREZ-SANTACREU‡ Abstract. In this work we prove that weak solutions constructed by a variational multiscale method are suitable in the sense of Scheffer. In order to prove this result, we consider a subgrid model that enforces orthogonality between subgrid and finite element components. Further, the subgrid component must be tracked in time. Since this type of schemes introduce pressure stabilization, we have proved the result for equal-order velocity and pressure finite element spaces that do not satisfy a discrete inf-sup condition. 2010 Mathematics Subject Classification: 35Q30; 65N30; 76N10. Keywords: Navier–Stokes equations; Suitable weak solutions; Stabilized finite element methods, Subgrid scales. Contents 1. Introduction 1 2. Statement of the problem 3 2.1. Notation 3 2.2. The Navier-Stokes equations 4 3. Finite element approximation 6 3.1. Hypotheses 6 3.2. The discrete problem 7 3.3. Discrete operators 8 4. Technical preliminary results 9 5. A priori energy estimates 13 6. Convergence towards weak and suitable weak solutions 17 Appendix A. Proof of the inverse inequalities (8) 21 Acknowledgment 21 References 22 1. Introduction Incompressible Newtonian fluids are governed by the Navier-Stokes equations. The existence of solutions is known from the works by Leray [31] and Hopf [27]. However, uniqueness is still an open Date: May 21, 2018. †Universitat Polit`ecnica de Catalunya, Jordi Girona1-3, Edifici C1, E-08034 Barcelona & Centre Internacional de M`etodes Num`erics en Enginyeria, Parc Mediterrani de la Tecnologia, Esteve Terrades 5, E-08860 Castelldefels, Spain Email: [email protected]. SB was partially supported by by the European Research Council under the FP7 Program Ideas through the Starting Grant No. 258443 - COMFUS: Computational Methods for Fusion Technology and the FP7 NUMEXAS project under grant agreement 611636. SB gratefully acknowledges the support received from the Catalan Government through the ICREA Acad`emia Research Program. ‡Dpto. de Matem´atica Aplicada I, E. T. S. I. Inform´atica, Universidad de Sevilla. Avda. Reina Mercedes, s/n. E41012 Sevilla, Spain. E-mail: [email protected]. JVGS was partially supported by the Spanish grant No. MTM2015-69875-P from Ministerio de Econom´ıa y Competitividad with the participation of FEDER. 1 2 S. BADIA AND J. V. GUTI´ ERREZ-SANTACREU question. The loss of regularity is related to turbulence [24], and Leray denoted weak solutions as turbulent solution. Scheffer defined the concept of suitable weak solutions in [36] and proved a bound for the Haussdorff dimension of the singular set of a weak suitable solution. This result was later improved by Cafarelli, Kohn, and Nirenberg [9], proving that this dimension is smaller than 1. This is the sharpest regularity result so far. Suitable weak solutions of the Navier-Stokes equations can be constructed by regularization (see, e.g., [33]). More recently, Guermond proved that inf-sup stable finite element (FE) approximations having a discrete commutator property also converge to suitable weak solutions, first for periodic boundary conditions in the three-dimensional torus [22], and next on general domains and no-slip boundary conditions [23]. The Fourier method does not satisfy the required assumptions, and it is still an open question whether it provides suitable solutions. The Navier-Stokes equations have a dissipative structure, due to the viscous term. The system has a singular limit in the assymptotic regime as the Reynolds (Re) number, which is the ratio of inertia forces to viscous forces, goes to infinity. The singular limit and the fact that the system is indefinite complicate its numerical approximation. The first property requires to introduce some kind of convection stabilization, whereas the second prevents the use of the same FE space for both the velocity and pressure unknowns, the discrete system is unstable. At the continuous level, the nonlinear convective term transfers energy from the largest to the smallest scales, till reaching the Kolmogorov scale, where energy is dissipated. In direct numerical simulations (DNS) the mesh needs to be fine enough to capture the smallest scales in the flow. However, this approach is unacceptable for industrial turbulent flows, due to the limits in computational resources. In real applications, under-resolved simulations are needed. The smallest scales that can be captured in these simulations are far from the Kolmogorov scale and dissipation is negligible. Thus, one has to add so-called large eddy simulation (LES) turbulent models that add artificial diffusion mechanisms. The concept of suitability and the fact that energy is dissipated at the mesh scale in a physically consistent way have been related in [24]. Otherwise, an energy pile-up occurs at the smallest grid scales, leading to instabilities. Convection stabilization and turbulence models are strongly related. In this sense, many authors have considered so-called implicit LES (ILES) methods that do not modify the original Navier-Stokes equations but introduce additional numerical artifacts when carrying out the discretization [7, 18]. In the frame of FE techniques, one approach is to consider variational multiscale (VMS) methods [28, 29]. The idea is to use a two-scale decomposition of the original problem and provide a numerically motivated closure for the fine scale (see, e.g., [21]). A similar stabilization procedure can be used for the convective term and the pressure term, leading to methods that do not require to satisfy a discrete inf-sup condition. An alternative to traditional residual-based methods is to consider subscales that are in some sense orthogonal to the FE space. This idea has been proposed by Codina [12], where L2(Ω) orthogonality was used. This method involves global projections, which has motivated the use of local projections (see, e.g., [5, 2]). The treatment of the time dimension in the subgrid model has also been object of active research. In particular, the use of dynamic subscales methods that track the subgrid scale in time have been proposed in [12]. Even though DNS is impractical in real applications, it is better understood than stabilized or ILES schemes. The groundbreaking works by Guermond have proved that the FE Galerkin method leads to weak suitable solutions in [22, 23]. However, the extension to ILES methods is not straightforward, due to the introduction of additional terms to the numerical formulation. The analysis of these methods has usually been restricted to a priori error estimates for smooth enough solution (see, e.g., [11]). Residualbased VMS schemes are not amenable for weak convergence analysis, due to the proliferation of terms, e.g., including new velocity-pressure coupling terms. However, enforcing the modelled subgrid scales to be orthogonal to the FE space and considering the dynamic formulation in [12], the authors have proved in [4] that the resulting scheme converges to weak (turbulent) solutions of the Navier-Stokes equations. For the same scheme, long-term stability estimates and existence of a global attractor have been proved in [3]. Further, a very detailed numerical experimentation of these methods for CONVERGENCE TO SUITABLE WEAK SOLUTIONS FOR A SUBGRID FEM MODEL 3 isotropic and wall-bounded turbulent flows can be found in [13], proving that these subgrid models act as accurate turbulence models. Theoretical analyses supporting these results can also be found in [20]. In this work, we want to analyze whether VMS-type FE ILES schemes converge to suitable weak solutions in the sense of Scheffer. We prove that subgrid closures that are orthogonal and dynamic converge in fact to suitable solutions for equal order FE pairs for the velocity and pressure unknowns. The outline of the work is the following. First, we state the problem and introduce the notation in Section 2. The FE approximation based on the VMS-type ILES scheme is introduced in Section 3. Section 4 includes some technical results in fractional Sobolev spaces. Energy estimates are proved in Section 5. Finally, the convergence towards weak and suitable solutions is proved in 6. 2. Statement of the problem Throughout this paper we follow faithfully the notation used in [26] and [23] so that the reader can trace with ease the main differences between these two works and the one presented herein. 2.1. Notation. Let Ω be an open subset of R3. For p∈[1,∞], we denote by Lp(Ω) the usual Lebesgue space, i.e., Lp(Ω) = {v: Ω →R:vLebesgue-measurable,ZΩ |v(x)|pdx<∞}, with the usual modification when p=∞. This space is a Banach space endowed with the norm kvkLp(Ω) = (RΩ|v(x)|pdx)1/p if p∈[1,∞) or kvkL∞(Ω) = ess supx∈Ω|v(x)|if p=∞. In particular, L2(Ω) is a Hilbert space. We shall use (u, v) = RΩu(x)v(x)dxfor its inner product and k · k for its norm. For m∈N, we denoted by Hm(Ω) the classical Sobolev-Hilbert spaces, i.e., Hm(Ω) = {v∈L2(Ω) : ∂kv∈L2(Ω) ∀ |k| ≤ m} associated to the norm kvkHm(Ω) = X 0≤|k|≤m k∂kvk2 L2(Ω)  1 2 , where k= (k1, ..., kd)∈Ndis a multi-index and |k|=Pd i=1 ki. Let D(Ω) be the space of infinitely times differentiable functions with compact support in Ω, i.e. the space of test functions on Ω. Thus Hm 0(Ω) is defined as the completion of D(Ω) with respect to the Hm(Ω)-norm. Fractional-order Hilbert-Sobolev spaces are defined by the real method or K-method of interpolation due to Peetre and Lions [1]. Thus, we consider two spaces: Hs(Ω) = [L2(Ω), H1(Ω)]s, for s∈(0,1), and ˜ Hs 0(Ω) = [L2(Ω), H1 0(Ω)]sfor s∈[0,1]. Moreover, for s∈(0,1), Hs 0(Ω) is the closure of D(Ω) with respect to the Hs(Ω)-norm. Note that the spaces Hs(Ω) and Hs 0(Ω) coincide for s∈[0,1 2], with uniform norms [25, Th 11.1], and the spaces Hs(Ω) and ˜ Hs 0(Ω) coincide with equivalent norms [34] for s∈[0,1 2). We also consider Hs(Ω) = [H1(Ω), H2(Ω)]sfor s∈(1,2] and ˜ Hs 0(Ω) = Hs(Ω) ∩H1 0(Ω) for s∈(1,2]. The dual space of D(Ω), the space of distributions, is denoted by D′(Ω). Moreover, for s < 0, ˜ Hs(Ω) is the dual of ˜ H−s 0(Ω) and the space H−s 0(Ω) is the complexion of D(Ω) under the norm kvkH−s(Ω) = sup w∈D(Ω)\{0} (v, w) kwkHs(Ω) , We use h·,·i to denote the duality pairing. For s∈[0,1 2)∪(1 2,3 2), H−s(Ω) coincides with ˜ H−s 0(Ω). We will use boldfaced letters for spaces of vector functions, e.g. L2(Ω) in place of L2(Ω)d. We will make use of the following space of vector fields: ϑ={v∈D(Ω) : ∇ · v= 0 in Ω}. Related to the space ϑ, we consider the closures in the L2(Ω) and H1(Ω)-norm, which are characterized by H={u∈L2(Ω) : ∇ · u= 0 in Ω,u·n= 0 on ∂Ω}, V={u∈H1(Ω) : ∇ · u= 0 in Ω,u=0on ∂Ω}, 4 S. BADIA AND J. V. GUTI´ ERREZ-SANTACREU where nis the outward normal to Ω on ∂Ω. This characterization is true for locally Lipschitz-continuous domains (see [38, Theorems 1.4 and 1.6] for a detailed proof). Furthermore, L2 R=0(Ω) (resp. H1 R=0(Ω)) is the space of zero-average L2(Ω)-functions (resp. zero-average H1(Ω)-functions ). Thus, by the real method of interpolation, Hs R=0(Ω) = [L2 R=0(Ω), H1 R=0(Ω)] for s∈(0,1) (see [25]). Let Xbe a Banach space. Thus, Lp(a, b;X) denotes the space of Bochner-measurable, X-valued functions on the interval (0, T ) such that RT 0kf(s)kp Xds < ∞if 1 ≤p < ∞or ess sups∈(0,T )kf(s)kX< ∞if p=∞. Moreover, W1,1(0, T ;X) is the space of functions f∈L1(0, T ;X) and d dsf∈L1(0, T ;X) such that RT 0(kf(s)kX+kd dsf(s)kX) ds < ∞and W1,1 0(0, T ;X) is the closure of D(0, T ;X) with respect to the W1,1(0, T ;X)-norm, with D(0, T ;X) being the space of infinitely times differentiable functions defined on (0, T ) having values into Xwith compact support in (0, T ). Additionally, the dual space of W1,1 0(0, T ;X) is denoted by W−1,∞(0, T ;X′) provided that Xis separable and reflexive. The Fourier transform of a function f∈L1(R;X) is denoted by Ff(ξ) := Z+∞ −∞ e−2πit·ξf(t)dt. Let Hbe a Hilbert space and let S′(R;H) be the space of tempered distributions taking value in H. Thus, for γ∈R, one defines Hγ(R;H) = {v∈ S′(R;H); ZR (1 + |ξ|)2γkFvk2 Hdξ}, where His a Hilbert space. Additionally, the space Hγ(0, T ;H) is made up of tempered distributions in S′(0, T ;H) with the norm kvkHγ(0,T ;H)= inf v∈S′(R;H)kvkHγ(R;H), where vis the extension of vby zero off (0, T ) belonging to S′(R;H). Note that throughout this paper we use the symbol C(with or without subscripts) to represent generic positive constants which can take different values at different places. 2.2. The Navier-Stokes equations. The Navier-Stokes equations for the motion of a viscous, incompressible, Newtonian fluid can be written as ∂tu−ν∆u+ (u· ∇)u+∇p=fin Ω ×(0, T ), ∇ · u= 0 in Ω ×(0, T ),(1) with Ω being a bounded, three-dimensional domain and with 0 < T < +∞. Here u: Ω ×(0, T )→R3 represents the incompressible fluid velocity and p: Ω ×(0, T )→Rrepresents the fluid pressure. Moreover, fis the external body force which acts on the system, and ν > 0 is the kinematic fluid viscosity. These equations are supplemented by the no-slip boundary condition u=0on ∂Ω×(0, T ),(2) and the initial condition u(0) = u0in Ω.(3) The first authors dealing with the concept of weak solutions for the Navier-Stokes equations were Leray [31] for the Cauchy problem in the whole space and later Hopf [27] for the initial-boundary value problem in bounded domains. Particularly, weak solutions were called turbulent by Leray due to the possible connection between the lack of regularity of weak solutions and turbulence. Definition 2.1. A function uis said to be a weak solution of problem (1)-(2) if: u∈L∞(0, T ;H)∩L2(0, T ;V) (4) CONVERGENCE TO SUITABLE WEAK SOLUTIONS FOR A SUBGRID FEM MODEL 5 and −ZT 0 (u(t), ∂tv(t)) dt+ZT 0 h(u(t)· ∇)u(t),v(t)idt+ZT 0 ν(∇u(t),∇v(t)) dt = (u0,v(0)) + ZT 0 hf(t),v(t)idt for all v∈W1,1(0, T ;V)with v(T) = 0. Moreover, the energy inequality 1 2ku(t)k2+νZt 0 k∇u(s)k2ds≤1 2ku0k2+Zt 0 hf(s),u(s)ids(5) holds a. e. in [0, T ]. An equivalent definition for weak solutions involving the pressure term is defined as follows. Definition 2.2. A pair (u, p)is said to be a weak solution of problem (1)-(2) if: u∈L∞(0, T ;H)∩L2(0, T ;V)and p∈W−1,∞(0, T, L2(Ω)/R) and ∂tu+ (u· ∇)u−ν∆u+∇p=fin W−1,∞(0, T ;H−1(Ω)), u(0) = u0in H. Moreover, the energy inequality 1 2ku(t)k2+νZt 0 k∇u(s)k2ds≤1 2ku0k2+Zt 0 hf(s),u(s)ids holds a. e. in [0, T ]. We refer the reader to [16, Th. 1.3, Ch. V] for a proof of the equivalence between Definitions 2.1 and 2.2 with p∈ D′((0, T )×Ω), that can easily be extended to p∈W−1,∞(0, T ;L2(Ω)/R), by using de Rham’s Lemma in [37, Lm. 2]. The two previous definitions of weak solutions can be proved for Ω being a bounded, Lipschitzian domain, and f∈L2(0, T ;H−1(Ω)) only. The weak solution that will be proved in this paper requires Ω to be, for instance, convex, and f∈L2(0, T + 1; H−1(Ω)) ∩Lp(0, T + 1; Lq(Ω)), with p∈[1,2] and q∈[1,3 2] satisfying 2 p+3 q= 4. Definition 2.3. A pair (u, p)is said to be a weak solution of problem (1)-(2) if: u∈L∞(0, T ;H)∩L2(0, T ;V)and p∈H−r(0, T, H1−s R=0(Ω)) with s∈(1 2,7 10 ]and r > ¯r=3 4−s 2, and (∂tu+ (u· ∇)u−ν∆u+∇p=fin H−r(0, T ;˜ H−s 0(Ω)), u(0) = u0in H. Moreover, the energy inequality 1 2ku(t)k2+νZt 0 k∇u(s)k2ds≤1 2ku0k2+Zt 0 hf(s),u(s)ids holds a. e. in [0, T ]. Scheffer [36] introduced the definition of suitable weak solutions so as to prove a partial regularity theorem. Afterwards, Caffarelli, Kohn, and Nirenberg [9] improved Scheffer’s results, and F.-H. Lin [32] simplified the proofs of the results in [9]. Definition 2.4. A weak solution (u, p)is said to be suitable if the local energy inequality ∂t(1 2u2) + ∇ · ((1 2u2+p)u)−ν∆(1 2u2) + ν(∇u)2−f·u≤0 holds in D′((0, T )×Ω; R+). 6 S. BADIA AND J. V. GUTI´ ERREZ-SANTACREU 3. Finite element approximation 3.1. Hypotheses. Throughout this paper we will assume the following hypotheses: (H1) Let Ω be a connected, bounded, open subset of R3having a polyhedral boundary such that there exist v∈V∩H2(Ω) and p∈H1 R=0(Ω) satisfying −∆v+∇p=gin Ω, ∇ · v= 0 in Ω. (H2) Consider {Th}h>0to be a shape-regular and quasi-uniform family of simplicial and conforming meshes of Ω such that Ω = ∪K∈ThKwith h= maxK∈ThhKwhere hK= diam K. (H3) Let {Wh}h>0and {Qh}h>0be two families of finite-element spaces associated with {Th}h>0 such that Wh⊂H1 0(Ω) and Qh⊂H1 R=0(Ω). Moreover, the finite-element spaces are required to satisfy the following conditions. Let πWh:L2(Ω) →Whand πQh:L2(Ω) →Qhbe the orthogonal projections onto Whand Qh, respectively, with respect to the L2(Ω)-inner product. Furthermore, we denote π⊥ Wh(·) := (·)−πWh(·) and π⊥ Qh(·) := (·)−πQh(·). (a) There exists a constant Cinv >0, independent of h, such that, for all wh∈Wh, kwhkL∞(Ω) ≤Cinvh−3 kkwhkLk(Ω) (6) and k∇whkLk(Ω) ≤Cinvh−1kwhkLk(Ω) (7) for k∈[2,∞], kwhkH1(Ω) ≤Cinvh−1+skwhk˜ Hs 0(Ω) (8) for each s∈[0,1], and kwhk˜ Hs 0(Ω) ≤Cinvh−skwhkand kwhk ≤ Cinvh−skwhk˜ Hs 0(Ω) (9) for s∈[0,1]. (b) There exists a constant Cst(s)>0, independent of h, such that, for s∈[0,3 2), kπWhwk˜ Hs 0(Ω) ≤Cst(s)kwk˜ Hs 0(Ω) for all w∈˜ Hs 0(Ω),(10) (c) There exists a constant Cint >0, independent of h, such that, for all land s, satisfying 0≤l≤min{1, s}and l≤s≤2, there holds kπ⊥ Whwk˜ Hl 0(Ω) ≤Cinths−lkwk˜ Hs 0(Ω) for all w∈˜ Hs 0(Ω),(11) and kπ⊥ QhqkHl(Ω) ≤Cinths−lkqkHs(Ω) for all q∈Hs R=0(Ω).(12) (d) There exists Ccom >0, independent of h, such that, for 0 ≤l≤m≤1 and ϕ∈W2,∞ 0(Ω), kπ⊥ Wh(ϕwh)kHl(Ω) ≤Ch1+m−lkwhkHm(Ω)kϕkWm+1,∞ 0(Ω) for all wh∈Wh,(13) and kπ⊥ Qh(ϕqh)kHl(Ω) ≤Ch1+m−lkqhkHm(Ω)kϕkWm+1,∞ 0(Ω) for all qh∈Qh.(14) (H4) Let u0∈Vand f∈L2(0, T + 1; H−1(Ω)) ∩Lp(0, T + 1; Lq(Ω)), with p∈[1,2] and q∈[1,3 2] satisfying 2 p+3 q= 4. Hypothesis (H1) is ensured for domains having a C1,1boundary or being a convex polygon (cf. [30] or [19]) or polyhedron (cf. [14] ), with continuous dependence on f. Hypothesis (H3) is extremely flexible and allows equal-order finite-element spaces for velocity and pressure. For instance, let Pk(K) be the set of piecewise polynomial functions of degree less than or CONVERGENCE TO SUITABLE WEAK SOLUTIONS FOR A SUBGRID FEM MODEL 7 equal to kon Kbeing a tetrahedra. Thus the space of continuous, piecewise polynomial functions of degree less than or equal to kon a mesh This denoted as Xh=vh∈C0(Ω) : vh|K∈ Pk(K),∀K∈ Th, We choose the following continuous finite-element spaces Wh=Xh∩H1 0(Ω) and Qh=Xh∩L2 R=0(Ω), for approximating velocity and pressure, respectively. The shape-regular and quasi-uniform properties of {Th}h>0assumed in (H2) suffice to ensure the properties of (H3)(a). We recommend the books [8, Sec. 4.5 ] and [15, Sec. 1.7] for a proof of (6) and (7), Appendix A for a proof of (8), and [17] for a proof of (9). Moreover, the error estimates stated in (H3) make use of (H2) as well (see [26, Lm A.3, Rm 2.1] for a proof). The local approximation properties for the orthogonal projection operators πWhand πQhguarantee hypothesis (H4). The reader is referred to [6]. Remark 3.1. Let pand qbe as in (H4). We know from Sobolev’s embeddings that ˜ Hs 0(Ω) is embedded in Lq′(Ω), where 1 q′+1 q= 1 and s= 3(1 q−1 2); hence Lq(Ω) is embedded in ˜ H−s 0(Ω). Moreover, Hr(R;H)is embedded in Lp′(R;H), where 1 p′+1 p= 1 and r > ¯r=1 p−1 2with Hbeing a Hilbert space; hence Lp(R;H)is embedded in H−r(R;H). Let fbe the extension of foutside [0, T ]as zero. Then, by Hausdorff-Young’s inequality for the Fourier transform, we have kFfkH−r(R;˜ H−s 0(Ω)) ≤CkFfkLp′(R;Lq(Ω)) ≤CkfkLp(R;Lq(Ω)) =CkfkLp(0,T ;Lq(Ω)).(15) Therefore, f∈H−r(0, T ;˜ H−s 0(Ω)) (16) As a reference for further development, it is well to point out, here, the conditions for p,q,sand ¯r: (C) Let s= 3(1 q−1 2)and ¯r=1 p−1 2be defined for pand qas in (H4). 3.2. The discrete problem. Find uh∈H1(0, T ;Wh), ph∈L2(0, T ;Qh) and ˜ uh∈H1(0, T ;˜ Wh) such that, for all (vh,˜ vh, qh)∈Whט Wh×Qh, (∂tuh,vh) + b(uh,uh,vh) + ν(∇uh,∇vh) −(ph,∇ · vh)−b(uh,vh,˜ uh) = (fh,vh),(17a) (uh,∇qh) + (˜ uh,∇qh) = 0,(17b) (∂t˜ uh,˜ vh) + b(uh,uh,˜ vh) +τ−1(˜ uh,˜ vh) + (∇ph,˜ vh) = 0,(17c) uh(0) = u0h,(17d) where τ=1 Csν h2+CckuhkL∞(Ω) h =h2 Csν+CchkuhkL∞(Ω) , with Csand Ccbeing algorithmic positive constants, and fh∈Whis defined by duality as (fh,wh) = hf,whi, for all wh∈Wh. Let us define b(uh,vh,wh) = hN (uh,vh),whi, where N(uh,vh) = (uh· ∇)vh+1 2(∇ · uh)vh. Let {ψi}i=1,...,nube a basis of Whand let {ψi}i=1,...,npbe a basis of Qh, where nuand npdenote the space dimension for Whand Qh, respectively. Thus, one defines ˜ Wh= span{π⊥ Wh(N(φi,φj)), π⊥ Wh(∇φk)}, and W⋆=Wh⊕˜ Wh. Moreover, one defines V⋆={v⋆∈W⋆: (vh,∇qh) + (˜ vh,∇qh) = 0 for all qh∈Qh}. 8 S. BADIA AND J. V. GUTI´ ERREZ-SANTACREU which is a non-conforming approximation space of V. The initialization of the discrete problem can be obtained by the following projection problem: find u0h∈Vh,˜ u0h∈˜ Vhand ξh∈Qhsuch that    (u0h,vh)−(ξh,∇ · vh) = (u0,vh),for all vh∈Vh, (˜ u0h,˜ v) + (∇ξh,˜ v) = (u0,˜ vh),for all ˜ vh∈˜ Vh, (∇ · u0h, qh)−(˜ u0h,∇qh) = 0,for all qh∈Qh. (18) 3.3. Discrete operators. This subsection is devoted to introducing the discrete operators that are used throughtout this paper. Firstly, we will consider a conforming and non-conforming approximation of the Laplace operator −∆ : ˜ H2 0(Ω) →L2(Ω). The non-conforming approximation is based on a stabilizing technique. Consider −∆h:H1 0(Ω) →Whto be the discrete Laplacian operator defined as: −(∆hw,¯ wh) = (∇wh,∇¯ wh) for all ¯ wh∈Wh. The restriction of this operator −∆hto Wh⊂H1 0(Ω) gives a self-adjoint, positive-definite operator. Therefore, we are allowed to define the fractional power of −∆h, say (−∆h)s, for all s∈R, by the Hilbert-Schmidt theorem. The domain of definition of (−∆h)sis D((−∆h)s)≡Whsince dim Wh< ∞. Hence, Ws hmakes reference to Whequipped with the Hilbert norm kwhkWs h= ((−∆h)s 2wh,(−∆h)s 2wh)1 2. The family {Ws h}s∈Ris a scale of Hilbert spaces with respect to the real method of interpolation. Analogously, consider −∆⋆:W⋆→W⋆to be the stabilized discrete Laplacian operator defined as −(∆⋆w⋆,¯ w⋆) = (∇πWhw⋆,∇πWh¯ w⋆) + h−2(π⊥ hw⋆, π⊥ h¯ w⋆) for all ¯ w⋆∈W⋆. It is easy to see that −∆⋆w⋆=−πWh∆⋆w⋆−π⊥ Wh∆⋆w⋆=−∆hπWhw⋆−h−2π⊥ Whw⋆. We have that −∆⋆is self-adjoint and positive-definite. Therefore, we are also allowed to define the fractional power of −∆⋆, say (−∆⋆)s, for all s∈R, by the Hilbert-Schmidt theorem. Thus, Ws ⋆is W⋆equipped with the Hilbert norm kw⋆kWs ⋆= ((−∆⋆)s 2w⋆,(−∆⋆)s 2w⋆). Secondly, we will consider a non-conforming approximation of the Stokes operator A:= P(−∆) : V∩H2(Ω) →Hwhere Pis the Leray-Helmholtz projector operator. Let A⋆:V⋆→V⋆be defined as (A⋆v⋆,¯ v⋆) = (∇πWhv⋆,∇πWh¯ v⋆) + h−2(π⊥ Whv⋆, π⊥ Wh¯ v⋆) for all ¯ v⋆∈V⋆. Equivalently, one can write A⋆=πWhA⋆+π⊥ WhA⋆:= Ah+˜ Ahsatisfying    (Ahv⋆,wh) + (∇rh,wh) = (∇πWhv⋆,∇wh) for all wh∈Wh, (Ahv⋆,∇qh) + ( ˜ Ahv⋆,∇qh) = 0 for all qh∈Qh, (˜ Ahv⋆,˜ wh) + (∇rh,˜ wh) = h−2(π⊥ Whv⋆,˜ wh) for all ˜ wh∈˜ Wh. (19) Again, A⋆is a self-adjoint, positive-definite operator. Therefore, the fractional power of A⋆, say As ⋆, is well-defined for all s∈R. Moreover, Vs ⋆denotes V⋆equipped with the Hilbert norm kv⋆kVs ⋆= (As 2 ⋆v⋆, As 2 ⋆v⋆)1 2. The family {Vs ⋆}s∈Ris a scale of Hilbert space with respect to the real method of interpolation. Next we will consider a non-conforming approximation of the Leray-Helmholtz projection operator. Let P⋆:L2(Ω) →V⋆be defined as (P⋆v,¯ v⋆) = (v,¯ v⋆) for all ¯ v⋆∈V⋆. Equivalently, one can write P⋆=πWhP⋆+π⊥ WhP⋆:= Ph+˜ Phsatisfying    (Phv,wh) + (∇rh,wh) = (πWhv,wh) for all wh∈Wh (Phv,∇qh) + ( ˜ Phv,∇qh) = 0 for all qh∈Qh, (˜ Phv,˜ wh) + (∇rh,˜ wh) = (π⊥ Whv,˜ wh) for all ˜ wh∈˜ Wh. (20) CONVERGENCE TO SUITABLE WEAK SOLUTIONS FOR A SUBGRID FEM MODEL 9 Finally, we define the stabilized Ritz projection operator onto V⋆. Let R⋆:H1 0(Ω) = πWhH1 0⊕ π⊥ WhH1 0(Ω) →V⋆be defined as (∇πWhR⋆v,∇πWhv⋆) + h−2(π⊥ WhR⋆v, π⊥ Wh¯ v⋆) = (∇πWhv,∇πWhv⋆) + h−2(π⊥ Whv,π⊥ Wh¯ v⋆), for all v⋆∈V⋆. Equivalently, one can write R⋆=πWhR⋆+π⊥ WhR⋆:= Rh+˜ Rhsatisfying    (∇Rhv,∇wh) + (∇rh,wh) = (∇πWhv,∇wh) for all wh∈Wh (Rhv,∇qh) + ( ˜ Rhv,∇qh) = 0 for all qh∈Qh, h−2(˜ Rhv,˜ wh) + (∇rh,˜ wh) = h−2(π⊥ Whv,˜ wh) for all ˜ wh∈˜ Wh. (21) 4. Technical preliminary results This section is mainly devote to some technical results concerning equivalence between norms and inf-sup conditions in fractional-order Sobolev spaces. Lemma 4.1. Suppose that conditions (H1)-(H3) hold. Then there exist two positive constants c, C such that, for all s∈R, c(kwhkWs h+h−sk˜ whk)≤ kw⋆kWs ⋆≤C(kwhkWs h+h−sk˜ whk),(22) for all w⋆=wh+˜ wh∈W⋆. Proof. The proof follows by observing that (−∆⋆w⋆)s= (−∆πWhw⋆)s+h−sπ⊥ Wh˜ w⋆for all w⋆.  Corollary 4.2. Suppose that conditions (H1)-(H3) hold. Then there exist two positive constants c, C such that, for all s∈(−3 2,3 2), c(kwhk˜ Hs 0(Ω) +h−sk˜ whk)≤ kw⋆kWs ⋆≤C(kwhk˜ Hs 0(Ω) +h−sk˜ whk),(23) for all w⋆=wh+˜ wh∈W⋆. Proof. The proof is based on the result of [26, Lemma 2.2]: ckwhk˜ Hs 0≤ kwhkWs h≤Ckwhk˜ Hs 0for all wh∈Wh. with s∈(−3 2,3 2).  In the next lemma, we prove the stability of the stabilized discrete Leray-Helmholtz operator P⋆= Ph+˜ Ph. Lemma 4.3. Assume that conditions (H1)-(H3) are satisfied. Then there exists a positive constant C, independent of h, such that, for all s∈[0,1 2), kPhvk˜ Hs 0(Ω) +h−sk˜ Phvk ≤ Ckvk˜ Hs 0(Ω) for all v∈˜ Hs 0(Ω),(24) where P⋆=Ph+˜ Phis the L2(Ω)-orthogonal projection operator onto V⋆. Proof. Let v∈˜ Hs 0(Ω). Then, by the Helmholtz-Hodge decomposition, there exists r∈H1 R=0(Ω) such that v=Pv+∇r, whose variational formulation reads as: (Pv,¯ v) + (∇r, ¯ v) = (v,¯ v) for all ¯ v∈L2(Ω), (v,∇q) = 0 for all H1 R=0(Ω),(25) Note that problem (20) is the stabilized discrete counterpart of (25). From [17, Chapter II, Theorem 1.1], we get kP⋆v−Pvk+k∇rh− ∇rk ≤ C( inf w⋆∈W⋆ kPv−w⋆k+ inf qh∈Qh k∇r− ∇qhk).(26) 16 S. BADIA AND J. V. GUTI´ ERREZ-SANTACREU where µ=2r 2−γ. Integrating over Rand using H¨older’s inequality and Plancherel’s equality gives ZR |ξ|2 2−γ−µkF˜ u⋆k2 V−α ⋆dξ ≤Ck˜ gk 2 2−γ H−r(R;W−s ⋆)k˜ u⋆k 2(1−γ) 2−γ L2(R;V⋆), which implies that ZR |ξ|2βkF˜ u⋆k2 V−α ⋆≤Ck˜ gk 2 2−γ H−r(0,T ;W−s ⋆)k˜ u⋆k 2(1−γ) 2−γ L2(0,T ;V⋆)),(50) for β < ¯ βwith ¯ β:= 1+α 1+s(1 −¯r) coming from the definition of γ,µ, ¯r, and α≤s≤1 + 2α. Next observe that we have, from (15) and (46) for s∈[0,3 2), k˜ gkH−r(R;W−s ⋆)≤C. (51) Inserting (43) and (51) into (50), we arrive at ZR |ξ|2βkF˜ u⋆k2 V−α ⋆≤C. For β≥0, we write ZR (1 + |ξ|)2βkF˜ u⋆k2 V−α ⋆dξ =Z|ξ|≤1 (1 + |ξ|)2βkF˜ u⋆k2 V−α ⋆dξ +Z|ξ|>1 (1 + |ξ|)2βkFu⋆k2 V−α ⋆dξ ≤CZ|ξ|≤1 kF ˜ u⋆k2 V−α ⋆dξ +CZ|ξ|>1 |ξ|2βkF˜ u⋆k2 V−α ⋆dξ. ≤CZR k˜ u⋆k2 V⋆dξ +CZR |ξ|2βkF˜ u⋆k2 V−α ⋆dξ, where Plancherel’s equality and the continuous embedding between V⋆and V−α ⋆were used in the last line. The above estimate also holds trivially for β < 0. Thus we get k∂t˜ u⋆kHβ−1(R;V−α ⋆)+k˜ u⋆kHβ(R;V−α ⋆)≤C. As a result of (33) for s∈[0,1 2), we obtain k∂t˜ u⋆kHβ−1(0,T ;W−α ⋆)+k˜ u⋆kHβ(0,T ;W−α ⋆)≤C, for all αsatisfying 0 ≤α≤s≤1 + 2α < 2, and for all βsatisfying β < ¯ β:= 1+α 1+s(s 2+1 4). This latter inequality leads to (47). Next, multiply (49) by A1−s ⋆F˜ u⋆and take the real part to get νkA⋆F˜ u⋆k2 V−s ⋆≤CkF ˜ gkW−s ⋆kA1−s ⋆F˜ u⋆kVs ⋆≤CkF ˜ gkW−s ⋆kA⋆F˜ u⋆kVs ⋆, where we have applied (33) for s∈[0,2). Thus, (1 + |ξ|)−2rkA⋆F˜ u⋆k2 V−s ⋆≤C(1 + |ξ|)−2rkF˜ gk2 W−s ⋆, and hence kA⋆˜ u⋆kH−r(R;V−s ⋆)≤Ck˜ gkH−r(R;W−s ⋆). It is not hard to see from (33) that k∆⋆v⋆kW−s ⋆≤CkA⋆v⋆kfor all v⋆∈V⋆and all s∈[0,3 2). Then, using (51) yields k∆⋆˜ u⋆kH−r(0,T ;W−s ⋆)≤C, for all r > ¯r, which implies (48).  Corollary 5.7. Suppose that assumptions (H1)-(H4) hold. For α∈[1 4,1 2)and β < ¯ β=2 5(1 + α), it follows that k∂tu⋆kHβ−1(0,T ;W−α ⋆)+ku⋆kHβ(0,T ;W−α ⋆)≤C, (52) where C > 0is a constante independent of h. CONVERGENCE TO SUITABLE WEAK SOLUTIONS FOR A SUBGRID FEM MODEL 17 Proof. From s∈[1 2,3 2) and 0 ≤α≤s≤1 + 2α < 2, we obtain α∈[1 4,1 2). Next note that 1 1+s(s 2+1 4) reaches its maximum 2 5at s=3 2. Therefore we can simplify the expression ¯ βin term of αonly as ¯ β=2 5(1 + α) in (47).  Using (23), one can also prove the following. Corollary 5.8. Assume that assumptions (H1)-(H4) hold. Then, for α∈[1 4,1 2)and β < ¯ β=2 5(1+α), it follows that there exists a constant C > 0, independent of h, such that k∂tuhkHβ−1(0,T ;˜ H−α 0(Ω)) +kuhkHβ(0,T ;˜ H−α 0(Ω)) ≤C. (53) Furthermore, for s∈[1 2,3 2]and rsuch that r > ¯r=3 4−s 2=1 p−1 2, it follows that k∆huhkH−r(0,T ;˜ H−s 0(Ω)) ≤C. (54) We now proceed to obtain an estimate for ph. Lemma 5.9. Suppose that conditions (H1)-(H4) hold. There exists a constant C > 0, independent of h, such that, for s∈[1 2,7 10 ]and r > ¯r=3 4−s 2, kphkH−r(0,T ;H1−s(Ω)) ≤C, (55) where C > 0is a constant independent of h. Proof. First we write (52) as k∂tu⋆kH−r(0,T ;W−α ⋆)≤C, where α∈[1 4,1 2) and r > ˜r:= 1 −¯ β=3 5−2 5α. As a result, we have that k∂tu⋆kH−r(0,T ;W−s ⋆)≤C(56) holds for α≤sand ˜r≤¯rprovided that s∈[1 2,7 10 ] and r > ¯r=3 4−s 2. From (45), we bound kphkH1−s(Ω) ≤sup v⋆∈W⋆\{0} (∇ph,v⋆) kv⋆kWs ⋆ ≤C(k∂tu⋆kW−s ⋆+k∆⋆u⋆kW−s ⋆+kN⋆(u⋆,u⋆)kW−s ⋆+kfhkW−s ⋆) ≤C(k∂tu⋆kW−α ⋆+k∆⋆u⋆kW−s ⋆+kN⋆(u⋆,u⋆)kW−s ⋆+kfhk˜ H−s 0(Ω)). The proof is completed via (56), (48), (46) and (16).  6. Convergence towards weak and suitable weak solutions In this section we will prove that the sequence of the approximate solutions provided by scheme (17) converges towards a weak solution in the sense of Definition 2.3 and towards a suitable weak solution in the sense of Definition 2.4. In order for these convergence results to hold, we will need to use the following compactness results `a la Aubin-Lions. The following compactness result is due to Lions [33]. Lemma 6.1. Let H0֒→H ֒→H1be three Hilbert spaces with dense and continuous embedding. Assume that the embedding H0֒→His compact. Then L2(0, T ;H0)∩Hγ(0, T ;H1)embeds compactly in L2(0, T ;H)for γ > 0. The proof of the two following compactness result can be found in [23, Ap. A.1, A.2]. Lemma 6.2. Let X ֒→Ybe two Hilbert spaces with compact embedding. Then Hβ(0, T ;X)embeds continuously and compactly in C0([0, T ]; Y)for β > 1 2. Lemma 6.3. Let H0֒→H1be two Hilbert spaces with compact embedding. Let γ > 0and γ > µ, then the injection Hγ(0, T ;H0)֒→Hµ(0, T H1)is compact. 18 S. BADIA AND J. V. GUTI´ ERREZ-SANTACREU Theorem 6.4. Assume that hypotheses (H1)-(H4) are satisfied. Then there exists a subsequence (denoted in the same way) of approximate solutions (uh, ph)converging toward a weak solution given in Definition (2.3) in the following sense as h→0: uh→uin L2(0, T ;H1 0(Ω)) −weak and in L2(0, T ;Hβ(Ω)) −strong for all β < 1 (57) and ph→pin H−r(0, T ;Hδ(Ω)) −weak for all δ∈[3 10,1 2]and r > 1 4+δ 2.(58) Proof. Let v∈Hr(0, T ;˜ Hs 0(Ω)), for s∈(1 2,7 10 ] and r > 3 4−s 2, and q∈L2(0, T ;H1 R=0(Ω)). From (11) and (12), we are allowed to construct three sequences {vh}h>0⊂Hr(0, T ;Wh), {˜ vh}h>0⊂ Hr(0, T ;˜ Wh) and {qh}h>0⊂L2(0, T ;Qh) such that vh→vin Hr(0, T ;˜ Hs 0(Ω))-strong, ˜ vh→0in L2(0, T ;L2(Ω))-strong and qh→qin L2(0, T ;H1 R=0(Ω))-strong as h→0. By virtue of (43), (53), (54) and (55), we know that there exist a subsequence of {vh}h>0and {ph}h>0, still denoted by itself, and a pair (u, p) such that uh→uin L∞(0, T ;L2(Ω)) −weak-⋆, uh→uin L2(0, T ;H1 0(Ω)) −weak, ∂tuh→∂tuin H−r(0, T ;˜ H−s 0(Ω)) −weak, ∆huh→∆uin H−r(0, T ;˜ H−s 0(Ω)) −weak, and ∇ph→ ∇pin H−r(0, T ;˜ H−s 0(Ω)) −weak, for all s∈(1 2,7 10 ] and r > 3 4−s 2. Observe that we have used that the fact that ˜ H−s(Ω) coincides with H−s 0(Ω) for s∈(1 4,7 10 ] for the pressure. We also have, from (43), that ˜ uh→0 in L2(0, T ;L2(Ω)) −strong,(59) since ν1 2 hk˜ uhkL2(0,T ;L2(Ω)) ≤ kτ−1 2˜ uhkL2(0,T ;L2(Ω)) ≤C. We can pass to the limit in (17b). Thus we find that ∇ · u= 0 in (0, T )×Ω, whence u∈ L∞(0, T ;H)∩L2(0, T ;V). For the trilinear terms, we proceed as follows. By Lemma 6.1, we have that uh→uin L2(0, T ;Hβ(Ω)) −strong for all β < 1, since {uh}h>0is bounded in L2(0, T ;H1 0(Ω))∩Hβ((0, T ); ˜ H−α 0(Ω)) for α∈[1 4,1 2) and 0 < β < 2 5(1+α) from (43) and (53). Therefore, N(uh,uh)→ N(u,u) in D′((0, T )×Ω). As a consequence of (46), we obtain N(uh,uh)→ N(u,u) in H−r(0, T ;˜ H−s 0(Ω)). On passing to the limit in (17b), we have had that ∇ · u= 0 in (0, T )×Ω, thereby N(uh,uh)→(u· ∇)uin H−r(0, T ;˜ H−s 0(Ω)). By an analogous argument, we find that ˜ N(uh,˜ uh)→0in H−r(0, T ;˜ H−s 0(Ω)), where h˜ N(uh,˜ wh),vhi=b(uh,vh,˜ wh) for all uh,wh∈Whand ˜ vh∈˜ Wh. From the above convergences, it is easy to see that ∂tuh+N(uh,uh)−ν∆huh+∇ph−˜ N(uh,˜ uh)−fh→∂tu+ (uh· ∇)uh−ν∆u+∇p−f. in H−r(0, T ;˜ H−s 0(Ω)) as h→0. CONVERGENCE TO SUITABLE WEAK SOLUTIONS FOR A SUBGRID FEM MODEL 19 For the initial condition, we have that uh→uin C0([0, T ]; ˜ H−α 0(Ω))-strong for α∈(1 4,1 2) by Lemma 6.2; therefore, uh(0) →u(0) in ˜ H−α 0(Ω). Furthermore, it follows from (18) and (27) that u0h→u0in ˜ H−α 0(Ω). We have thus shown that u(0) = u0. The energy inequality can be verified by the lower semicontinuity of the norm for the weak topology; for complete details, see [4].  Theorem 6.5. Under hypotheses (H1)-(H4), the sequence of approximate solutions (uh, ph)converges, up to a subsequence, to a suitable weak solution given in Definition 2.4 as h→0. Proof. Let φ∈ D((0, T )×Ω; R+) and substitute vh=πWh(uhφ) into (17a) to get ZT 0 {(∂tuh, πWh(uhφ)) + b(uh,uh, πWh(uhφ)) + ν(∇uh,∇πWh(uhφ)) −(ph,∇ · πWh(uhφ)) −b(uh, πWh(uhφ),˜ uh)−(fh, πWh(uhφ))}dt= 0. (60) We are ready to take the limit in (60) as h→0 so as to prove that the weak solution (u, p) found in Theorem 6.4 is suitable. We will only focus on passing to the limit in the terms of (60) involving the subscale velocity ˜ uhand the pressure term. The remaining terms appear in a rudimentary finite element formulation so that a proof can be found in [25]. In particular, from (13) and (9) and in virtue of Lemma 6.3, it follows that lim h→0ZT 0 (∂tuh,uhφ) dt=−1 2ZT 0 (|u|2, ∂tφ) dt, lim h→0ZT 0 b(uh,uh, πWh(uhφ)) dt=−1 2ZT 0 (|u|2u,∇φ)dt, lim inf h→0ZT 0 ν(∇uh,∇πWh(uhφ)) dt≥ZT 0 (|∇u|2, φ) dt−ZT 0 (1 2|u|2,∆φ) dt, and lim h→0−ZT 0 hfh,uhφidt=−ZT 0 hf,uφidt. To begin with, we first turn our attention to passing to the limit in the convective term. b(uh, πWh(uhφ),˜ uh) dt=b(uh,uhφ, ˜ uh) dt+b(uh, πWh(uhφ)−uhφ, ˜ uh) dt =b(uh,uh,˜ uhφ) + (uh· ∇φuh,˜ uh) + b(uh, π⊥ Wh(uhφ),˜ uh) = (π⊥ Wh(N(uh,uh)φ),˜ uh) + (π⊥ Wh(uh· ∇φuh),˜ uh) +b(uh, π⊥ Wh(uhφ),˜ uh). (61) From (13) and (6), we have: ZT 0 (π⊥ Wh(N(uh,uh)φ),˜ uh) dt≤ZT 0 kπ⊥ Wh(N(uh,uh)φ),˜ uh)kk˜ uhkdt≤CZT 0 hkuhkL∞(Ω)k∇uhkk˜ uhkdt ≤C ZT 0 τh2kuhk2 L∞(Ω)k∇uhk2dt!1 2 ZT 0 τ−1k˜ uhk2dt!1 2 ≤Ch3 4kuhk1 2 L∞(0,T ;L2(Ω))kuhkL2(0,T ;H1 0(Ω))k˜ uhkτ−1 2L2(0,T ;L2(Ω)) and hence lim h→0ZT 0 (π⊥ Wh(N(uh,uh)φ),˜ uh) dt= 0. Analogously, we bound ZT 0 (uh· ∇φ, uh·˜ uh) dt≤CT h3 4kuhkL∞(0,T ;L2(Ω))k˜ uhkτ−1 2L(0,T ;L2(Ω)) 20 S. BADIA AND J. V. GUTI´ ERREZ-SANTACREU and ZT 0 b(uh, πWh(uhφ)−uhφ, ˜ uh) dt≤Ch3 4kuhk1 2 L∞(0,T ;L2(Ω))kuhkL2(0,T ;H1 0(Ω))k˜ uhkτ−1 2L2(0,T ;L2(Ω)). Thus lim h→0ZT 0 (uh· ∇φ, uh·˜ uh) dt= 0, and lim h→0ZT 0 b(uh, πWh(uhφ)−uhφ, ˜ uh)dt = 0. For the “viscous” term, it is not hard to see that lim inf h→0ZT 0 τ−1(|˜ uh|2, φ) dt≥0. For the pressure terms, we write ZT 0 (ph,∇·πWh(uhφ)) dt=ZT 0 (phuh,∇φ)dt+ZT 0 (ph,∇·(πWh(uhφ)−(uhφ))) dt+ZT 0 (φph,∇·uh) dt. It was proved in [23] that lim h→0ZT 0 (phuh,∇φ) dt=ZT 0 (pu, ∇φ) dt and lim h→0=ZT 0 (ph,∇ · (πWh(uhφ)−(uhφ))) dt= 0. For the remaining pressure terms, we use (17b) with qh=πQh(φph) to obtain ZT 0 (φph,∇ · uh) dt+ZT 0 (∇ph,˜ uhφ) = ZT 0 (phφ−πQh(phφ),∇ · uh) dt +ZT 0 (∇(phφ)− ∇πQh(phφ),˜ uh) dt−ZT 0 (ph∇φ, ˜ uh) dt. We know from [23] that lim h→0ZT 0 (phφ−πQh(phφ),∇ · uh) dt= 0. Let ε > 0 and set s=1 2+16 9εand ¯r=3 4−s 2=1 4−4 9ε. Now choose r=1 2−4 9ε. Moreover, set α=1 4−5 9εand ¯ β=2 5(1 + α) = 1 5−2 9ε. Thus we have 1 −s > α and ¯ β > r since 1−s=1 2−16 9ε > 1 2(1 2−16 9ε) = 1 4−4 9ε > 1 4−5 9ε=α and r=1 2−4 9ε < 1 2−2 9ε=2 5(5 4−5 9ε) = 2 5(1 + 1 4−5 9ε) = 2 5(1 + α) = ¯ β. From the a priori energy estimates (52) and (55) and the commutator property (14), our choice of parameters yields ZT 0 (∇(pφ)− ∇πQh(phφ),˜ uh) dt≤ k∇(pφ)− ∇πQh(phφ)kH−r(0,T ;L2(Ω))k˜ uhkHr(0,T ;L2(Ω)) ≤Ch1−s−αkphkH−r(0,T ;H1−s(Ω))k˜ uhkhαHr(0,T ;L2(Ω)) and hence lim h→0ZT 0 (∇(φ)− ∇πQh(phφ),˜ uh) dt= 0. CONVERGENCE TO SUITABLE WEAK SOLUTIONS FOR A SUBGRID FEM MODEL 21 Finally, it is easy to see in a similar fashion that lim h→0ZT 0 (φph,∇ · uh) dt= lim h→ZT 0 (φph−πWh(φph),∇ · uh) dt= 0.  Appendix A. Proof of the inverse inequalities (8) To prove inequalities (8), we follow very closely the arguments developed in [8, Thm. 4.5.11]. We first need to introduce an equivalent norm for fractional order Hilbert spaces as follows. Let s∈(0,1). Then kuk2 Hs(Ω) =kuk2+|u|2 Hs(Ω), where |u|2 Hs(Ω) =ZΩZΩ |u(x)−u(y)|2 |x−y|3+2sdxdy. Given (K, P,Σ), we define ( ˜ K, ˜ P, ˜ Σ) where ˆ K={(1/hK)x:x∈K}. Thus, if uhis a function defined on K, then ˆuhis defined on ˜ Kby ˆu(ˆ x) = u(h−1 Kx) for all ˆ x∈ˆ K. Thus we can write k∇uhkL2(K)=h1 2 Kkˆ ∇ˆuhkL2(K). As ˆ ∇ˆuhbelongs to a space of finite and fixed dimension on ˆ K, on which all norms are equivalent, it is not hard to see that there is a constant Cˆ T>0 such that kˆ ∇ˆuhkL2(K)≤Cˆ T|ˆuh|Hs(ˆ K). Reverting to K, this leads to kˆ ∇ˆukL2≤Cˆ Th−3 2+s K|uh|Hs(K) and hence k∇uhkL2(K)≤Cˆ Th−1+s KkuhkHs(Ω). An argument in the proof of [8, Prop. 4.4.11 ] shows that if ( ˜ K, ˜ P,˜ Σ) is a referent element, we have that there exists a constant C˜ T>0 such that Cˆ T≤C˜ T. Summing over all elements Kand using the quasi-uniformity of the mesh leads to k∇uhk ≤ Ch−1+s X K∈Th kuhk2 Hs(Ω)!. Then (8) follows because the sum of the fractional norms over all elements is smaller than the fractional norm over the union of the elements. Acknowledgment The authors are very grateful to Professor Vivette Girault who provided a proof of a particular case of inequality (8). 22 S. BADIA AND J. V. GUTI´ ERREZ-SANTACREU References [1] Adams, R. A.; Fournier, J. J. F. Sobolev spaces, volume 140 of Pure and Applied Mathematics (Amsterdam). Elsevier/Academic Press, Amsterdam, second edition, 2003. [2] Badia, S. On stabilized finite element methods based on the Scott-Zhang projector. Circumventing the inf-sup condition for the Stokes problem. Computer Methods in Applied Mechanics and Engineering, 247-248(0):65–72, 2012. [3] Badia, S.; Codina, R.; Guti´ errez-Santacreu, J. V. Long-term stability estimates and existence of a global attractor in a finite element approximation of the Navier-Stokes equations with numerical subgrid scale modeling. SIAM J. Numer. Anal., 48(3):1013–1037, 2010. [4] Badia, S; Guti´ errez-Santacreu, J. V. Convergence towards weak solutions of the Navier-Stokes equations for a finite element approximation with numerical subgrid-scale modelling. IMA J. Numer. Anal., 34(3):1193–1221, 2014. [5] Becker, R.; Braack, M. A finite element pressure gradient stabilization for the Stokes equations based on local projections. Calcolo, 38(4):173–199, 2001. [6] Bertoluzza, S. The discrete commutator property of approximation spaces. C. R. Acad. Sci. Paris S´er. I Math., 329(12):1097–1102, 1999. [7] Boris, J. P.; Grinstein, F. F.; Oran, E. S.; Kolbe, R. L. New insights into large eddy simulation. Fluid Dynamics Research, 10(4–6):199, 1992. [8] Brenner, S. C.; Scott, L. R. The mathematical theory of finite element methods, volume 15 of Texts in Applied Mathematics. Springer, New York, third edition, 2008. [9] Caffarelli, L.; Kohn, R.; Nirenberg, L. Partial regularity of suitable weak solutions of the Navier-Stokes equations. Comm. Pure Appl. Math., 35(6):771–831, 1982. [10] Codina, R. Analysis of a stabilized finite element approximation of the Oseen equations using orthogonal subscales. Applied Numerical Mathematics, 58(3):264–283, 2008. [11] Codina, R.; Blasco, J. Stabilized finite element method for the transient Navier-Stokes equations based on a pressure gradient projection. Computer Methods in Applied Mechanics and Engineering, 182:277–300, 2000. [12] Codina, R.; Principe, J.; Guasch, O. Badia, S. Time dependent subscales in the stabilized finite element approximation of incompressible flow problems. Computer Methods in Applied Mechanics and Engineering, 196(21– 24):2413–2430, 2007. [13] Colom´ es, O; Badia, S.; Codina, R.; Principe, J. Assessment of variational multiscale models for the large eddy simulation of turbulent incompressible flows. Computer Methods in Applied Mechanics and Engineering, 285:32–63, March 2015. [14] Dauge, M. Stationary Stokes and Navier-Stokes systems on twoor three-dimensional domains with corners. I. Linearized equations. SIAM J. Math. Anal., 20(1):74–97, 1989. [15] Ern, A.; and Guermond, J.-L. Theory and practice of finite elements, volume 159 of Applied Mathematical Sciences. Springer-Verlag, New York, 2004. [16] Girault, V.; Raviart, P.-A. Finite element approximation of the Navier-Stokes equations, volume 749 of Lecture Notes in Mathematics. Springer-Verlag, Berlin-New York, 1979. [17] Girault, V.; Raviart, P.-A. Finite element methods for Navier-Stokes equations, volume 5 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 1986. Theory and algorithms. [18] Grinstein, F. F.; Margollin, L. G.; Rider, W. J. Implicit large eddy simulation: computing turbulent fluid dynamics. Cambridge university press, 2007. [19] Grisvard, P. Elliptic problems in nonsmooth domains, volume 24 of Monographs and Studies in Mathematics. Pitman (Advanced Publishing Program), Boston, MA, 1985. [20] Guasch, O.; Codina, R. Statistical behavior of the orthogonal subgrid scale stabilization terms in the finite element large eddy simulation of turbulent flows. Computer Methods in Applied Mechanics and Engineering, 261–262:154– 166, July 2013. [21] Guermond, J.-L. Stabilization of Galerkin approximations of transport equations by subgrid modeling. ESAIM: Mathematical Modelling and Numerical Analysis - Mod´elisation Math´ematique et Analyse Num´erique, 33(6):1293– 1316, 1999. [22] Guermond, J.-L. Finite-element-based Faedo–Galerkin weak solutions to the Navier–Stokes equations in the threedimensional torus are suitable. Journal de Math´ematiques Pures et Appliqu´ees, 85(3):451–464, March 2006. [23] Guermond, J.-L. Faedo-Galerkin weak solutions of the Navier-Stokes equations with Dirichlet boundary conditions are suitable. J. Math. Pures Appl. (9), 88(1):87–106, 2007. [24] Guermond, J.-L. On the use of the notion of suitable weak solutions in CFD. International Journal for Numerical Methods in Fluids, 57(9):1153–1170, July 2008. [25] Guermond, J.-L. The LBB condition in fractional Sobolev spaces and applications. IMA J. Numer. Anal., 29(3):790–805, 2009. [26] Guermond, J.-L.; Pasciak, J. E. Stability of discrete Stokes operators in fractional Sobolev spaces. J. Math. Fluid Mech., 10(4):588–610, 2008. [27] Hopf, E. ¨ Uber die Anfangswertaufgabe f¨ur die hydrodynamischen Grundgleichungen. Math. Nachr., 4:213–231, 1951. CONVERGENCE TO SUITABLE WEAK SOLUTIONS FOR A SUBGRID FEM MODEL 23 [28] Hughes, T. J. R.; Feij´ oo, G. R.; Mazzei, L.; Quincy, J.-B. The variational multiscale method - A paradigm for computational mechanics. Computer Methods in Applied Mechanics and Engineering, 166(1–2):3–24, 1998. [29] Hughes, T. J. R.; Mazzei, L.; Jansen, K. E. Large eddy simulation and the variational multiscale method. Computing and Visualization in Science, 3:47–59, 2000. [30] Kellogg, R. B.; Osborn, J. E. A regularity result for the Stokes problem in a convex polygon. J. Functional Analysis, 21(4):397–431, 1976. [31] Leray, J. Sur le mouvement d’un liquide visqueux emplissant l’espace. Acta Math., 63(1):193–248, 1934. [32] Lin, F. H. A new proof of the Caffarelli-Kohn-Nirenberg theorem. Comm. Pure Appl. Math., 51(3):241–257, 1998. [33] Lions, J.-L. Quelques m´ethodes de r´esolution des probl`emes aux limites non lin´eaires. Dunod; Gauthier-Villars, Paris, 1969. [34] Lions, J.-L.; Magenes, E. Probl`emes aux limites non homog`enes et applications. Vol. 1. Travaux et Recherches Math´ematiques, No. 17. Dunod, Paris, 1968. [35] Lions, J.-L.; Magenes, E. Non-homogeneous boundary value problems and applications. Vol. I. Springer-Verlag, New York-Heidelberg, 1972. Translated from the French by P. Kenneth, Die Grundlehren der mathematischen Wissenschaften, Band 181. [36] Scheffer, V. Hausdorff measure and the Navier-Stokes equations. Comm. Math. Phys., 55(2):97–112, 1977. [37] Simon, J. Nonhomogeneous viscous incompressible fluids: existence of velocity, density, and pressure. SIAM J. Math. Anal., 21(5):1093–1117, 1990. [38] Temam, R. Navier-Stokes equations. AMS Chelsea Publishing, Providence, RI, 2001. Theory and numerical analysis, Reprint of the 1984 edition.