On the convergence order of the finite element error in the kinetic energy for high Reynolds number incompressible flows

Available online at www.sciencedirect.com ScienceDirect Comput. Methods Appl. Mech. Engrg. 385 (2021) 114032 www.elsevier.com/locate/cma On the convergence order of the finite element error in the kinetic energy for high Reynolds number incompressible flows Bosco García-Archillaa, Volker Johnb,c, Julia Novod,∗ aDepartamento de Matemática Aplicada II, Universidad de Sevilla, Sevilla, Spain bWeierstrass Institute for Applied Analysis and Stochastics (WIAS), Mohrenstr. 39, 10117 Berlin, Germany cFreie Universität Berlin, Department of Mathematics and Computer Science, Arnimallee 6, 14195 Berlin, Germany dDepartamento de Matemáticas, Universidad Autónoma de Madrid, Spain Received 23 December 2020; received in revised form 30 June 2021; accepted 1 July 2021 Available online 24 July 2021 Abstract The kinetic energy of a flow is proportional to the square of the L2(Ω) norm of the velocity. Given a sufficient regular velocity field and a velocity finite element space with polynomials of degree r, then the best approximation error in L2(Ω) is of order r+1. In this survey, the available finite element error analysis for the velocity error in L∞(0,T;L2(Ω)) is reviewed, where Tis a final time. Since in practice the case of small viscosity coefficients or dominant convection is of particular interest, which may result in turbulent flows, robust error estimates are considered, i.e., estimates where the constant in the error bound does not depend on inverse powers of the viscosity coefficient. Methods for which robust estimates can be derived enable stable flow simulations for small viscosity coefficients on comparatively coarse grids, which is often the situation encountered in practice. To introduce stabilization techniques for the convection-dominated regime and tools used in the error analysis, evolutionary linear convection–diffusion equations are studied at the beginning. The main part of this survey considers robust finite element methods for the incompressible Navier–Stokes equations of order r−1, r, and r+1/2 for the velocity error in L∞(0,T;L2(Ω)). All these methods are discussed in detail. In particular, a sketch of the proof for the error bound is given that explains the estimate of important terms which determine finally the order of convergence. Among them, there are methods for inf–sup stable pairs of finite element spaces as well as for pressure-stabilized discretizations. Numerical studies support the analytic results for several of these methods. In addition, methods are surveyed that behave in a robust way but for which only a non-robust error analysis is available. The conclusion of this survey is that the problem of whether or not there is a robust method with optimal convergence order for the kinetic energy is still open. c 2021TheAuthor(s).PublishedbyElsevier B.V.ThisisanopenaccessarticleundertheCCBY-NC-NDlicense (http://creativecommons.org/licenses/by-nc-nd/4.0/). MSC: 65M60; 76D05 Keywords: Incompressible Navier–Stokes equations; Convection–diffusion equations; Convection-dominated regime; Finite element methods; Convergence of the error of the kinetic energy; Robust error bounds ∗Corresponding author. E-mail addresses: [email protected] (B. Garc´ ıa-Archilla), [email protected] (V. John), [email protected] (J. Novo). https://doi.org/10.1016/j.cma.2021.114032 0045-7825/ c 2021 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 1. Introduction Incompressible flows are modeled by considering the conservation of linear momentum and the conservation of mass, leading to the incompressible Navier–Stokes equations, given here already in dimensionless form, ∂tu−ν∆u+(u·∇)u+∇p=fin (0,T]×Ω, ∇ ·u=0 in (0,T]×Ω, u=0in [0,T]×∂Ω, u(0,·)=u0(·) in Ω, (1) where Ωis a bounded domain in Rd,d∈ {2,3}, with polyhedral and Lipschitz continuous boundary ∂Ω,T∈Ris a final time, uis the velocity field, pthe pressure, ν > 0 the kinematic viscosity coefficient, u0a given initial velocity, and frepresents the external body accelerations acting on the fluid. Eq. (1) is already given as an equation without physical dimensions, as it is used in numerical analysis and simulations. Consider now the underlying physical problem. Let ρ[kg/m3] be the density of the fluid, µ[kg/m s] be its dynamic viscosity and define a characteristic length scale L[m] and a characteristic velocity scale U[m/s]. Then, the dimensionless number Re =ρU L µ is called Reynolds number. Utilizing the characteristic time scale L/U[s], then one finds the relation ν=Re−1. Thus, a high Reynolds number flow leads to a small coefficient νin (1). An important physical quantity of incompressible flows is the contained energy, which is the sum of the kinetic energy 1 2∥u(t)∥2 0at time t∈(0,T], where ∥·∥0denotes the norm in L2(Ω), and the dissipative energy. It is known that a distributional solution exists, satisfying the energy inequality 1 2∥u(t)∥2 0+ν∥∇u∥2 L2((0,t),L2(Ω)) ≤1 2∥u(0)∥2 0+∫t 0∫Ω f·udxdt, for all t∈(0,T], e.g., see [1, Chapter V, Theorem 3.1.1] or [2, Lemma 7.21]. In the present paper, we will assume sufficient regularity so that weak solutions are indeed strong solutions and, hence, unique. More precisely, Leray–Hopf solutions are considered, e.g., see [3]. The unique weak solution in two dimensions satisfies even an energy equality. The kinetic energy of a flow belongs to the quantities of interest in many simulations. This paper surveys available analytic results for the order of convergence of the finite element error in the kinetic energy, more precisely, of the convergence of the velocity in L∞(L2), where Lp(Ls):= Lp(0,T;Ls(Ω)) =⎧ ⎨ ⎩ v:(∫T 0(∫Ω|v|sdx)p/s dt)1/p <∞,p,s∈[1,∞)⎫ ⎬ ⎭ , with the application to each component of a vector-valued function and the usual replacement of the integrals with the essential supremum if p= ∞ or s= ∞. Since in applications often the case of small viscosity coefficients is of importance, which leads eventually to turbulent flows, especially error estimates are considered where the constant of the error bound does not explicitly depend on inverse powers of ν, such that it does not blow up for ν→0. The case where νis small (relative to the product of a typical velocity scale and a typical length) is often referred to as convection-dominated, since, if ∇uand ∆uhave norms of similar sizes, the first order term (u·∇)udominates the second order term −ν∆u. Our aim is to investigate finite element approximations of the Navier–Stokes equations that are suitable for the convection-dominated regime. This requires to analyze how the different terms in the Navier–Stokes equations affect their finite element approximation at high Reynolds numbers. To this end, we begin this survey with a simpler case, time-dependent linear convection–diffusion–reaction equations, in order to study a linear convective term. Then, we consider the incompressibility condition and finally we also add the nonlinear term of the Navier–Stokes equations. It will be analyzed how the different terms affect the rate of convergence of finite element methods for a regular solution and ν→0. Typically, in the error bounds, it appears a product of an error constant and of norms of the solution in Sobolev spaces. Our error analysis aims for error bounds with constants that are independent of inverse powers of ν. Of course, one must be aware that in practical situations the norms of the solutions will grow as ν→0, but error constants that do not explode as ν→0 can be considered as the minimal requirement or 2 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 necessary condition for a robust performance of a method in practical convection-dominated problems. Hence, this survey provides both, some insight into which methods are expected to fail in practice, because they do not meet the minimum requirement, and in particular a survey of those methods which are expected to be robust. Finally, from the theoretical point of view, the question of deriving local error bounds for the methods is still open. Since for realistic problems, a high global regularity of the solution cannot be expected, such local error estimates would be more appropriate. That means, to prove bounds for the methods in a subset of the domain where the solution is smooth, and to get optimal bounds for the error depending on norms of the solution in this subdomain. To the best of our knowledge, there is almost no error analysis in this direction and the techniques needed probably have not been developed yet. For evolutionary problems, the only references we know are [4] and [5], where local error bounds for the SUPG method applied to linear evolutionary convection–reaction–diffusion equations are obtained. Accordingly, we first consider the following time-dependent linear convection–diffusion–reaction problem ∂tu−ν∆u+b·∇u+cu =fin (0,T]×Ω, u=0 on [0,T]×∂Ω, u(0,x)=u0(x) in Ω, (2) where b(t,x), c(t,x), and f(t,x) are given functions, ν > 0 is a constant diffusion coefficient, and u0(x) is the initial data. Problem (2) describes the passive transport of a scalar quantity, like temperature or concentration, by molecular motion, modeled with the diffusive term −ν∆u, and by a convective transport via a flow field b. One speaks of a convection-dominated regime if ν≪ ∥b∥L∞(L∞)L, where Lis a characteristic length scale. To provide an idea of a typical range for diffusion and convection, in the temperature balance of the crystallization process studied in [6], the diffusion coefficient is around 1.5·10−7m2/sand the maximal velocity is of order 0.1m/s. If L=1 m is used for deriving non-dimensionalized equations, the convection in these equations is stronger than the diffusion by around six orders of magnitude. A common feature of the solutions of both the Navier–Stokes equations (1) and the linear convection–diffusion– reaction equation (2) in the convection-dominated regime is the appearance of (very) small structures. The smallest physically important structures or scales of turbulent incompressible flows are the so-called Kolmogorov scales whose size, away from walls, is proportional to Re−3/4L. The solution of (2) typically possesses layers, i.e., regions of a small width of O(ν) to O(ν1/2) with very steep slopes of the gradient. Standard spatial discretizations of (1) and (2), like central finite differences or the Galerkin finite element method, try to predict the behavior of all important scales of the solution. However, in the convection-dominated regime, the smallest scales are usually much smaller than the affordable mesh width, such that most of the small scales even cannot be represented on such grids. Of course, scales that cannot be represented cannot be simulated. Standard discretizations cannot cope with this situation and they lead usually either to a blow up of the simulations or to meaningless numerical solutions that are globally polluted with spurious oscillations, compare also the numerical results presented in Section 4.1.3. The remedy consists in using so-called stabilized discretizations, which are usually denoted as turbulence models in the case of turbulent flows. Such discretizations introduce in some sophisticated way numerical (artificial) diffusion/viscosity in the discrete equations. A mile stone in this direction was the proposal of the Streamline- Upwind Petrov–Galerkin (SUPG) method in [7,8] for steady-state convection–diffusion–reaction equations. The first finite element error analysis for the steady-state as well the time-dependent convection–diffusion equations can be found in [4,9]. The analysis was extended to the steady-state Navier–Stokes equations in [10,11] and the timedependent equations in [12,13]. The SUPG method is the core of a variational multiscale turbulence model proposed in [14]. For other stabilized discretizations, like the Continuous Interior Penalty (CIP) method or Local Projection Stabilization (LPS) methods, one can observe a similar development: they were first proposed for the steadystate convection–diffusion–reaction equation and then extended to other equations, including (2) and (1). Since the same ideas can be utilized for discretizing scalar convection–diffusion–reaction equations and incompressible flow equations, these classes of equations are sometimes jointly presented in the literature, like in the most often cited monograph in the field [15]. For completeness, it should be noted that, besides the purely mathematical approaches for turbulence modeling mentioned above, there are many turbulence models whose derivation is based on physical insight in turbulent flows, e.g., see [16], which do not have a counterpart for the scalar equation (2). Often, a numerical analysis for such turbulence models is not available. Finite element convergence analysis of partial differential equations usually assumes that the solution of these equations is sufficiently regular, i.e., that it belongs to certain Sobolev spaces. Then, the goal consists in deriving 3 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 so-called optimal estimates, where the order of convergence equals the order of the best approximation error in the corresponding norm. The techniques for proving such estimates differ for different norms. For convection-dominated problems, there is an additional important aspect: one has to track the dependency of the error bound on the data of the problem, in particular on the coefficient ν. The goal is to obtain error bounds where the constant does not contain inverse powers of νand therefore it does not blow up as ν→0. Such estimates are called robust or semi-robust. The latter notion is inspired from the fact that the error bounds contain norms of the solution of the partial differential equation in Sobolev spaces, which are usually of higher order with respect to the spatial variable, these norms might depend implicitly on ν, and they might grow as ν→0. Another aspect for convection-dominated problems is that the convergence of the method, in the sense that the mesh width tends to zero, is not of interest, because if the mesh is sufficiently fine such that all important scales are resolved, a stabilization is not longer necessary. Instead, the behavior on coarse grids is of importance, especially in practice. To distinguish these aspects, one speaks often of ‘order of error reduction’ or ‘order of error decay’ instead of ‘order of convergence’ in the latter situation. The main topic of this survey is the provable order of error reduction of the velocity error, where with respect to the spatial variable the space L2(Ω)dis considered. Given a finite element velocity space with continuous piecewise polynomials of order rand assuming for the velocity solution of (1) that u(t)∈Hr+1(Ω)d, then the best approximation error in L2(Ω)dis of order r+1. Since for the kinetic energy error, applying the Cauchy–Schwarz inequality, one has ∥u(t)∥2 0−∥uh(t)∥2 0≤ ∥u(t)−uh(t)∥0∥u(t)+uh(t)∥0, this error is proportional to the L2(Ω) norm of the velocity error and the optimal order of error reduction for the kinetic energy is (r+1). The finite element error analysis of time-dependent problems leads usually to estimates where a sum of errors in different norms is finally bounded. One of these terms, emerging from the temporal derivative, is the velocity error in L∞(L2). To facilitate the comparison with the literature, also in this survey the error analysis will be presented in terms of the velocity error in L∞(L2) instead of the kinetic energy. One of our motivations for compiling this survey is an open problem stated recently in [17]: It is open whether optimal and semi-robust L∞(0,T;L2(Ω)d)error bounds for the velocity can be proved for some method. The present paper provides a survey on the state of the art of error estimates in L∞(L2)1for the convectiondominated regime of problems (1) and (2). The available analysis was revised for this survey such that the formulas are dimensionally correct, i.e., if all quantities in a formula are equipped with physical units, then the operations in this formula are well defined. In this way, the numerical analysis becomes physically consistent. The complete error analysis of finite element discretizations of (1) and (2) needs usually several pages of estimates. In this survey, we restrict the presentation only to the key steps of the analysis, which lead to the dominating term in the error bounds. With this presentation, also important differences of the analysis between different methods are illustrated. Further aspects of the realization of the methods, e.g., fully implicit approaches vs. implicit–explicit (IMEX) methods, solution as a coupled velocity–pressure problem vs. projection methods, etc. will not be discussed here to focus the presentation on the main topic. Several results from the finite element error analysis are supported with numerical examples. Here, the expectations are illustrated for an academic problem with smooth solution, because only with such a setup, it is possible to support the rates of convergence predicted by the analysis. Numerical simulations for realistic problems, i.e., turbulent flows, possess usually quantities of interest that are not covered by the numerical analysis and such simulations are outside the scope of this survey. Section 2sets up weak formulations of (1) and (2), introduces some notations, and provides some tools that are used in the numerical analysis. Discretizations of the convection–diffusion–reaction equation (2) are studied in Section 3. The bulk of the paper, Section 4, is devoted to finite element discretizations of the incompressible Navier–Stokes equations (1). The current state of the art is summarized in Section 5. 2. Weak formulations and notations Throughout the paper, standard notations will be used for Lebesgue and Sobolev spaces. The dimension is not indicated in symbols of time–space function spaces and the corresponding norms. 1When surveying the literature, error bounds will be discussed that were obtained for fully discrete problems. Strictly speaking, these bounds are valid only in the discrete time instants and not in the whole time interval. For not unnecessarily complicating the notation, we will use L∞(L2) also in the case of fully discrete problems. 4 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 For the convection–diffusion–reaction equation (2), it will be assumed that band care sufficiently smooth functions with respect to xand that µ0≤µ(t,x)=(c−1 2∇ ·b)(t,x)∀(t,x)∈[0,T]×Ω(3) is satisfied. Let V=H1 0(Ω). A variational form of (2) reads as follows: Given f∈H−1(Ω), find u:(0,T]→V such that (∂tu, v)+ν(∇u,∇v)+(b·∇u, v)+(cu, v)= ⟨f, v⟩ ∀ v∈V,(4) and u(0,x)=u0(x). Here, (·,·) denotes the inner product in L2(Ω)dand ⟨·,·⟩ the dual pairing of H1 0(Ω)dand H−1(Ω)d,d∈ {1,2,3}. Inner products in L2(ω) with ω= Ωwill be denoted with (·,·)ωand a similar convention is used for the notation of norms. Using for the incompressible Navier–Stokes equations (1) the function spaces V=H1 0(Ω)d,Q=L2 0(Ω)={q∈L2(Ω):(q,1) =0}, a weak formulation of problem (1) is as follows: Given f∈H−1(Ω)d, find (u,p):(0,T]→V×Qsuch that for all (v,q)∈V×Q, (∂tu,v)+ν(∇u,∇v)+((u·∇)u,v)−(∇ ·v,p)+(∇ ·u,q)= ⟨f,v⟩,(5) and u(0,·)=u0(·)∈V. An important subspace of Vis the space of weakly divergence-free functions Vdiv ={v∈V:(∇ ·v,q)=0∀q∈Q}.(6) Let {Th}be a family of triangulations of Ω, where h>0 indicates the fineness of Th. The mesh cells of Th are denoted by Kwith hKbeing the diameter of Kand h=maxK∈ThhK. It will be assumed that the meshes are conforming and quasi-uniform. In practice, anisotropic grids are of considerable importance. However, the analysis of finite element discretizations on such grids has many open questions, in particular there are no contributions with respect to the topic of this survey. In numerical simulations, the infinite-dimensional spaces are replaced by finite-dimensional ones, e.g., Vby Vh,r, Vby Vh,r, and Qby Qh,lwhere r,l∈Ndenote the degrees of the local finite element polynomials. Most of the available finite element analysis is for the case of conforming finite element spaces, i.e., the finite element spaces are subspaces of the corresponding infinite-dimensional spaces, e.g., Vh,r⊂V. For the Navier–Stokes equations, the space of discretely divergence-free functions is defined by Vdiv h={vh∈Vh,r:(∇ ·vh,qh)=0∀qh∈Qh,l}.(7) Clearly, there are less conditions for a function to be contained in Vdiv hthan to belong to Vdiv. In fact, for many pairs Vh,r/Qh,lof finite element spaces, it turns out that Vdiv h⊂ Vdiv. The following inverse inequality holds for each vh∈Vh,r, e.g., see [18, Theorem 3.2.6], ∥vh∥Wm,q(K)≤cinvh l−m−d(1 q′−1 q) K∥vh∥Wl,q′(K),(8) where 0 ≤l≤m≤1, 1 ≤q′≤q≤ ∞, and ∥·∥Wm,q(K)is the norm in Wm,q(K). We will denote by πhu∈Vh,rthe elliptic projection defined by (∇πhu,∇ϕh)=(∇u,∇ϕh)∀ϕh∈Vh,r.(9) The following bound holds for any u∈H1 0(Ω)∩Hr+1(Ω), e.g., see [18], ∥u−πhu∥0+h∥u−πhu∥1≤Chr+1∥u∥r+1.(10) Generic constants that do not depend neither on the mesh width hnor on the diffusion/viscosity coefficient ν will be denoted by C. To avoid technicalities that are not significant for the topic of this survey, the actual regularity assumptions on the weak solution will not be stated explicitly in the formulation of the theorems, but it will be assumed that the solution is sufficiently regular such that all terms appearing in the error bounds are well defined. 5 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 3. Convection–diffusion–reaction equations The first part of this section discusses briefly the Galerkin finite element discretization of the convection– diffusion–reaction equation (4). If the best approximation error in L2(Ω) is r+1, then an error bound in L∞(L2) with constants independent of inverse powers of the diffusion coefficient of order rcan be proved. Stabilized methods, with the emphasis on the SUPG method, are the topic of the second part of the section. For these methods, robust estimates of order r+1/2 can be derived. 3.1. The Galerkin method The time-continuous finite element problem aims to find a function uh∈Vh,rthat fulfills a problem of form (4) for all test functions from Vh,r. More precisely, uh:(0,T]→Vh,ris called Galerkin approximation if (∂tuh, vh)+ν(∇uh,∇vh)+(b·∇uh, vh)+(cuh, vh)= ⟨f, vh⟩ ∀ vh∈Vh,r,(11) with uh(0,x) being an appropriate approximation of u0(x). It is well known that the standard Galerkin method leads to numerical solutions with non-physical spurious oscillations in the convection-dominated regime when the solution of (4) possesses layers. However, if the solution of (4) does not exhibit layers and if it is sufficiently smooth, the Galerkin method produces meaningful approximations that tend to the solution of (4) with an order of error reduction rin L∞(L2), independently of the value of ν, which will be proved in Theorem 3.1 and illustrated numerically in Example 3.2. In most works on convection-dominated problems, the case µ0>0, the constant in (3), is assumed; see [15]. In Theorem 3.1, we consider both cases µ0>0 and µ0≤0. Theorem 3.1 (Error Estimate for the Galerkin Method).Let u be the sufficiently smooth solution of (2) and uhbe the finite element approximation defined in (11). Assume b∈L∞(L∞)and c ∈L∞(L∞), then the following bound holds for 0<t≤T ∥u(t)−uh(t)∥0≤e−µ0t∥uh(0) −πhu0∥2 0+ChrKu⏐⏐1−e−µ0t⏐⏐ |µ0|+Chr+1∥u∥L∞(Hr+1),(12) the bound being also valid for µ0=0if ⏐⏐1−e−µ0t⏐⏐/|µ0|is replaced by t. The constant is Ku= ∥∂tu∥L∞(Hr)+(∥b∥L∞(L∞)+h∥c∥L∞(L∞))∥u∥L∞(Hr+1).(13) Proof. The error is split into e=u−uh=(u−πhu)−(uh−πhu)=η−eh, where πhuis the elliptic projection defined in (9). Subtracting (11) from (2) gives the error equation (∂teh, vh)+ν(∇eh,∇vh)+(b·∇eh, vh)+(ceh, vh)=(τ1+τ2, vh),(14) with the truncation errors τ1=(∂tu−πh∂tu)+cηand τ2=b·∇η. Here, the commutation of elliptic projection and temporal differentiation ∂t(πhu)=πh∂tuwas used. Choosing vh=ehin the error equation (14) and integrating by parts yields 1 2 d dt ∥eh∥2 0+µ0∥eh∥2 0≤(τ1+τ2,eh).(15) Because of using the elliptic projection, the right-hand side of (15) does not depend on ν. The non-negative term with νon the left-hand side was neglected. Altogether, estimate (15) does not contain νsuch that the arising error bound will also not depend on ν. Applying the Cauchy–Schwarz inequality gives (τ1+τ2,eh)≤(∥τ1∥0+∥τ2∥0)∥eh∥0, and then, arguing as in [19, Lemma 3.1] leads to ∥eh∥0 d dt ∥eh∥0+µ0∥eh∥2 0=1 2 d dt ∥eh∥2 0+µ0∥eh∥2 0≤(∥τ1∥0+∥τ2∥0)∥eh∥0, 6 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 which yields d dt ∥eh∥0+µ0∥eh∥0≤ ∥τ1∥0+∥τ2∥0. Multiplying by eµ0tand integrating with respect to time, one obtains ∥eh(t)∥0≤e−µ0t∥eh(0)∥0+∫t 0 e−µ0(t−s)(∥τ1(s)∥0+∥τ2(s)∥0)ds ≤e−µ0t∥eh(0)∥0+(1 −e−µ0t) µ0 max 0≤s≤T(∥τ1(s)∥0+∥τ2(s)∥0),(16) where, in the last term, the factor (1−e−µ0t)/µ0must be replaced by tif µ0=0. Utilizing (10) and the assumptions on the regularity of the coefficient functions leads to max 0≤s≤T(∥τ1(s)∥0+∥τ2(s)∥0)≤ChrKu,(17) where Kuis the constant in (13). Inserting (17) in (16) and applying the triangle inequality together with (10) gives (12).□ A similar result to that of Theorem 3.1 can be found in [20], where the error of the time derivative ∥∂tu−∂tuh∥0 is proved to be O(hr) with a constant independent of inverse powers of ν, assuming sufficient regularity for the solution. Example 3.2 (Galerkin Method).Let Ω=(0,1)2and consider the convection–diffusion–reaction equation (2) with the prescribed solution u(x,y)=sin(πt) sin(2πx) sin(2πy), the convection field b=(1 −y,−2/3+x)T, the reaction field c=1, the final time T=2, and Dirichlet boundary conditions on ∂Ω. Fig. 1 displays the initial grid (level 1). This grid was refined uniformly using the standard regular refinement by dividing each triangle into four triangles. Piecewise linear P1and piecewise quadratic P2finite elements were used, see Table 1 for information on the numbers of degrees of freedom. To mimic the continuous-in-time situation, a second order time stepping scheme, the Crank–Nicolson scheme, was applied with the very small time step ∆t=10−5. We checked that almost the same results were obtained for ∆t=10−4, such that the temporal error can be considered to be negligible. If not mentioned otherwise, the linear systems of equations in all examples of this paper were solved with the sparse direct solver UMFPACK [21]. The computational results obtained with the Galerkin method for four small values of the diffusion coefficient are presented in Fig. 2. First of all, one can see that the errors in L∞(L2) are small, which indicates that for this smooth problem even the standard Galerkin method computes meaningful numerical solutions. For P1finite elements, the optimal second order of convergence can be observed in Fig. 2. In contrast, for P2finite elements, the error reduction on coarse grids is not optimal. In the cases ν=10−6and ν=10−8, it is only of second order, exactly as predicted by the analysis. For ν=10−4, it can be already seen that the order of convergence improves to third order when the grids become sufficiently fine. Finally, if the diffusion is sufficiently large, ν=10−2, the optimal order is observed already on very coarse grids. Since similar observations were made in all other numerical studies for convection–diffusion–reaction equations and the emphasis of this survey is on the convection-dominated regime, only numerical results for ν≤10−4will be presented in the remainder of this section. □ 3.2. Stabilized methods of order r +1/2in L∞(L2) As mentioned before, the Galerkin method produces numerical solutions in the convection-dominated regime that are globally polluted with spurious oscillations if the solution of (4) possesses layers. In this situation, one has to introduce a stabilizing mechanism in the discretization. In the ideal case, this mechanism removes all spurious oscillations and gives numerical solutions with sharp layers. However, numerical assessments in [22,23] show that there are only very few finite element methods with these properties, which are methods that use a so-called algebraic stabilization. Most stabilized finite element methods localize the spurious oscillations in a vicinity of the layers and reduce them considerably compared with the Galerkin method. 7 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 Fig. 1. Initial grid (level 1). Fig. 2. Example 3.2 Galerkin method for the convection–diffusion–reaction equation. Table 1 Numerical examples for convection–diffusion–reaction equation, number of degrees of freedom (including Dirichlet nodes). Level P1P2 3 41 145 4 145 545 5 545 2113 6 2113 8321 7 8321 33025 8 33025 131585 3.2.1. The SUPG method The most popular method of this kind is the streamline-upwind Petrov–Galerkin (SUPG) method. This residualbased stabilization adds artificial diffusion along the streamlines of the solution. It reads as follows: Given f∈ L2(Ω), find uh:(0,T]→Vh,rsuch that (∂tuh, vh)+aSUPG(uh, vh)+∑ K∈Th δK(∂tuh,b·∇vh)K=(f, vh)+∑ K∈Th δK(f,b·∇vh)K∀vh∈Vh,r,(18) with uh(0,x) being an appropriate approximation of u0(x) and aSUPG(uh, vh)=ν(∇uh,∇vh)+(b·∇uh, vh)+(cuh, vh)+∑ K∈Th δK(−ε∆uh+b·∇uh+cuh,b·∇vh)K. 8 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 Here, {δK}are local parameters that have to be chosen appropriately. Let uhave the physical unit unit(u), then the term (∂tuh, vh) possesses the unit unit(u)2md/s, where the unit md comes from integration and the term 1/sfrom the differentiation with respect to time. For all other terms in (18), one finds the same physical unit unit(u)2md/sby taking into account the units for the diffusion coefficient m2/s, the convection field m/s, and the reaction field 1/s. The stabilization term of the SUPG method has the physical unit unit(δK)unit(u)2md/s2. In order that adding this term to the other terms is correct from the dimensional point of view, it has to have the same unit as the other terms. Consequently, it follows that unit(δK)=s. Thus, the stabilization parameter is a time scale. The numerical analysis for a stabilized method is performed usually in an associated norm. For uh∈Vh,r, this norm for the SUPG method is as follows ∥uh∥SUPG := ⎛ ⎝ν∥∇uh∥2 0+∑ K∈Th δK∥b·∇uh∥2 0,K+∥µ1/2uh∥2 0⎞ ⎠ 1/2 . The SUPG method for evolutionary convection–diffusion equations is analyzed in [24]. A main topic in [24], for a full discretization in time and space, is the discussion whether the stabilization parameters {δK}should depend on the length of the time step or on the mesh width. The latter situation can be studied by considering the continuous-in-time method, because this method does not contain a time step. Assume that b(t,x)=b(x),∇ ·b(x)=0,c(t,x)=c(x),(19) µ0>0, the mesh is uniform with mesh width h, and the stabilization parameters are the same for all mesh cells, i.e., δK=δ. Also, considered only the convection-dominated regime, i.e., it is assumed that ν≤∥b∥L∞h. The stabilization parameter is defined to be δ=min {h 4cinv∥b∥L∞min {1 2,µ0 4∥c∥L∞ ,µ1/2 0 ∥c∥1/2 L∞ ,∥b∥L∞h 4νcinv },1 µ0 ,1 ∥c∥L∞}.(20) Note that the physical unit of δis seconds. Formula (20) for the stabilization parameter can be obtained with a slight revision of the analysis of [24] by carefully taking into account the dimensional issue. Then, the following error estimate is derived in [24, Theorem 5.2]. Theorem 3.3 (Error Estimate for the SUPG Method).Let t ≤T<∞, let ν < ∥b∥L∞h, let the solution of (2) be sufficiently regular, and let (19) and (20) be satisfied. Then, there exists a positive constant Mudepending on ∥u0∥r+1,∥u∥L∞(Hr+1),∥∂tu(0)∥L∞(Hr+1),∥∂tt u∥L2(Hr+1), and the coefficients of the problem, but not on inverse powers of ν, such that the following error estimate holds: ∥(u−uh)(t)∥0+∥u−uh∥L2(SUPG) ≤Muhr+1/2.(21) Proof (Sketch).The proof of (21) in [24] covers several pages. Here, only a sketch is given that concentrates on explaining the reason why estimate (21) leads to the order r+1/2 in L∞(L2). As usual, the error is decomposed into an interpolation or projection error and a discrete remainder. The analysis in [24] defines for this decomposition a function Πhu(t) which is the solution of a formally steady-state problem. The idea is to consider the solution u(t) of the evolutionary problem (2) as the solution of a steady-state problem with right-hand side f(t)−∂tu(t). Then, given t∈(0,T], Πhu(t)∈Vh,rsatisfies aSUPG (Πhu(t), vh)=aSUPG (u(t), vh)∀vh∈Vh,r,(22) 9 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 Using shas test function in (32) gives ∥∇sh∥0≤ ∥∇u∥0. Then, one can conclude with (34) and the Riesz–Thorin interpolation theorem that ∥∇sh∥Lp≤C∥∇u∥Lpfor p∈(2,∞).(35) Assuming the necessary smoothness in time, a Stokes projection like (32) can be defined for ∂tu. Then, an error bound of form (33) can be derived also for ∂t(u−sh). In this section, the Galerkin method for the Oseen equations (29) without any added stabilization is analyzed. Let uh:(0,T]→Vh,r,ph:(0,T]→Qh,r−1be the Galerkin approximation satisfying (∂tuh,vh)+ν(∇uh,vh)+((b·∇)uh,vh)−(∇ ·vh,ph)= ⟨f,vh⟩ ∀ vh∈Vh,r, (∇ ·uh,qh)=0∀qh∈Qh,r−1.(36) Since the pair Vh,r/Qh,r−1satisfies the discrete inf–sup condition (27), this problem is well posed. The following theorem states an error bound for the Galerkin approximation to the velocity in L∞(L2) in which the constant in the error bound does not depend on inverse powers of the viscosity. Theorem 4.1 (Error Estimate for the Galerkin Method).Let (u,p)be the sufficiently smooth solution of (29) and assume b∈L∞(L∞). Let uhbe the velocity finite element approximation defined in (36). Then, the following bound holds for 0<t≤T ∥u(t)−uh(t)∥0≤ ∥uh(0) −sh(0)∥0+Cthr−1Mu,p+Chr+1∥u∥L∞(Hr+1), where Mu,p=h(∥∂tu∥L∞(Hr)+∥b∥L∞(L∞)∥u∥L∞(Hr+1))+∥p∥L∞(Hr).(37) Proof. The proof is similar to the proof of Theorem 3.1. The error is split into an interpolation error between uand its modified Stokes projection shand the discrete remainder φh=uh−sh∈Vdiv h. A straightforward calculation, taking into account the definition of the Stokes projection, leads to the error equation (∂tφh,vh)+ν(∇φh,∇vh)+((b·∇)φh,vh)=(τ1,vh)+(∇ ·vh, τ2)∀vh∈Vdiv h,(38) where τ1=(∂tu−∂tsh)+((b·∇)(u−sh)), τ2=p−Πr−1 hp. Taking vh=φhin (38), neglecting a non-negative term on the left-hand side, and using the skew-symmetry (30), which gives in particular ((b·∇)φh,φh)=0, yields ∥φh∥0 d dt ∥φh∥0≤(τ1,φh)+(∇ ·φh, τ2).(39) This inequality does not contain the viscosity νany more, hence the resulting estimate will not depend explicitly on ν. The key step for estimating the last term consists in utilizing integration by parts (∇ ·φh, τ2)=(∇ ·φh,p−Πr−1 hp)= −(φh,∇(p−Πr−1 hp)).(40) Now, the estimate of the right-hand side of (39) continues by applying the Cauchy–Schwarz inequality, (33), and (31) leading to (τ1,φh)+(∇ ·φh, τ2)≤C Mu,phr−1∥φh∥0, where Mu,pis the constant in (37). To conclude proof, one argues as in the proof of Theorem 3.1, applies the triangle inequality and again (33).□ In view of Theorem 4.1, one observes that for the Galerkin approximation to the velocity in the Oseen equations, the suboptimal order of convergence r−1 is proved with constants independent of inverse powers of νwhile, as Theorem 3.1 shows, for the convection–reaction–diffusion model an order of convergence rcan be proved. The reduction of one order comes from the error bound of the second term in (39), which is the truncation error of the pressure, see (40). Alternatively, an inverse estimate could be used, which leads finally to the same error bound. 16 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 Table 2 Numerical examples for the Stokes and Navier–Stokes equations, number of degrees of freedom (including Dirichlet nodes) for inf–sup stable pairs of finite element spaces on triangular grids. P2/P1P2/Pdisk 1, barycentric P3/P2 Level Velocity Pressure Velocity Pressure Velocity Pressure 2 210 144 170 41 3 290 41 802 576 626 145 4 1090 145 3138 2304 2402 545 5 4226 545 12418 9216 9410 2113 6 16642 2113 49410 36864 37250 8321 7 66050 8321 Fig. 5. Example 4.2 Galerkin method for the evolutionary Stokes equations, P2/P1(left) and P3/P2(right). A first way to improve the situation is the use of inf–sup stable pairs of finite element spaces where the term (∇ · φh, τ2) does not appear. This is the case if discretely divergence-free functions are weakly divergence-free, see Section 4.2.1. Another way consists in adding appropriate stabilization terms to the Galerkin method, see Section 4.2.2. Example 4.2 (Galerkin Discretization of the Evolutionary Stokes Equations).The order of convergence predicted by Theorem 4.1 is determined by the pressure term in Mu,p, the last term in (37). As already mentioned, the term (∇ · φh, τ2), a coupling term of the finite element velocity and pressure that does not vanish if the finite element velocity solution is not weakly divergence-free, is responsible for the appearance of the pressure in Mu,p. Dominating convection does not play any role in this respect. To support this observation, a numerical study for the evolutionary Stokes equations, i.e., for (28) with b=0, will be presented here. The considered example is defined in Ω=(0,1)2, the final time is T=2, and the prescribed solution is given by u(t;x,y)=2πsin(πt)(sin2(πx) sin(πy) cos(πy) −sin2(πy) sin(πx) cos(πx)),p(t,x,y)=20 sin(πt)(x2y−1 6).(41) The same spatial meshes and the same temporal discretization as in Example 3.2 were used. For the spatial discretization, the Taylor–Hood pairs P2/P1and P3/P2of finite element spaces were utilized. The numbers of degrees of freedom are given in Table 2. A quadrature rule was applied that is exact for polynomials of degree eight on triangles. This order of the quadrature rule is not necessary for assembling the matrices but it keeps the quadrature error small for the calculation of the right-hand side and the L2(Ω) error of the velocity. Simulations were performed for ν∈ {10−2,10−4,10−6,10−8}. Fig. 5 presents the numerical results. There are two important observations. First, if the mesh is sufficiently fine with respect to the viscosity coefficient, one can see numerically the optimal order of convergence, as predicted by the classical finite element convergence theory. But second, in the other situation, there is a clear decrease of 17 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 the order of error reduction for ν≤10−4. Note that coarse meshes with respect to the viscosity coefficient are the typical situation for high Reynolds number flows. In this case, the order of error reduction is described obviously better by robust error bounds. From the results for the P3/P2pair of spaces, one can conclude that the order r−1 is a sharp prediction. □ 4.1.3. On the impact of the discretization of the nonlinear term In this section, the Galerkin approximation to the Navier–Stokes equations (5) is studied. Again, the situation Vdiv h⊂ Vdiv will be considered and for fixing ideas, as in previous section, the analysis will be presented for the inf–sup stable Taylor–Hood pair of mixed finite element spaces. Given f∈L2(Ω)d, let uh:(0,T]→Vh,r,ph:(0,T]→Qh,r−1be the solution of the following Galerkin approximation of (5) (∂tuh,vh)+ν(∇uh,∇vh)+((uh·∇)uh,vh)−(∇ ·vh,ph)=(f,vh)∀vh∈Vh,r, (∇ ·uh,qh)=0∀qh∈Qh,r−1.(42) The nonlinear term in (42) is called convective form of the nonlinear term in [2, Chapter 6]. It is not skew-symmetric for Taylor–Hood pairs of finite element spaces. For (uh,ph) satisfying (42), one cannot even prove boundedness of the approximations, since in the stability analysis, the term ((uh·∇)uh,uh) cannot be bounded. Next, in addition to (42), a Galerkin approximation with a skew-symmetric bilinear form of the nonlinear term is considered: Find uh:(0,T]→Vh,r,ph:(0,T]→Qh,r−1such that (∂tuh,vh)+ν(∇uh,∇vh)+b(uh,uh,vh)−(∇ ·vh,ph)=(f,vh)∀vh∈Vh,r, (∇ ·uh,qh)=0∀qh∈Qh,r−1,(43) where, the so-called divergence form of the nonlinear term is given by b(u,v,w)=((u·∇)v,w)+1 2((∇ ·u)v,w)∀u,v,w∈V.(44) Notice the well-known property b(u,v,w)= −b(u,w,v)∀u,v,w∈V,(45) such that, in particular, b(u,v,v)=0 for all u,v∈V. There are several ways for estimating the convective term (44) in the error analysis, depending on the smoothness assumptions on its arguments. Assuming somewhat more than the minimal regularity, it holds that b(u,v,w)≤ ∥u∥0∥∇v∥L∞∥w∥0+1 2∥∇ ·u∥0∥v∥L∞∥w∥0∀u,w∈V,v∈W1,∞(Ω)d.(46) A standard stability analysis proves that the velocity approximation using the skew-symmetric form (44) is bounded in L∞(L2) with a constant that does not depend on inverse powers of ν. Theorem 4.3 (Stability for Skew-Symmetric Form of the Nonlinear Term).Let uhbe the velocity approximation defined in (43). Then, the following stability estimate is valid for all t ∈[0,T] ∥uh(t)∥0≤ ∥uh(0)∥0+∫t 0∥f(s)∥0ds.(47) Proof. From the theory of the Galerkin method developed in [3], it is known that (43) has a unique velocity solution uh. Taking vh=uhin (43), using the skew-symmetric property (45) and the Cauchy–Schwarz inequality gives ∥uh∥0 d dt ∥uh∥0+ν∥∇uh∥2 0=1 2 d dt ∥uh∥2 0+ν∥∇uh∥2 0≤ ∥f∥0∥uh∥0. Neglecting the second term on the left-hand side to get rid of ν, dividing by ∥uh∥0and integrating in (0,T) lead to (47).□ Despite the stability, it is not possible to derive error bounds with constants that are independent of inverse powers of the viscosity. The important terms of the error equation are 1 2 d dt ∥φh∥2 0+ν∥∇φh∥2 0= −b(φh,sh,φh)−b(uh,φh,φh)+··· = −b(φh,sh,φh)+··· ,(48) 18 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 Fig. 6. Example 4.4 Galerkin method for the Navier–Stokes equations, ν=10−6, level 5, solution at t=1.3409, left: convective form of the convective term (42), right: divergence form of the convective term (44); note the different scalings. where φh=uh−sh, with shbeing the Stokes projection defined in (32). The second term on the right-hand side of (48) vanishes due to (45). Then, using the estimate (46) for the convective term, Young’s inequality, and bounding the norm of the divergence by the norm of the gradient yields 1 2 d dt ∥φh∥2 0+ν∥∇φh∥2 0≤ ∥∇sh∥L∞∥φh∥2 0+1 2∥sh∥L∞∥∇ ·φh∥0∥φh∥0+··· ≤ν 2∥∇φh∥2 0+(∥∇sh∥L∞+1 8ν∥sh∥2 L∞)∥φh∥2 0+··· .(49) Hence, the application of Gronwall’s lemma leads to the appearance of ν−1in the exponential factor. Note that with minimal regularity assumptions on u, one even gets ν−3in the exponential factor, e.g., see [2, Theorem 7.35]. An alternative estimate that can be considered uses the inverse inequality (8) ∥∇ ·φh∥0≤ ∥∇φh∥0≤cinv h∥φh∥0, which, however, leads finally to an inverse power of the mesh width in the exponential factor of the error bound. Example 4.4 (Galerkin Method: Impact of the Discretization of the Nonlinear Term: Convective and Divergence Form).For this example, in principle the same configuration is considered as in Example 4.2. Only, the right-hand side of the problem is different since the prescribed solution (41) is now inserted in the evolutionary Navier–Stokes equations (1) instead of the evolutionary Stokes equations. In each time instant, a nonlinear problem has to be solved. For this purpose, a standard fixed point iteration (Picard iteration) was employed. This iteration was stopped if the Euclidean norm of the residual vector was smaller than 10−9. Simulations were performed again on levels 3 to 7 and for ν∈ {10−4,10−6,10−8}. Using the convective form of the nonlinear term (42), usually the numerical simulations blew up. Only if the spatial resolution was sufficiently fine with respect to the viscosity, there was no blow-up: here for ν=10−4on levels 5–7 and for ν=10−6on level 7. In contrast, we could not observe any blow-up if the divergence form of the convective term (43) was applied. Fig. 6 presents a comparison of numerical solutions obtained with both discretizations of the convective term, shortly before the blow-up of the simulation with the convective form of the convective term. Even if we could not observe a blow-up if the divergence form was used, the numerical solution is globally polluted with spurious oscillations. The errors of the velocity in L∞(L2) for the simulations with the divergence form are presented in Fig. 7. One can observe that the errors are quite large, of order 1, which fits the theory since both ∥u(t)∥0and ∥uh(t)∥0are bounded for a finite time interval, Theorem 4.3, so that also the maximal L2(Ω) error is bounded. Apart of this fact, the Galerkin method with the divergence form of the nonlinear term cannot be used without stabilization in the convection-dominated regime on coarse grids, in view of Figs. 6 and 7, where again coarse is to be understood with respect to the size of the viscosity coefficient. □ 19 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 Fig. 7. Example 4.4 Galerkin method for the Navier–Stokes equations, P2/P1, divergence form of the convective term. Recently, a new discrete form of the nonlinear term has been proposed in [52], the so-called EMAC form c(u,v,w)=(2D(u)v,w)+((∇ ·u)v,w), where D(u)=(∇u+(∇u)T)/2 is the velocity deformation tensor. The EMAC form is the only discrete form that conserves kinetic energy (for ν=0,f=0), momentum (for fwith zero linear momentum), and angular momentum (for fwith zero angular momentum) in the case that uis not weakly divergence-free. Using the EMAC form requires a modification of the pressure. In [53], the Galerkin discretization of the Navier–Stokes equations (5) with the EMAC form of the nonlinear term is analyzed. It is shown that under the assumption of a sufficiently smooth solution, the nonlinear term can be bounded in a robust way. Altogether, a non-robust estimate of order r is presented in [53], where the dependency on ν−1appears in the pressure term of the error bound. However, it is possible to derive a robust estimate, but only of order r−1. The numerical study in Example 4.6 will demonstrate that this estimate is sharp. Theorem 4.5 (Error Estimate for the Galerkin Method with EMAC Form of the Nonlinear Term).Let the solution of the Navier–Stokes equations (5) be sufficiently smooth and define L(T)=C∫T 0∥∇u(s)∥L∞ds, Mu,p=C(∫T 0(h∥∂tu(s)∥r+h(∥u(s)∥1∥u(s)∥2)1/2∥u(s)∥r+1+∥p(s)∥r)ds). Then, the following error estimate holds for the Galerkin method with EMAC form of the nonlinear term ∥u(t)−uh(t)∥0≤eL(T)(∥uh(0) −sh(0)∥0+Mu,phr−1)+Chr+1∥u∥L∞(Hr+1), with shbeing the Stokes projection. The constant Mu,pis dimensionally correct in three dimensions. For Mu,pto be dimensionally correct in two-dimensional problems, the term (∥u(s)∥1∥∥u(s)∥2)1/2in Mu,pmust be replaced by (h∥u(s)∥1∥u(s)∥2)1/2+(∥u(s)∥0∥u(s)∥2)1/2.(50) Proof (Sketch).The error is decomposed in the form u−uh=(u−sh)−(uh−sh)=η−φh, where shis the Stokes projection defined in (32). Then, deriving the error equation in the usual way and using φhas test function lead to ∥φh∥0 d dt ∥φh∥0≤⏐⏐c(u,u,φh)−c(uh,uh,φh)⏐⏐+|(∇ ·φh,P−qh)|+··· (51) for arbitrary qh∈Qh,r−1. The estimate of the difference of the convective terms in (51) follows [53]. A direct calculation yields c(u,u,φh)−c(uh,uh,φh) 20 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 =c(η,u,φh)+c(sh,η,φh)−c(sh,φh,φh)−c(φh,sh,φh)−c(φh,φh,φh).(52) The last term on the right-hand side vanishes due to a property of the EMAC form. Now, a part of the fourth term on the right-hand side is reformulated. Using the definition of the deformation tensor, integration by parts, the skew symmetry of the divergence form of the nonlinear term, and again the definition of the deformation tensor yields 2(D(φh)sh,φh) =((∇φh)sh,φh)+(sh,(∇φh)φh) =((sh·∇)φh,φh)+((φh·∇)φh,sh) = −((∇ ·sh)φh,φh)−((sh·∇)φh,φh)−((∇ ·φh)φh,sh)−((φh·∇)sh,φh) = −1 2((∇ ·sh)φh,φh)−((∇ ·φh)φh,sh)−((φh·∇)sh,φh) = −1 2((∇ ·sh)φh,φh)−((∇ ·φh)φh,sh)−(D(sh)φh,φh). Inserting this expression in (52) and collecting terms give c(u,u,φh)−c(uh,uh,φh)=c(sh,η,φh)+c(η,u,φh)−1 2((∇ ·sh)φh,φh)−(D(sh)φh,φh).(53) Applying H¨ older’s inequality, the stability (35) of the Stokes projection, a Sobolev embedding, and an interpolation in Sobolev spaces yields for the first term c(sh,η,φh)=(2D(sh)η,φh)+((∇ ·sh)η,φh) ≤2∥D(sh)∥L3∥η∥L6∥φh∥0+∥∇ ·sh∥L3∥η∥L6∥φh∥0 ≤C∥∇sh∥L3∥η∥L6∥φh∥0≤C∥∇u∥L3∥η∥1∥φh∥0≤C∥∇u∥1/2∥η∥1∥φh∥0 ≤C(∥u∥1∥u∥2)1/2∥η∥1∥φh∥0. To estimate the second term, H¨ older’s and Agmon’s inequality are used to get c(η,u,φh)≤C∥∇η∥0∥u∥L∞∥φh∥0≤C(∥u∥1∥u∥2)1/2∥η∥1∥φh∥0. The third and the fourth term in (53) are bounded just with H¨ older’s inequality. One observes that there is no inverse power of the viscosity in the error bound of the discrete convective term. The arguments above are dimensionally correct in three dimensions. In two dimensions, both ∥·∥L3and ∥·∥L6 should be replaced by ∥·∥L4, and arguing similarly, one concludes that the term (∥u(s)∥1∥∥u(s)∥2)1/2in Mu,p should be replaced by (50). The pressure term on the right-hand side of (51) can be bounded robustly in the same way as in the proof of Theorem 4.1, which gives a term of order r−1. □ Example 4.6 (Galerkin Method: EMAC Form of the Nonlinear Term).The same problem as in Example 4.4 is considered, with the same setup. Only the EMAC form was used as discrete nonlinear form. The results are presented in Fig. 8. Comparing with Fig. 7, it can be seen that there is a much better error reduction if the EMAC form is used instead of the divergence form. For small viscosity coefficients and on coarse grids, one can observe for the P2/P1 pair of spaces an order of error reduction of 1.5, which is the same observation as for the evolutionary Stokes problem in Example 4.2. In case of the P3/P2pair of spaces, the order of error reduction is 2, which is exactly the order that is proposed by the robust error bound from Theorem 4.5.□ Theorem 4.5 and Example 4.6 allow to address another open problem stated in [17]: From the theoretical point of view, the new EMAC form offers several attractive features. Further numerical studies are necessary to clarify whether this form should be preferred also in simulations. With respect to robust velocity estimates in L∞(L2) for the evolutionary Navier–Stokes equations, the answer is twofold. On the one hand, for the Galerkin method and non weakly divergence-free and inf–sup stable pairs of finite element spaces, the EMAC form outperforms all other 21 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 Fig. 8. Example 4.6 Galerkin method for the Navier–Stokes equations, EMAC form of the discrete convective term. discrete nonlinear forms. It is the only form for which a robust estimate is possible, which is of order r−1. But on the other hand, as it will be shown in Section 4.2, using weakly divergence-free pairs of spaces or applying a rather simple stabilization (grad–div stabilization) for non weakly divergence-free pairs with arbitrary form of the discrete convective term leads to robust estimates of order r, i.e., one order higher. 4.2. Methods of order r in L∞(L2) This section presents discretizations whose order of robust error reduction for the velocity error in L∞(L2) is r. First, the Galerkin discretization with inf–sup stable pairs of finite element spaces and weakly divergence-free velocity is studied, in Section 4.2.1. Second, Section 4.2.2 considers non weakly divergence-free and inf–sup stable pairs of finite element spaces, where the Galerkin method is augmented either with a grad–div stabilization or with a LPS stabilization for the velocity gradient. The final part, Section 4.2.3, discusses three methods for pairs of finite element spaces that are not inf–sup stable and thus require a pressure stabilization. Two of these methods achieve stability with respect to dominating convection with the same mechanisms as the methods from Section 4.2.2. For the last method, it is shown that an appropriate pressure stabilization might be sufficient for stabilizing also dominant convection. 4.2.1. Inf–sup stable and weakly divergence-free pairs of finite element spaces This section considers the Galerkin method (42) for conforming inf–sup stable pairs of finite element spaces with Vdiv h⊂Vdiv. There are a number of pairs that satisfy this property. The most popular one is the Scott–Vogelius pair [54], Pr/Pdisk r−1,r≥d. However, the discrete inf–sup condition for the case r=dcan be proved only for special grids, so-called barycentric-refined grids, see [55,56]. Barycentric-refined grids are constructed from a given simplicial grid by connecting the vertices of each mesh cell with the barycenter of that mesh cell. Notice that although Vdiv h⊂Vdiv, one should not use the notion that the approximation uhis pointwise divergence-free, since ∇ ·uhis not continuous across boundaries of mesh cells. The use of weakly divergence-free pairs of finite element spaces is attractive from the point of view of the physical consistency of numerical solutions. First, there is an exact conservation of mass in the sense that ∥∇·uh(t)∥0=0 for all times. Second, provided that possible consistency errors do not depend on the pressure, the numerical solution satisfies a so-called fundamental invariance property: gradient forces on the right-hand side possess an impact only on the discrete pressure. Closely connected to this property is the so-called pressure robustness of the method, i.e., error bounds for the velocity do not depend on the pressure, compare (58). Thus, possibly large norms of the pressure are not of importance for velocity errors. All these properties are in contrast to pairs of finite element spaces that are not weakly divergence-free, see [47] for a comprehensive discussion of this topic. If Vdiv h⊂Vdiv, the finite element velocity solution uhis weakly divergence-free, hence the second term of the divergence form (44) of the nonlinear term vanishes for uhand the convective form is skew-symmetric. Consequently, one can use the convective form of the convective term in the numerical analysis as well as 22 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 in simulations. There are less terms to be estimated with the convective form compared with the divergence form. The error analysis of the Galerkin method can be found in [57]. Error bounds of order rare proved for the L∞(L2) norm of the velocity. Here, a simplified proof is presented. In this proof, the velocity approximation is compared with the Stokes projection shdefined in (32) and the decomposition u−uh=(u−sh)−(uh−sh)=(u−sh)−φhis utilized. From the error equation, one obtains by using the definition of the Stokes projection and that φhis weakly divergence-free ∥φh∥0 d dt ∥φh∥0≤⏐⏐((uh·∇)uh,φh)−((sh·∇)sh,φh)⏐⏐+⏐⏐(τ1,φh)⏐⏐,(54) where τ1=(∂tu−∂tsh)+((u·∇)u−(sh·∇)sh).Note that in contrast to the case where the finite element velocity is not weakly divergence-free, there is no pressure contribution in the error equation, i.e., a term like the last term on the right-hand side of (39) does not appear. It was shown in Section 4.1.2 that exactly this term is responsible for obtaining only order r−1 for the velocity error in L∞(L2). Using the skew-symmetric property of the nonlinear term for weakly divergence-free functions and H¨ older’s inequality leads to ⏐⏐((uh·∇)uh,φh)−((sh·∇)sh,φh)⏐⏐≤⏐⏐((φh·∇)sh,φh)⏐⏐+⏐⏐((uh·∇)φh,φh)⏐⏐(55) =⏐⏐((φh·∇)sh,φh)⏐⏐≤ ∥∇sh∥L∞∥φh∥2 0. The first term of the truncation error is bounded by using the Cauchy–Schwarz inequality and (33) for the temporal derivative (∂tu−∂tsh,φh)≤ ∥∂tu−∂tsh∥0∥φh∥0≤Chr∥∂tu∥r∥φh∥0. For estimating the second term, H¨ older’s inequality, Sobolev embedding and interpolations, (34), and (33) are utilized ((u·∇)u−(sh·∇)sh,φh)=(((u−sh)·∇)u−(sh·∇)(sh−u),φh) ≤(∥u−sh∥L6∥∇u∥L3+∥sh∥L∞∥∇(sh−u)∥0)∥φh∥0 ≤C(∥∇u∥1/2+∥sh∥L∞)∥sh−u∥1∥φh∥0 ≤C((∥u∥1∥u∥2)1/2+∥sh∥L∞)∥sh−u∥1∥φh∥0 ≤Chr(∥u∥1∥u∥2)1/2∥u∥r+1∥φh∥0.(56) Inserting the bounds for the truncation error and (55) in the error equation (54) and applying (34) yield d dt ∥φh∥0≤ ∥∇sh∥L∞∥φh∥0+Chr(∥∂tu∥r+(∥u∥1∥u∥2)1/2∥u∥r+1) ≤C∥∇u∥L∞∥φh∥0+Chr(∥∂tu∥r+(∥u∥1∥u∥2)1/2∥u∥r+1). Defining L(T)=C∫T 0∥∇u(s)∥L∞ds,Mu=C(∫T 0(∥∂tu(s)∥r+(∥u∥1∥u∥2)1/2∥u(s)∥r+1)ds),(57) and applying Gronwall’s lemma give ∥φh(t)∥0≤eL(T)∥φh(0)∥0+eL(T)Muhr. 23 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 The proof is finished by applying the triangle inequality. Note that all steps of the proof are dimensionally correct for d=3, with the unit md+2/s3on both sides of all estimates.2 Theorem 4.7 (Error Estimate for the Galerkin Method and Weakly Divergence-Free Pairs of Finite Element Spaces). For T >0, assume that the solution of (5) is sufficiently regular and choose uh(0) =sh(0). Then, it holds for t∈(0,T] ∥u(t)−uh(t)∥0≤eL(T)Muhr+Chr+1∥u∥L∞(Hr+1),(58) where L(T)and Muare defined in (57). It should be emphasized that the error bound (58) for the velocity error does not depend on the pressure, since, as already noted above, the pressure terms in the error equation vanish. This kind of error estimate is called pressure-robust. In [58], a fully discrete finite element method for the incompressible Navier–Stokes equations is analyzed where Scott–Vogelius pairs of spaces Pr/Pdisk r−1,r≥d, are utilized for the spatial discretization and the so-called rotational form of the nonlinear convective term is used. The rotational form goes along with a modified pressure, the socalled Bernoulli pressure. Thanks to the use of weakly divergence-free finite elements for the velocity, the method is pressure robust. The velocity error in L∞(L2) is proved to be of order r, but the corresponding error bound is not robust for small values of the viscosity coefficient. Example 4.8 (Galerkin Method with a Weakly Divergence-Free Pair of Finite Element Spaces).The same problem as in Example 4.4 is considered but with the Scott–Vogelius pair of spaces P2/Pdisk 1. A barycentric refinement needs to be applied to the grids arising from refining the initial grid from Fig. 1 leading to the number of degrees of freedom given in Table 2. The errors obtained for the Galerkin method (42) are shown in Fig. 9. One can see that the errors are small, where the order of error reduction is two on coarse grids and for small viscosity coefficients. This observation corresponds to the analytic prediction. Third order convergence can be observed for the large viscosity coefficient ν=10−2already on coarse grids and for smaller viscosity coefficients on finer grids. □ The optimal behavior for large viscosity coefficients on coarse grids was observed also for other methods and it will not be reported in further numerical studies. Remark 4.9 (Summary of the Results for the Plain Galerkin Method).At this point, we can summarize the changes in the rate of error reduction of the L∞(L2) error of the velocity for the plain Galerkin method applied to the different equations, which are discretized with H1-conforming velocity finite elements. For linear convection– reaction–diffusion equations, using piecewise polynomials of degree r, while the optimal rate r+1 is achieved for νlarge relative to the coefficients of the first order terms, only rate rcan be obtained when ν→0. For the Stokes and Oseen equations only the rate r−1 is achieved, while for the Navier–Stokes equations, one cannot even prove robust error reduction unless the EMAC form of the nonlinear term is utilized (in that case the rate r−1 is recovered) or weakly divergence-free pairs of finite elements are used (in that case the rate is increased to r). By far most of the finite element analysis for evolutionary incompressible flow problems has been performed for H1-conforming velocity finite element spaces. However, many pairs with such velocity spaces do not lead to weakly 2In order not to overload the presentation, dimensional correctness will be taken care of in the following only for d=3. Only for this method, the changes for d=2 are addressed in some detail once more. Instead of (56), the following estimate is performed, applying the same tools, ((u·∇)u−(sh·∇)sh,φh) =(((u−sh)·∇)u−(sh·∇)(sh−u),φh) ≤(∥u−sh∥L4∥∇u∥L4+∥sh∥L∞∥∇(sh−u)∥0)∥φh∥0 ≤C(∥u−sh∥1/2∥∇u∥1/2+∥sh∥L∞∥∇(sh−u)∥0)∥φh∥0 ≤C((∥u−sh∥0∥u−sh∥1)1/2(∥u∥1∥u∥2)1/2+(∥u∥0∥u∥2)1/2∥u−sh∥1)∥φh∥0 ≤Chr(h1/2(∥u∥1∥u∥2)1/2+(∥u∥0∥u∥2)1/2)∥u∥r+1∥φh∥0. Then, the expression for Muchanges accordingly. 24 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 Fig. 9. Example 4.8 Galerkin method for the Navier–Stokes equation, P2/Pdisk 1on a barycentric-refined grid, convective form of the convective term. divergence-free discrete velocity solutions. Exceptions require usually a high polynomial degree and sometimes special grids, like the Scott–Vogelius pair of spaces for r=d. An attractive alternative approach consists in using H(div)-conforming velocity spaces. Such spaces are actually used for simulations of flows with small viscosity, e.g., in [59], and one can find numerical analysis for stationary problems, e.g., in [60,61]. To the best of our knowledge, there are only few contributions that study the finite element analysis of weakly divergence-free H(div)- conforming finite elements for evolutionary incompressible flow problems. In [62], this kind of discrete velocity spaces is considered together with pressure spaces such that a discrete inf–sup condition is satisfied. Concrete examples for such pairs are Raviart–Thomas (RT) spaces of order rfor the velocity and discontinuous pressure spaces of order ron simplicial grids. Also, Brezzi–Douglas–Marini (BDM) spaces can be used for approximating the velocity. It is referred to [47, Sec. 4.4] for more details concerning these spaces. The analysis in [62] is performed for the evolutionary Oseen equations (29). Because H(div)-conforming functions are generally not H1-conforming, special care has to be taken for discretizing the viscous term such that the consistency error becomes sufficiently small. In [62], the symmetric interior penalty (SIP) form, which is well known from discontinuous Galerkin methods, is applied. Since weakly divergence-free velocity finite element functions are considered, the convective form of the convective term can be applied. It is augmented with an upwind term defined on the facets of the mesh cells. An error bound of order rwith constants independent of inverse powers of the viscosity is proved for the L∞(L2) norm of the velocity. Like in the H1-conforming case, (58), the bound is independent of the pressure. The restriction to order ris caused from the bound of the convective term. The results from [62] are extended to the Navier–Stokes equations in [63]. In [64], the authors improve the error bounds of [62], obtaining order r+1/2 for the Navier–Stokes equations. This result fits into Section 4.3 and a sketch of the proof will be presented there, see Theorem 4.23. 4.2.2. Non weakly divergence-free inf–sup stable pairs of finite element spaces This section presents two methods that achieve robust error reduction of order rfor the velocity error in L∞(L2) if the standard skew-symmetric form (44) of the nonlinear term is utilized. For one of these methods, even convergence for the convective form of the nonlinear term is shown. The Galerkin method with grad–div stabilization and inf–sup stable pairs of mixed finite elements is analyzed in [50]. More precisely, the following discretization is considered: Find uh:(0,T]→Vh,r,ph:(0,T]→Qh,r−1 satisfying (∂tuh,vh)+ν(∇uh,∇vh)+b(uh,uh,vh)+µ(∇ ·uh,∇ ·vh)−(∇ ·vh,ph)=(f,vh)∀vh∈Vh,r, (∇ ·uh,qh)=0∀qh∈Qh,r−1,(59) where the nonlinear term is defined in (44) and µ > 0 is the grad–div stabilization parameter. Note that the grad–div term vanishes for weakly divergence-free finite element velocities, which gives the Galerkin method considered in the previous subsection. The additional term in (59) can be derived by adding to the strong form of the momentum balance of the Navier–Stokes equations (1) the vanishing term −µ∇(∇ · u), integrating the corresponding weak 25 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 (∂tuh,vh)+ν(∇uh,∇vh)+b(uh,uh,vh)−(∇ ·vh,ph)+µ(∇ ·uh,∇ ·vh)=(f,vh)∀vh∈Vh,r, (∇ ·uh,qh)+spres(ph,qh)=0∀qh∈Qh,r,(73) where the nonlinear term is defined in (44), the LPS pressure stabilization is given by spres(ph,qh)=∑ K∈Th τp,K(σ∗ h(∇ph), σ∗ h(∇qh))K,(74) µand {τp,K}are the grad–div and pressure stabilization parameters, respectively, and σ∗ h=Id −σr−1 his the fluctuation operator, where Id is the identity operator and σr−1 his a locally stable projection or interpolation operator from L2(Ω)donto Yh,r−1. This operator can be chosen as the Bernardi–Girault [75] or the Scott–Zhang [76], [77, § 4.8], interpolation operator. The spaces Yh,lare defined as follows: Yh,l={vh∈C0(Ω)|vh|K∈Pl(K),∀K∈Th},l≥1,Yh,l=(Yh,l)d, where Pl(K) is the space of polynomials on Kof degree less than or equal to l. The following approximation property holds for the fluctuations ∥σ∗ h(v)∥0= ∥(Id −σj h)(v)∥0≤Chs|v|s,1≤s≤j+1.(75) The physical unit of the pressure stabilization parameters is second, hence it is a time scale. The proof of the following theorem can be obtained arguing as in [72, Theorem 3.5]. Theorem 4.14 (Non inf–sup Stable Pairs of Finite Element Spaces with grad–div Stabilization and LPS Pressure Stabilization).Consider T >0and assume sufficient regularity for the solution of (5). Let r ≥2and let α1h2 K≤τp,K≤α2h2 K,(76) for some positive constants α1, α2independent of h and ν. Then, there exists a positive constant Mu,pdepending on ∥u0∥rand (∫T 0((α2+µ−1)∥p(s)∥2 Hr/R+T∥∂tu(s)∥2 r +(ν+µ+α−1 1+T|Ω|4−d d∥u(s)∥2 2)∥u(s)∥2 r+1)ds)1/2 such that for t ∈[0,T]the following bound holds ∥u(t)−uh(t)∥0≤Mu,pexp (L(T) 2)hr+Chr+1∥u∥L∞(Hr+1),(77) where L(T)is L(T)=1+C∫T 0|Ω|4−d 2d(2∥u(s)∥3+|Ω|4−d 2d∥u(s)∥2 2 2µ)ds. Proof (Sketch).This sketch stresses only the main new aspect if non inf–sup stable pairs of finite element spaces are studied instead of inf–sup stable pairs. For the proof of Theorem 4.14, see [72, Theorem 3.5], one compares the approximation uhwith ˆuh=Rhu∈Vh,r satisfying (u− ˆuh,vh)=0,∀vh∈Yh,r−1.(78) The interpolation operator Rhuexists and satisfies optimal approximation properties, see [73]. In the proof of Theorem 4.14, the nonlinear term is bounded as in (62) with shreplaced by ˆuh. In [72, (3.27), (3.28)] it is shown that ∥ˆuh∥L∞and ∥∇ˆuh∥L∞can be bounded in terms of ∥u∥2and ∥u∥3, respectively, but with constants that are not scale-invariant. Careful dimensional analysis shows that ∥ˆuh∥L∞≤C|Ω|4−d 2d∥u∥2,∥∇ˆuh∥L∞≤C|Ω|4−d 2d∥u∥3. 32 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 A new aspect for non inf–sup stable pairs of finite element spaces is that one cannot apply test functions that are discretely divergence-free, since such finite element functions do not necessarily exist. Consequently, the error equation contains a term that couples the discrete pressure ph, which does not appear in the error equation for inf–sup stable pairs of finite element spaces like (38), and the divergence of the velocity test function. The finite element pressure phis compared with ˆph, which is defined as the Lagrange interpolant of the continuous pressure pto which one subtracts its mean. Property (78) is required for bounding the term (∇ · ˆuh, λh), which appears in the error equation, where λh= ˆph−ph. Using that uis divergence-free, integration by parts, the range of σr−1 h, and (78) yields |(∇ · ˆuh, λh)|=|(∇ ·(ˆuh−u), λh)|=|(ˆuh−u,∇λh)|=|(ˆuh−u, σ∗ h(∇λh))|. From this equation, one obtains with (76) and the approximation properties of ˆuh |(∇ · ˆuh, λh)| ≤ Cα−1 1h2r∥u∥2 r+1+1 4spres(λh, λh), such that the second term on the right-hand side can be absorbed into the left-hand side of the error equation. □ From the error bound (77), one observes the order of convergence r,r≥2, for the L∞(L2) error of the velocity. The error bound (77) is the analog to (61) for the Galerkin method with grad–div stabilization and inf–sup stable elements. As for inf–sup stable pairs of finite element spaces, an asymptotic optimal choice of the grad–div stabilization parameter with respect to the mesh width is µ∼1, due to the appearance of both µand µ−1on the right-hand side of (77). It is noted in [72] that under assuming additional regularity of the pressure, the choice µ∼his likewise asymptotically optimal. From the point of view of the physical dimension, the parameters {τp,K}are time scales. Hence, the inverse of α1and α2are viscosity parameters. In fact, for the steady-state Stokes equations with viscosity, the parameter for LPS pressure stabilizations, and also several other pressure stabilizations, has the form Ch2 K/ν, e.g., compare the review on pressure stabilizations [71]. But for these equations, one cannot derive robust error bounds, because the dependency on inverse powers of the viscosity enters from estimating the pressure term. Inspecting the analysis from [72], which leads to the robust estimate (77), shows that the choice τp,K=Ch2 K/ν does not enable to derive a robust estimate since τp,Khas to be bounded from above without having the factor νavailable in the corresponding expression. Another method analyzed in [72] uses LPS stabilization of the velocity gradient in the discrete momentum equation instead of grad–div stabilization: Find uh:(0,T]→Vh,r,ph:(0,T]→Qh,rsuch that (∂tuh,vh)+ν(∇uh,∇vh)+b(uh,uh,vh)−(∇ ·vh,ph) +∑ K∈Th τν,K(σ∗ h(∇uh), σ∗ h(∇vh))K=(f,vh)∀vh∈Vh,r,(79) (∇ ·uh,qh)+spres(ph,qh)=0∀qh∈Qh,r, where the nonlinear term is defined in (44) and the LPS pressure term in (74). The LPS velocity term is a viscous term, hence the parameters {τν,K}represent viscosity parameters. Theorem 4.15 (Non inf–sup Stable Pairs of Finite Element Spaces with Velocity Gradient LPS Term).Let T >0, assume that the solution of (5) is sufficiently smooth, let r ≥2, and let the stabilization parameters satisfy α1h2 K≤τp,K≤α2h2 K, τ1≤τν,K≤τ2,(80) for some positive constants α1, α2,τ1, and τ2, independent of h. Then, there is a positive constant Mu,p, which depends on ∥u0∥rand (∫T 0((α2+τ−1 1)∥p(s)∥2 Hr/R+T∥∂tu(s)∥2 r+(ν+τ1+T|Ω|4−d d∥u(s)∥2 2)∥u(s)∥2 r+1)ds)1/2 , such that for t ∈[0,T]the error bound ∥u(t)−uh(t)∥0≤Mu,pexp (L(T) 2)hr+Chr+1∥u∥L∞(Hr+1) 33 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 holds, where L(T)is given by L(T)=1+C∫T 0|Ω|4−d 2d(∥u(s)∥3+|Ω|4−d 2d∥u(s)∥2 2 τ2)ds. Proof (Sketch).Again, the presentation will be restricted to the main issues in the proof of Theorem 4.15, see the proof of [72, Theorem 4.5] for more details. As in the proof of Theorem 4.14,uhwill be compared with ˆuhsatisfying (78). To bound error terms containing the pressure appropriately, phis compared with ˆph=Rhp, where ˆphsatisfies a property analogous to (78) for scalar functions (p− ˆph,qh)=0,∀qh∈Yh,r−1.(81) Denote by φh= ˆuh−uhand by λh= ˆph−ph. The main step in the proof is the bound of the nonlinear convective term in the error equation. It will be shown that this term can be bounded with a constant independent of inverse powers of the viscosity using only the pressure stabilization. Arguing as in (46), the first step of bounding the nonlinear term is |b(uh,uh,φh)−b(ˆuh,ˆuh,φh)| ≤ ∥∇ˆuh∥L∞∥φh∥2 0+1 2((∇ ·φh)ˆuh,φh). For the second term on the right-hand side, the key observation is that it has a structure that appears in the discrete continuity equation. A decomposition of this term will be applied that allows to utilize the error equation with respect to the continuity equation, introducing in this way spres(·,·). More precisely, the decomposition has the form ((∇ ·φh)ˆuh,φh)=(∇ ·φh, σr h(ˆuh·φh))+(∇ ·φh,(I−σr h)(ˆuh·φh)).(82) For the first term on the right-hand side of (82), one can use (0, σr h(ˆuh·φh))as test function in the error equation, since by the definition of the projection operator, it is an admissible choice. In this way, a link between the pressure stabilization and the estimate of the nonlinear convective term is established (∇ ·φh, σr h(ˆuh·φh))=spres(ph, σr h(ˆuh·φh)) +(∇ · ˆuh, σr h(ˆun+1 h·φh)).(83) Using the approximation properties of ˆuhand the L2(Ω) stability of the operator σr h, one obtains for the second term on the right-hand side of (83) (∇ · ˆuh, σr h(ˆun+1 h·φh)) ≤Ch2rτ2∥u∥2 r+1+Cτ−1 2∥ˆuh∥2 L∞∥φh∥2 0. For the first term on the right-hand side of (83), one adds and subtracts ˆphin the first argument of spres(·,·). Since spres(·,·) is a continuous, symmetric, positive semi-definite bilinear form, compare Assumptions A1 and A2 for a more general pressure stabilization, there holds in particular a Cauchy–Schwarz inequality. Utilizing this Cauchy– Schwarz inequality, the bound of the pressure stabilization parameters (80), the L2(Ω) stability of σr h, the inverse inequality (8), the approximation properties of ˆph, and the approximation properties (75) of the fluctuation operator gives spres(ph, σr h(ˆuh·φh)) ≤1 4spres(ˆph,ˆph)+3spres(σr h(ˆuh·φh), σr h(ˆuh·φh)) +1 8spres(λh, λh) ≤Cα2h2(∥σ∗ h(∇(p− ˆph))∥2 0+∥σ∗ h(∇p)∥2 0+∥σ∗ h(∇(σr h(ˆuh·φh)))∥2 0)+1 8spres(λh, λh) ≤Cα2h2r∥p∥2 r+Cα2∥ˆuh∥2 L∞∥φh∥2 0+1 8spres(λh, λh). To conclude the estimate of the nonlinear term, one needs to bound the second term on the right-hand side of (82). To this end, the following result from [78, Theorem 2.2] is applied: for v∈W1,∞(Ω)dand vh∈Yh,r, it holds ∥(I−σr h)(v·vh)∥0≤Ch∥∇v∥L∞∥vh∥0.(84) Applying (84) to v= ˆuhand vh=φhtogether with the inverse inequality (8) yields ((∇ ·φh),(I−σr h)(ˆuh·φh))≤C∥∇ ·φh∥0h∥∇ˆuh∥L∞∥φh∥0≤C∥∇ˆuh∥L∞∥φh∥2 0. 34 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 Collecting all estimates leads to ⏐⏐b(uh,uh,φh)−b(ˆuh,ˆuh,φh)⏐⏐≤C(∥∇ˆuh∥∞+(α2+τ−1 2)∥ˆuh∥2 ∞)∥φh∥2 0. +Ch2r(α2∥p∥2 r+τ2∥u∥2 r+1)+1 8spres(λh, λh).(85) The first term on the right-hand side is handled as in the previous theorems by applying Gronwall’s lemma, the second term becomes part of the error bound, and the last term is absorbed into the left-hand side of the error equation. We want to emphasize that in the above proof, it is essential to have velocity and pressure finite element spaces containing polynomials of the same degree. This fact comes from the application of (84) to a function in the velocity space (vh∈Yh,r) together with the use of σr h(ˆuh·φh) as a test function for the pressure. To finish with the sketch of the proof of Theorem 4.15, it shall be indicated in which step property (81) is required. Since there is no grad–div stabilization in the method, there is no direct control on the divergence of the error and, as a consequence, the following term, which appears in the error equation, cannot be bounded in the same way as in Theorem 4.14. Using the orthogonality property (81),∥∇ ·φh∥0,K≤√d∥∇φh∥0,K,τν,K> τ1, and the approximation properties of ˆphgives (p− ˆph,∇ ·φh)=(p− ˆph, σ∗ h(∇ ·φh)) ≤⎛ ⎝∑ K∈Th τ−1 ν,K∥p− ˆph∥2 0,K⎞ ⎠ 1/2 √d⎛ ⎝∑ K∈Th τν,K∥σ∗ h(∇φh)∥2 0,K⎞ ⎠ 1/2 ≤Ch2rτ−1 1∥p∥2 r+1 4∑ K∈Th τν,K∥σ∗ h(∇φh)∥2 0,K.(86) The second term on the right-hand side is now absorbed into the left-hand side of the error equation. □ Hence, also method (79) provides convergence of order rfor the velocity in L∞(L2) with constants independent of inverse powers of ν. We like to mention that velocity fluctuations are also used to define turbulence models for the simulation of incompressible turbulent flows. In so-called three-scale algebraic VMS-multigrid methods, proposed in [79,80], the deformation tensor of the velocity fluctuations is the main part of the turbulence model and in the three-scale coarse space projection-based VMS method, proposed in [81], fluctuations of the deformation tensor are the main part, e.g., see [82] for a survey of VMS methods. The parameter used in these methods depends typically on the numerical solution and therefore it is nonlinear. Inspecting the proof of Theorem 4.15, it can be observed that the LPS term with the velocity gradient can be replaced with a respective LPS term for the divergence ∑ K∈Th τµ,K(σr h(∇ ·uh), σr h(∇ ·vh))K, with τµ,K∼1, and the same order of convergence as in Theorem 4.15 can be proved, compare [72, Theorem 4.6]. There is another LPS method for equal order pairs of finite element spaces for which a finite element error analysis can be found in the literature. This method includes term-by-term LPS stabilization for the convective term, the divergence and the pressure. However, the available error bound for this method is not robust. For this reason, its discussion is postponed to Section 4.4.2. It was already emphasized in the proof of Theorem 4.15 that the convective term in the error equation can be bounded with a constant independent of inverse powers of the viscosity using only the pressure stabilization. Now, the natural question is whether using only an appropriate pressure stabilization can lead to a method where error bounds of this kind can be derived. This question is answered positively in [83]. While finishing the present review, we became aware of Ref. [84], where equal order pairs of finite elements with symmetric pressure stabilization are considered for discretizing the evolutionary Navier–Stokes equations. Robust error bounds of order rare proved for the L∞(L2) error of the velocity. We found that essentially the same idea explained in the proof of Theorem 4.15 is applied to obtain the bounds in [84]. A fully discrete version of the following method with symmetric pressure stabilization, see Assumption A1 below for justifying this notion, using the implicit Euler scheme, is studied in [83]: Find uh:(0,T]→Vh,r, 35 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 ph:(0,T]→Qh,rsuch that (∂tuh,vh)+ν(∇uh,∇vh)+b(uh,uh,vh)−(∇ ·vh,ph)=(f,vh), (∇ ·uh,qh)+αspres(ph,qh)=0,(87) where the nonlinear term is defined in (44),αis a parameter whose units are time over length squared (i.e., the units of the inverse of the kinematic viscosity) and the pressure stabilization term should satisfy a number of assumptions. Remark 4.16 (Assumptions on the Pressure Stabilization Term). A1 The pressure stabilization spres :Qh,r×Qh,r→Ris a symmetric, positive semi-definite bilinear form. A2 Continuity. For all qh∈Qh,r, it holds spres(qh,qh)≤C∥qh∥2 0with a constant C>0 independent of hand ν. A3 Weak consistency. There exists a projection Πh:Q→Qh,rsuch that ∥q−Πhq∥0≤Chs∥q∥s,q∈Q∩Hs,0≤s≤r+1, and spres(Πhq,Πhq)1/2≤Chs∥q∥s,q∈Q∩Hs,1≤s≤r. Typically, Πhis the L2(Ω) projection or some interpolant for non-smooth functions like the Scott–Zhang interpolant. A4 There is an interpolation operator Ir h:V→Vh,rsuch that for 1 ≤s≤rand for all v∈Vand qh∈Qh,r (∇ ·(v−Ir hv),qh)≤C⎛ ⎝∑ K∈Th h−2 K∥v−Il hv∥2 0,K+∥v−Ir hv∥2 1⎞ ⎠ 1 2 spres(qh,qh)1 2, and ∥v−Ir hv∥m≤Chs+1−m|v|s+1∀v∈Hs(Ω)d, for m∈ {0,1}, 1 ≤s≤r, with C>0 independent of hand ν. Assumptions A1–A4 are quite similar to those in [85], where finite element methods for the time-dependent Stokes equations are analyzed. Pressure stabilizations can be divided into two classes, e.g., see the survey [71]: residual-based stabilizations and stabilizations that use only the pressure. It turns out that practically all stabilizations of the latter class satisfy Assumptions A1–A4: the stabilization of Brezzi and Pitk¨ aranta [86], the stabilization of Dohrmann and Bochev [87], the orthogonal subscale stabilization proposed in [88], the classical LPS stabilization [89], the term by term local projection stabilization from [73], and the continuous interior penalty (CIP) stabilization proposed in [85,90]. The following error estimate is derived in [83, Theorem 2]. Theorem 4.17 (Non inf–sup Stable Pairs of Finite Element Spaces with Symmetric Pressure Stabilization).Let, for T >0, the solution of (5) be sufficiently regular and let r ≥2. Let Assumptions A1–A4 for the pressure stabilizations be satisfied, and let α∼1. Then, there holds for t ∈[0,T]that ∥u(t)−uh(t)∥0≤Mu,pexp (L(T) 2)hr+Chr+1∥u∥L∞(Hr+1), where Mu,p>0depends on (∥u0∥2 r+∫T 0[α∥p(s)∥2 Hr/R+T∥p(s)∥2 Hr+1/R+T∥∂tu(s)∥2 r +(ν+α−1+|Ω|4−d d∥u(s)∥2 2)∥u(s)∥2 r+1]ds)1/2 and L(T)is defined by L(T)=1+C∫T 0|Ω|4−d 2d(∥u(s)∥3+α|Ω|4−d 2d∥u(s)∥2 2)ds. 36 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 Table 4 Numerical examples for the Navier–Stokes equations, number of degrees of freedom (including Dirichlet nodes) for non inf–sup stable pairs of finite element spaces on triangular grids. P2/P2 Level Velocity Pressure 3 290 145 4 1090 545 5 4226 2113 6 16642 8321 7 66050 33025 Proof (Sketch).For proving Theorem 4.17,uhis compared with shdefined as the velocity component of the pair (sh,lh)∈(Vh,r,Qh,r) solving (∇sh,∇vh)−α(lh,∇ ·vh)=(∇u,∇vh)∀vh∈Vh,r, (∇ ·sh,qh)+αspres(lh,qh)=0∀qh∈Qh,r. A standard finite element error analysis shows that ∥sh−u∥0+h∥∇(sh−u)∥0+αhspres(lh,lh)1/2≤Chr+1∥u∥r+1. The finite element pressure solution phis compared with ˆph=Πhp, where Πhis the projection in Assumption A3. The same arguments as in the proof of Theorem 4.15 are applied to bound the nonlinear term. Estimating the term (p− ˆph,∇ · φh), with φh=sh−uh, is different than in (86) and it requires to increase the regularity assumptions on the pressure from Hr(Ω) to Hr+1(Ω), compared with Theorem 4.15. Integrating by parts yields (p− ˆph,∇ ·φh)= −(∇(p− ˆph),φh)≤Ch2rT∥p∥r+1+1 T∥φh∥2 0, such that the second term on the right-hand side can be handled with Gronwall’s lemma. □ Theorem 4.17 is formulated for quadratic and higher order finite elements. However, as stated in [83, Remark 3], an analogous error bound can be proved also for linear elements. Concerning low order methods, in [91] a weakly divergence-free method using continuous linear elements for the velocity and piecewise constants for the pressure is considered. A stabilization term for the pressure is added to the method that is enhanced with a discrete divergencefree convective field. A priori bounds of order hare obtained for the L∞(L2) error of the velocity with constants independent of inverse powers of the viscosity. Example 4.18 (Non inf–sup Stable Pairs of Finite Element Spaces).Numerical results for method (73) with grad–div and LPS pressure stabilization and for method (87) with only LPS pressure stabilization are presented in Fig. 13. Exactly the same problem and the same meshes are considered as in Example 4.4. Simulations were performed for the P2/P2pair of finite element spaces, see Table 4 for the corresponding numbers of degrees of freedom. For method (73), the grad–div stabilization parameter was set to be µ=0.1 and the pressure stabilization parameter to be τp,K=1.0. The pressure stabilization parameter in method (87) was chosen to be τp,K=0.1. The numerical results in Fig. 13 show for small viscosity coefficients and on coarse grids second order error reduction for the velocity error in L∞(L2), as predicted by numerical analysis. If the spatial grids become sufficiently fine with respect to the viscosity coefficient, an increase of the order can be observed in this example for both methods. □ 4.3. Methods of order r +1/2in L∞(L2) This section surveys methods where the L∞(L2) error of the velocity can be proved to be of order r+1/2 with a robust error bound. These methods are a special realization of the SUPG method, a CIP method, a modification of the method with LPS stabilization of the velocity gradient (79), a similar method with a different type of LPS stabilization, and a method based on H(div)-conforming discontinuous Galerkin elements. Except for the last one, all the methods are formulated for equal order pairs of finite element spaces. 37 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 Fig. 13. Non inf–sup stable pair of finite element spaces, method (73) (left) and method (87) (right). Since the SUPG or streamline diffusion method is one of the most popular finite element discretizations for scalar convection–diffusion equations, it is natural that it is utilized also in the context of the incompressible Navier– Stokes equations (5). In fact, probably the first paper that analyzes a finite element method for the time-dependent Navier–Stokes equations in the standard velocity–pressure formulation for high Reynolds number flows (small ν), paper [13], uses the SUPG method as stabilization. A main topic of the analysis of the SUPG method for scalar convection–diffusion equations presented in [24] is the choice of the stabilization parameter. By construction, this parameter is a time scale, see Section 3.2. However, the major difficulties of the numerical solution of these equations arise from very small spatial structures, the layers, such that it seems to be reasonable that the stabilization parameter should depend on the local mesh width and not on the length of the time step. The first numerical analysis for stabilization parameters depending on the mesh width was provided in [24], e.g., by considering the continuous-in-time situation, where a time step is not present. However, for the SUPG method in the framework of the incompressible Navier–Stokes equations, an analysis of the continuous-in-time situation does not yet seem to be available, to the best of our knowledge. For this reason, unlike to the other methods discussed in this paper, a fully discrete version of the SUPG stabilization will be discussed here. In [13], a SUPG method is analyzed for space–time finite elements, where velocity and pressure spaces consist of continuous piecewise linear functions (P1) in space and discontinuous piecewise linear functions in time. To the best of our knowledge, there seems to be no extension of this analysis available to higher order finite elements in space. Hence, the presentation below will be restricted to this case. A main assumption of the analysis in [13] is that ∆t∼h. Since it suffices for the discussion of the estimate of the velocity error in L∞(L2) and it simplifies the presentation considerably, the SUPG method will be considered here only in a form which was analyzed in [13] for |∇ ·uh| ≥ C>0, where Cis an arbitrary constant not depending on νand the mesh width. Let 0 =t0<t1<··· <tN=T, denote by Ωn=Ω×(tn−1,tn), Vn h,1=Vh,1×P1(tn−1,tn), and Qn h,1=Qh,1×P1(tn−1,tn). Then, the SUPG method reads as follows: Find (uh,ph)∈Vn h,1×Qn h,1such that for n=0,1,...,Nand all (vh,qh)∈Vn h,1×Qn h,1 (∂tuh,vh)Ωn+aSUPG(uh;(uh,ph),(vh,qh)) +µ(∇ ·uh,∇ ·vh)Ωn +∑ K∈Th δK(∂tuh, ∂tvh+(uh·∇)vh+∇qh)Kn+(un h,+−un h,−,vn h,+) =(f,vh)Ωn+∑ K∈Th δK(f, ∂tvh+(uh·∇)vh+∇qh)Kn,(88) 38 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 where Kn=K×(tn−1,tn), aSUPG(wh;(uh,ph),(vh,qh)) =ν(∇uh,∇vh)Ωn+((wh·∇)uh,vh)Ωn+1 2(∇ ·wh,uh·vh)Ωn−(∇ ·vh,ph)Ωn+(∇ ·uh,qh)Ωn +∑ K∈Th δK((wh·∇)uh+∇ph, ∂tvh+(wh·∇)vh+∇qh)Kn, and vn ±=lim s→0±v(tn+s),u0 −=u0. Note that in contrast to formulation (18) for the scalar equations, the temporal derivative of the velocity test function appears. The following error estimate is proved in [13]. Theorem 4.19 (Error Estimate for the SUPG Method).Assume that the solution of the Navier–Stokes equations (5) is sufficiently smooth for T >0. Let δK=δ∼h, ∆t∼h, µ∼h, and ν≤h. Then, for sufficiently small mesh width h, there is a constant C that does not depend on h and νsuch that the velocity solution of the SUPG method (88) satisfies ∥u(tn)−un h∥0≤Ch3/2(89) for all time instants tn, n =0,...,N. Proof (Sketch).The error analysis in [13] is performed in the norm |||(v,q)|||2 Sn=1 2(∥vn −∥2 0+ n−1 ∑ l=1∥vl +−vl −∥2 0+∥v0 +∥2 0)+ν∥∇v∥2 0,Sn +δ∥∂tv+(uh·∇)v+∇q∥2 0,Sn+µ∥∇ ·v∥2 0,Sn,(90) with Sn=Ω×(0,tn). It is clear that the bound (89) is obtained as a consequence of bounding the error in the norm (90). First, a stability estimate for the discrete solution is proved in [13]. This estimate is derived in the usual way by utilizing the discrete solution as test function. As final step, a discrete Gronwall inequality is applied, which requires that the length of the time step, and since ∆t∼his assumed, the mesh width, is sufficiently small. The proof of the error bound starts with considering the error equation, which is obtained by subtracting the discrete equation (88) from the continuous equation. Splitting the errors u−uh=(u−Ihu)−(Ihu−uh)=η−φh and p−ph=(p−Jhp)−(Jhp−ph)=(p−Jhp)−λh, where Ihand Jhare interpolation operators with the usual interpolation properties, and rearranging terms of the error equation yields |||(φh, λh)|||2 Sn=(∂tη+(uh·∇)η,φh)Sn+(∇ ·η, λh)Sn+ n−1 ∑ l=0∑ K∈Th δν(∆η, ∂tφh+(uh·∇)φh+∇λh)Kl+··· . Applying integration by parts gives |||(φh, λh)|||2 Sn= −(η, ∂tφh)Sn− n−1 ∑ l=0 (ηl+1 −,φl+1 h,−)−(ηl +,φl h,+) −(η,(uh·∇)φh)Sn−(∇ ·uh,η·φh)Sn−(η,∇λh)Sn + n−1 ∑ l=0∑ K∈Th δν(∆η, ∂tφh+(uh·∇)φh+∇λh)Kl+··· . 39 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 Concentrating on some important terms, using the Cauchy–Schwarz inequality and Young’s inequality lead to |||(φh, λh)|||2 Sn≤2 δ∥η∥2 0,Sn+δ 8∥∂tφh+(uh·∇)φh+∇λh∥2 0,Sn +2δν2 n−1 ∑ l=0∑ K∈Th∥∆η∥2 0,Kl+δ 8∥∂tφh+(uh·∇)φh+∇λh∥2 0,Sn+··· .(91) The second and fourth terms on the right-hand side of (91) are absorbed by the contribution of the SUPG term on the left-hand side. Using the interpolation property of ∥η∥0,Snand the choice of δgives that the first term is of order h3, with respect to the square of the error. The asymptotic of δand the condition that ν≤hlead also for the third term of (91) to third order. There are even more terms leading to order 3/2, but in this sketch of the proof, it is concentrated to terms that are related to the SUPG stabilization. The final step of the proof is again the application of a discrete Gronwall lemma. □ An error analysis of the SUPG method applied to the non-conforming Crouzeix–Raviart finite element of lowest order and combined with the backward Euler scheme as temporal discretization is presented in [92]. The SUPG stabilization parameter depends on the length of the time step. Robust estimates are derived. For a certain relation of the mesh width, the length of the time step, and a stabilization parameter for a velocity jump term, the velocity error in L∞(L2) is proved to be of order 4/3. The stabilization terms of the SUPG method are part of a two-scale residual-based variational multiscale (VMS) method for the simulation of turbulent incompressible flows proposed in [14]. This VMS method possesses two more stabilization terms. Recent studies of a two-dimensional flow at high Reynolds number and at three-dimensional turbulent channel flows in [93,94] show that the results obtained with the SUPG method and the residual-based VMS method are similarly accurate. In [90], a CIP finite element method to approximate the solution of the Navier–Stokes equations (5) is considered. Equal order approximations of the velocity and pressure (uh,ph), belonging to the spaces Vh,rand Qh,r, respectively, are taken and the boundary conditions are weakly imposed in the method. Incompressibility, convective dominance, and pressure are stabilized by adding the following stabilization terms to the Galerkin formulation: j(uh,vh):= ∑ K∈Th∑ F∈∂K∫F h2 F[|∇uh|]F:[|∇vh|]Fds, juh(uh,vh):= ∑ K∈Th∑ F∈∂K∫F h2 F|I1 huh·nF|2[|∇uh|]F:[|∇vh|]Fds,(92) ι(ph,qh)=∑ K∈Th∑ F∈∂K∫F h2 F[|∇ph|]F·[|∇qh|]Fds, where the same notations as for convection–diffusion equations are used. In addition, I1 huhdenotes the interpolation of uhonto the space Vh,1of continuous piecewise linear polynomials. The robust, with respect to ν, stability of the method relies on the fact that the stabilization terms (92) control some interpolation errors of the streamline derivative, the divergence, and the pressure gradient. More precisely, let rhvdenote the Oswald quasi-interpolant, compare [90, Definition 1], then the following bounds hold for any wh∈Vh,1,vh∈Vh,r,qh∈Qh,r, see [90, Lemma 3]) ∥h1/2((wh·∇)vh−rh((wh·∇)vh))∥2 0≤C1jwh(vh,vh), ∥h1/2(∇ ·vh−rh(∇ ·vh))∥2 0≤C2j(vh,vh),(93) ∥h1/2(∇qh−rh(∇qh))∥2 0≤C3ι(qh,qh), where C1,C2,C3are positive constants.4 4We remark that the error analysis of the CIP method presented here follows closely [90], although, contrary to most methods in the present paper, the analysis and the corresponding estimates are not necessarily dimensionally correct. 40 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 Concerning the viscosity-independent control of the divergence, the following result is applied in the error analysis of the method, see [90, Lemma 7], C∥h1/2∇ ·uh∥2 0≤j(uh,uh)+ι(ph,ph)+∥uh·n∥2 0,∂Ω.(94) Let Πr hdenote both, the L2(Ω) orthogonal projection onto Vh,rand Qh,r, then the following bound, analogous to (84), is also essential in the error analysis, see [90, Corollary 1] and [78, Theorem 2.2], ∥Πr h(u·vh)−u·vh∥0≤Ch∥∇u∥L∞∥vh∥0,∀u∈W1,∞(Ω)d,vh∈Vh,r.(95) In the following theorem, whose proof can be found in [90, Theorem 2, Corollary 3], the rate of convergence of the velocity is stated. Theorem 4.20 (Error Estimate for the CIP Method).Let the solution of (5) be sufficiently smooth and assume ν < h. Then, the following bound holds for the velocity approximation uhof the CIP method with stabilization terms (92) ∥u−uh∥L∞(L2)≤C(T,u,p)hr+1/2, where C(T,u,p)does not depend on inverse powers of ν. Proof (Sketch).For proving the statement of Theorem 4.20,uhis compared with Πr huand phwith Πr hp. Denote φh=Πr hu−uh,η=u−Πr hu, and λh=Πr hp−ph. Some details concerning the bounds of the more involved terms appearing in the error equation will be presented. First, terms that come from a decomposition of the nonlinear terms in the error equations are considered. One of them is (η,(uh·∇)φh). Adding and subtracting I1 huhand using the orthogonality properties of the L2(Ω) projection with respect to the range of rhyield (η,(uh·∇)φh)=(η,((uh−I1 huh)·∇)φh)+(η,(I1 huh·∇)φh−rh(I1 huh·∇φh)).(96) For the first term on the right-hand side, one obtains after some manipulations, including inverse inequalities and the L∞(Ω) stability of the L2(Ω) projection, (η,((uh−I1 huh)·∇)φh)≤C∥∇u∥L∞(∥φh∥2 0+h−1/2∥η∥2 0). The second term on the right-hand side of (96) is a truncation error, even of higher order than 2r+1. For this term, applying (93), one can prove (η,I1 huh·∇φh−rh(I1 huh·∇φh)) ≤Cε−1 1h−1∥η∥2 0+ε1C1juh(φh,φh), where ε1is any positive parameter. Note that the square root of the first term on the right-hand side is of order r+1/2. For the term (∇ ·uh,η·φh), appearing also in the error equation, the control of the divergence given by (94) is applied. Bounding the term (∇ ·uh,u·φh) covers three pages in the proof of [90, Theorem 2]. The first observation is that one can take vh=0and qh=Πr h(u·φh) in the error equation, which gives (∇ ·uh,Πr h(u·φh)) =(uh·n,Πr h(u·φh))∂Ω−ι(ph,Πr h(u·φh)). Then, adding and subtracting Πr h(u·φh) and using this equality yield (∇ ·uh,u·φh)=(∇ ·uh,u·φh−Πr h(u·φh)) +(uh·n,Πr h(u·φh))∂Ω−ι(ph,Πr h(u·φh)). After some work and with the help of (95), one finally arrives at (∇ ·uh,u·φh)≤Cε2(ι(Πr hp,Πr hp)+ι(λh, λh))+Cε−1 2h∥∇u∥2 L∞∥φh∥2 0+Cε−1 2∥u∥2 L∞j(φh,φh), for any positive parameter ε2. Then, it can be observed that, apart from applying (95), the stabilization terms included to achieve the control over the divergence and the gradient of the pressure are essential to finish the estimate of this term. 41 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 4.4.3. Other methods The grad–div stabilization introduced in Section 4.2.2 results in a mutual coupling between all components of the velocity. There are much more non-zero matrix entries than for other terms, e.g., the standard gradient form of the viscous term or the convective, skew-symmetric, and divergence formulations of the convective term. In [116], a sparse grad–div stabilization is proposed and analyzed that neglects several couplings of the standard grad–div stabilization. The natural expectation is that the sparse approach possesses weaker stability properties than the standard method. In [116], only the splitting error for a projection scheme is analyzed and the corresponding error bound is not robust. Another method to facilitate the application of the grad–div stabilization in practice is proposed in [117]. This modular method computes in the first step an intermediate solution (ˆuh,ph) by solving the Galerkin discretization of the Navier–Stokes equations. In the second step, a kind of regularization of the intermediate velocity field is applied, which is based on the grad–div operator. Several variants of this step are investigated in [117], where the simplest one is the solution of the variational problem µ∆t(∇ ·uh,∇ ·vh)+(uh,vh)=(ˆuh,vh)∀vh∈Vh,r, with an appropriate stabilization parameter µ. For a fully discretized problem, a non-robust estimate for the velocity error in L∞(L2) of order ris derived for different variants of the modular grad–div stabilization. An extension of this idea to the BDF2 temporal discretization is studied in [118], where a non-robust velocity estimate in L∞(L2) of order ris proved. Inserting the differential filter (103) in the momentum equation gives the additional term (wh−wh,vh)/(∆t). This term has in some sense a similar structure like (104): it is the difference of some term containing the finite element solution and a smoothed term with this solution, where smooth has to be understood in this context in the sense of less oscillations or variations of values. Approaches of this kind has been explored also with respect to the temporal evolution of the solution, leading to so-called temporal regularizations, which add, e.g., terms of the form γ(u−u)or −γ∆(u−u)(106) to the momentum equation. If fand uare smooth and γis sufficiently large, then both terms in (106) enforce that uis also smooth. In temporal regularizations, uis chosen as the solution from a previous time instant, or a combination of solutions from previous time instants, starting with an already smooth or smoothed initial condition. In [119], a temporal regularization with −γh∆(un+1−un) is analyzed for inf–sup stable pairs of finite element spaces. A non-robust bound of the velocity error in L∞(L2) of order ris proved. A so-called stabilization of the curvature is considered in [120], which adds a temporal regularization of the form γ ν(−ν∆Un+1+(un+1·∇)Un+1+∇Pn+1) with Un+1=un+1−2un+un−1and Pn+1=pn+1−2pn+pn−1. It is shown in [120], for inf–sup stable pairs of finite element spaces, that the velocity error in L∞(L2) is bounded of order rin a non-robust way. 5. Summary In the convergence theory of finite element methods for evolutionary convection–diffusion equations and incompressible Navier–Stokes equations, one can find two types of estimates: robust ones, where the constant in the error bound does not depend on inverse powers of the diffusion or viscosity, and non-robust ones. The classical convergence theory for the Galerkin discretization of the Navier–Stokes equations with inf–sup stable pairs of finite element spaces predicts an optimal order of convergence for the velocity error in L∞(L2), and with that for the kinetic energy, in a non-robust way. However, numerical simulations with non weakly divergence-free pairs of finite element spaces reveal that in the situation of small viscosity coefficients and coarse grids, the solution obtained with the Galerkin method is polluted with spurious oscillations, even if the analytic solution is smooth, compare Example 4.4. This is not the case for methods for which robust estimates can be derived. For such methods, one obtains meaningful numerical solutions in this situation and even a certain order of error reduction, which is the 48 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 Table 5 Discretizations of the incompressible Navier–Stokes equations for which robust estimates for the velocity error in L∞(L2) are proved. Method Finite element Order Analysis Section Galerkin inf–sup stable, EMAC form of the nonlinear term r−1 [53] + Theorem 4.1 Section 4.1.3 Galerkin inf–sup stable, weakly divergence-free r[57] Section 4.2.1 Galerkin H(div)-conforming, inf–sup stable, weakly divergence-free r[62,63] Section 4.2.1 grad–div stabilization inf–sup stable r[50] Section 4.2.2 LPS for velocity gradient inf–sup stable r[70] Section 4.2.2 grad–div stabilization, LPS for pressure non inf–sup stable r[72] Section 4.2.3 LPS for velocity gradient, LPS for pressure non inf–sup stable r[72] Section 4.2.3 symmetric pressure stabilization non inf–sup stable r[83,84] Section 4.2.3 SUPG P1/P1, space–time finite elements 1.5 [13] Section 4.3 CIP non inf–sup stable r+1/2 [90] Section 4.3 LPS for velocity gradient, LPS for pressure, modified parameter (two types of methods) non inf–sup stable r+1/2 [72,96] Section 4.3 H(div)-conforming DG with upwind inf–sup stable r+1/2 [64] Section 4.3 notion used in this survey for the behavior of a method on coarse grids, can be observed. Thus, in the case of small viscosity coefficients and coarse grids, only robust estimates provide useful information about the behavior of a numerical method on coarse grids if the analytic solution is smooth. In this respect, robust error bounds can be considered as the necessary condition (although perhaps not sufficient condition) for a good practical performance in more realistic problems, where the solutions are unlikely to have bounded derivatives as ν→0. For scalar linear convection–diffusion equations, robust estimates of order r+1/2 can be proved for many popular stabilized methods. Sometimes, there are restrictions in the analysis, e.g., that the coefficients of the equation do not depend on time. There are special situations, like in Example 3.4 for the SUPG method and P1finite elements, where an optimal order of convergence can be observed in numerical simulations. The situation is more complex for the incompressible Navier–Stokes equations due to the nonlinear convective term and the interaction of velocity and pressure. Robust estimates could be proved so far only for methods whose structure is relatively simple compared with methods that are actually used for turbulent flow simulations. These more complicated methods might contain stabilization mechanisms already used in the simple methods, like the grad–div and SUPG stabilizations in the residual-based VMS method. An overview of discretizations for which robust estimates for the velocity error in L∞(L2) are known is provided in Table 5. For the Galerkin method and inf–sup stable pairs of finite element spaces that are not weakly divergence-free, only the EMAC form of the nonlinear term allows a robust estimate whose order is r−1. The analysis carried out in the present paper reveals how the different terms in the Navier–Stokes equations affect this order reduction. In particular, the interaction between velocity and pressure is responsible for the order reduction to r−1, unless stabilization terms are added, while the nonlinear term is responsible for the lack of convergence, unless the EMAC form is used or stabilization terms are added. For other methods using inf–sup stable pairs of finite element spaces, robust estimates of order rfor the velocity error in L∞(L2) are proved. Note that for many pairs of such finite element spaces, the velocity space contains piecewise polynomials of degree rand the pressure space of degree r−1. From the point of view of physical consistency (mass conservation, fundamental invariance property, pressure robustness), the use of weakly divergence-free pairs is attractive. However, their implementation might be somewhat involved. The question about the possibility of yielding r+1/2 convergence for these methods adding some kind of stabilization remains open. For non weakly divergence-free pairs, the use of the grad–div stabilization is an appealing approach, because grad– div stabilization is easily to implement. Implementing the LPS method for the velocity gradient is certainly more 49 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 involved. Altogether, the best order of convergence of the velocity error in L∞(L2) that could be proved so far with a robust estimate for a method with inf–sup stable pairs of finite element spaces and using continuous velocity elements is r. This is one order less than the best approximation error for the velocity with respect to L2(Ω)d. Consequently, the convergence order of the finite element error with respect to the kinetic energy is reduced by one compared with the optimal order. Concerning inf–sup stable pairs of finite element spaces with H(div)-conforming discontinuous velocity finite elements, a robust estimate of order r+1/2 could be proved for an upwind discretization of the convective term. There is a larger variety of methods where robust error bounds have been proved for non inf–sup stable pairs of finite element spaces. Such pairs are often so-called equal order pairs, where both the finite element velocity and pressure consist of piecewise polynomials of degree r. The stabilization mechanisms that lead for inf–sup stable pairs to methods of order rgive also methods of at least order rfor equal order pairs. The parameter choice in the LPS method for the velocity gradient can be modified such that even order r+1/2 is reached, which can be proven even for two different types of LPS approaches. There are two more methods with this order, the CIP method and the SUPG method, where for the latter only analysis for r=1 is available. Thus, with equal order pairs of finite element spaces it is possible to achieve order r+1/2 for the convergence of the velocity error in L∞(L2) with robust error bounds, which is half an order less than the optimal order, and so it is the case of the kinetic energy. By far most of the finite element error analysis that can be found in the literature is for H1-conforming methods. Nevertheless, the analysis is not complete, e.g., for the SUPG method only a special case is covered and there are many methods, discussed in Section 4.4, where robust error bounds are not yet available. A generalization of H(div)-conforming discontinuous methods are so-called hybrid or hybridizable discontinuous Galerkin (HDG) methods. Different formulations for such methods have been derived. However, a finite element analysis for the time-dependent Navier–Stokes equations does not seem to be available, although these methods are used for the simulation of turbulent flows, e.g., in [121]. Numerical analysis of HDG methods for the Stokes equations can be found, e.g., in [61,122]. The state of the art is that there is no method for which robust estimates for the velocity error in L∞(L2), and with that for the error in the kinetic energy, of optimal order can be proved so far. This survey presented numerical results for several methods from Table 5 and all of them showed that the analytic prediction of the order of error reduction is sharp. An optimal order could not be observed for any of these methods. We performed with these methods also simulations on a family of triangulations of Friedrichs–Keller type, like for the SUPG method for convection–diffusion equations in Fig. 4. For brevity, the results are not presented here. Even on these regular grids, none of the methods showed optimal error reduction of the velocity error in L∞(L2). Thus, one can conclude that among the methods for which numerical studies are presented in this survey, there is no method with optimal error reduction, even in special situations like regular grids. In summary, whether or not there is a method with an optimal robust error bound for the velocity error in L∞(L2), perhaps at least in special situations, remains an open problem. Another widely open problem is the question of local error estimates, that is, to obtain (robust) error bounds in subsets of the domain away from layers. For time-dependent problems, only in [4,5] such estimates can be found, for the SUPG method applied to linear evolutionary convection–reaction–diffusion equations. Further open problems, like the analysis of the considered equations on anisotropic grids or numerical analysis for temporal or spatio-temporal averages of errors are already stated in [17]. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments The research of Bosco Garc´ ıa-Archilla has been supported by Spanish MICIU under grant PGC2018-096265- B-I00 and MICI under grant PID2019-104141GB-I00. The Research of Julia Novo has been supported by Spanish MICI under grants PID2019-104141GB-I00 and VA169P20 (Junta de Castilla y Leon, ES) cofinanced by FEDER funds. 50 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 References [1] Hermann Sohr, The navier-stokes equations, in: Birkhäuser Advanced Texts: Basler Lehrbücher. [Birkhäuser Advanced Texts: Basel Textbooks], Birkhäuser Verlag, Basel, 2001, p. x+367, An elementary functional analytic approach. [2] Volker John, Finite Element Methods for Incompressible Flow Problems, in: Springer Series in Computational Mathematics, vol. 51, Springer, Cham, 2016, p. xiii+812. [3] Eberhard Hopf, Über die Anfangswertaufgabe für die hydrodynamischen grundgleichungen, Math. Nachr. 4 (1951) 213–231. [4] U. Nävert, A Finite Element Method for Convection-Diffusion Problems (Ph.D. thesis), Chalmers University of Technology, Göteborg, 1982. [5] Javier de Frutos, Bosco García-Archilla, Julia Novo, Local error estimates for the SUPG method applied to evolutionary convection-reaction-diffusion equations, J. Sci. Comput. 66 (2) (2016) 528–554. [6] Clemens Bartsch, Viktoria Wiedmeyer, Zahra Lakdawala, Robert I.A. Patterson, Andreas Voigt, Kai Sundmacher, Volker John, Stochastic-deterministic population balance modeling and simulation of a fluidized bed crystallizer experiment, Chem. Eng. Sci. 208 (2019) 115102. [7] T.J.R. Hughes, A. Brooks, A multidimensional upwind scheme with no crosswind diffusion, in: Finite Element Methods for Convection Dominated Flows (Papers, Winter Ann. Meeting Amer. Soc. Mech. Engrs., New York, 1979), in: AMD, vol. 34, Amer. Soc. Mech. Engrs. (ASME), New York, 1979, pp. 19–35. [8] Alexander N. Brooks, Thomas J.R. Hughes, Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible navier-Stokes equations, in: FENOMECH ’81, Part I (Stuttgart, 1981), Comput. Methods Appl. Mech. Engrg. 32 (1–3) (1982) 199–259. [9] Claes Johnson, Uno Nävert, An analysis of some finite element methods for advection-diffusion problems, in: Analytical and Numerical Approaches to Asymptotic Problems in Analysis (Proc. Conf., Univ. Nijmegen, Nijmegen, 1980), in: North-Holland Math. Stud, vol. 47, North-Holland, Amsterdam-New York, 1981, pp. 99–116. [10] G. Lube, L. Tobiska, A nonconforming finite element method of streamline diffusion type for the incompressible Navier-Stokes equations, J. Comput. Math. 8 (2) (1990) 147–158. [11] Lutz Tobiska, Gert Lube, A modified streamline diffusion method for solving the stationary Navier-Stokes equation, Numer. Math. 59 (1) (1991) 13–29. [12] Claes Johnson, Jukka Saranen, Streamline diffusion methods for the incompressible euler and Navier-Stokes equations, Math. Comp. 47 (175) (1986) 1–18. [13] Peter Hansbo, Anders Szepessy, A velocity-pressure streamline diffusion finite element method for the incompressible Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 84 (2) (1990) 175–192. [14] Y. Bazilevs, V.M. Calo, J.A. Cottrell, T.J.R. Hughes, A. Reali, G. Scovazzi, Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows, Comput. Methods Appl. Mech. Engrg. 197 (1–4) (2007) 173–201. [15] Hans-Görg Roos, Martin Stynes, Lutz Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations, second ed., in: Springer Series in Computational Mathematics, vol. 24, Springer-Verlag, Berlin, 2008, p. xiv+604, Convection-diffusion-reaction and flow problems. [16] Stephen B. Pope, Turbulent Flows, Cambridge University Press, Cambridge, 2000, p. xxxiv+771. [17] Volker John, Petr Knobloch, Julia Novo, Finite elements for scalar convection-dominated equations and incompressible flow problems: A never ending story? Comput. Vis. Sci. 19 (5–6) (2018) 47–63. [18] Philippe G. Ciarlet, The Finite Element Method for Elliptic Problems, in: Studies in Mathematics and its Applications, vol. 4, North-Holland Publishing Co., Amsterdam-New York-Oxford, 1978, p. xix+530. [19] Daniel Arndt, Helene Dallmann, Gert Lube, Local projection FEM stabilization for the time-dependent incompressible Navier-Stokes problem, Numer. Methods Partial Differential Equations 31 (4) (2015) 1224–1250. [20] Javier de Frutos, Bosco García-Archilla, Julia Novo, Stabilization of Galerkin finite element approximations to transient convection-diffusion problems, SIAM J. Numer. Anal. 48 (3) (2010) 953–979. [21] Timothy A. Davis, Algorithm 832: UMFPACK V4.3—An unsymmetric-pattern multifrontal method, ACM Trans. Math. Softw. 30 (2) (2004) 196–199. [22] Volker John, Ellen Schmeyer, Finite element methods for time-dependent convection-diffusion-reaction equations with small diffusion, Comput. Methods Appl. Mech. Engrg. 198 (3–4) (2008) 475–494. [23] Volker John, Julia Novo, On (essentially) non-oscillatory discretizations of evolutionary convection-diffusion equations, J. Comput. Phys. 231 (4) (2012) 1570–1586. [24] Volker John, Julia Novo, Error analysis of the SUPG finite element discretization of evolutionary convection-diffusion-reaction equations, SIAM J. Numer. Anal. 49 (3) (2011) 1149–1176. [25] Erik Burman, Giles Smith, Analysis of the space semi-discretized SUPG method for transient convection-diffusion equations, Math. Models Methods Appl. Sci. 21 (10) (2011) 2049–2068. [26] G. Lube, D. Weiss, Stabilized finite element methods for singularly perturbed parabolic problems, Appl. Numer. Math. 17 (4) (1995) 431–459. [27] Ramon Codina, Jordi Blasco, Analysis of a stabilized finite element approximation of the transient convection-diffusion-reaction equation using orthogonal subscales, Comput. Vis. Sci. 4 (3) (2002) 167–174. [28] N. Ahmed, G. Matthies, L. Tobiska, H. Xie, Discontinuous Galerkin time stepping with local projection stabilization for transient convection-diffusion-reaction problems, Comput. Methods Appl. Mech. Engrg. 200 (21–22) (2011) 1747–1756. [29] Naveed Ahmed, Gunar Matthies, Higher order continuous Galerkin-Petrov time stepping schemes for transient convection-diffusion- reaction equations, ESAIM Math. Model. Numer. Anal. 49 (5) (2015) 1429–1450. 51 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 [30] Gabriel R. Barrenechea, Volker John, Petr Knobloch, A local projection stabilization finite element method with nonlinear crosswind diffusion for convection-diffusion-reaction equations, ESAIM Math. Model. Numer. Anal. 47 (5) (2013) 1335–1366. [31] Erik Burman, Miguel A. Fernández, Finite element methods with symmetric stabilization for the transient convection-diffusion-reaction equation, Comput. Methods Appl. Mech. Engrg. 198 (33–36) (2009) 2508–2519. [32] Jean-Luc Guermond, Subgrid stabilization of Galerkin approximations of linear contraction semi-groups of class C0in Hilbert spaces, Numer. Methods Partial Differential Equations 17 (1) (2001) 1–25. [33] Thomas J.R. Hughes, Multiscale phenomena: Green’s functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods, Comput. Methods Appl. Mech. Engrg. 127 (1–4) (1995) 387–401. [34] M. Feistauer, K. Švadlenka, Discontinuous Galerkin method of lines for solving nonstationary singularly perturbed linear problems, J. Numer. Math. 12 (2) (2004) 97–117. [35] Vít Dolejší, Miloslav Feistauer, Discontinuous Galerkin Method, in: Springer Series in Computational Mathematics, vol. 48, Springer, Cham, 2015, p. xiv+572, Analysis and applications to compressible flow. [36] J. Gopalakrishnan, G. Kanschat, A multilevel discontinuous Galerkin method, Numer. Math. 95 (3) (2003) 527–550. [37] Erik Burman, Consistent SUPG-method for transient transport problems: Stability and convergence, Comput. Methods Appl. Mech. Engrg. 199 (17–20) (2010) 1114–1123. [38] Santiago Badia, Alba Hierro, On monotonicity-preserving stabilized finite element approximations of transport problems, SIAM J. Sci. Comput. 36 (6) (2014) A2673–A2697. [39] Hans-Görg Roos, Martin Stynes, Some open questions in the numerical analysis of singularly perturbed differential equations, Comput. Methods Appl. Math. 15 (4) (2015) 531–550. [40] Guohui Zhou, How accurate is the streamline diffusion finite element method? Math. Comp. 66 (217) (1997) 31–44. [41] John G. Heywood, Rolf Rannacher, Finite element approximation of the nonstationary Navier-Stokes problem. I. Regularity of solutions and second-order error estimates for spatial discretization, SIAM J. Numer. Anal. 19 (2) (1982) 275–311. [42] John G. Heywood, Rolf Rannacher, Finite element approximation of the nonstationary Navier-Stokes problem. II. Stability of solutions and error estimates uniform in time, SIAM J. Numer. Anal. 23 (4) (1986) 750–777. [43] John G. Heywood, Rolf Rannacher, Finite element approximation of the nonstationary Navier-Stokes problem. III. Smoothing property and higher order error estimates for spatial discretization, SIAM J. Numer. Anal. 25 (3) (1988) 489–512. [44] John G. Heywood, Rolf Rannacher, Finite-element approximation of the nonstationary Navier-Stokes problem. IV. Error analysis for second-order time discretization, SIAM J. Numer. Anal. 27 (2) (1990) 353–384. [45] D.N. Arnold, F. Brezzi, M. Fortin, A stable finite element for the Stokes equations, Calcolo 21 (4) (1984) 337–344 (1985). [46] P. Hood, C. Taylor, Navier–Stokes equations using mixed interpolation, in: J.T. Oden, R.H. Gallagher, O.C. Zienkiewicz, C. Taylor (Eds.), Finite Element Methods in Flow Problems, University of Alabama in Huntsville Press, 1974, pp. 121–132. [47] Volker John, Alexander Linke, Christian Merdon, Michael Neilan, Leo G. Rebholz, On the divergence constraint in mixed finite element methods for incompressible flows, SIAM Rev. 59 (3) (2017) 492–544. [48] Javier de Frutos, Bosco García-Archilla, Volker John, Julia Novo, Grad-div stabilization for the evolutionary Oseen problem with inf-sup stable finite elements, J. Sci. Comput. 66 (3) (2016) 991–1024. [49] Bosco García-Archilla, Julia Novo, Edriss S. Titi, Uniform in time error estimates for a finite element method applied to a downscaling data assimilation algorithm for the Navier-Stokes equations, SIAM J. Numer. Anal. 58 (1) (2020) 410–429. [50] Javier de Frutos, Bosco García-Archilla, Volker John, Julia Novo, Analysis of the grad-div stabilization for the time-dependent Navier-Stokes equations with inf-sup stable finite elements, Adv. Comput. Math. 44 (1) (2018) 195–225. [51] Hongsen Chen, Pointwise error estimates for finite element solutions of the Stokes problem, SIAM J. Numer. Anal. 44 (1) (2006) 1–28. [52] Sergey Charnyi, Timo Heister, Maxim A. Olshanskii, Leo G. Rebholz, On conservation laws of Navier-Stokes Galerkin discretizations, J. Comput. Phys. 337 (2017) 289–308. [53] Maxim A. Olshanskii, Leo G. Rebholz, Longer time accuracy for incompressible Navier-Stokes simulations with the EMAC formulation, Comput. Methods Appl. Mech. Engrg. 372 (2020) 113369, 17. [54] L.R. Scott, M. Vogelius, Conforming finite element methods for incompressible and nearly incompressible continua, in: Large-Scale Computations in Fluid Mechanics, Part 2 (La Jolla, Calif., 1983), in: Lectures in Appl. Math, vol. 22, Amer. Math. Soc., Providence, RI, 1985, pp. 221–244. [55] Jinshui Qin, On the Convergence of Some Low Order Mixed Finite Elements for Incompressible Fluids (Ph.D. thesis), Department of Mathematics, Pennsylvania State University, 1994. [56] Shangyou Zhang, A new family of stable mixed finite elements for the 3D Stokes equations, Math. Comp. 74 (250) (2005) 543–554. [57] Philipp W. Schroeder, Gert Lube, Pressure-robust analysis of divergence-free and conforming FEM for evolutionary incompressible Navier-Stokes flows, J. Numer. Math. 25 (4) (2017) 249–276. [58] Benjamin R. Cousins, Leo G. Rebholz, Nicholas E. Wilson, Enforcing energy, helicity and strong mass conservation in finite element computations for incompressible Navier-Stokes simulations, Appl. Math. Comput. 218 (4) (2011) 1208–1221. [59] Philipp W. Schroeder, Volker John, Philip L. Lederer, Christoph Lehrenfeld, Gert Lube, Joachim Schöberl, On reference solutions and the sensitivity of the 2D Kelvin-Helmholtz instability problem, Comput. Math. Appl. 77 (4) (2019) 1010–1028. [60] Bernardo Cockburn, Guido Kanschat, Dominik Schötzau, A note on discontinuous Galerkin divergence-free solutions of the Navier-Stokes equations, J. Sci. Comput. 31 (1–2) (2007) 61–73. [61] Philip L. Lederer, Christoph Lehrenfeld, Joachim Schöberl, Hybrid discontinuous Galerkin methods with relaxed H(div)-conformity for incompressible flows. Part II, ESAIM Math. Model. Numer. Anal. 53 (2) (2019) 503–522. [62] Philipp W. Schroeder, Gert Lube, Divergence-free H(div)-FEM for time-dependent incompressible flows with applications to high Reynolds number vortex dynamics, J. Sci. Comput. 75 (2) (2018) 830–858. 52 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 [63] Philipp W. Schroeder, Christoph Lehrenfeld, Alexander Linke, Gert Lube, Towards computable flows and robust estimates for inf-sup stable FEM applied to the time-dependent incompressible Navier-Stokes equations, SeMA J. 75 (4) (2018) 629–653. [64] Yongbin Han, Yanren Hou, Robust error analysis of H(div)-conforming DG method for the time-dependent incompressible Navier-Stokes equations, J. Comput. Appl. Math. 390 (2021) 113365, 13. [65] Eleanor W. Jenkins, Volker John, Alexander Linke, Leo G. Rebholz, On the parameter choice in grad-div stabilization for the Stokes equations, Adv. Comput. Math. 40 (2) (2014) 491–516. [66] Leopoldo P. Franca, Thomas J.R. Hughes, Two classes of mixed finite element methods, Comput. Methods Appl. Mech. Engrg. 69 (1) (1988) 89–129. [67] Maxim Olshanskii, Gert Lube, Timo Heister, Johannes Löwe, Grad-div stabilization and subgrid pressure models for the incompressible Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 198 (49–52) (2009) 3975–3988. [68] Blanca Ayuso, Bosco García-Archilla, Julia Novo, The postprocessed mixed finite-element method for the Navier-Stokes equations, SIAM J. Numer. Anal. 43 (3) (2005) 1091–1111. [69] Michael A. Case, Vincent J. Ervin, Alexander Linke, Leo G. Rebholz, A connection between Scott-Vogelius and grad-div stabilized Taylor-Hood FE approximations of the Navier-Stokes equations, SIAM J. Numer. Anal. 49 (4) (2011) 1461–1481. [70] Naveed Ahmed, Gunar Matthies, Higher-order discontinuous Galerkin time discretizations for the evolutionary Navier–Stokes equations, IMA J. Numer. Anal. (2021) in press. [71] Volker John, Petr Knobloch, Ulrich Wilbrandt, Finite element pressure stabilizations for incompressible flow problems, in: T. Bodnár, G.P. Galdi, Š. Neˇ casová (Eds.), Fluids Under Pressure, in: Adv. Math. Fluid Mech, Birkhäuser, 2020, pp. 483–573. [72] Javier de Frutos, Bosco García-Archilla, Volker John, Julia Novo, Error analysis of non inf-sup stable discretizations of the time-dependent Navier-Stokes equations with local projection stabilization, IMA J. Numer. Anal. 39 (4) (2019) 1747–1786. [73] T. Chacón Rebollo, M. Gómez Mármol, V. Girault, I. Sánchez Muñoz, A high order term-by-term stabilization solver for incompressible flow problems, IMA J. Numer. Anal. 33 (3) (2013) 974–1007. [74] T. Chacón Rebollo, M. Gómez Mármol, M. Restelli, Numerical analysis of penalty stabilized finite element discretizations of evolution Navier-Stokes equations, J. Sci. Comput. 63 (3) (2015) 885–912. [75] C. Bernardi, V. Girault, A local regularization operator for triangular and quadrilateral finite elements, SIAM J. Numer. Anal. 35 (5) (1998) 1893–1916. [76] L. Ridgway Scott, Shangyou Zhang, Finite element interpolation of nonsmooth functions satisfying boundary conditions, Math. Comp. 54 (190) (1990) 483–493. [77] Susanne C. Brenner, L. Ridgway Scott, The Mathematical Theory of Finite Element Methods, third ed., Texts in Applied Mathematics, vol. 15, Springer, New York, 2008, p. xviii+397. [78] Silvia Bertoluzza, The discrete commutator property of approximation spaces, C. R. Acad. Sci., Paris 329 (12) (1999) 1097–1102. [79] Volker Gravemeier, Michael W. Gee, Wolfgang A. Wall, An algebraic variational multiscale-multigrid method based on plain aggregation for convection-diffusion problems, Comput. Methods Appl. Mech. Engrg. 198 (47–48) (2009) 3821–3835. [80] Volker Gravemeier, Michael W. Gee, Martin Kronbichler, Wolfgang A. Wall, An algebraic variational multiscale-multigrid method for large eddy simulation of turbulent flow, Comput. Methods Appl. Mech. Engrg. 199 (13–16) (2010) 853–864. [81] Volker John, Songul Kaya, A finite element variational multiscale method for the Navier-Stokes equations, SIAM J. Sci. Comput. 26 (5) (2005) 1485–1503 (electronic). [82] Naveed Ahmed, Tomás Chacón Rebollo, Volker John, Samuele Rubino, A review of variational multiscale methods for the simulation of turbulent incompressible flows, Arch. Comput. Methods Eng. 24 (1) (2017) 115–164. [83] Bosco García-Archilla, Volker John, Julia Novo, Symmetric pressure stabilization for equal-order finite element approximations to the time-dependent Navier–Stokes equations, IMA J. Numer. Anal. 41 (2) (2021) 1093–1129. [84] Gang Chen, Minfu Feng, Analysis of solving Galerkin finite element methods with symmetric pressure stabilization for the unsteady Navier-Stokes equations using conforming equal order interpolation, Adv. Appl. Math. Mech. 9 (2) (2017) 362–377. [85] Erik Burman, Miguel A. Fernández, Galerkin finite element methods with symmetric pressure stabilization for the transient stokes equations: Stability and convergence analysis, SIAM J. Numer. Anal. 47 (1) (2008/09) 409–439. [86] F. Brezzi, J. Pitkäranta, On the stabilization of finite element approximations of the Stokes equations, in: Efficient Solutions of Elliptic Systems (Kiel, 1984), in: Notes Numer. Fluid Mech, vol. 10, Friedr. Vieweg, Braunschweig, 1984, pp. 11–19. [87] Clark R. Dohrmann, Pavel B. Bochev, A stabilized finite element method for the Stokes problem based on polynomial pressure projections, Internat. J. Numer. Methods Fluids 46 (2) (2004) 183–201. [88] Ramon Codina, Jordi Blasco, A finite element formulation for the Stokes problem allowing equal velocity-pressure interpolation, Comput. Methods Appl. Mech. Engrg. 143 (3–4) (1997) 373–391. [89] R. Becker, M. Braack, A finite element pressure gradient stabilization for the Stokes equations based on local projections, Calcolo 38 (4) (2001) 173–199. [90] Erik Burman, Miguel A. Fernández, Continuous interior penalty finite element method for the time-dependent Navier-stokes equations: Space discretization and convergence, Numer. Math. 107 (1) (2007) 39–77. [91] Alejandro Allendes, Gabriel Barrenechea, Julia Novo, A divergence-free low-order stabilized finite element method for the evolutionary Navier-Stokes equations, SIAM J. Sci. Comput. (2021) in press. [92] Gang Chen, Min-fu Feng, Yin-nian He, Finite difference streamline diffusion method using nonconforming space for incompressible time-dependent Navier-Stokes equations, Appl. Math. Mech. (English Ed.) 34 (9) (2013) 1083–1096. [93] Naveed Ahmed, Samuele Rubino, Numerical comparisons of finite element stabilized methods for a 2D vortex dynamics simulation at high Reynolds number, Comput. Methods Appl. Mech. Engrg. 349 (2019) 191–212. [94] Naveed Ahmed, Volker John, An assessment of two classes of variational multiscale methods for the simulation of incompressible turbulent flows, Comput. Methods Appl. Mech. Engrg. 365 (2020) 112997. 53 B. Garc´ ıa-Archilla, V. John and J. Novo Computer Methods in Applied Mechanics and Engineering 385 (2021) 114032 [95] Erik Burman, Alexandre Ern, Miguel A. Fernández, Fractional-step methods and finite elements with symmetric stabilization for the transient Oseen problem, ESAIM Math. Model. Numer. Anal. 51 (2) (2017) 487–507. [96] Gang Chen, Minfu Feng, Hong Zhou, Local projection stabilized method on unsteady Navier-Stokes equations with high Reynolds number using equal order interpolation, Appl. Math. Comput. 243 (2014) 465–481. [97] J. Smagorinsky, General circulation experiments with the primitive equations, Mon. Weather Rev. 91 (1963) 99–164. [98] O.A. Ladyženskaja, New equations for the description of the motions of viscous incompressible fluids, and global solvability for their boundary value problems, Tr. Mat. Inst. Steklova 102 (1967) 85–104. [99] V. John, W.J. Layton, Analysis of numerical errors in large eddy simulation, SIAM J. Numer. Anal. 40 (3) (2002) 995–1020. [100] Tomás Chacón Rebollo, Roger Lewandowski, Mathematical and numerical foundations of turbulence models and applications, in: Modeling and Simulation in Science, Engineering and Technology, Birkhäuser/Springer, New York, 2014, p. xviii+517. [101] Erik Burman, Peter Hansbo, Mats G. Larson, Error estimates for the Smagorinsky turbulence model: Enhanced stability through scale separation and numerical stabilization, 2021, arXiv:2102.00043v1. [102] Snorre H. Christiansen, Kaibo Hu, Generalized finite element systems for smooth differential forms and Stokes’ problem, Numer. Math. 140 (2) (2018) 327–371. [103] Vincent J. Ervin, William J. Layton, Monika Neda, Numerical analysis of filter-based stabilization for evolution equations, SIAM J. Numer. Anal. 50 (5) (2012) 2307–2335. [104] William Layton, Leo G. Rebholz, Catalin Trenchea, Modular nonlinear filter stabilization of methods for higher Reynolds numbers flow, J. Math. Fluid Mech. 14 (2) (2012) 325–354. [105] William Layton, Nathaniel Mays, Monika Neda, Catalin Trenchea, Numerical analysis of modular regularization methods for the BDF2 time discretization of the Navier-Stokes equations, ESAIM Math. Model. Numer. Anal. 48 (3) (2014) 765–793. [106] Jean-Luc Guermond, Stabilization of Galerkin approximations of transport equations by subgrid modeling, M2AN Math. Model. Numer. Anal. 33 (6) (1999) 1293–1316. [107] Ramon Codina, Santiago Badia, Joan Baiges, Javier Principe, Variational multiscale methods in computational fluid dynamics, in: E. Stein, R. Borst, T.J.R. Hughes (Eds.), Encyclopedia of Computational Mechanics, second ed., American Cancer Society, 2017, pp. 1–28. [108] Ursula Rasthofer, Volker Gravemeier, Recent developments in variational multiscale methods for large-eddy simulation of turbulent flow, Arch. Comput. Methods Eng. 25 (3) (2018) 647–690. [109] M. Braack, E. Burman, Local projection stabilization for the Oseen problem and its interpretation as a variational multiscale method, SIAM J. Numer. Anal. 43 (6) (2006) 2544–2566 (electronic). [110] W. Layton, A connection between subgrid scale eddy viscosity and mixed methods, Appl. Math. Comput. 133 (1) (2002) 147–157. [111] Volker John, Songul Kaya, Finite element error analysis of a variational multiscale method for the Navier-Stokes equations, Adv. Comput. Math. 28 (1) (2008) 43–61. [112] Lars Röhe, Gert Lube, Analysis of a variational multiscale method for large-eddy simulation and its application to homogeneous isotropic turbulence, Comput. Methods Appl. Mech. Engrg. 199 (37–40) (2010) 2331–2342. [113] William Layton, Lars Röhe, Hoang Tran, Explicitly uncoupled VMS stabilization of fluid flow, Comput. Methods Appl. Mech. Engrg. 200 (45–46) (2011) 3183–3199. [114] Li Shan, Yanren Hou, Haibiao Zheng, Variational multiscale method based on the Crank-Nicolson extrapolation scheme for the non-stationary Navier-Stokes equations, Int. J. Comput. Math. 89 (16) (2012) 2198–2223. [115] Naveed Ahmed, Tomás Chacón Rebollo, Volker John, Samuele Rubino, Analysis of a full space-time discretization of the navier-Stokes equations by a local projection stabilization method, IMA J. Numer. Anal. 37 (3) (2017) 1437–1467. [116] Abigail L. Bowers, Sabine Le Borne, Leo G. Rebholz, Error analysis and iterative solvers for Navier-Stokes projection methods with standard and sparse grad-div stabilization, Comput. Methods Appl. Mech. Engrg. 275 (2014) 1–19. [117] J.A. Fiordilino, W. Layton, Y. Rong, An efficient and modular grad-div stabilization, Comput. Methods Appl. Mech. Engrg. 335 (2018) 327–346. [118] Y. Rong, J.A. Fiordilino, Numerical analysis of a BDF2 modular grad–div stabilization method for the Navier–Stokes equations, J. Sci. Comput. 82 (3) (2020). [119] Alexandr Labovsky, William J. Layton, Carolina C. Manica, Monika Neda, Leo G. Rebholz, The stabilized extrapolated trapezoidal finite-element method for the Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 198 (9–12) (2009) 958–974. [120] Nan Jiang, Muhammad Mohebujjaman, Leo G. Rebholz, Catalin Trenchea, An optimally accurate discrete regularization for second order timestepping methods for Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 310 (2016) 388–405. [121] Niklas Fehn, Martin Kronbichler, Christoph Lehrenfeld, Gert Lube, Philipp W. Schroeder, High-order DG solvers for underresolved turbulent incompressible flows: A comparison of L2and H(div) methods, Internat. J. Numer. Methods Fluids 91 (11) (2019) 533–556. [122] Sander Rhebergen, Garth N. Wells, Analysis of a hybridized/interface stabilized finite element method for the Stokes equations, SIAM J. Numer. Anal. 55 (4) (2017) 1982–2003. 54