# Grad-div stabilized discretizations on S-type meshes for the Oseen problem

Grad-div stabilized discretizations on S-type meshes for the Oseen problem Abstract We consider discretizations of the singularly perturbed Oseen equations on properly layer-adapted meshes. Using a suitable solution decomposition, we are able to prove optimal convergence orders in the associated energy norm for grad-div stabilized finite element methods in a general setting. Two families of pairs of discrete function spaces, namely $$Q_k\times Q_{k-1}$$ and $$Q_k\times P_{k-1}^{\text{disc}}$$, $$k\ge 2$$, are investigated in detail. The usage of a standard nonstabilized Galerkin method reduces the order by 1 while stabilization outside the layers is enough to regain the full optimal order. 1. Introduction The numerical solution of the time-dependent Navier–Stokes equations, especially for high Reynolds numbers, is still a great challenge of computational fluid dynamics. Using time discretization and linearization, the Oseen problem arises as an important subproblem that is for high Reynolds numbers singularly perturbed. It is well known that standard finite element methods applied to the Oseen problem may suffer from two difficulties: the generally dominating convection and the incompatibility between the approximation spaces for velocity and pressure. The streamline-upwind Petrov–Galerkin method (SUPG), introduced in Brooks & Hughes (1982), and the pressure-stabilized Petrov–Galerkin method (PSPG), introduced in Hughes et al. (1986), Johnson & Saranen (1986), Tezduyar et al. (1992), provide a unique frame to tackle both problems simultaneously. An additional elementwise stabilization of the divergence constraint, abbreviated as grad-div stabilization, increases the robustness; see Franca & Frey (1992), Hansbo & Szepessy (1990), Tobiska & Lube (1991). In order to ensure the strong consistency of these residual-based stabilization methods, various terms have to be added to the weak formulation. The PSPG term can be skipped if inf–sup stable pairs of finite element spaces for approximating velocity and pressure are used. This results in reduced stabilized schemes; see Gelhard et al. (2005), Matthies et al. (2009). The coupling between velocity and pressure that is caused by the SUPG term complicates the analysis and the grad-div stabilization seems to be even more important; compare Burman & Linke (2008), Gelhard et al. (2005), Linke & Rebholz (2013), Jenkins et al. (2014), Matthies et al. (2009), Olshanskii & Reusken (2004). On general meshes, the use of stabilization techniques provides semirobust error estimates in the sense of Roos (2012). This means that the constants in the error bounds are independent of the small viscosity parameter while the higher-order norms of the solution that appear strongly increase with decreasing viscosity. The role of grad-div stabilization has been investigated by several authors. The combination of SUPG and grad-div stabilization applied to the stationary Navier–Stokes equations was considered in Olshanskii (2002). The grad-div term enhances the accuracy of the numerical solution and the convergence properties of iterative solvers if the parameter of the grad-div stabilization is chosen properly. The influence of grad-div stabilization for the Stokes equations has been studied in Olshanskii & Reusken (2004), Jenkins et al. (2014). Within the scope of variational multiscale methods, the grad-div stabilization as a pressure subgrid model has been investigated in Olshanskii et al. (2009). The importance of the grad-div stabilization inside local projection stabilizations was stated in Dallmann et al. (2016). Standard methods for stationary singularly perturbed convection–diffusion problems do also lack stability; see, e.g. the book Roos et al. (2008). For these problems, knowledge of the layer structure of their solutions allows creation and usage of properly layer-adapted meshes; see the survey book Linß (2010). Such special meshes have quite a long history. First ideas were proposed by Bakhvalov (1969) using exponentially fitted boundary meshes. The piecewise uniform Shishkin meshes (Miller et al., 1996; Stynes & O’Riordan, 1997) are easier in structure and in the analysis of methods based on them, but their approximation quality is inferior to Bakhvalov meshes. The combination of Bakhvalov’s graded fine meshes with Shishkin’s choice of the transition point (Linß, 1999, 2000) led to general S-type meshes (Roos & Linß, 1999) where the mesh inside the layers is defined with the help of a mesh-generating function. The use of layer-adapted meshes allows for robust or uniform error estimates in the sense of Roos (2012). Hence, the error bounds are completely independent of the viscosity parameter. On such adapted meshes, the oscillations induced by boundary conditions are strongly reduced. However, the lack of stability is still there. Therefore, stabilization methods like SUPG (also called SDFEM), local projection stabilization or interior penalty stabilization are used in addition. Some results for those combinations can be found in Roos et al. (2008), Franz et al. (2008), Franz (2008), Franz & Matthies (2010), Stynes & Tobiska (2008). This article now combines these two approaches. We apply a grad-div stabilization on layer-adapted meshes. Using suitable pairs of function spaces for approximating velocity and pressure, we prove error bounds that are uniform in $$\varepsilon$$ and yield an optimal convergence rate. Stabilized finite element methods for the Oseen problem on anisotropic meshes have been analysed in Apel et al. (2008). The authors considered equal-order interpolation for velocity and pressure, but S-type meshes were explicitly excluded. We consider the stationary Oseen equations   −εΔu−(b⋅∇)u+cu+∇p =f in Ω=(0,1)2,divu =0 in Ω,u =0 on Γ=∂Ω, (1.1) where the perturbation parameter $$\varepsilon$$ with $$0<\varepsilon\ll 1$$ corresponds to the viscosity, the convection field $${\boldsymbol{b}}=(b_1,b_2)\in W^{1,\infty}({\it{\Omega}})^2$$ with $${\text{div}}{\boldsymbol{b}}=0$$ fulfils $$b_1\geq \beta_1 >0,\ b_2\geq \beta_2 >0$$ with some positive constants $$\beta_1,\beta_2$$ and $$0<c_0\leq c\in L^\infty({\it{\Omega}})$$ is assumed. Then $${\boldsymbol{u}}$$ describes the velocity of our problem and $$p$$ the pressure. Since $$0<\varepsilon\ll 1$$, we expect layers at $$x=0$$ and $$y=0$$. A typical solution of (1.1) can be seen in Fig. 1. Fig. 1. View largeDownload slide Typical solution of (1.1) with exponential layers in the velocity. Fig. 1. View largeDownload slide Typical solution of (1.1) with exponential layers in the velocity. As the solution exhibits boundary layers, the use of quasi-uniform meshes gives accurate approximations of $$(1.1)$$ only if the mesh size is of the order of the perturbation parameter $$\varepsilon$$. Since this gives a prohibitive restriction, layer-adapted meshes will be used for the discretizations. We are aware that the considered problem class is not a prototype for Oseen problems coming from applications as we utilize a solution decomposition (see Assumption 2.1) which is usually unavailable in real-world problems. However, this article is intended as a starting point for the robust numerical analysis of such problems on layer-adapted meshes. The article is organized as follows. Section 2 introduces a weak formulation, a solution decomposition and layer-adapted meshes. In Section 3, we discretize and analyse our problem in an abstract setting without specifying the finite element spaces. We prove optimal convergence orders for the fully grad-div stabilized method and for stabilization only outside the layers while the order is reduced by 1 for the standard Galerkin method. We present in Section 4, two pairs of spaces that will suit our purpose. Finally, Section 5 presents numerical simulations confirming the theoretical results. Notation: On a domain $$D,$$ we denote the $$L^2$$-norm by $$\|\cdot\|_D$$ and the $$L^\infty$$-norm by $$\|\cdot\|_{\infty,D}$$. In the case of $$D={\it{\Omega}},$$ we drop the reference to the domain from the notation. Moreover, $$C$$ will denote a generic constant independent of the perturbation parameter $${\varepsilon}$$ and the mesh parameter $$N$$. 2. Solution decomposition and layer-adapted meshes Let us start with a weak formulation of (1.1). We define $$V:=H_0^1({\it{\Omega}})^2$$ as the space of vector functions with components in the Sobolev space $$H_0^1({\it{\Omega}})$$ of $$L^2({\it{\Omega}})$$-functions with weak derivatives in $$L^2({\it{\Omega}})$$ and vanishing boundary traces, and $$Q:=L_0^2({\it{\Omega}})$$ as the space of $$L^2({\it{\Omega}})$$-functions with vanishing integral mean value. A weak formulation then reads as follows: Find $$({\boldsymbol{u}},p)\in V\times Q$$ such that for all $$({\boldsymbol{v}},q)\in V\times Q$$,   ε(∇u,∇v)−((b⋅∇)u,v)+(cu,v)−(p,divv) =(f,v),(q,divu) =0. (2.1) Note that $${\text{div}} {\boldsymbol{b}}=0$$ implies   ((b⋅∇)u,v)=−((b⋅∇)v,u)∀u,v∈V, (2.2) in particular,   ((b⋅∇)v,v)=0∀v∈V. (2.3) Thus, applying the Lax–Milgram lemma in the subspace of divergence-free functions, we establish the unique velocity field $${\boldsymbol{u}}$$. The unique pressure $$p\in Q$$ such that $$({\boldsymbol{u}},p)$$ solves (2.1) follows from the Babuška–Brezzi condition for the pair $$(V,Q)$$; see Girault & Raviart (1986, Ch. I). We will make some assumptions about a solution decomposition of $${\boldsymbol{u}}=(u_1,\,u_2)$$ and $$p$$ in order to carry out our analysis. We assume that $$u_1$$ and $$u_2$$ will consist of layer and regular components while the pressure $$p$$ will not show any layer behaviour. To be more precise, we assume the following bounds. Assumption 2.1 The solution $${\boldsymbol{u}}=(u_1,u_2)$$ can be decomposed as   u1 =g1+v1+w1+z1andu2=g2+v2+w2+z2, where we have for fixed $$k\in \mathbb{N}$$, all $$x,y\in [0,1]$$ and $$0\leq i+j \leq k+1$$ the pointwise estimates   |∂xi∂yjg1(x,y)| ≤C,|∂xi∂yjg2(x,y)| ≤C,|∂xi∂yjv1(x,y)| ≤Cε1−ie−β1x/ε,|∂xi∂yjv2(x,y)| ≤Cε−ie−β1x/ε,|∂xi∂yjw1(x,y)| ≤Cε−ie−β2y/ε,|∂xi∂yjw2(x,y)| ≤Cε1−ie−β2y/ε,|∂xi∂yjz1(x,y)| ≤Cε1−(i+j)e−β1x/εe−β2y/ε,|∂xi∂yjz2(x,y)| ≤Cε1−(i+j)e−β1x/εe−β2y/ε, where $$g_1,\,g_2$$ are the regular parts, $$v_1,\,v_2$$ the exponential boundary layers in the $$x$$-direction, $$w_1,\,w_2$$ the exponential boundary layers in the $$y$$-direction and $$z_1,\,z_2$$ the corner layers. The bounds   |∂xi∂yjp(x,y)|≤C,0≤i+j≤k, are assumed for the pressure $$p$$. Remark 2.2 The divergence constraint guarantees the existence of a scalar stream function $$\psi$$ such that $${\boldsymbol{u}}=(\partial_y \psi, -\partial_x \psi)$$. Inserting this into the Oseen equations and applying the (two-dimensional) $${\text{curl}}$$ operator leads to a fourth-order problem. In Franz et al. (2016), it was established that its solution has a principle structure that allows us to derive the assumptions stated above providing that enough compatibility conditions on the data of our problem are fulfilled. The nature of these additional assumptions is similar to those of convection–diffusion problems (Kellogg & Stynes, 2005, 2007). Furthermore, our numerical studies confirm the assumed decomposition. Concerning the pressure, our analysis does also work if the pressure has strong boundary layers similar to those of $$u_1$$ or $$u_2$$. Then we would estimate $$p$$ similarly to $${\boldsymbol{u}}$$ which is possible if we use the same layer-adapted mesh for all solution components. With Assumption 2.1 in mind, we are able to construct a priori adapted meshes. The solution $${\boldsymbol{u}}$$ has weak and strong layers. In order to resolve the strong layers, we use fine meshes along $$y=0$$ for $$u_1$$ and along $$x=0$$ for $$u_2$$. Weak layers are not visible at first sight. In the process of a numerical analysis, their unpleasant $${\varepsilon}$$-dependence comes into play via derivatives. Thus, the numerical analysis works best if the weak layers are also resolved. Therefore, we will use layer-adapted meshes that are fine at $$x=0$$ and at $$y=0$$ for both components $$u_1$$ and $$u_2$$. For simplicity, these meshes are also used for $$p$$ although we did assume that $$p$$ has no boundary layers. We will use so-called S-type meshes; see Roos & Linß (1999). They are constructed using transition points where the mesh changes between fine and coarse. To be more precise, let us assume the number $$N$$ of mesh intervals in each coordinate direction is even. Then we define the transition parameters   λx:=min{12,σεβ1ln⁡N}andλy:=min{12,σεβ2ln⁡N} with some user-chosen positive parameter $$\sigma$$. For the sake of mere simplicity in our subsequent analysis, we assume that   ε<min{β1,β2}2σln⁡N⇔N<exp⁡min{β1,β2}2σε such that $$\lambda_x$$ and $$\lambda_y$$ are smaller than $$1/2$$. This is the typical case for (1.1) as otherwise a standard analysis on an equidistant mesh suffices. Using the transition parameters $$\lambda_x$$ and $$\lambda_y$$, the domain $${\it{\Omega}}$$ is divided into the subdomains $${\it{\Omega}}_{11},{\it{\Omega}}_{21},{\it{\Omega}}_{12}$$ and $${\it{\Omega}}_{22}$$ according to   Ω11 :=[λx,1]×[λy,1],Ω12:=[0,λx]×[λy,1],Ω21 :=[λx,1]×[0,λy]andΩ22:=[0,λx]×[0,λy]; see Fig. 2. Fig. 2. View largeDownload slide Decomposition of $${\it{\Omega}}$$ into subregions. Fig. 2. View largeDownload slide Decomposition of $${\it{\Omega}}$$ into subregions. Note that $${\it{\Omega}}_{11}$$ covers the nonlayer region, $${\it{\Omega}}_{22}$$ the corner layers and $${\it{\Omega}}_{12}$$ and $${\it{\Omega}}_{21}$$ the two boundary layer regions. With the help of a mesh-generating function$$\phi$$ which is monotonically increasing with $$\phi(0)=0$$ and $$\phi(1/2)=\ln N$$, we construct the layer-adapted mesh as a tensor product mesh with nodes according to   xi :={σεβ1ϕ(iN),i=0,…,N/2,1−2(1−λx)(1−iN),i=N/2,…,N,  and   yj :={σεβ2ϕ(jN),j=0,…,N/2,1−2(1−λy)(1−jN),j=N/2,…,N.  Thus, the dimensions of cells in the layer regions depend on $$\phi$$ while the mesh is uniform outside. Related to the mesh-generating function $$\phi,$$ we define the mesh-characterizing function$$\psi$$ by   ϕ=−ln⁡ψwhich impliesψ=e−ϕ. The derivative of this function will be used to characterize the convergence behaviour on general S-type meshes. Some examples of S-type meshes from Roos & Linß (1999) are presented in Table 1. Table 1 Some examples of mesh-generating and mesh-characterizing functions of S-type meshes Name  $$\phi(t)$$  $$\max\phi'$$  $$\psi(t)$$  $$\max|\psi'|$$  Shishkin mesh  $$2t\ln N$$  $$2\ln N$$  $$N^{-2t}$$  $$2\ln N$$  Bakhvalov–Shishkin mesh  $$-\!\ln(1\!-\!2t(1\!-\!N^{-1}))$$  $$2N$$  $$1\!-\!2t(1\!-\!N^{-1})$$  2  Modified Bakhvalov–Shishkin mesh  $$\frac{t}{q-t},\,q=\frac{1}{2}(1+\frac{1}{\ln N})$$  $$3(\ln N)^2$$  $${\,\mathit{e}}^{-{t}/{q-t}}$$  $$3/(2q)\leq 3$$  Name  $$\phi(t)$$  $$\max\phi'$$  $$\psi(t)$$  $$\max|\psi'|$$  Shishkin mesh  $$2t\ln N$$  $$2\ln N$$  $$N^{-2t}$$  $$2\ln N$$  Bakhvalov–Shishkin mesh  $$-\!\ln(1\!-\!2t(1\!-\!N^{-1}))$$  $$2N$$  $$1\!-\!2t(1\!-\!N^{-1})$$  2  Modified Bakhvalov–Shishkin mesh  $$\frac{t}{q-t},\,q=\frac{1}{2}(1+\frac{1}{\ln N})$$  $$3(\ln N)^2$$  $${\,\mathit{e}}^{-{t}/{q-t}}$$  $$3/(2q)\leq 3$$  In order to analyse the convergence of numerical methods on these meshes, the mesh-generating function $$\phi$$ has to fulfil an additional assumption; see Roos & Linß (1999). Assumption 2.3 Let the mesh-generating function $$\phi$$ be piecewise differentiable such that   maxt∈[0,1/2]ϕ′(t)≤CN or equivalently maxt∈[0,1/2]|ψ′(t)|ψ(t)≤CN (2.4) is fulfilled. Note that (2.4) is satisfied for all meshes given in Table 1. We will now use (2.4) to bound the mesh width from above inside the layers. Let $$h_i:=x_i-x_{i-1}$$ and $$t_i=i/N$$ for $$i=1,\dots,N/2$$. Then it holds for $$t\in[t_{i-1},t_i]$$ (with $$\max \phi'$$ taken over $$t\in[t_{i-1},t_i]$$),   ψ(ti)=e−ϕ(ti)=e−(ϕ(ti)−ϕ(t))e−ϕ(t) ≥e−(ϕ(ti)−ϕ(ti−1))ψ(t)≥e−N−1maxϕ′ψ(t)≥Cψ(t), (2.5) where we used (2.4) for the last estimate. Furthermore, we have   x=σεβ1ϕ(t)=−σεβ1ln⁡ψ(t)which givesψ(t)=e−β1x/(σε). Using this, the monotonicity of $$\psi$$ and (2.5), we obtain for $$i=1,\ldots,N/2$$ and $$x\in[x_{i-1},x_i]$$,   hi =σεβ1(ϕ(ti)−ϕ(ti−1))≤σβ1εN−1maxt∈[ti−1,ti]ϕ′(t)≤σβ1εN−1(maxt∈[ti−1,ti]|ψ′(t)|)/ψ(ti) ≤CεN−1(maxt∈[ti−1,ti]|ψ′(t)|)/ψ(t)≤CεN−1max|ψ′|eβ1x/(σε) (2.6) with $$\max|\psi'|:=\max\limits_{t\in[0,1/2]}|\psi'(t)|$$. Similarly, we get for $$j=1,\dots,N/2$$,   ℓj:=yj−yj−1 ≤CεN−1max|ψ′|eβ2y/(σε) (2.7) with $$y\in[y_{j-1},y_j]$$. Of course, the simpler bounds   hi ≤CεN−1maxϕ′ ≤Cε,i=1,…,N/2,ℓj ≤CεN−1maxϕ′ ≤Cε,j=1,…,N/2, follow also from (2.4). Let us denote by   h:=maxi=1,…,N/2hiandℓ:=maxj=1,…,N/2ℓj the maximal dimensions of cells in the layer region perpendicular to the boundaries. We will see later on that the errors will be bounded by powers of the quantity   E(N,ε):=h+ℓ+N−1max|ψ′|, (2.8) which depends on $$h$$, $$\ell$$, $$N$$ and $$\psi$$. By the above considerations, we have   E(N,ε)≤{Cε+N−1max|ψ′|for a general S-type mesh,CN−1ln⁡N,for a Shishkin mesh.  Thus, for $${\varepsilon}\leq C N^{-1}$$ this quantity is (almost) of order 1 and especially for the Bakhvalov–Shishkin and the modified Bakhvalov–Shishkin meshes we obtain   E(N,ε)≤CN−1. We denote in the following by $$\tau_{ij}=[x_{i-1},x_i]\times[y_{j-1},y_j]$$ a specific element with dimensions $$h_i$$ and $$\ell_j$$ and by $$\tau$$ a generic mesh rectangle. Note that the mesh cells are assumed to be closed. 3. Abstract problem and its convergence analysis We are interested in conforming discretizations of the Oseen problem, i.e., we choose finite element spaces $$V_N\subset V$$ and $$Q_N\subset Q$$. For this section, we will not specify these discrete spaces and consider examples in Section 4. The standard Galerkin formulation of (2.1) reads as follows: Find $$({\boldsymbol{u}}_N,p_N)\in V_N\times Q_N$$ such that   ε(∇uN,∇vN)−((b⋅∇)uN,vN)+(cuN,vN)−(pN,divvN) =(f,vN),(qN,divuN) =0 (3.1) for all $$({\boldsymbol{v}}_N,q_N)\in V_N\times Q_N$$. In order to ensure the unique solvability of the discrete problem (3.1), the discrete spaces $$V_N$$ and $$Q_N$$ have to fulfil a compatibility condition. Assumption 3.1 For any given $$N$$ and $${\varepsilon}$$, there exists a positive number $$\beta_0(N,{\varepsilon})$$ such that   infqN∈QNsupvN∈VN(qN,divvN)‖∇vN‖‖qN‖≥β0(N,ε)>0 (3.2) holds true. Note that in the assumption above, we assume only a lower positive bound for fixed $$N$$ and $${\varepsilon}$$, not a uniform bound in $$N$$ and $${\varepsilon}$$ that would be an inf–sup or Babuška–Brezzi condition. The $${\varepsilon}$$-dependence of $$\beta_0(N,{\varepsilon})$$ can be caused only via the underlying mesh which depends on $${\varepsilon}$$. If $${\boldsymbol{f}}$$ contains a large irrotational part in the sense of the Helmholtz decomposition, the numerical accuracy of (3.1) is known to be poor; see Linke (2014). In order to overcome this, the grad-div stabilization term $$(\gamma{\text{div}}{\boldsymbol{u}}_N,\gamma{\text{div}}{\boldsymbol{v}}_N)$$ will be added where $$\gamma\in L^{\infty}({\it{\Omega}})$$ fulfils $$\gamma\ge\gamma_0>0$$ in $$\overline{{\it{\Omega}}}$$ unless specified otherwise. Our stabilized discrete formulation looks as follows: Find $$({\boldsymbol{u}}_N,p_N)\in V_N\times Q_N$$ such that   ε(∇uN,∇vN)−((b⋅∇)uN,vN)+(cuN,vN)+(γdivuN,γdivvN)−(pN,divvN) =(f,vN),(qN,divuN) =0 (3.3) for all $$({\boldsymbol{v}}_N,q_N)\in V_N\times Q_N$$. We define the bilinear form   A((u,p),(v,q)):=ε(∇u,∇v)−((b⋅∇)u,v)+(cu,v)+(γdivu,γdivv)−(p,divv)+(q,divu) (3.4) on the product space $$V\times Q$$ and the norms   |||v|||:=(ε‖∇v‖2+c0‖v‖2+‖γdivv‖2)1/2,|||(v,q)|||:=(|||v|||2+α‖q‖2)1/2 (3.5) on $$V$$ and $$V\times Q$$ with $$\alpha=\alpha(\varepsilon,c_0,c_\infty,b_\infty)$$ given by   α=1max{ε+c∞2c0CF2,b∞2}, where $$C_{\rm F}$$ denotes the constant in the Friedrichs’ inequality for the domain $${\it{\Omega}}$$, $$c_\infty:=\|c\|_\infty$$ and $$b_{\infty}:=\|{\boldsymbol{b}}\|_\infty$$. Note that the term $$c_\infty^2/c_0$$ is bounded since we assumed $$0 <c_0\le c$$ in $${\it{\Omega}}$$ and $$c\in L^{\infty}({\it{\Omega}})$$ in the introduction. Lemma 3.2 Let $$({\boldsymbol{u}},p)$$ and $$({\boldsymbol{u}}_N,p_N)$$ denote the solutions of (2.1) and (3.3), respectively. Then, the Galerkin orthogonality   A((u−uN,p−pN),(vN,qN))=0 (3.6) holds true for all $$({\boldsymbol{v}}_N,q_N)\in V_N\times Q_N$$. Proof. Subtract (3.3) from (2.1), use $$V_N\subset V$$, $$Q_N\subset Q$$ and $${\text{div}}{\boldsymbol{u}} =0$$. □ Our analysis relies on an upper bound for the parameter $$\gamma$$ of the grad-div stabilization. Assumption 3.3 Let $$\gamma$$ be chosen such that   γ∞2≤Cγ2α holds true with some positive constant $$C_\gamma$$ where $$\gamma_\infty:=\|\gamma\|_\infty$$. Let us define for any $$q_N\in Q_N$$ the positive number   B0(N,ε,qN):={1‖qN‖supvN∈VN(qN,divvN)‖∇vN‖,qN≠0,1,qN=0.  (3.7) Note that by definition and Assumption 3.1 it follows that   B0(N,ε,qN)≥β0(N,ε). The stability of the bilinear form $$A$$ can be shown similarly to Matthies & Tobiska (2015, Lemma 3.1). Since we have to take into consideration that here the constant $$\alpha$$ is defined differently, $$c$$ and $$\gamma$$ are no longer constant and we will use $$B_0(N,{\varepsilon},q_N)$$ instead of $$\beta_0(N,{\varepsilon})$$; the proof will be given here. Lemma 3.4 Under Assumptions 3.1 and 3.3, there exists for each $$N>0$$, $${\varepsilon}>0$$ and each pair $$({\boldsymbol{v}}_N,q_N)\in V_N\times Q_N$$ a positive number $$B(N,{\varepsilon},q_N)$$ such that   |||(vN,qN)|||≤1B(N,ε,qN)sup(wN,rN)A((vN,qN),(wN,rN))|||(wN,rN)||| (3.8) holds true. Proof. Let $$({\boldsymbol{v}}_N,q_N)$$ be an arbitrary element of $$V_N \times Q_N$$ and $$B_0=B_0(N,{\varepsilon},q_N)$$. We obtain immediately   A((vN,qN),(vN,qN)) ≥ε‖∇vN‖2+c0‖vN‖2+‖γdivvN‖2=|||vN|||2. (3.9) It follows from the definition (3.7) of $$B_0(N,{\varepsilon},q_N)$$ that for $$q_n\in Q_n$$ there exists a $${\boldsymbol{z}}_N={\boldsymbol{z}}_N(q_N) \in V_N$$ such that   (divzN,qN)=−α‖qN‖2,‖∇zN‖=αB0‖qN‖, (3.10) where we already have used the stronger estimate with $$B_0$$ instead of $$\beta_0(N,{\varepsilon})$$. Hence, we have   A((vN,qN),(zN,0))=ε(∇vN,∇zN)−((b⋅∇)vN,zN)+(cvN,zN)+(γdivvN,γdivzN)+α‖qN‖2. Now, we estimate the first four terms of that equation. We obtain directly using Friedrichs’ inequality and (3.10),   |ε(∇vN,∇zN)+(cvN,zN)| ≤(ε‖∇vN‖2+c0‖vN‖2)1/2(ε‖∇zN‖2+c∞2c0‖zN‖2)1/2 ≤(ε+c∞2c0CF2)1/2αB0‖qN‖|||vN||| ≤α4‖qN‖2+1B02|||vN|||2 and   |−((b⋅∇)vN,zN)|=|((b⋅∇)zN,vN)| ≤b∞‖∇zN‖‖vN‖≤α4‖qN‖2+αb∞2B02‖vN‖2 ≤α4‖qN‖2+1B02|||vN|||2, where the definition of $$\alpha$$ was used. For the stabilization term, we are using (3.10) and obtain   (γdivvN,γdivzN) ≤‖γdivvN‖‖γdivzN‖ ≤|||vN|||2γ∞‖∇zN‖≤|||vN|||Cγ1/2α1/2αB0‖qN‖ ≤α4‖qN‖2+CγB02|||vN|||2 with the constant $$C_\gamma$$ from Assumption 3.3. Combining above estimates we get   A((vN,qN),(zN,0))≥α4‖qN‖2−C~4|||vN|||2, (3.11) where $$\tilde C:=4(2+C_\gamma)/B_0^2$$. Thus, we define for $$({\boldsymbol{v}}_N,q_N) \in V_N\times Q_N$$ the pair $$({\boldsymbol{w}}_N,r_N)\in V_N\times Q_N$$ by   (wN,rN):=(vN,qN)+41+C~(zN,0) and obtain with (3.9) and (3.11),   A((vN,qN),(wN,rN)) ≥α1+C~‖qN‖2+(1−C~1+C~)|||vN|||2=11+C~|||(vN,qN)|||2. It remains to show that $$\left|\!\!\;\left|\!\!\;\left| {({\boldsymbol{w}}_N,r_N)} \right|\!\!\;\right|\!\!\;\right|\leq C \left|\!\!\;\left|\!\!\;\left| {({\boldsymbol{v}}_N,q_N)} \right|\!\!\;\right|\!\!\;\right|$$. To this end, we estimate   |||(wN,rN)|||≤|||(vN,qN)|||+41+C~|||zN||| and with Friedrichs’ inequality and (3.10),   |||zN|||2=ε‖∇zN‖2+c0‖zN‖2+‖γdivzN‖2≤(ε+c∞2c0CF2)‖∇zN‖2+CγαB02‖qN‖2≤C~4α‖qN‖2. Taking into account that $$\sup\nolimits_{0 <t <\infty}\frac{4}{1+t}\sqrt{\frac{t}{4}}=1$$ we get   |||(wN,rN)|||≤2|||(vN,qN)|||. Hence, with   B(N,ε,qN):=12(1+C~)=B0(N,ε,qN)22B0(N,ε,qN)2+16+8Cγ, the stated condition holds. □ 3.1 Interpolation error Assumption 3.5 There exists an interpolation operator $$I=(I_1,I_2):(C({\it{\Omega}}))^2\to V_N$$ such that the anisotropic interpolation error bounds   ‖wm−Imwm‖Lp(τij) ≤C(hik+1‖∂xk+1wm‖Lp(τij)+ℓjk+1‖∂yk+1wm‖Lp(τij)),m=1,2, (3.12)  ‖(wm−Imwm)x‖Lp(τij) ≤C(hik‖∂xk+1wm‖Lp(τij)+ℓjk‖∂xk∂ywm‖Lp(τij)),m=1,2, (3.13) and similarly for the $$y$$-derivative, hold true for $$p\in [1,\infty]$$ and arbitrary $${\boldsymbol{w}}=(w_1,w_2)\in (W^{k+1,p}({\it{\Omega}}))^2$$. Furthermore, we have the stability estimates   ‖Imwm‖∞,τ≤C‖wm‖∞,τ,m=1,2, (3.14) and   ‖∂xImwm‖∞,τ≤C‖∂xwm‖∞,τand‖∂yImwm‖∞,τ≤C‖∂ywm‖∞,τ,m=1,2. (3.15) Assumption 3.6 There exists an interpolation operator $$J:C({\it{\Omega}})\to Q_N$$ such that   ‖q−Jq‖∞≤CE(N,ε)k hold true for arbitrary $$q\in H^{k}({\it{\Omega}})$$ where $${E(N,{\varepsilon})}$$ is given by (2.8). Using these assumptions on the existence of interpolation operators, we estimate the interpolation errors and their derivatives in the next two lemmas. Lemma 3.7 Let $$\sigma>k+1$$, $$r=v_1+w_1+z_1$$ and $$s=v_2+w_2+z_2$$ from Assumption 2.1, $$I=(I_1,I_2)$$ the interpolation operator from Assumption 3.5, $$J$$ the interpolation operator from Assumption 3.6 and $${E(N,{\varepsilon})}$$ be given by (2.8). Then, the interpolation errors in the $$L^2$$-norm can be bounded by   ‖r−I1r‖Ω11+‖s−I2s‖Ω11 ≤CN−σ(ε1/2+N−1/2), (3.16)  ‖r−I1r‖Ω∖Ω11+‖s−I2s‖Ω∖Ω11 ≤Cε1/2E(N,ε)k+1. (3.17) Moreover, the estimates   ‖g1−I1g1‖Ω11+‖g2−I2g2‖Ω11 ≤CN−(k+1),‖g1−I1g1‖Ω∖Ω11+‖g2−I2g2‖Ω∖Ω11 ≤CE(N,ε)k+1 hold true for the smooth parts $$g_1$$ and $$g_2$$. Proof. Since the proof of the error bounds for $$s$$ follows the same line as the proof for $$r$$, we restrict ourselves to the latter one. We start with $${\it{\Omega}}_{11}$$. Clearly, we have   ‖r−I1r‖Ω11≤‖r‖Ω11+‖I1r‖Ω11. (3.18) With Assumption 2.1, we get for $$\|r\|_{{\it{\Omega}}_{11}}$$ directly,   ‖r‖Ω11 ≤‖v1‖Ω11+‖w1‖Ω11+‖z1‖Ω11≤Cε1/2N−σ, (3.19) where we used integration and the definition of $$\lambda_x$$ and $$\lambda_y$$. Let $$\overline{h}$$ and $$\overline{\ell}$$ denote the mesh width and height in $${\it{\Omega}}_{11}$$. In order to estimate $$I_1 r,$$ we use an idea from Stynes & Tobiska (2003) and divide the domain $${\it{\Omega}}_{11}$$ into $$S=[\lambda_x+\overline{h},1]\times [\lambda_y+\overline{\ell},1]$$ and $${\it{\Omega}}_{11}\setminus S$$ with the mesh transition points $$\lambda_x$$ and $$\lambda_y$$. Note that $${\it{\Omega}}_{11}\setminus S$$ consists of only a single ply of $$N-1$$ mesh elements and $$\overline{h},\overline{\ell}\leq CN^{-1}$$. Thus, we obtain   ‖I1r‖Ω11∖S2≤∑τ⊂Ω11∖Sh¯ℓ¯‖I1r‖∞,τ2≤CN−1‖r‖∞,Ω112≤CN−2σ−1 (3.20) using the stability (3.14) and in $$S$$,   ‖I1r‖S2≤C{∫Ω11(ε2e−2β1x/ε+e−2β2y/ε(1+ε2e−2β1x/ε))}≤CεN−2σ. (3.21) Combining (3.18)–(3.21) completes the proof for (3.16). On $${\it{\Omega}}_{12}$$, we get by using (3.12), (2.6) and (2.7),   ‖v1−I1v1‖Ω122 ≤C(‖hik+1∂xk+1v1‖Ω122+‖ℓjk+1∂yk+1v1‖Ω122)≤Cε3(N−1max|ψ′|)2(k+1). Using the stability (3.14), we obtain for $$w_1$$ and $$z_1$$,   ‖(w1+z1)−I1(w1+z1)‖Ω12≤C(measΩ12)1/2‖w1+z1‖∞,Ω12≤Cε1/2N−σ(ln⁡N)1/2. Similarly we have on $${\it{\Omega}}_{21}$$,   ‖w1−I1w1‖Ω212≤Cε(N−1max|ψ′|)2(k+1) and   ‖(v1+z1)−I1(v1+z1)‖Ω21≤Cε3/2N−σ(ln⁡N)1/2. Using the anisotropic interpolation error bounds (3.12), we get on $${\it{\Omega}}_{22}$$ analogously,   ‖v1−I1v1‖Ω222 ≤Cε4ln⁡N(N−1max|ψ′|+ℓ)2(k+1),‖w1−I1w1‖Ω222 ≤Cε2ln⁡N(h+N−1max|ψ′|)2(k+1),‖z1−I1z1‖Ω222 ≤Cε4(N−1max|ψ′|)2(k+1). Together with $$\sigma>k+1$$ and the definition of $${E(N,{\varepsilon})}$$ we get (3.17). For the smooth part, we present only the proofs for $$g_1$$ since $$g_2$$ can be estimated similarly. We have   ‖g1−I1g1‖Ω11 ≤C(measΩ11)1/2‖g1−I1g1‖∞,Ω11≤C(h¯k+1+ℓ¯k+1)≤CN−(k+1),‖g1−I1g1‖Ω∖Ω11 ≤C(measΩ∖Ω11)1/2‖g1−I1g1‖∞,Ω∖Ω11≤C(h+ℓ+N−1)k+1, where (3.12) in the $$L^\infty$$-norm was used. □ Lemma 3.8 Let $$\sigma>k+1$$ and $${E(N,{\varepsilon})}$$ be given by (2.8). Then the interpolation errors in the $$H^1$$-seminorm for $${\boldsymbol{u}}=(u_1,u_2)$$ and $$I=(I_1,I_2)$$ from Assumption 3.5 can be bounded by   ‖∇(u−Iu)‖Ω11 ≤Cε−1/2N−k,‖∇(u−Iu)‖Ω∖Ω11 ≤Cε−1/2E(N,ε)k while the estimates   ‖div(u−Iu)‖Ω11 ≤CN−k,‖div(u−Iu)‖Ω∖Ω11 ≤CE(N,ε)k for the divergence of the interpolation error hold true. Proof. Since the proofs of the estimates in the $$H^1$$-seminorm are similar for $$u_1$$ and $$u_2$$, we show only the proofs for $$u_1$$. With the anisotropic interpolation error bounds (3.13) and its $$y$$-derivative counterpart we get for the smooth part,   ‖∇(g1−I1g1)‖Ω11 ≤C(‖(g1−I1g1)x‖∞,Ω11+‖(g1−I1g1)y‖∞,Ω11)≤CN−k and similarly   ‖∇(g1−I1g1)‖Ω∖Ω11 ≤C(h+ℓ+N−1)k. Now, we estimate the layer parts. Consider first the term $$v_1$$. Here we have by the $$W^{1,\infty}$$-stability (3.15),   ‖(v1−I1v1)x‖Ω11∪Ω21 ≤C(‖v1,x‖Ω11∪Ω21+(measΩ11∪Ω21)1/2‖(I1v1)x‖∞,Ω11∪Ω21) ≤C(‖v1,x‖Ω11∪Ω21+‖v1,x‖∞,Ω11∪Ω21) ≤CN−σ. In the remaining regions, we use (3.13) of $$I_1$$ to deduce   ‖(v1−I1v1)x‖Ω12∪Ω22 ≤C{‖hik∂xk+1v1‖Ω12∪Ω22+‖ℓjk∂x∂ykv1‖Ω12∪Ω22} ≤Cε1/2(ℓ+N−1max|ψ′|)k. The techniques for the $$y$$-derivatives are similar. Hence, we obtain   ‖(v1−I1v1)y‖Ω11∪Ω21 ≤C(‖v1,y‖Ω11∪Ω21+‖I1v1,y‖∞,Ω11∪Ω21) ≤C(‖v1,y‖Ω11∪Ω21+‖v1,y‖∞,Ω11∪Ω21) ≤CεN−σ and   ‖(v1−I1v1)y‖Ω12∪Ω22 ≤C{‖hik∂xk∂yv1‖Ω12∪Ω22+‖ℓjk∂yk+1v1‖Ω12∪Ω22} ≤Cε3/2(ℓ+N−1max|ψ′|)k. For the second layer part $$w_1$$, we get analogously   ‖(w1−I1w1)x‖Ω11∪Ω12 ≤CN−σ,‖(w1−I1w1)x‖Ω21Ω22 ≤Cε1/2(h+N−1max|ψ′|)k,‖(w1−I1w1)y‖Ω11∪Ω12 ≤C(ε−1/2+N)N−σ,‖(w1−I1w1)y‖Ω21∪Ω22 ≤Cε−1/2(h+N−1max|ψ′|)k. Now, we estimate the error of the function $$z_1$$ using the $$W^{1,\infty}$$-stability (3.15) of $$I_1$$,   ‖(z1−I1z1)x‖Ω∖Ω22 ≤C(‖z1,x‖Ω∖Ω22+‖z1,x‖∞,Ω∖Ω22)≤CN−σ, and in $${\it{\Omega}}_{22}$$ with (3.13),   ‖(z1−I1z1)x‖Ω22 ≤C{‖hik∂xk+1z1‖Ω22+‖ℓjk∂x∂ykz1‖Ω22}≤Cε(N−1max|ψ′|)k. The bounds for the $$y$$-derivatives follow analogously:   ‖(z1−I1z1)y‖Ω∖Ω22 ≤CN−σ,‖(z1−I1z1)y‖Ω22 ≤Cε(N−1max|ψ′|)k. Putting all these estimates together, we obtain with the definition of $${E(N,{\varepsilon})}$$ the bounds for the $$H^1$$-seminorms. We remark that the bounds of the $$x$$-derivatives of $$u_1$$ are better than those of the $$y$$-derivatives. For $$u_2$$ and $$I_2,$$ we obtain similar results with the difference that the $$x$$-derivatives are worse than the $$y$$-derivatives. Now combining the $$x$$-derivatives of $$u_1$$ and the $$y$$-derivatives of $$u_2$$, we get   ‖div(u−Iu)‖Ω11 ≤‖(u1−I1u1)x‖Ω11+‖(u2−I2u2)y‖Ω11≤CN−k,‖div(u−Iu)‖Ω∖Ω11 ≤C(εln⁡N)1/2(h+ℓ+N−1max|ψ′|)k≤CE(N,ε)k, which are the last two statements of this lemma. □ Lemma 3.9 Let $$\sigma>k+1$$, $$I=(I_1,I_2),$$ the interpolation operator from Assumption 3.5 and $$J$$ the interpolation operator from Assumption 3.6. Then the estimate   |||(u−Iu,p−Jp)|||≤CE(N,ε)k (3.22) holds true where $${E(N,{\varepsilon})}$$ is given by (2.8). Proof. For this proof, we use the previous interpolation error estimates. The definition of the energy norm gives   |||(u−Iu,p−Jp)|||2=ε‖∇(u−Iu)|02+c0‖u−Iu‖2+α‖p−Jp‖2+‖γdiv(u−Iu)‖2. (3.23) For the first term on the right-hand side of (3.23), we get with Lemma (3.8),   ε1/2‖∇(u−Iu)‖ ≤ε1/2(‖∇(u1−I1u1)‖+‖∇(u2−I2u2)‖)≤CE(N,ε)k. The second term of (3.23) yields with Lemma 3.7,   ‖u−Iu‖ ≤C(‖u1−I1u1‖+‖u2−I2u2‖)≤CE(N,ε)k+1. The bound   ‖p−Jp‖≤CE(N,ε)k for the pressure term follows from Assumption 3.6. The stabilization term can be bounded by the divergence bounds of Lemma (3.8),   ‖γdiv(u−Iu)‖≤γ∞‖div(u−Iu)‖≤CE(N,ε)k. Summarizing all these estimates finishes the proof. □ 3.2 A priori error estimates Lemma 3.10 Let $$\gamma\geq \gamma_0>0$$ be independent of $$N$$ and $${\varepsilon}$$, and $$\sigma>k+1$$. The solutions of 2.1 and 3.3 are denoted by $$({\boldsymbol{u}},p)$$ and $$({\boldsymbol{u}}_N,p_N)$$, respectively. Then, the estimate   |||(Iu−uN,Jp−pN)|||≤C1B(N,ε,Jp−pN)E(N,ε)k holds true where $$C\to\infty$$ for $$\gamma_0\to 0$$. Proof Abbreviating $$B=B(N,{\varepsilon},J p-p_N)$$ in this proof, the estimate 3.8 gives   |||(Iu−uN,Jp−pN)||| ≤1Bsup(wN,rN)∈VN×QNA((Iu−uN,Jp−pN),(wN,rN))|||(wN,rN)||| =1Bsup(wN,rN)∈VN×QNA((Iu−u,Jp−p),(wN,rN))|||(wN,rN)||| due to the Galerkin orthogonality of Lemma 3.2. The expression   A((Iu−u,Jp−p),(wN,rN))  =ε(∇(Iu−u),∇wN)−(b⋅∇(Iu−u),wN)+(c(Iu−u),wN)  −(Jp−p,divwN)+(rN,div(Iu−u))+(γdiv(Iu−u),γdiv(wN)) will be estimated term by term. We get with the interpolation error bounds of Lemmas 3.7 and 3.8    |ε(∇(Iu−u),∇wN)+(c(Iu−u),wN)|  ≤(ε‖∇(Iu−u)‖2+c∞2c0‖Iu−u‖2)1/2(ε‖∇wN‖2+c0‖wN‖2)1/2  ≤CE(N,ε)k|||(wN,rN)|||, by Lemma 3.7 for the pressure,   |(Jp−p,divwN)| =|(Jp−p,divwN)| ≤‖Jp−p‖γ0−1|||(wN,rN)||| ≤Cγ0−1E(N,ε)k|||(wN,rN)||| and with the divergence bounds of Lemma 3.8,   |(rN,div(Iu−u))+(γdiv(Iu−u),γdivwN))| ≤‖div(Iu−u)‖(‖rN‖+γ∞‖γdivwN‖) ≤(α−1/2+γ∞)‖div(Iu−u)‖|||(wN,rN)||| ≤C(1+γ∞)E(N,ε)k|||(wN,rN)|||. The last two estimates show that both $$\gamma_0^{-1}$$ and $$\gamma_{\infty}$$ should be bounded independent of $$N$$ and $${\varepsilon}$$. Hence, $$\gamma$$ should be contained in the interval $$[\gamma_0,\gamma_{\infty}]$$ with end points independent of $$N$$ and $${\varepsilon}$$. For the convective term, we estimate the terms of   |(b⋅∇(Iu−u),wN)| |(b1(I1u1−u1)x,wN,1)|+|(b2(I1u1−u1)y,wN,1)|+ |(b1(I2u2−u2)x,wN,2)|+|(b2(I2u2−u2)y,wN,2)| individually. Here we have directly,   |(b1(I1u1−u1)x,wN,1)| ≤C‖(I1u1−u1)x‖‖wN,1‖≤CE(N,ε)k|||(wN,rN)|||,|(b2(I2u2−u2)y,wN,2)| ≤C‖(I2u2−u2)y‖‖wN,2‖≤CE(N,ε)k|||(wN,rN)|||. For the other terms, we use integration by parts to obtain   |(b2(I1u1−u1)y,wN,1)| ≤C(‖I1u1−u1‖‖wN,1‖+|(b2(I1u1−u1),wN,1,y)|) ≤C(‖I1u1−u1‖‖wN,1‖+‖I1u1−u1‖Ω11‖wN,1,y‖Ω11 +‖I1u1−u1‖Ω∖Ω11‖wN,1,y‖Ω∖Ω11) ≤C(‖u1−u1‖+N‖I1u1−u1‖Ω11+ε−1/2‖I1u1−u1‖Ω∖Ω11)|||(wN,rN)||| ≤CE(N,ε)k|||(wN,rN)|||. The remaining term can be estimated analogously. Upon summarizing the proof is done. □ Theorem 3.11 Let $$\gamma\geq\gamma_0>0$$ be independent of $$N$$ and $${\varepsilon}$$, and $$\sigma>k+1$$. The error between the solution $$({\boldsymbol{u}},p)$$ of 2.1 and the discrete solution $$({\boldsymbol{u}}_N,p_N)$$ of 3.3 can be bounded by   |||(u−uN,p−pN)|||≤C(1+1B(N,ε,Jp−pN))E(N,ε)k, where $${E(N,{\varepsilon})}$$ is given by 2.8 and $$J$$ is an interpolation operator fulfilling Assumption 3.6. Proof It follows directly from Lemmas 3.9 and 3.10, and the triangle inequality. □ Corollary 3.12 If $$\gamma=0$$ is chosen on the whole domain, then we have the standard Galerkin method. For the solution of this method we have the estimate   |||(u−uN,p−pN)|||≤C(1+1B(N,ε,Jp−pN))(N+(ln⁡N)1/2)E(N,ε)k, i.e., the convergence order reduces by 1 compared to the stabilized case. Proof The only difference to the estimation before is the bounding of the term $$|(J p-p,{\text{div}} {\boldsymbol{w}}_N)|$$. This term cannot be estimated against $$\|\gamma{\text{div}} {\boldsymbol{w}}_N\|$$ due to $$\gamma=0$$. Alternatively, we get   |(Jp−p,divwN)| =|(Jp−p,divwN)| ≤C(‖Jp−p‖Ω11‖∇wN‖Ω11+‖Jp−p‖Ω∖Ω11‖∇wN‖Ω∖Ω11) ≤C(N‖Jp−p‖Ω11‖wN‖Ω11+(εln⁡N)1/2‖Jp−p‖∞,Ω∖Ω11‖∇wN‖Ω∖Ω11) ≤C(N‖Jp−p‖+(ln⁡N)1/2‖Jp−p‖∞)|||(wN,rN)||| ≤C(N+(ln⁡N)1/2)E(N,ε)k|||(wN,rN)||| and the reduced order can be seen. □ Corollary 3.13 If we choose $$\gamma=0$$ only in the layer region, we can recover the convergence rate $$k$$ compared with $$\gamma=0$$ on the whole domain. Proof See the proof before. In the layer region, we obtain an additional factor of only $$(\ln N)^{1/2}$$ and not $$N$$. □ 4 Examples for suitable spaces We will consider two families of finite element pairs on the introduced meshes. The first one is the Taylor–Hood family $$Q_k\times Q_{k-1}$$, $$k\ge 2$$, and the second one is $$Q_k\times P_{k-1}^{\text{disc}}$$, $$k\ge 2$$. On isotropic meshes, the pairs $$Q_k\times Q_{k-1}$$, $$k\ge 2$$, and $$Q_k\times P_{k-1}^{\text{disc}}$$, $$k\ge 2$$, are known to fulfil Assumption 3.1 even with a constant $$\beta_0>0$$ that is independent of $$N$$; see Girault & Raviart (1986, Ch. II), and Assumption 3.1 becomes just the well-known inf–sup or Babuška–Brezzi condition. Since the $${\varepsilon}$$-dependence of $$\beta_0(N,{\varepsilon})$$ comes from the underlying mesh only, there is also no $${\varepsilon}$$-dependence on general isotropic meshes. On anisotropic meshes, the situation is much more complicated. The existence of an inf–sup constant that is independent of the maximal aspect ratio of the underlying mesh depends strongly on the precise mesh structure and the chosen spaces for approximating velocity and pressure. Proofs for the existence of inf–sup constants are mainly available for situations where the difference between the polynomial orders of velocity and pressure is $$2$$; see Sch¨otzau & Schwab (1998), Sch¨otzau et al. (1999) for the cases $$Q_k\times Q_{k-2}$$ on two-dimensional domains. A review of results on inf–sup constants for the Stokes problem on several types of anisotropic meshes can be found in Apel & Maharavo Randrianarivony (2003). Nonconforming discretizations of the Stokes problem on anisotropic rectangular meshes were investigated in Apel & Matthies (2008). Although the four nonconforming velocity approximation spaces in Apel & Matthies (2008) differ only slightly, just two of them provide an inf–sup constant independent of the aspect ratio. Another approach for generating inf–sup-stable pairs by eliminating certain pressure modes is presented in Ainsworth et al. (2015). In Section 5.1.1, we investigate the values of $$\beta_0=\beta_0(N,{\varepsilon})$$ and $$B_0=B_0(N,{\varepsilon},q_N)$$ for both element pairs. It can be observed that for the continuous pressure space, $$\beta_0$$ is bounded independent of $${\varepsilon}$$ but not independent of $$N$$, while for the discontinuous pressure space it goes to zero for $${\varepsilon}\to 0$$. In contrast to this negative observation, the values of $$B_0$$ stay bounded away from zero, independent of $$N$$ and $${\varepsilon}$$, for both cases if $$q_N$$ is the difference between the interpolation of $$p$$ and the discrete solution $$p_N$$. This shows that $$B_0$$ might be more appropriate to estimate the errors although we do not have a theoretical insight into the boundedness of $$B_0$$. In order to define the interpolation operators, let $$I_k$$ be the scalar nodal interpolator of order $$k$$. Then we set $$I=(I_k,I_k)$$. In Apel (1999, Theorem 2.7) the anisotropic interpolation error bounds for Assumption 3.5 are given. The stability estimates are standard and can be shown directly. For the interpolation operator into the discrete pressure space $$Q_N$$, we have to distinguish two cases. In the case $$Q_N=P_{k-1}^{\text{disc}}$$, we can use the $$L^2$$-projection into $$Q_N$$ as $$J$$. The global $$L^2$$-projection furthermore localizes to $$L^2$$-projections on each mesh cell. Hence, the proof of Assumption 3.6 can exploit the anisotropic error estimates for the $$L^2$$-projection that follow the same reasoning as in Apel (1999). In the case $$Q_N=Q_{k-1}$$, the construction of $$J$$ is two stage. Let $$I_{k-1}$$ denote the scalar nodal interpolation into the space of continuous, piecewise $$Q_{k-1}$$ functions. Then we set   Jq=Ik−1q−c¯,wherec¯=1|Ω|∫ΩIk−1q to ensure that $$J$$ maps into $$L^2_0({\it{\Omega}})$$. Using again the anisotropic interpolation error estimates by Apel (1999), we have   ‖p−Ik−1p‖Ω11≤CN−kand‖p−Ik−1p‖Ω∖Ω11≤CE(N,ε)k. From the definition of $$J$$ and the fact that constants are $$L^2({\it{\Omega}})$$-orthogonal to functions from $$L^2_0({\it{\Omega}})$$, we get   ‖p−Ik−1p‖2=‖p−Jp−c¯‖2=‖p−Jp‖2+‖c¯‖2 and therefore   ‖p−Jp‖≤‖p−Ik−1p‖≤CE(N,ε)k. Furthermore, we obtain in $$L^\infty$$,   ‖p−Jp‖∞ ≤‖p−Ik−1p‖∞+|c¯|≤‖p−Ik−1p‖∞+|∫ΩIk−1p−p|≤CE(N,ε)k which provides Assumption 3.6. Corollary 4.1 For both choices of the discrete pressure space, Corollary 3.12 can be improved slightly. If $$Q_N=P_{k-1}^{\rm disc}$$ then $$\|J p-p\|_{{\it{\Omega}}_{11}}\leq CN^{-k}$$ and the final lines in the proof of Corollary 3.12 become   |(Jp−p,divwN)| ≤C(N‖Jp−p‖Ω11+(ln⁡N)1/2‖Jp−p‖∞,Ω∖Ω11)|||(wN,rN)||| ≤C(N−k+1+(ln⁡N)1/2E(N,ε)k)|||(wN,rN)||| ≤CE(N,ε)k−1|||(wN,rN)|||. If $$Q_N=Q_{k-1}$$ then $$(J p-p,{\text{div}} {\boldsymbol{w}}_N)=(I_{k-1}p-p,{\text{div}}{\boldsymbol{w}}_N)$$ and the proof of Corollary 3.12 can be done with $$I_{k-1}$$ instead of $$J$$. Using $$\|I_{k-1} p-p\|_{{\it{\Omega}}_{11}}\leq CN^{-k}$$ leads to   |(Jp−p,divwN)| =|(Ik−1p−p,divwN)| ≤C(N−k+1+(ln⁡N)1/2E(N,ε)k)|||(wN,rN)||| ≤CE(N,ε)k−1|||(wN,rN)|||. 5 Numerical results In the numerical section, we show the results of computations supporting our theoretical findings. We start with a problem where the exact solution is known and end with some considerations on problems with unknown exact solutions where our theory does not hold to its full extent. 5.1 Problem with known solution We consider the singularly perturbed Oseen equations   −εΔu−((2+3x3+2y2)⋅∇)u+(1+2xy)u+∇p =f in Ω=(0,1)2,divu =0 in Ω,u =0 on ∂Ω, (5.1) where the right-hand side $${\boldsymbol{f}}$$ of (5.1) was chosen such that   u1(x,y)=∂Ψ(x,y)∂y,u2(x,y)=−∂Ψ(x,y)∂x,p(x,y)=2cos⁡(x)sin⁡(y)−2sin⁡(1)(1−cos⁡(1)) with $$\Psi(x,y):=F(x)G(y)$$ and   F(x) :=c3(1−x)3+c2(1−x)2+c1(1−x)−ε(e−2x/ε−e−2/ε)−sin⁡(1−x),G(y) :=d3(1−y)3+d2(1−y)2+d1(1−y)−1+ε(e−3y/ε−e−3/ε)+cos⁡(1−y) and constants   c1 :=(1+2e−2/ε),c2:=−((4+3ε)e−2/ε+4+cos⁡(1)−3sin⁡(1)−3ε),c3 :=((2+2ε)e−2/ε+3+cos⁡(1)−2ε−2sin⁡(1)),d1 :=−3e−3/ε,d2:=((6+3ε)e−3/ε+6−3cos⁡(1)−sin⁡(1)−3ε),d3 :=−((3+2ε)e−3/ε+5−2cos⁡(1)−sin⁡(1)−2ε) is the solution of (5.1). The function $$u_1$$ shows a strong exponential boundary layer at $$y=0$$ and a weak exponential boundary layer at $$x=0$$, and vice versa for $$u_2$$. The pressure shows no layer behaviour. Moreover, Assumption 2.1 is satisfied. For a visualization of the solution, see Fig. 1. All calculations were done in MATLAB using the finite element suite $$\mathbb{SOFE}$$ developed by Lars Ludwig; see Ludwig (2016). In the following, ‘order’ will always denote the exponent $$r$$ in a convergence order of form $$\mathcal{O}(N^{-r})$$ while ‘$$\ln$$-order’ corresponds to the exponent $$r$$ in a convergence order of form $$\mathcal{O}((N^{-1}\ln N)^r)$$. If not stated otherwise then the results on the two finest meshes were used to calculate the convergence order. We have chosen $$\gamma=1$$ on the whole domain, unless otherwise stated. The parameter $$\sigma$$ was chosen to be $$\sigma=k+2$$, where $$k$$ is the order of the elements for the velocity $${\boldsymbol{u}}$$. The errors will be measured in the norm   |||(u,p)|||ε2:=ε‖∇u‖2+‖u‖2+‖p‖2+‖γdivu‖2; (5.2) this means that $$\alpha$$ is set to $$1$$. 5.1.1 Dependence of $${\beta_0}$$ and $${B_0}$$ on $${N}$$ and $${{\varepsilon}}$$ First we will consider the dependence of $$\beta_0=\beta_0(N,{\varepsilon})$$ and $$B_0=B_0(N,{\varepsilon},Jp-p_N)$$ on $$N$$ and $${\varepsilon}$$. The pictures in Fig. 3 show the dependence of $$\beta_0$$ on $$N$$ and $${\varepsilon}$$ for the discretizations $$Q_2\times Q_1$$ and $$Q_3\times Q_2$$ on Shishkin meshes. Clearly, there is an $${\varepsilon}$$-uniform bound for any fixed $$N$$ while the value of $$\beta_0$$ will decrease with increasing $$N$$ if $${\varepsilon}$$ is fixed. If $$N$$ and $${\varepsilon}$$ are fixed, the constant $$\beta_0$$ decreases with increasing approximation order $$k$$. This effect appears also on isotropic meshes and spectral methods; see Bernardi & Maday (1997). The behaviour of $$\beta_0$$ on Bakhvalov–Shishkin meshes is presented in Fig. 4. We observe that $$\beta_0$$ takes much larger values and the dependence on $$N$$ is much weaker compared to Shishkin meshes. Numerical experiments for discretizations $$Q_k\times P_{k-1}^{\text{disc}}$$ show that $$\beta_0$$ approaches zero if $${\varepsilon}$$ tends to zero. Fig. 3. View largeDownload slide Problem (5.1): The inf–sup constants for different values of $$N$$ as a function of $$\varepsilon$$ on Shishkin meshes: $$Q_2\times Q_1$$ (left) and $$Q_3\times Q_2$$ (right). Fig. 3. View largeDownload slide Problem (5.1): The inf–sup constants for different values of $$N$$ as a function of $$\varepsilon$$ on Shishkin meshes: $$Q_2\times Q_1$$ (left) and $$Q_3\times Q_2$$ (right). Fig. 4. View largeDownload slide Problem (5.1): The inf–sup constants for different values of $$N$$ as a function of $$\varepsilon$$ on Bakhvalov–Shishkin meshes: $$Q_2\times Q_1$$ (left) and $$Q_3\times Q_2$$ (right). Fig. 4. View largeDownload slide Problem (5.1): The inf–sup constants for different values of $$N$$ as a function of $$\varepsilon$$ on Bakhvalov–Shishkin meshes: $$Q_2\times Q_1$$ (left) and $$Q_3\times Q_2$$ (right). Figure 5 shows the dependence of $$\beta_0$$ and $$B_0$$ on $$N$$ for $${\varepsilon}=10^{-8}$$. We observe that the values of $$B_0$$ are much larger than the corresponding values of $$\beta_0$$. Moreover, the value of $$B_0$$ is bounded away from zero for both continuous and discontinuous pressure spaces and is nondecreasing in $$N$$. Additional studies showed that there is also an $${\varepsilon}$$-uniform lower bound for $$B_0$$. Hence, Theorem 3.11 provides an $${\varepsilon}$$-uniform error bound. Fig. 5. View largeDownload slide Problem (5.1): The behaviour of $$\beta_0$$ and $$B_0$$ with $${\varepsilon}=10^{-8}$$ as a function of $$N$$ for $$Q_3\times Q_2$$ (left) and $$Q_3\times\,P_2^{\text{disc}}$$ (right), respectively. Fig. 5. View largeDownload slide Problem (5.1): The behaviour of $$\beta_0$$ and $$B_0$$ with $${\varepsilon}=10^{-8}$$ as a function of $$N$$ for $$Q_3\times Q_2$$ (left) and $$Q_3\times\,P_2^{\text{disc}}$$ (right), respectively. 5.1.2 Uniformity in $${{\varepsilon}}$$ For this study, we have chosen $$\gamma=1$$ on the whole domain in the computations. Table 2 presents the numerical results on Shishkin meshes (S-mesh) and Bakhvalov–Shishkin meshes (B–S-mesh) with $$N=32$$ and $${\varepsilon}$$ varying from $$10^{-1}$$ to $$10^{-10}$$. We show second-order approximations by $$Q_2\times Q_1$$ and $$Q_2\times P_1^{\text{disc}}$$ as well as third-order approximations by $$Q_3\times Q_2$$ and $$Q_3\times P_2^{\text{disc}}$$. Note that for large values of $${\varepsilon}$$ the mesh is equidistant while the characteristic meshes are obtained for smaller values of $${\varepsilon}$$. We clearly see from Table 2 that the error for small $${\varepsilon}$$ is independent of $${\varepsilon}$$, as predicted by our theory and the $${\varepsilon}$$-uniform bound for $$B_0$$. Furthermore, the errors for discretizations of the same order and on the same mesh family differ only slightly. Moreover, the errors on Bakhvalov–Shishkin meshes are smaller than those on Shishkin meshes. Table 2 Problem (5.1): Error in energy norm $${\left|\!\!\;\left|\!\!\;\left| {(\cdot,\cdot)} \right|\!\!\;\right|\!\!\;\right|}_{\varepsilon}$$ for $$N=32$$    $$Q_2\times Q_1$$  $$Q_2\times P^{\text{disc}}_1$$  $$Q_3\times Q_2$$  $$Q_3\times P^{\text{disc}}_2$$  $${\varepsilon}$$  S-mesh  B–S-mesh  S-mesh  B–S-mesh  S-mesh  B–S-mesh  S-mesh  B–S-mesh  $$10^{-1}$$  1.7439-02  8.6659-03  1.7439-02  8.6661-03  1.4779-03  1.4779-03  1.4779-03  1.4779-03  $$10^{-2}$$  2.4899-02  3.5717-03  2.4901-02  3.5810-03  3.4299-03  1.2697-04  3.4299-03  1.2702-04  $$10^{-3}$$  2.5197-02  3.7450-03  2.5198-02  3.7565-03  3.4685-03  1.2840-04  3.4685-03  1.2850-04  $$10^{-4}$$  2.5230-02  3.8047-03  2.5233-02  3.8166-03  3.4723-03  1.2854-04  3.4723-03  1.2865-04  $$10^{-5}$$  2.5234-02  3.8146-03  2.5237-02  3.8273-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04  $$10^{-6}$$  2.5235-02  3.8156-03  2.5237-02  3.8292-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04  $$10^{-7}$$  2.5235-02  3.8157-03  2.5237-02  3.8282-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04  $$10^{-8}$$  2.5235-02  3.8157-03  2.5237-02  3.8272-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04  $$10^{-9}$$  2.5235-02  3.8157-03  2.5237-02  3.8271-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04  $$10^{-10}$$  2.5235-02  3.8157-03  2.5237-02  3.8271-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04     $$Q_2\times Q_1$$  $$Q_2\times P^{\text{disc}}_1$$  $$Q_3\times Q_2$$  $$Q_3\times P^{\text{disc}}_2$$  $${\varepsilon}$$  S-mesh  B–S-mesh  S-mesh  B–S-mesh  S-mesh  B–S-mesh  S-mesh  B–S-mesh  $$10^{-1}$$  1.7439-02  8.6659-03  1.7439-02  8.6661-03  1.4779-03  1.4779-03  1.4779-03  1.4779-03  $$10^{-2}$$  2.4899-02  3.5717-03  2.4901-02  3.5810-03  3.4299-03  1.2697-04  3.4299-03  1.2702-04  $$10^{-3}$$  2.5197-02  3.7450-03  2.5198-02  3.7565-03  3.4685-03  1.2840-04  3.4685-03  1.2850-04  $$10^{-4}$$  2.5230-02  3.8047-03  2.5233-02  3.8166-03  3.4723-03  1.2854-04  3.4723-03  1.2865-04  $$10^{-5}$$  2.5234-02  3.8146-03  2.5237-02  3.8273-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04  $$10^{-6}$$  2.5235-02  3.8156-03  2.5237-02  3.8292-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04  $$10^{-7}$$  2.5235-02  3.8157-03  2.5237-02  3.8282-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04  $$10^{-8}$$  2.5235-02  3.8157-03  2.5237-02  3.8272-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04  $$10^{-9}$$  2.5235-02  3.8157-03  2.5237-02  3.8271-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04  $$10^{-10}$$  2.5235-02  3.8157-03  2.5237-02  3.8271-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04  5.1.3 Convergence studies Table 3 shows for Shishkin meshes and $${\varepsilon}=10^{-8}$$, which is sufficiently small to have a solution with layers, the errors in the energy norm for pairs of second-, third-, and fourth-order with continuous and discontinuous pressure approximations. It can be seen that the theoretically predicted convergence orders are achieved. The differences in the energy norm between the two different pressure approximation spaces are again very small. Here the different pressure errors are similar and in both cases dominated by the velocity errors. Table 3 Problem (5.1): Error in energy norm $${\left|\!\!\;\left|\!\!\;\left| {(\cdot,\cdot)} \right|\!\!\;\right|\!\!\;\right|}_{\varepsilon}$$ on Shishkin meshes with $${\varepsilon}=10^{-8}$$ $$N$$  $$Q_2\times Q_1$$  $$Q_2\times P_1^{\text{disc}}$$  $$Q_3\times Q_2$$  $$Q_3\times P_2^{\text{disc}}$$  $$Q_4\times Q_3$$  $$Q_4\times P_3^{\text{disc}}$$  4  2.7153-01  2.5880-01  7.2169-02  7.2220-02  2.5802-02  2.5803-02  8  1.3134-01  1.3131-01  3.6858-02  3.6859-02  1.0660-02  1.0660-02  16  6.1849-02  6.1859-02  1.2918-02  1.2918-02  2.6828-03  2.6828-03  32  2.5235-02  2.5237-02  3.4727-03  3.4727-03  4.7297-04  4.7297-04  64  9.2647-03  9.2650-03  7.8324-04  7.8324-04  6.5527-05  6.5527-05  $$\ln$$-order  1.96  1.96  2.92  2.92  3.87  3.87  Theory  2  2  3  3  4  4  $$N$$  $$Q_2\times Q_1$$  $$Q_2\times P_1^{\text{disc}}$$  $$Q_3\times Q_2$$  $$Q_3\times P_2^{\text{disc}}$$  $$Q_4\times Q_3$$  $$Q_4\times P_3^{\text{disc}}$$  4  2.7153-01  2.5880-01  7.2169-02  7.2220-02  2.5802-02  2.5803-02  8  1.3134-01  1.3131-01  3.6858-02  3.6859-02  1.0660-02  1.0660-02  16  6.1849-02  6.1859-02  1.2918-02  1.2918-02  2.6828-03  2.6828-03  32  2.5235-02  2.5237-02  3.4727-03  3.4727-03  4.7297-04  4.7297-04  64  9.2647-03  9.2650-03  7.8324-04  7.8324-04  6.5527-05  6.5527-05  $$\ln$$-order  1.96  1.96  2.92  2.92  3.87  3.87  Theory  2  2  3  3  4  4  The results on Bakhvalov–Shishkin meshes for $${\varepsilon}=10^{-8}$$ are shown in Table 4. The same approximation spaces as on Shishkin meshes have been used. The convergence orders predicted by our theory are obtained. Moreover, one sees that the errors on Bakhvalov–Shishkin meshes are smaller than those on the corresponding Shishkin meshes. This is due to the logarithmic factor which is present only for Shishkin meshes. The difference between Shishkin and Bakhvalov–Shishkin meshes becomes larger with increasing approximation order and already reaches two orders of magnitude for the fourth-order pairs $$Q_4\times Q_3$$ and $$Q_4\times P_3^{\text{disc}}$$. Table 4 Problem (5.1): Error in energy norm $${\left|\!\!\;\left|\!\!\;\left| {(\cdot,\cdot)} \right|\!\!\;\right|\!\!\;\right|}_{\varepsilon}$$ on Bakhvalov–Shishkin meshes with $${\varepsilon}=10^{-8}$$ $$N$$  $$Q_2\times Q_1$$  $$Q_2\times P_1^{\text{disc}}$$  $$Q_3\times Q_2$$  $$Q_3\times P_2^{\text{disc}}$$  $$Q_4\times Q_3$$  $$Q_4\times P_3^{\text{disc}}$$  4  2.2054-01  2.0379-01  3.1063-02  3.1183-02  8.3515-03  8.3526-03  8  5.6584-02  5.6509-02  6.1054-03  6.1147-03  9.5490-04  9.5494-04  16  1.4800-02  1.4845-02  9.3378-04  9.3470-04  7.8247-05  7.8249-05  32  3.8157-03  3.8272-03  1.2856-04  1.2866-04  5.5671-06  5.5672-06  64  9.7042-04  9.7364-04  1.6851-05  1.6864-05  3.7078-07  3.7079-07  Order  1.98  1.97  2.93  2.93  3.91  3.91  Theory  2  2  3  3  4  4  $$N$$  $$Q_2\times Q_1$$  $$Q_2\times P_1^{\text{disc}}$$  $$Q_3\times Q_2$$  $$Q_3\times P_2^{\text{disc}}$$  $$Q_4\times Q_3$$  $$Q_4\times P_3^{\text{disc}}$$  4  2.2054-01  2.0379-01  3.1063-02  3.1183-02  8.3515-03  8.3526-03  8  5.6584-02  5.6509-02  6.1054-03  6.1147-03  9.5490-04  9.5494-04  16  1.4800-02  1.4845-02  9.3378-04  9.3470-04  7.8247-05  7.8249-05  32  3.8157-03  3.8272-03  1.2856-04  1.2866-04  5.5671-06  5.5672-06  64  9.7042-04  9.7364-04  1.6851-05  1.6864-05  3.7078-07  3.7079-07  Order  1.98  1.97  2.93  2.93  3.91  3.91  Theory  2  2  3  3  4  4  5.1.4 Necessity of stabilization The size of the stabilization parameter $$\gamma$$ is arbitrary and our analysis gives only the restriction that it should be constant with respect to $$N$$ and $${\varepsilon}$$. Figure 6 shows for the energy norm of the solution and the $$L^2$$-norm of the pressure alone the dependency on $$\gamma$$ for a Bakhvalov–Shishkin mesh with fixed $$N=64$$, $${\varepsilon}=10^{-8}$$, $$Q_3\times Q_2$$ and stabilization in the full domain. It can clearly be seen that a moderate value of $$\gamma$$ gives best results. Furthermore, we observe almost constant error values on a wide range of values for $$\gamma$$ that include $$\gamma=1$$. This is the reason for choosing this particular value in the numerical studies. Fig. 6. View largeDownload slide Problem (5.1): Dependency of errors on the stabilization parameter $$\gamma$$. Fig. 6. View largeDownload slide Problem (5.1): Dependency of errors on the stabilization parameter $$\gamma$$. Table 5 shows the influence of the grad-div stabilization for $${\varepsilon}=10^{-8}$$ and $$Q_3\times Q_2$$. We consider three choices for $$\gamma$$: first $$\gamma=1$$ in $${\it{\Omega}}$$, second $$\gamma=1$$ in $${\it{\Omega}}_{11}$$ only, and third $$\gamma=0$$ in $${\it{\Omega}}$$. The last choice gives the standard Galerkin method, thus no stabilization. The errors are measured in the corresponding energy norm $$\left|\!\!\;\left|\!\!\;\left| {(\cdot,\cdot)} \right|\!\!\;\right|\!\!\;\right|_{\varepsilon}$$. We observe that the errors for discrete problems without grad-div stabilization are larger than those for both stabilized discrete problems. Moreover, the convergence order for problems without stabilization is less than for problems with stabilization. This is in agreement with Corollary 3.12. Furthermore, we observe that the errors of both stabilized methods are almost indistinguishable, confirming Corollary 3.13. Table 5 Problem (5.1): Error in energy norm $${\left|\!\!\;\left|\!\!\;\left| {(\cdot,\cdot)} \right|\!\!\;\right|\!\!\;\right|}_{\varepsilon}$$ for $$Q_3\times Q_2$$ and several choices of $$\gamma$$ on Bakhvalov–Shishkin meshes with $${\varepsilon}=10^{-8}$$    $$\gamma=1$$ in $${\it{\Omega}}$$  $$\gamma=1$$ in $${\it{\Omega}}_{11}$$  $$\gamma=0$$ in $${\it{\Omega}}$$  $$N$$  Error  Order  Error  Order  Error  Order  4  3.1063-02     3.1075-02     3.1302-02     8  6.1054-03  2.35  6.1061-02  2.35  6.1727-03  2.34  16  9.3378-04  2.71  9.3381-04  2.71  9.6048-04  2.68  32  1.2856-04  2.86  1.2856-04  2.86  1.4026-04  2.78  64  1.6851-05  2.93  1.6851-05  2.93  2.1901-05  2.68  128  2.1567-06  2.97  2.1567-06  2.97  4.1048-06  2.42  196  6.0560-07  2.98  6.0560-07  2.98  1.6071-06  2.20     $$\gamma=1$$ in $${\it{\Omega}}$$  $$\gamma=1$$ in $${\it{\Omega}}_{11}$$  $$\gamma=0$$ in $${\it{\Omega}}$$  $$N$$  Error  Order  Error  Order  Error  Order  4  3.1063-02     3.1075-02     3.1302-02     8  6.1054-03  2.35  6.1061-02  2.35  6.1727-03  2.34  16  9.3378-04  2.71  9.3381-04  2.71  9.6048-04  2.68  32  1.2856-04  2.86  1.2856-04  2.86  1.4026-04  2.78  64  1.6851-05  2.93  1.6851-05  2.93  2.1901-05  2.68  128  2.1567-06  2.97  2.1567-06  2.97  4.1048-06  2.42  196  6.0560-07  2.98  6.0560-07  2.98  1.6071-06  2.20  5.2 Problems with unknown solution In this part of our numerical study, we do not prescribe a solution but give only the data for the problem. Therefore it is not clear whether Assumption 2.1 holds. In order to calculate the errors and the value of $$B_0$$, we replace the exact solution with a reference solution computed on a Bakhvalov–Shishkin mesh with $$N=128$$ and $$Q_4\times Q_3$$. This approach is a variant of the double-mesh principle; see Farrell et al. (2000). Let us start with   −εΔu−((11)⋅∇)u+u+∇p =f in Ω=(0,1)2,divu =0 in Ω,u =0 on ∂Ω, (5.3) where the right-hand side $${\boldsymbol{f}}=(1,x)^{\rm T}$$ was chosen to be simple but not a gradient. First we look at the convergence of our method with full stabilization for $${\varepsilon}=10^{-8}$$ on a Bakhvalov–Shishkin mesh. We present the results for $$Q_2\times Q_1$$ and $$Q_3\times Q_2$$ only as the counterparts with discontinuous pressure approximations behave similarly. The left part of Table 6 shows the error behaviour measured in the energy norm. Clearly the convergence order is far from the theoretically predicted values and does not improve upon increasing the polynomial order. A possible explanation is that Assumption 2.1 is not fulfilled and some corner singularities arise. This is similar to convection–diffusion problems where additional compatibility conditions of the right-hand side have to be fulfilled; see, e.g., Linß & Stynes (2001), Kellogg & Stynes (2005, 2007). Thus, we modify our problem by changing the right-hand side to   f2:=(1−1)(1−x−y)(x−y). Table 6 Problem (5.4): Error in energy norm $${\left|\!\!\;\left|\!\!\;\left| {(\cdot,\cdot)} \right|\!\!\;\right|\!\!\;\right|}_{\varepsilon}$$ on Bakhvalov–Shishkin meshes with $${\varepsilon}=10^{-8}$$ $${\boldsymbol{f}}=(1,x)^{\rm T}$$  $$N$$  $$Q_2\times Q_1$$  $$Q_3\times Q_2$$  8  3.178-02  2.244-02  16  1.787-02  1.258-02  32  1.121-02  8.035-03  64  6.780-03  4.891-03  Order  0.73  0.72  $${\boldsymbol{f}}=(1,x)^{\rm T}$$  $$N$$  $$Q_2\times Q_1$$  $$Q_3\times Q_2$$  8  3.178-02  2.244-02  16  1.787-02  1.258-02  32  1.121-02  8.035-03  64  6.780-03  4.891-03  Order  0.73  0.72  $${\boldsymbol{f}}=(1-x-y)(x-y)(1,-1)^{\rm T}$$  $$N$$  $$Q_2\times Q_1$$  $$Q_3\times Q_2$$  8  9.403-03  1.300-03  16  2.488-03  2.495-04  32  6.322-04  4.794-05  64  1.571-04  9.191-06  Order  2.00  2.38  $${\boldsymbol{f}}=(1-x-y)(x-y)(1,-1)^{\rm T}$$  $$N$$  $$Q_2\times Q_1$$  $$Q_3\times Q_2$$  8  9.403-03  1.300-03  16  2.488-03  2.495-04  32  6.322-04  4.794-05  64  1.571-04  9.191-06  Order  2.00  2.38  This function $${\boldsymbol{f}}_2$$ is zero in the corners of $${\it{\Omega}}$$ which is a minimal compatibility condition in convection–diffusion problems. The convergence results are presented in the right part of Table 6. We do observe the optimal convergence order for $$Q_2\times Q_1$$ but still a reduced order for the higher-order pair. Thus again, not enough compatibility conditions are fulfilled to ensure that Assumption 2.1 holds for the higher derivatives. It remains an open question which conditions are needed. Using $${\boldsymbol{f}}_2$$ we can also look into the constant $$B_0$$. The obtained results are shown in Fig. 7. Again we observe for all values of $$N$$ and $${\varepsilon}$$ that $$B_0$$ is bounded away from zero. In order to get more insight into the behaviour of $$B_0$$, we now change the vector $${\boldsymbol{b}}$$ in (5.3) to   b=(cos⁡(φ)sin⁡(φ)) for $$\varphi\in(0,\pi/4]$$. For fixed values of $$N=32$$ and $${\varepsilon}=10^{-8}$$, we vary in Fig. 8 the angle $$\varphi$$. Note that $$\varphi=0$$ would give characteristic layers that we will investigate in a moment. Although the value of $$B_0$$ depends on $$\varphi$$, it is bounded away from zero over the full range of $$\varphi$$. Fig. 7. View largeDownload slide Problem (5.3): The value of $$B_0$$ for different values of $$N$$ as a function of $$\varepsilon$$ on Bakhvalov–Shishkin meshes: $$Q_2\times Q_1$$ (left) and $$Q_3\times Q_2$$ (right). Fig. 7. View largeDownload slide Problem (5.3): The value of $$B_0$$ for different values of $$N$$ as a function of $$\varepsilon$$ on Bakhvalov–Shishkin meshes: $$Q_2\times Q_1$$ (left) and $$Q_3\times Q_2$$ (right). Fig. 8. View largeDownload slide Problem (5.3): The value of $$B_0$$ as function of $$\varphi$$ on Bakhvalov–Shishkin meshes. Fig. 8. View largeDownload slide Problem (5.3): The value of $$B_0$$ as function of $$\varphi$$ on Bakhvalov–Shishkin meshes. Let us now consider the case of characteristic layers, e.g.,   −εΔu−((10)⋅∇)u+u+∇p =(1−1)(1−x−y)(x−y) in Ω=(0,1)2,divu =0 in Ω,u =0 on ∂Ω. (5.4) Here our theory about the decomposition of $${\boldsymbol{u}}$$ does not hold at all. Nevertheless, we can look at the corresponding convection–diffusion case where a decomposition is known; see Kellogg & Stynes (2005, 2007), and use the analogy of the layer structure. Thus we expect $$u_1$$ to have characteristic layers at $$y=0$$ and $$y=1$$ and $$u_2$$ to have an exponential layer at $$x=0$$. An approximation to the solution is shown in Fig. 9. With this assumption on the layer structure, a suitable layer-adapted mesh can be constructed in the usual way and the analysis presented could also be done. We present only the convergence results in Table 7 that are similar to those of only exponential layers with reduced regularity. Also the constant $$B_0$$ is bounded away from zero for all $${\varepsilon}$$ and $$N$$. Table 7 Problem (5.4): Error in energy norm $${\left|\!\!\;\left|\!\!\;\left| {{(\cdot,\cdot)}} \right|\!\!\;\right|\!\!\;\right|}_{\varepsilon}$$ on Bakhvalov–Shishkin meshes with $${\varepsilon}=10^{-8}$$ $$N$$  $$Q_2\times Q_1$$  $$Q_3\times Q_2$$  8  9.403-03  1.300-03  16  2.488-03  2.495-04  32  6.322-04  4.794-05  64  1.571-04  9.191-06  Order  2.00  2.38  $$N$$  $$Q_2\times Q_1$$  $$Q_3\times Q_2$$  8  9.403-03  1.300-03  16  2.488-03  2.495-04  32  6.322-04  4.794-05  64  1.571-04  9.191-06  Order  2.00  2.38  Fig. 9. View largeDownload slide Problem (5.4): Solution of (5.4) with exponential and characteristic layers in the velocity. Fig. 9. View largeDownload slide Problem (5.4): Solution of (5.4) with exponential and characteristic layers in the velocity. 6. Outlook This article is a starting point for the analysis of more complex domains and problems. The analysis relies on a solution decomposition and therefore a precise knowledge of the layer structure of $${\boldsymbol{u}}$$. Having such a structure, it is relatively easy to extend the techniques of this article to the three-dimensional case or more general domains. Numerical experiments do actually indicate that the solution structure in three dimensions is similar to the two-dimensional situation. In the case of second-order reaction–diffusion problems, more general but smooth domains were investigated in, e.g., Xenophontos & Fulton (2003), Kopteva (2009). Their techniques could be incorporated into the system case. For more general domains and problems, the type of layers can change from exponential to characteristic (see Roos et al., 2008 for these types of layers), but the exact transition from one type to another is not obvious. A technique to handle characteristic layers is the use of problem-specific enriched function spaces where the enrichment could be based on the law of wall; see Krank & Wall (2016) for a recent development. For flows in channels and around obstacles in unbounded domains, one tries to prevent strong boundary layers at the outflow. The idea is to apply suitable boundary conditions that allow us to restrict the domains without significant changes of the flow solutions. Such methods can be found, e.g., in Griffiths (1997), Braack & Mucha (2014). Another topic for further research is to show that $$B_0$$ is actually bounded away from zero, independent of $$N$$ and $${\varepsilon}$$. Having this would provide a tool for proving the uniform error estimation. As indicated in Fig. 8, the value of $$B_0$$ stays away from zero even for the transition from exponential to characteristic layers. Assuming a suitable layer structure for the solution of the Oseen problem with characteristic layers, a convergence behaviour similar to the case with only exponential layers is obtained; compare Tables 6 (right) and 7. Funding K.H. was funded by the grant FR 3052/2-1 of the German Research Foundation (DFG). References Ainsworth M., Barrenechea G. R. & Wachtel A. ( 2015) Stabilization of high aspect ratio mixed finite elements for incompressible flow. SIAM J. Numer. Anal. , 53, 1107– 1120. Google Scholar CrossRef Search ADS   Apel T. ( 1999) Anisotropic Finite Elements: Local Estimates and Applications . Advances in Numerical Mathematics. Stuttgart: B. G. Teubner. Apel T., Knopp T. & Lube G. ( 2008) Stabilized finite element methods with anisotropic mesh refinement for the Oseen problem. Appl. Numer. Math. , 58, 1830– 1843. Google Scholar CrossRef Search ADS   Apel T. & Maharavo Randrianarivony H. ( 2003) Stability of discretizations of the Stokes problem on anisotropic meshes. Math. Comput. Simul. , 61, 437– 447. Google Scholar CrossRef Search ADS   Apel T. & Matthies G. ( 2008) Nonconforming, anisotropic, rectangular finite elements of arbitrary order for the Stokes problem. SIAM J. Numer. Anal. , 46, 1867– 1891. Google Scholar CrossRef Search ADS   Bakhvalov N. S. ( 1969) The optimization of methods of solving boundary value problems with a boundary layer. U.S.S.R. Comput. Math. Math. Phys. , 9, 139– 166. Google Scholar CrossRef Search ADS   Bernardi C. & Maday Y. ( 1997) Spectral methods. Handbook of Numerical Analysis , vol. V, pp. 209– 485. Amsterdam: North-Holland. Google Scholar CrossRef Search ADS   Braack M. & Mucha P. B. ( 2014) Directional do-nothing condition for the Navier-Stokes equations. J. Comput. Math. , 32, 507– 521. Google Scholar CrossRef Search ADS   Brooks A. N. & Hughes T. J. R. ( 1982) Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations. Comput. Methods Appl. Mech. Eng. , 32, 199– 259. Google Scholar CrossRef Search ADS   Burman E. & Linke A. ( 2008) Stabilized finite element schemes for incompressible flow using Scott–Vogelius elements. Appl. Numer. Math. , 58, 1704– 1719. Google Scholar CrossRef Search ADS   Dallmann H., Arndt D. & Lube G. ( 2016) Local projection stabilization for the Oseen problem. IMA J. Numer. Anal. , 36, 796– 823. Google Scholar CrossRef Search ADS   Farrell P. A., Hegarty A. F., Miller J. J. H., O’Riordan E. & Shishkin G. I. ( 2000) Robust Computational Techniques for Boundary Layers . Applied Mathematics (Boca Raton), vol. 16. Boca Raton, FL: Chapman & Hall/CRC. Franca L. P. & Frey S. L. ( 1992) Stabilized finite element methods: II. The incompressible Navier–Stokes equations. Comput. Methods Appl. Mech. Eng. , 99, 209– 233. Google Scholar CrossRef Search ADS   Franz S. ( 2008) Continuous interior penalty method on a Shishkin mesh for convection-diffusion problems with characteristic boundary layers. Comput. Meth. Appl. Mech. Eng. , 197, 3679– 3686. Google Scholar CrossRef Search ADS   Franz S., Höhne K. & Waurick M. ( 2016) A solution decomposition for a singularly perturbed 4th order problem (submitted). Franz S., Linß T. & Roos H.-G. ( 2008) Superconvergence analysis of the SDFEM for elliptic problems with characteristic layers. Appl. Numer. Math. , 58, 1818– 1829. Google Scholar CrossRef Search ADS   Franz S. & Matthies G. ( 2010) Local projection stabilisation on S-type meshes for convection-diffusion problems with characteristic layers. Computing , 87, 135– 167. Google Scholar CrossRef Search ADS   Gelhard T., Lube G., Olshanskii M. A. & Starcke J.-H. ( 2005) Stabilized finite element schemes with LBB-stable elements for incompressible flows. J. Comput. Appl. Math. , 177, 243– 267. Google Scholar CrossRef Search ADS   Girault V. & Raviart P.-A. ( 1986) Finite Element Methods for Navier–Stokes Equations . Berlin: Springer. Google Scholar CrossRef Search ADS   Griffiths D. F. ( 1997) The “no boundary condition” outflow boundary condition. Int. J. Numer. Methods Fluids , 24, 393– 411. Google Scholar CrossRef Search ADS   Hansbo P. & Szepessy A. ( 1990) A velocity-pressure streamline diffusion method for the incompressible Navier–Stokes equations. Comput. Methods Appl. Mech. Eng. , 84, 175– 192. Google Scholar CrossRef Search ADS   Hughes T. J. R., Franca L. P. & Balestra M. ( 1986) A new finite element formulation for computational fluid dynamics. V: circumventing the Babuska–Brezzi condition: a stable Petrov–Galerkin formulation of the Stokes problem accommodating equal-order interpolations. Comput. Methods Appl. Mech. Eng. , 59, 85– 99. Google Scholar CrossRef Search ADS   Jenkins E. W., John V., Linke A. & Rebholz L. G. ( 2014) On the parameter choice in grad-div stabilization for the Stokes equations. Adv. Comput. Math. , 40, 491– 516. Google Scholar CrossRef Search ADS   Johnson C. & Saranen J. ( 1986) Streamline diffusion methods for the incompressible Euler and Navier–Stokes equations. Math. Comput. , 47, 1– 18. Google Scholar CrossRef Search ADS   Kellogg R. B. & Stynes M. ( 2005) Sharpened and corrected version of: corner singularities and boundary layers in a simple convection-diffusion problem. J. Differ. Equ. , 213, 81– 120. Google Scholar CrossRef Search ADS   Kellogg R. B. & Stynes M. ( 2007) Sharpened bounds for corner singularities and boundary layers in a simple convection-diffusion problem. Appl. Math. Lett. , 20, 539– 544. Google Scholar CrossRef Search ADS   Kopteva N. ( 2009) Numerical analysis of a 2d singularly perturbed semilinear reaction-diffusion problem. Numerical Analysis and Its Applications  ( Margenov S. Vulkov L. G. & Waśniewski J. eds). Lecture Notes in Computer Science, vol. 5434. Berlin: Springer, pp. 80– 91. Google Scholar CrossRef Search ADS   Krank B. & Wall W. A. ( 2016) A new approach to wall modeling in LES of incompressible flow via function enrichment. J. Comput. Phys. , 316, 94– 116. Google Scholar CrossRef Search ADS   Linke A. ( 2014) On the role of the Helmholtz decomposition in mixed methods for incompressible flows and a new variational crime. Comput. Methods Appl. Mech. Eng. , 268, 782– 800. Google Scholar CrossRef Search ADS   Linke A. & Rebholz L. G. ( 2013) On a reduced sparsity stabilization of grad-div type for incompressible flow problems. Comput. Methods Appl. Mech. Eng. , 261/262, 142– 153. Google Scholar CrossRef Search ADS   Linß T. ( 1999) An upwind difference scheme on a novel Shishkin-type mesh for a linear convection-diffusion problem. J. Comput. Appl. Math. , 110, 93– 104. Google Scholar CrossRef Search ADS   Linß T. ( 2000) Analysis of a Galerkin finite element method on a Bakhvalov-Shishkin mesh for a linear convection-diffusion problem. IMA J. Numer. Anal. , 20, 621– 632. Google Scholar CrossRef Search ADS   Linß T. ( 2010) Layer-Adapted Meshes for Reaction-Convection-Diffusion Problems . Lecture Notes in Mathematics, vol. 1985. Berlin: Springer. Google Scholar CrossRef Search ADS   Linß T. & Stynes M. ( 2001) Asymptotic analysis and Shishkin-type decomposition for an elliptic convection-diffusion problem. J. Math. Anal. Appl. , 261, 604– 632. Google Scholar CrossRef Search ADS   Ludwig L. ( 2017) https://github.com/SOFE-Developers/SOFE, last accessed 17 March 2017. Matthies G., Lube G. & Röhe L. ( 2009) Some remarks on residual-based stabilisation of inf-sup stable discretisations of the generalised Oseen problem. Comput. Methods Appl. Math. , 9, 368– 390. Google Scholar CrossRef Search ADS   Matthies G. & Tobiska L. ( 2015) Local projection type stabilization applied to inf–sup stable discretizations of the Oseen problem. IMA J. Numer. Anal. , 35, 239– 269. Google Scholar CrossRef Search ADS   Miller J. J. H., O’Riordan E. & Shishkin G. I. ( 1996) Fitted Numerical Methods for Singular Perturbation Problems: Error Estimates in the Maximum Norm for Linear Problems in One and Two Dimensions.  River Edge, NJ: World Scientific Publishing Co. Google Scholar CrossRef Search ADS   Olshanskii M. A. ( 2002) A low order Galerkin finite element method for the Navier-Stokes equations of steady incompressible flow: a stabilization issue and iterative methods. Comput. Methods Appl. Mech. Eng. , 191, 5515– 5536. Google Scholar CrossRef Search ADS   Olshanskii M. A., Lube G., Heister T. & Löwe J. ( 2009) Grad-div stabilization and subgrid pressure models for the incompressible Navier-Stokes equations. Comput. Methods Appl. Mech. Eng. , 198, 3975– 3988. Google Scholar CrossRef Search ADS   Olshanskii M. A. & Reusken A. ( 2004) Grad-div stabilization for the Stokes equations. Math. Comput. , 73, 1699– 1718. Google Scholar CrossRef Search ADS   Roos H.-G. ( 2012) Robust numerical methods for singularly perturbed differential equations: a survey covering 2008–2012. ISRN Appl. Math. , 2012, Article ID 379547. Roos H.-G. & Linß T. ( 1999) Sufficient conditions for uniform convergence on layer-adapted grids. Computing , 63, 27– 45. Google Scholar CrossRef Search ADS   Roos H.-G., Stynes M. & Tobiska L. ( 2008) Robust Numerical Methods for Singularly Perturbed Differential Equations . Springer Series in Computational Mathematics, vol. 24, 2nd edn. Berlin: Springer. Schötzau D. & Schwab C. ( 1998) Mixed $$hp$$-FEM on anisotropic meshes. Math. Models Methods Appl. Sci. , 8, 787– 820. Google Scholar CrossRef Search ADS   Schötzau D., Schwab C. & Stenberg R. ( 1999) Mixed $$hp$$-FEM on anisotropic meshes. II. Hanging nodes and tensor products of boundary layer meshes. Numer. Math. , 83, 667– 697. Google Scholar CrossRef Search ADS   Stynes M. & O’Riordan E. ( 1997) A uniformly convergent Galerkin method on a Shishkin mesh for a convection-diffusion problem. J. Math. Anal. Appl. , 214, 36– 54. Google Scholar CrossRef Search ADS   Stynes M. & Tobiska L. ( 2003) The SDFEM for a convection-diffusion problem with a boundary layer: optimal error analysis and enhancement of accuracy. SIAM J. Numer. Anal. , 41, 1620– 1642. Google Scholar CrossRef Search ADS   Stynes M. & Tobiska L. ( 2008) Using rectangular $$Q_p$$ elements in the SDFEM for a convection-diffusion problem with a boundary layer. Appl. Numer. Math. , 58, 1709– 1802. Google Scholar CrossRef Search ADS   Tezduyar T. E., Mittal S., Ray S. E. & Shih R. ( 1992) Incompressible flow computations with stabilized bilinear and linear equal order interpolation velocity pressure elements. Comput. Methods Appl. Mech. Eng. , 95, 221– 242. Google Scholar CrossRef Search ADS   Tobiska L. & Lube G. ( 1991) A modified streamline diffusion method for solving the stationary Navier–Stokes equations. Numer. Math. , 59, 13– 29. Google Scholar CrossRef Search ADS   Xenophontos C. & Fulton S. R. ( 2003) Uniform approximation of singularly perturbed reaction-diffusion problems by the finite element method on a Shishkin mesh. Numer. Methods Partial Differ. Equ. , 19, 89– 111. Google Scholar CrossRef Search ADS   © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

# Grad-div stabilized discretizations on S-type meshes for the Oseen problem

, Volume 38 (1) – Jan 1, 2018
31 pages

Publisher
Oxford University Press
ISSN
0272-4979
eISSN
1464-3642
D.O.I.
10.1093/imanum/drw069
Publisher site
See Article on Publisher Site

### Abstract

Abstract We consider discretizations of the singularly perturbed Oseen equations on properly layer-adapted meshes. Using a suitable solution decomposition, we are able to prove optimal convergence orders in the associated energy norm for grad-div stabilized finite element methods in a general setting. Two families of pairs of discrete function spaces, namely $$Q_k\times Q_{k-1}$$ and $$Q_k\times P_{k-1}^{\text{disc}}$$, $$k\ge 2$$, are investigated in detail. The usage of a standard nonstabilized Galerkin method reduces the order by 1 while stabilization outside the layers is enough to regain the full optimal order. 1. Introduction The numerical solution of the time-dependent Navier–Stokes equations, especially for high Reynolds numbers, is still a great challenge of computational fluid dynamics. Using time discretization and linearization, the Oseen problem arises as an important subproblem that is for high Reynolds numbers singularly perturbed. It is well known that standard finite element methods applied to the Oseen problem may suffer from two difficulties: the generally dominating convection and the incompatibility between the approximation spaces for velocity and pressure. The streamline-upwind Petrov–Galerkin method (SUPG), introduced in Brooks & Hughes (1982), and the pressure-stabilized Petrov–Galerkin method (PSPG), introduced in Hughes et al. (1986), Johnson & Saranen (1986), Tezduyar et al. (1992), provide a unique frame to tackle both problems simultaneously. An additional elementwise stabilization of the divergence constraint, abbreviated as grad-div stabilization, increases the robustness; see Franca & Frey (1992), Hansbo & Szepessy (1990), Tobiska & Lube (1991). In order to ensure the strong consistency of these residual-based stabilization methods, various terms have to be added to the weak formulation. The PSPG term can be skipped if inf–sup stable pairs of finite element spaces for approximating velocity and pressure are used. This results in reduced stabilized schemes; see Gelhard et al. (2005), Matthies et al. (2009). The coupling between velocity and pressure that is caused by the SUPG term complicates the analysis and the grad-div stabilization seems to be even more important; compare Burman & Linke (2008), Gelhard et al. (2005), Linke & Rebholz (2013), Jenkins et al. (2014), Matthies et al. (2009), Olshanskii & Reusken (2004). On general meshes, the use of stabilization techniques provides semirobust error estimates in the sense of Roos (2012). This means that the constants in the error bounds are independent of the small viscosity parameter while the higher-order norms of the solution that appear strongly increase with decreasing viscosity. The role of grad-div stabilization has been investigated by several authors. The combination of SUPG and grad-div stabilization applied to the stationary Navier–Stokes equations was considered in Olshanskii (2002). The grad-div term enhances the accuracy of the numerical solution and the convergence properties of iterative solvers if the parameter of the grad-div stabilization is chosen properly. The influence of grad-div stabilization for the Stokes equations has been studied in Olshanskii & Reusken (2004), Jenkins et al. (2014). Within the scope of variational multiscale methods, the grad-div stabilization as a pressure subgrid model has been investigated in Olshanskii et al. (2009). The importance of the grad-div stabilization inside local projection stabilizations was stated in Dallmann et al. (2016). Standard methods for stationary singularly perturbed convection–diffusion problems do also lack stability; see, e.g. the book Roos et al. (2008). For these problems, knowledge of the layer structure of their solutions allows creation and usage of properly layer-adapted meshes; see the survey book Linß (2010). Such special meshes have quite a long history. First ideas were proposed by Bakhvalov (1969) using exponentially fitted boundary meshes. The piecewise uniform Shishkin meshes (Miller et al., 1996; Stynes & O’Riordan, 1997) are easier in structure and in the analysis of methods based on them, but their approximation quality is inferior to Bakhvalov meshes. The combination of Bakhvalov’s graded fine meshes with Shishkin’s choice of the transition point (Linß, 1999, 2000) led to general S-type meshes (Roos & Linß, 1999) where the mesh inside the layers is defined with the help of a mesh-generating function. The use of layer-adapted meshes allows for robust or uniform error estimates in the sense of Roos (2012). Hence, the error bounds are completely independent of the viscosity parameter. On such adapted meshes, the oscillations induced by boundary conditions are strongly reduced. However, the lack of stability is still there. Therefore, stabilization methods like SUPG (also called SDFEM), local projection stabilization or interior penalty stabilization are used in addition. Some results for those combinations can be found in Roos et al. (2008), Franz et al. (2008), Franz (2008), Franz & Matthies (2010), Stynes & Tobiska (2008). This article now combines these two approaches. We apply a grad-div stabilization on layer-adapted meshes. Using suitable pairs of function spaces for approximating velocity and pressure, we prove error bounds that are uniform in $$\varepsilon$$ and yield an optimal convergence rate. Stabilized finite element methods for the Oseen problem on anisotropic meshes have been analysed in Apel et al. (2008). The authors considered equal-order interpolation for velocity and pressure, but S-type meshes were explicitly excluded. We consider the stationary Oseen equations   −εΔu−(b⋅∇)u+cu+∇p =f in Ω=(0,1)2,divu =0 in Ω,u =0 on Γ=∂Ω, (1.1) where the perturbation parameter $$\varepsilon$$ with $$0<\varepsilon\ll 1$$ corresponds to the viscosity, the convection field $${\boldsymbol{b}}=(b_1,b_2)\in W^{1,\infty}({\it{\Omega}})^2$$ with $${\text{div}}{\boldsymbol{b}}=0$$ fulfils $$b_1\geq \beta_1 >0,\ b_2\geq \beta_2 >0$$ with some positive constants $$\beta_1,\beta_2$$ and $$0<c_0\leq c\in L^\infty({\it{\Omega}})$$ is assumed. Then $${\boldsymbol{u}}$$ describes the velocity of our problem and $$p$$ the pressure. Since $$0<\varepsilon\ll 1$$, we expect layers at $$x=0$$ and $$y=0$$. A typical solution of (1.1) can be seen in Fig. 1. Fig. 1. View largeDownload slide Typical solution of (1.1) with exponential layers in the velocity. Fig. 1. View largeDownload slide Typical solution of (1.1) with exponential layers in the velocity. As the solution exhibits boundary layers, the use of quasi-uniform meshes gives accurate approximations of $$(1.1)$$ only if the mesh size is of the order of the perturbation parameter $$\varepsilon$$. Since this gives a prohibitive restriction, layer-adapted meshes will be used for the discretizations. We are aware that the considered problem class is not a prototype for Oseen problems coming from applications as we utilize a solution decomposition (see Assumption 2.1) which is usually unavailable in real-world problems. However, this article is intended as a starting point for the robust numerical analysis of such problems on layer-adapted meshes. The article is organized as follows. Section 2 introduces a weak formulation, a solution decomposition and layer-adapted meshes. In Section 3, we discretize and analyse our problem in an abstract setting without specifying the finite element spaces. We prove optimal convergence orders for the fully grad-div stabilized method and for stabilization only outside the layers while the order is reduced by 1 for the standard Galerkin method. We present in Section 4, two pairs of spaces that will suit our purpose. Finally, Section 5 presents numerical simulations confirming the theoretical results. Notation: On a domain $$D,$$ we denote the $$L^2$$-norm by $$\|\cdot\|_D$$ and the $$L^\infty$$-norm by $$\|\cdot\|_{\infty,D}$$. In the case of $$D={\it{\Omega}},$$ we drop the reference to the domain from the notation. Moreover, $$C$$ will denote a generic constant independent of the perturbation parameter $${\varepsilon}$$ and the mesh parameter $$N$$. 2. Solution decomposition and layer-adapted meshes Let us start with a weak formulation of (1.1). We define $$V:=H_0^1({\it{\Omega}})^2$$ as the space of vector functions with components in the Sobolev space $$H_0^1({\it{\Omega}})$$ of $$L^2({\it{\Omega}})$$-functions with weak derivatives in $$L^2({\it{\Omega}})$$ and vanishing boundary traces, and $$Q:=L_0^2({\it{\Omega}})$$ as the space of $$L^2({\it{\Omega}})$$-functions with vanishing integral mean value. A weak formulation then reads as follows: Find $$({\boldsymbol{u}},p)\in V\times Q$$ such that for all $$({\boldsymbol{v}},q)\in V\times Q$$,   ε(∇u,∇v)−((b⋅∇)u,v)+(cu,v)−(p,divv) =(f,v),(q,divu) =0. (2.1) Note that $${\text{div}} {\boldsymbol{b}}=0$$ implies   ((b⋅∇)u,v)=−((b⋅∇)v,u)∀u,v∈V, (2.2) in particular,   ((b⋅∇)v,v)=0∀v∈V. (2.3) Thus, applying the Lax–Milgram lemma in the subspace of divergence-free functions, we establish the unique velocity field $${\boldsymbol{u}}$$. The unique pressure $$p\in Q$$ such that $$({\boldsymbol{u}},p)$$ solves (2.1) follows from the Babuška–Brezzi condition for the pair $$(V,Q)$$; see Girault & Raviart (1986, Ch. I). We will make some assumptions about a solution decomposition of $${\boldsymbol{u}}=(u_1,\,u_2)$$ and $$p$$ in order to carry out our analysis. We assume that $$u_1$$ and $$u_2$$ will consist of layer and regular components while the pressure $$p$$ will not show any layer behaviour. To be more precise, we assume the following bounds. Assumption 2.1 The solution $${\boldsymbol{u}}=(u_1,u_2)$$ can be decomposed as   u1 =g1+v1+w1+z1andu2=g2+v2+w2+z2, where we have for fixed $$k\in \mathbb{N}$$, all $$x,y\in [0,1]$$ and $$0\leq i+j \leq k+1$$ the pointwise estimates   |∂xi∂yjg1(x,y)| ≤C,|∂xi∂yjg2(x,y)| ≤C,|∂xi∂yjv1(x,y)| ≤Cε1−ie−β1x/ε,|∂xi∂yjv2(x,y)| ≤Cε−ie−β1x/ε,|∂xi∂yjw1(x,y)| ≤Cε−ie−β2y/ε,|∂xi∂yjw2(x,y)| ≤Cε1−ie−β2y/ε,|∂xi∂yjz1(x,y)| ≤Cε1−(i+j)e−β1x/εe−β2y/ε,|∂xi∂yjz2(x,y)| ≤Cε1−(i+j)e−β1x/εe−β2y/ε, where $$g_1,\,g_2$$ are the regular parts, $$v_1,\,v_2$$ the exponential boundary layers in the $$x$$-direction, $$w_1,\,w_2$$ the exponential boundary layers in the $$y$$-direction and $$z_1,\,z_2$$ the corner layers. The bounds   |∂xi∂yjp(x,y)|≤C,0≤i+j≤k, are assumed for the pressure $$p$$. Remark 2.2 The divergence constraint guarantees the existence of a scalar stream function $$\psi$$ such that $${\boldsymbol{u}}=(\partial_y \psi, -\partial_x \psi)$$. Inserting this into the Oseen equations and applying the (two-dimensional) $${\text{curl}}$$ operator leads to a fourth-order problem. In Franz et al. (2016), it was established that its solution has a principle structure that allows us to derive the assumptions stated above providing that enough compatibility conditions on the data of our problem are fulfilled. The nature of these additional assumptions is similar to those of convection–diffusion problems (Kellogg & Stynes, 2005, 2007). Furthermore, our numerical studies confirm the assumed decomposition. Concerning the pressure, our analysis does also work if the pressure has strong boundary layers similar to those of $$u_1$$ or $$u_2$$. Then we would estimate $$p$$ similarly to $${\boldsymbol{u}}$$ which is possible if we use the same layer-adapted mesh for all solution components. With Assumption 2.1 in mind, we are able to construct a priori adapted meshes. The solution $${\boldsymbol{u}}$$ has weak and strong layers. In order to resolve the strong layers, we use fine meshes along $$y=0$$ for $$u_1$$ and along $$x=0$$ for $$u_2$$. Weak layers are not visible at first sight. In the process of a numerical analysis, their unpleasant $${\varepsilon}$$-dependence comes into play via derivatives. Thus, the numerical analysis works best if the weak layers are also resolved. Therefore, we will use layer-adapted meshes that are fine at $$x=0$$ and at $$y=0$$ for both components $$u_1$$ and $$u_2$$. For simplicity, these meshes are also used for $$p$$ although we did assume that $$p$$ has no boundary layers. We will use so-called S-type meshes; see Roos & Linß (1999). They are constructed using transition points where the mesh changes between fine and coarse. To be more precise, let us assume the number $$N$$ of mesh intervals in each coordinate direction is even. Then we define the transition parameters   λx:=min{12,σεβ1ln⁡N}andλy:=min{12,σεβ2ln⁡N} with some user-chosen positive parameter $$\sigma$$. For the sake of mere simplicity in our subsequent analysis, we assume that   ε<min{β1,β2}2σln⁡N⇔N<exp⁡min{β1,β2}2σε such that $$\lambda_x$$ and $$\lambda_y$$ are smaller than $$1/2$$. This is the typical case for (1.1) as otherwise a standard analysis on an equidistant mesh suffices. Using the transition parameters $$\lambda_x$$ and $$\lambda_y$$, the domain $${\it{\Omega}}$$ is divided into the subdomains $${\it{\Omega}}_{11},{\it{\Omega}}_{21},{\it{\Omega}}_{12}$$ and $${\it{\Omega}}_{22}$$ according to   Ω11 :=[λx,1]×[λy,1],Ω12:=[0,λx]×[λy,1],Ω21 :=[λx,1]×[0,λy]andΩ22:=[0,λx]×[0,λy]; see Fig. 2. Fig. 2. View largeDownload slide Decomposition of $${\it{\Omega}}$$ into subregions. Fig. 2. View largeDownload slide Decomposition of $${\it{\Omega}}$$ into subregions. Note that $${\it{\Omega}}_{11}$$ covers the nonlayer region, $${\it{\Omega}}_{22}$$ the corner layers and $${\it{\Omega}}_{12}$$ and $${\it{\Omega}}_{21}$$ the two boundary layer regions. With the help of a mesh-generating function$$\phi$$ which is monotonically increasing with $$\phi(0)=0$$ and $$\phi(1/2)=\ln N$$, we construct the layer-adapted mesh as a tensor product mesh with nodes according to   xi :={σεβ1ϕ(iN),i=0,…,N/2,1−2(1−λx)(1−iN),i=N/2,…,N,  and   yj :={σεβ2ϕ(jN),j=0,…,N/2,1−2(1−λy)(1−jN),j=N/2,…,N.  Thus, the dimensions of cells in the layer regions depend on $$\phi$$ while the mesh is uniform outside. Related to the mesh-generating function $$\phi,$$ we define the mesh-characterizing function$$\psi$$ by   ϕ=−ln⁡ψwhich impliesψ=e−ϕ. The derivative of this function will be used to characterize the convergence behaviour on general S-type meshes. Some examples of S-type meshes from Roos & Linß (1999) are presented in Table 1. Table 1 Some examples of mesh-generating and mesh-characterizing functions of S-type meshes Name  $$\phi(t)$$  $$\max\phi'$$  $$\psi(t)$$  $$\max|\psi'|$$  Shishkin mesh  $$2t\ln N$$  $$2\ln N$$  $$N^{-2t}$$  $$2\ln N$$  Bakhvalov–Shishkin mesh  $$-\!\ln(1\!-\!2t(1\!-\!N^{-1}))$$  $$2N$$  $$1\!-\!2t(1\!-\!N^{-1})$$  2  Modified Bakhvalov–Shishkin mesh  $$\frac{t}{q-t},\,q=\frac{1}{2}(1+\frac{1}{\ln N})$$  $$3(\ln N)^2$$  $${\,\mathit{e}}^{-{t}/{q-t}}$$  $$3/(2q)\leq 3$$  Name  $$\phi(t)$$  $$\max\phi'$$  $$\psi(t)$$  $$\max|\psi'|$$  Shishkin mesh  $$2t\ln N$$  $$2\ln N$$  $$N^{-2t}$$  $$2\ln N$$  Bakhvalov–Shishkin mesh  $$-\!\ln(1\!-\!2t(1\!-\!N^{-1}))$$  $$2N$$  $$1\!-\!2t(1\!-\!N^{-1})$$  2  Modified Bakhvalov–Shishkin mesh  $$\frac{t}{q-t},\,q=\frac{1}{2}(1+\frac{1}{\ln N})$$  $$3(\ln N)^2$$  $${\,\mathit{e}}^{-{t}/{q-t}}$$  $$3/(2q)\leq 3$$  In order to analyse the convergence of numerical methods on these meshes, the mesh-generating function $$\phi$$ has to fulfil an additional assumption; see Roos & Linß (1999). Assumption 2.3 Let the mesh-generating function $$\phi$$ be piecewise differentiable such that   maxt∈[0,1/2]ϕ′(t)≤CN or equivalently maxt∈[0,1/2]|ψ′(t)|ψ(t)≤CN (2.4) is fulfilled. Note that (2.4) is satisfied for all meshes given in Table 1. We will now use (2.4) to bound the mesh width from above inside the layers. Let $$h_i:=x_i-x_{i-1}$$ and $$t_i=i/N$$ for $$i=1,\dots,N/2$$. Then it holds for $$t\in[t_{i-1},t_i]$$ (with $$\max \phi'$$ taken over $$t\in[t_{i-1},t_i]$$),   ψ(ti)=e−ϕ(ti)=e−(ϕ(ti)−ϕ(t))e−ϕ(t) ≥e−(ϕ(ti)−ϕ(ti−1))ψ(t)≥e−N−1maxϕ′ψ(t)≥Cψ(t), (2.5) where we used (2.4) for the last estimate. Furthermore, we have   x=σεβ1ϕ(t)=−σεβ1ln⁡ψ(t)which givesψ(t)=e−β1x/(σε). Using this, the monotonicity of $$\psi$$ and (2.5), we obtain for $$i=1,\ldots,N/2$$ and $$x\in[x_{i-1},x_i]$$,   hi =σεβ1(ϕ(ti)−ϕ(ti−1))≤σβ1εN−1maxt∈[ti−1,ti]ϕ′(t)≤σβ1εN−1(maxt∈[ti−1,ti]|ψ′(t)|)/ψ(ti) ≤CεN−1(maxt∈[ti−1,ti]|ψ′(t)|)/ψ(t)≤CεN−1max|ψ′|eβ1x/(σε) (2.6) with $$\max|\psi'|:=\max\limits_{t\in[0,1/2]}|\psi'(t)|$$. Similarly, we get for $$j=1,\dots,N/2$$,   ℓj:=yj−yj−1 ≤CεN−1max|ψ′|eβ2y/(σε) (2.7) with $$y\in[y_{j-1},y_j]$$. Of course, the simpler bounds   hi ≤CεN−1maxϕ′ ≤Cε,i=1,…,N/2,ℓj ≤CεN−1maxϕ′ ≤Cε,j=1,…,N/2, follow also from (2.4). Let us denote by   h:=maxi=1,…,N/2hiandℓ:=maxj=1,…,N/2ℓj the maximal dimensions of cells in the layer region perpendicular to the boundaries. We will see later on that the errors will be bounded by powers of the quantity   E(N,ε):=h+ℓ+N−1max|ψ′|, (2.8) which depends on $$h$$, $$\ell$$, $$N$$ and $$\psi$$. By the above considerations, we have   E(N,ε)≤{Cε+N−1max|ψ′|for a general S-type mesh,CN−1ln⁡N,for a Shishkin mesh.  Thus, for $${\varepsilon}\leq C N^{-1}$$ this quantity is (almost) of order 1 and especially for the Bakhvalov–Shishkin and the modified Bakhvalov–Shishkin meshes we obtain   E(N,ε)≤CN−1. We denote in the following by $$\tau_{ij}=[x_{i-1},x_i]\times[y_{j-1},y_j]$$ a specific element with dimensions $$h_i$$ and $$\ell_j$$ and by $$\tau$$ a generic mesh rectangle. Note that the mesh cells are assumed to be closed. 3. Abstract problem and its convergence analysis We are interested in conforming discretizations of the Oseen problem, i.e., we choose finite element spaces $$V_N\subset V$$ and $$Q_N\subset Q$$. For this section, we will not specify these discrete spaces and consider examples in Section 4. The standard Galerkin formulation of (2.1) reads as follows: Find $$({\boldsymbol{u}}_N,p_N)\in V_N\times Q_N$$ such that   ε(∇uN,∇vN)−((b⋅∇)uN,vN)+(cuN,vN)−(pN,divvN) =(f,vN),(qN,divuN) =0 (3.1) for all $$({\boldsymbol{v}}_N,q_N)\in V_N\times Q_N$$. In order to ensure the unique solvability of the discrete problem (3.1), the discrete spaces $$V_N$$ and $$Q_N$$ have to fulfil a compatibility condition. Assumption 3.1 For any given $$N$$ and $${\varepsilon}$$, there exists a positive number $$\beta_0(N,{\varepsilon})$$ such that   infqN∈QNsupvN∈VN(qN,divvN)‖∇vN‖‖qN‖≥β0(N,ε)>0 (3.2) holds true. Note that in the assumption above, we assume only a lower positive bound for fixed $$N$$ and $${\varepsilon}$$, not a uniform bound in $$N$$ and $${\varepsilon}$$ that would be an inf–sup or Babuška–Brezzi condition. The $${\varepsilon}$$-dependence of $$\beta_0(N,{\varepsilon})$$ can be caused only via the underlying mesh which depends on $${\varepsilon}$$. If $${\boldsymbol{f}}$$ contains a large irrotational part in the sense of the Helmholtz decomposition, the numerical accuracy of (3.1) is known to be poor; see Linke (2014). In order to overcome this, the grad-div stabilization term $$(\gamma{\text{div}}{\boldsymbol{u}}_N,\gamma{\text{div}}{\boldsymbol{v}}_N)$$ will be added where $$\gamma\in L^{\infty}({\it{\Omega}})$$ fulfils $$\gamma\ge\gamma_0>0$$ in $$\overline{{\it{\Omega}}}$$ unless specified otherwise. Our stabilized discrete formulation looks as follows: Find $$({\boldsymbol{u}}_N,p_N)\in V_N\times Q_N$$ such that   ε(∇uN,∇vN)−((b⋅∇)uN,vN)+(cuN,vN)+(γdivuN,γdivvN)−(pN,divvN) =(f,vN),(qN,divuN) =0 (3.3) for all $$({\boldsymbol{v}}_N,q_N)\in V_N\times Q_N$$. We define the bilinear form   A((u,p),(v,q)):=ε(∇u,∇v)−((b⋅∇)u,v)+(cu,v)+(γdivu,γdivv)−(p,divv)+(q,divu) (3.4) on the product space $$V\times Q$$ and the norms   |||v|||:=(ε‖∇v‖2+c0‖v‖2+‖γdivv‖2)1/2,|||(v,q)|||:=(|||v|||2+α‖q‖2)1/2 (3.5) on $$V$$ and $$V\times Q$$ with $$\alpha=\alpha(\varepsilon,c_0,c_\infty,b_\infty)$$ given by   α=1max{ε+c∞2c0CF2,b∞2}, where $$C_{\rm F}$$ denotes the constant in the Friedrichs’ inequality for the domain $${\it{\Omega}}$$, $$c_\infty:=\|c\|_\infty$$ and $$b_{\infty}:=\|{\boldsymbol{b}}\|_\infty$$. Note that the term $$c_\infty^2/c_0$$ is bounded since we assumed $$0 <c_0\le c$$ in $${\it{\Omega}}$$ and $$c\in L^{\infty}({\it{\Omega}})$$ in the introduction. Lemma 3.2 Let $$({\boldsymbol{u}},p)$$ and $$({\boldsymbol{u}}_N,p_N)$$ denote the solutions of (2.1) and (3.3), respectively. Then, the Galerkin orthogonality   A((u−uN,p−pN),(vN,qN))=0 (3.6) holds true for all $$({\boldsymbol{v}}_N,q_N)\in V_N\times Q_N$$. Proof. Subtract (3.3) from (2.1), use $$V_N\subset V$$, $$Q_N\subset Q$$ and $${\text{div}}{\boldsymbol{u}} =0$$. □ Our analysis relies on an upper bound for the parameter $$\gamma$$ of the grad-div stabilization. Assumption 3.3 Let $$\gamma$$ be chosen such that   γ∞2≤Cγ2α holds true with some positive constant $$C_\gamma$$ where $$\gamma_\infty:=\|\gamma\|_\infty$$. Let us define for any $$q_N\in Q_N$$ the positive number   B0(N,ε,qN):={1‖qN‖supvN∈VN(qN,divvN)‖∇vN‖,qN≠0,1,qN=0.  (3.7) Note that by definition and Assumption 3.1 it follows that   B0(N,ε,qN)≥β0(N,ε). The stability of the bilinear form $$A$$ can be shown similarly to Matthies & Tobiska (2015, Lemma 3.1). Since we have to take into consideration that here the constant $$\alpha$$ is defined differently, $$c$$ and $$\gamma$$ are no longer constant and we will use $$B_0(N,{\varepsilon},q_N)$$ instead of $$\beta_0(N,{\varepsilon})$$; the proof will be given here. Lemma 3.4 Under Assumptions 3.1 and 3.3, there exists for each $$N>0$$, $${\varepsilon}>0$$ and each pair $$({\boldsymbol{v}}_N,q_N)\in V_N\times Q_N$$ a positive number $$B(N,{\varepsilon},q_N)$$ such that   |||(vN,qN)|||≤1B(N,ε,qN)sup(wN,rN)A((vN,qN),(wN,rN))|||(wN,rN)||| (3.8) holds true. Proof. Let $$({\boldsymbol{v}}_N,q_N)$$ be an arbitrary element of $$V_N \times Q_N$$ and $$B_0=B_0(N,{\varepsilon},q_N)$$. We obtain immediately   A((vN,qN),(vN,qN)) ≥ε‖∇vN‖2+c0‖vN‖2+‖γdivvN‖2=|||vN|||2. (3.9) It follows from the definition (3.7) of $$B_0(N,{\varepsilon},q_N)$$ that for $$q_n\in Q_n$$ there exists a $${\boldsymbol{z}}_N={\boldsymbol{z}}_N(q_N) \in V_N$$ such that   (divzN,qN)=−α‖qN‖2,‖∇zN‖=αB0‖qN‖, (3.10) where we already have used the stronger estimate with $$B_0$$ instead of $$\beta_0(N,{\varepsilon})$$. Hence, we have   A((vN,qN),(zN,0))=ε(∇vN,∇zN)−((b⋅∇)vN,zN)+(cvN,zN)+(γdivvN,γdivzN)+α‖qN‖2. Now, we estimate the first four terms of that equation. We obtain directly using Friedrichs’ inequality and (3.10),   |ε(∇vN,∇zN)+(cvN,zN)| ≤(ε‖∇vN‖2+c0‖vN‖2)1/2(ε‖∇zN‖2+c∞2c0‖zN‖2)1/2 ≤(ε+c∞2c0CF2)1/2αB0‖qN‖|||vN||| ≤α4‖qN‖2+1B02|||vN|||2 and   |−((b⋅∇)vN,zN)|=|((b⋅∇)zN,vN)| ≤b∞‖∇zN‖‖vN‖≤α4‖qN‖2+αb∞2B02‖vN‖2 ≤α4‖qN‖2+1B02|||vN|||2, where the definition of $$\alpha$$ was used. For the stabilization term, we are using (3.10) and obtain   (γdivvN,γdivzN) ≤‖γdivvN‖‖γdivzN‖ ≤|||vN|||2γ∞‖∇zN‖≤|||vN|||Cγ1/2α1/2αB0‖qN‖ ≤α4‖qN‖2+CγB02|||vN|||2 with the constant $$C_\gamma$$ from Assumption 3.3. Combining above estimates we get   A((vN,qN),(zN,0))≥α4‖qN‖2−C~4|||vN|||2, (3.11) where $$\tilde C:=4(2+C_\gamma)/B_0^2$$. Thus, we define for $$({\boldsymbol{v}}_N,q_N) \in V_N\times Q_N$$ the pair $$({\boldsymbol{w}}_N,r_N)\in V_N\times Q_N$$ by   (wN,rN):=(vN,qN)+41+C~(zN,0) and obtain with (3.9) and (3.11),   A((vN,qN),(wN,rN)) ≥α1+C~‖qN‖2+(1−C~1+C~)|||vN|||2=11+C~|||(vN,qN)|||2. It remains to show that $$\left|\!\!\;\left|\!\!\;\left| {({\boldsymbol{w}}_N,r_N)} \right|\!\!\;\right|\!\!\;\right|\leq C \left|\!\!\;\left|\!\!\;\left| {({\boldsymbol{v}}_N,q_N)} \right|\!\!\;\right|\!\!\;\right|$$. To this end, we estimate   |||(wN,rN)|||≤|||(vN,qN)|||+41+C~|||zN||| and with Friedrichs’ inequality and (3.10),   |||zN|||2=ε‖∇zN‖2+c0‖zN‖2+‖γdivzN‖2≤(ε+c∞2c0CF2)‖∇zN‖2+CγαB02‖qN‖2≤C~4α‖qN‖2. Taking into account that $$\sup\nolimits_{0 <t <\infty}\frac{4}{1+t}\sqrt{\frac{t}{4}}=1$$ we get   |||(wN,rN)|||≤2|||(vN,qN)|||. Hence, with   B(N,ε,qN):=12(1+C~)=B0(N,ε,qN)22B0(N,ε,qN)2+16+8Cγ, the stated condition holds. □ 3.1 Interpolation error Assumption 3.5 There exists an interpolation operator $$I=(I_1,I_2):(C({\it{\Omega}}))^2\to V_N$$ such that the anisotropic interpolation error bounds   ‖wm−Imwm‖Lp(τij) ≤C(hik+1‖∂xk+1wm‖Lp(τij)+ℓjk+1‖∂yk+1wm‖Lp(τij)),m=1,2, (3.12)  ‖(wm−Imwm)x‖Lp(τij) ≤C(hik‖∂xk+1wm‖Lp(τij)+ℓjk‖∂xk∂ywm‖Lp(τij)),m=1,2, (3.13) and similarly for the $$y$$-derivative, hold true for $$p\in [1,\infty]$$ and arbitrary $${\boldsymbol{w}}=(w_1,w_2)\in (W^{k+1,p}({\it{\Omega}}))^2$$. Furthermore, we have the stability estimates   ‖Imwm‖∞,τ≤C‖wm‖∞,τ,m=1,2, (3.14) and   ‖∂xImwm‖∞,τ≤C‖∂xwm‖∞,τand‖∂yImwm‖∞,τ≤C‖∂ywm‖∞,τ,m=1,2. (3.15) Assumption 3.6 There exists an interpolation operator $$J:C({\it{\Omega}})\to Q_N$$ such that   ‖q−Jq‖∞≤CE(N,ε)k hold true for arbitrary $$q\in H^{k}({\it{\Omega}})$$ where $${E(N,{\varepsilon})}$$ is given by (2.8). Using these assumptions on the existence of interpolation operators, we estimate the interpolation errors and their derivatives in the next two lemmas. Lemma 3.7 Let $$\sigma>k+1$$, $$r=v_1+w_1+z_1$$ and $$s=v_2+w_2+z_2$$ from Assumption 2.1, $$I=(I_1,I_2)$$ the interpolation operator from Assumption 3.5, $$J$$ the interpolation operator from Assumption 3.6 and $${E(N,{\varepsilon})}$$ be given by (2.8). Then, the interpolation errors in the $$L^2$$-norm can be bounded by   ‖r−I1r‖Ω11+‖s−I2s‖Ω11 ≤CN−σ(ε1/2+N−1/2), (3.16)  ‖r−I1r‖Ω∖Ω11+‖s−I2s‖Ω∖Ω11 ≤Cε1/2E(N,ε)k+1. (3.17) Moreover, the estimates   ‖g1−I1g1‖Ω11+‖g2−I2g2‖Ω11 ≤CN−(k+1),‖g1−I1g1‖Ω∖Ω11+‖g2−I2g2‖Ω∖Ω11 ≤CE(N,ε)k+1 hold true for the smooth parts $$g_1$$ and $$g_2$$. Proof. Since the proof of the error bounds for $$s$$ follows the same line as the proof for $$r$$, we restrict ourselves to the latter one. We start with $${\it{\Omega}}_{11}$$. Clearly, we have   ‖r−I1r‖Ω11≤‖r‖Ω11+‖I1r‖Ω11. (3.18) With Assumption 2.1, we get for $$\|r\|_{{\it{\Omega}}_{11}}$$ directly,   ‖r‖Ω11 ≤‖v1‖Ω11+‖w1‖Ω11+‖z1‖Ω11≤Cε1/2N−σ, (3.19) where we used integration and the definition of $$\lambda_x$$ and $$\lambda_y$$. Let $$\overline{h}$$ and $$\overline{\ell}$$ denote the mesh width and height in $${\it{\Omega}}_{11}$$. In order to estimate $$I_1 r,$$ we use an idea from Stynes & Tobiska (2003) and divide the domain $${\it{\Omega}}_{11}$$ into $$S=[\lambda_x+\overline{h},1]\times [\lambda_y+\overline{\ell},1]$$ and $${\it{\Omega}}_{11}\setminus S$$ with the mesh transition points $$\lambda_x$$ and $$\lambda_y$$. Note that $${\it{\Omega}}_{11}\setminus S$$ consists of only a single ply of $$N-1$$ mesh elements and $$\overline{h},\overline{\ell}\leq CN^{-1}$$. Thus, we obtain   ‖I1r‖Ω11∖S2≤∑τ⊂Ω11∖Sh¯ℓ¯‖I1r‖∞,τ2≤CN−1‖r‖∞,Ω112≤CN−2σ−1 (3.20) using the stability (3.14) and in $$S$$,   ‖I1r‖S2≤C{∫Ω11(ε2e−2β1x/ε+e−2β2y/ε(1+ε2e−2β1x/ε))}≤CεN−2σ. (3.21) Combining (3.18)–(3.21) completes the proof for (3.16). On $${\it{\Omega}}_{12}$$, we get by using (3.12), (2.6) and (2.7),   ‖v1−I1v1‖Ω122 ≤C(‖hik+1∂xk+1v1‖Ω122+‖ℓjk+1∂yk+1v1‖Ω122)≤Cε3(N−1max|ψ′|)2(k+1). Using the stability (3.14), we obtain for $$w_1$$ and $$z_1$$,   ‖(w1+z1)−I1(w1+z1)‖Ω12≤C(measΩ12)1/2‖w1+z1‖∞,Ω12≤Cε1/2N−σ(ln⁡N)1/2. Similarly we have on $${\it{\Omega}}_{21}$$,   ‖w1−I1w1‖Ω212≤Cε(N−1max|ψ′|)2(k+1) and   ‖(v1+z1)−I1(v1+z1)‖Ω21≤Cε3/2N−σ(ln⁡N)1/2. Using the anisotropic interpolation error bounds (3.12), we get on $${\it{\Omega}}_{22}$$ analogously,   ‖v1−I1v1‖Ω222 ≤Cε4ln⁡N(N−1max|ψ′|+ℓ)2(k+1),‖w1−I1w1‖Ω222 ≤Cε2ln⁡N(h+N−1max|ψ′|)2(k+1),‖z1−I1z1‖Ω222 ≤Cε4(N−1max|ψ′|)2(k+1). Together with $$\sigma>k+1$$ and the definition of $${E(N,{\varepsilon})}$$ we get (3.17). For the smooth part, we present only the proofs for $$g_1$$ since $$g_2$$ can be estimated similarly. We have   ‖g1−I1g1‖Ω11 ≤C(measΩ11)1/2‖g1−I1g1‖∞,Ω11≤C(h¯k+1+ℓ¯k+1)≤CN−(k+1),‖g1−I1g1‖Ω∖Ω11 ≤C(measΩ∖Ω11)1/2‖g1−I1g1‖∞,Ω∖Ω11≤C(h+ℓ+N−1)k+1, where (3.12) in the $$L^\infty$$-norm was used. □ Lemma 3.8 Let $$\sigma>k+1$$ and $${E(N,{\varepsilon})}$$ be given by (2.8). Then the interpolation errors in the $$H^1$$-seminorm for $${\boldsymbol{u}}=(u_1,u_2)$$ and $$I=(I_1,I_2)$$ from Assumption 3.5 can be bounded by   ‖∇(u−Iu)‖Ω11 ≤Cε−1/2N−k,‖∇(u−Iu)‖Ω∖Ω11 ≤Cε−1/2E(N,ε)k while the estimates   ‖div(u−Iu)‖Ω11 ≤CN−k,‖div(u−Iu)‖Ω∖Ω11 ≤CE(N,ε)k for the divergence of the interpolation error hold true. Proof. Since the proofs of the estimates in the $$H^1$$-seminorm are similar for $$u_1$$ and $$u_2$$, we show only the proofs for $$u_1$$. With the anisotropic interpolation error bounds (3.13) and its $$y$$-derivative counterpart we get for the smooth part,   ‖∇(g1−I1g1)‖Ω11 ≤C(‖(g1−I1g1)x‖∞,Ω11+‖(g1−I1g1)y‖∞,Ω11)≤CN−k and similarly   ‖∇(g1−I1g1)‖Ω∖Ω11 ≤C(h+ℓ+N−1)k. Now, we estimate the layer parts. Consider first the term $$v_1$$. Here we have by the $$W^{1,\infty}$$-stability (3.15),   ‖(v1−I1v1)x‖Ω11∪Ω21 ≤C(‖v1,x‖Ω11∪Ω21+(measΩ11∪Ω21)1/2‖(I1v1)x‖∞,Ω11∪Ω21) ≤C(‖v1,x‖Ω11∪Ω21+‖v1,x‖∞,Ω11∪Ω21) ≤CN−σ. In the remaining regions, we use (3.13) of $$I_1$$ to deduce   ‖(v1−I1v1)x‖Ω12∪Ω22 ≤C{‖hik∂xk+1v1‖Ω12∪Ω22+‖ℓjk∂x∂ykv1‖Ω12∪Ω22} ≤Cε1/2(ℓ+N−1max|ψ′|)k. The techniques for the $$y$$-derivatives are similar. Hence, we obtain   ‖(v1−I1v1)y‖Ω11∪Ω21 ≤C(‖v1,y‖Ω11∪Ω21+‖I1v1,y‖∞,Ω11∪Ω21) ≤C(‖v1,y‖Ω11∪Ω21+‖v1,y‖∞,Ω11∪Ω21) ≤CεN−σ and   ‖(v1−I1v1)y‖Ω12∪Ω22 ≤C{‖hik∂xk∂yv1‖Ω12∪Ω22+‖ℓjk∂yk+1v1‖Ω12∪Ω22} ≤Cε3/2(ℓ+N−1max|ψ′|)k. For the second layer part $$w_1$$, we get analogously   ‖(w1−I1w1)x‖Ω11∪Ω12 ≤CN−σ,‖(w1−I1w1)x‖Ω21Ω22 ≤Cε1/2(h+N−1max|ψ′|)k,‖(w1−I1w1)y‖Ω11∪Ω12 ≤C(ε−1/2+N)N−σ,‖(w1−I1w1)y‖Ω21∪Ω22 ≤Cε−1/2(h+N−1max|ψ′|)k. Now, we estimate the error of the function $$z_1$$ using the $$W^{1,\infty}$$-stability (3.15) of $$I_1$$,   ‖(z1−I1z1)x‖Ω∖Ω22 ≤C(‖z1,x‖Ω∖Ω22+‖z1,x‖∞,Ω∖Ω22)≤CN−σ, and in $${\it{\Omega}}_{22}$$ with (3.13),   ‖(z1−I1z1)x‖Ω22 ≤C{‖hik∂xk+1z1‖Ω22+‖ℓjk∂x∂ykz1‖Ω22}≤Cε(N−1max|ψ′|)k. The bounds for the $$y$$-derivatives follow analogously:   ‖(z1−I1z1)y‖Ω∖Ω22 ≤CN−σ,‖(z1−I1z1)y‖Ω22 ≤Cε(N−1max|ψ′|)k. Putting all these estimates together, we obtain with the definition of $${E(N,{\varepsilon})}$$ the bounds for the $$H^1$$-seminorms. We remark that the bounds of the $$x$$-derivatives of $$u_1$$ are better than those of the $$y$$-derivatives. For $$u_2$$ and $$I_2,$$ we obtain similar results with the difference that the $$x$$-derivatives are worse than the $$y$$-derivatives. Now combining the $$x$$-derivatives of $$u_1$$ and the $$y$$-derivatives of $$u_2$$, we get   ‖div(u−Iu)‖Ω11 ≤‖(u1−I1u1)x‖Ω11+‖(u2−I2u2)y‖Ω11≤CN−k,‖div(u−Iu)‖Ω∖Ω11 ≤C(εln⁡N)1/2(h+ℓ+N−1max|ψ′|)k≤CE(N,ε)k, which are the last two statements of this lemma. □ Lemma 3.9 Let $$\sigma>k+1$$, $$I=(I_1,I_2),$$ the interpolation operator from Assumption 3.5 and $$J$$ the interpolation operator from Assumption 3.6. Then the estimate   |||(u−Iu,p−Jp)|||≤CE(N,ε)k (3.22) holds true where $${E(N,{\varepsilon})}$$ is given by (2.8). Proof. For this proof, we use the previous interpolation error estimates. The definition of the energy norm gives   |||(u−Iu,p−Jp)|||2=ε‖∇(u−Iu)|02+c0‖u−Iu‖2+α‖p−Jp‖2+‖γdiv(u−Iu)‖2. (3.23) For the first term on the right-hand side of (3.23), we get with Lemma (3.8),   ε1/2‖∇(u−Iu)‖ ≤ε1/2(‖∇(u1−I1u1)‖+‖∇(u2−I2u2)‖)≤CE(N,ε)k. The second term of (3.23) yields with Lemma 3.7,   ‖u−Iu‖ ≤C(‖u1−I1u1‖+‖u2−I2u2‖)≤CE(N,ε)k+1. The bound   ‖p−Jp‖≤CE(N,ε)k for the pressure term follows from Assumption 3.6. The stabilization term can be bounded by the divergence bounds of Lemma (3.8),   ‖γdiv(u−Iu)‖≤γ∞‖div(u−Iu)‖≤CE(N,ε)k. Summarizing all these estimates finishes the proof. □ 3.2 A priori error estimates Lemma 3.10 Let $$\gamma\geq \gamma_0>0$$ be independent of $$N$$ and $${\varepsilon}$$, and $$\sigma>k+1$$. The solutions of 2.1 and 3.3 are denoted by $$({\boldsymbol{u}},p)$$ and $$({\boldsymbol{u}}_N,p_N)$$, respectively. Then, the estimate   |||(Iu−uN,Jp−pN)|||≤C1B(N,ε,Jp−pN)E(N,ε)k holds true where $$C\to\infty$$ for $$\gamma_0\to 0$$. Proof Abbreviating $$B=B(N,{\varepsilon},J p-p_N)$$ in this proof, the estimate 3.8 gives   |||(Iu−uN,Jp−pN)||| ≤1Bsup(wN,rN)∈VN×QNA((Iu−uN,Jp−pN),(wN,rN))|||(wN,rN)||| =1Bsup(wN,rN)∈VN×QNA((Iu−u,Jp−p),(wN,rN))|||(wN,rN)||| due to the Galerkin orthogonality of Lemma 3.2. The expression   A((Iu−u,Jp−p),(wN,rN))  =ε(∇(Iu−u),∇wN)−(b⋅∇(Iu−u),wN)+(c(Iu−u),wN)  −(Jp−p,divwN)+(rN,div(Iu−u))+(γdiv(Iu−u),γdiv(wN)) will be estimated term by term. We get with the interpolation error bounds of Lemmas 3.7 and 3.8    |ε(∇(Iu−u),∇wN)+(c(Iu−u),wN)|  ≤(ε‖∇(Iu−u)‖2+c∞2c0‖Iu−u‖2)1/2(ε‖∇wN‖2+c0‖wN‖2)1/2  ≤CE(N,ε)k|||(wN,rN)|||, by Lemma 3.7 for the pressure,   |(Jp−p,divwN)| =|(Jp−p,divwN)| ≤‖Jp−p‖γ0−1|||(wN,rN)||| ≤Cγ0−1E(N,ε)k|||(wN,rN)||| and with the divergence bounds of Lemma 3.8,   |(rN,div(Iu−u))+(γdiv(Iu−u),γdivwN))| ≤‖div(Iu−u)‖(‖rN‖+γ∞‖γdivwN‖) ≤(α−1/2+γ∞)‖div(Iu−u)‖|||(wN,rN)||| ≤C(1+γ∞)E(N,ε)k|||(wN,rN)|||. The last two estimates show that both $$\gamma_0^{-1}$$ and $$\gamma_{\infty}$$ should be bounded independent of $$N$$ and $${\varepsilon}$$. Hence, $$\gamma$$ should be contained in the interval $$[\gamma_0,\gamma_{\infty}]$$ with end points independent of $$N$$ and $${\varepsilon}$$. For the convective term, we estimate the terms of   |(b⋅∇(Iu−u),wN)| |(b1(I1u1−u1)x,wN,1)|+|(b2(I1u1−u1)y,wN,1)|+ |(b1(I2u2−u2)x,wN,2)|+|(b2(I2u2−u2)y,wN,2)| individually. Here we have directly,   |(b1(I1u1−u1)x,wN,1)| ≤C‖(I1u1−u1)x‖‖wN,1‖≤CE(N,ε)k|||(wN,rN)|||,|(b2(I2u2−u2)y,wN,2)| ≤C‖(I2u2−u2)y‖‖wN,2‖≤CE(N,ε)k|||(wN,rN)|||. For the other terms, we use integration by parts to obtain   |(b2(I1u1−u1)y,wN,1)| ≤C(‖I1u1−u1‖‖wN,1‖+|(b2(I1u1−u1),wN,1,y)|) ≤C(‖I1u1−u1‖‖wN,1‖+‖I1u1−u1‖Ω11‖wN,1,y‖Ω11 +‖I1u1−u1‖Ω∖Ω11‖wN,1,y‖Ω∖Ω11) ≤C(‖u1−u1‖+N‖I1u1−u1‖Ω11+ε−1/2‖I1u1−u1‖Ω∖Ω11)|||(wN,rN)||| ≤CE(N,ε)k|||(wN,rN)|||. The remaining term can be estimated analogously. Upon summarizing the proof is done. □ Theorem 3.11 Let $$\gamma\geq\gamma_0>0$$ be independent of $$N$$ and $${\varepsilon}$$, and $$\sigma>k+1$$. The error between the solution $$({\boldsymbol{u}},p)$$ of 2.1 and the discrete solution $$({\boldsymbol{u}}_N,p_N)$$ of 3.3 can be bounded by   |||(u−uN,p−pN)|||≤C(1+1B(N,ε,Jp−pN))E(N,ε)k, where $${E(N,{\varepsilon})}$$ is given by 2.8 and $$J$$ is an interpolation operator fulfilling Assumption 3.6. Proof It follows directly from Lemmas 3.9 and 3.10, and the triangle inequality. □ Corollary 3.12 If $$\gamma=0$$ is chosen on the whole domain, then we have the standard Galerkin method. For the solution of this method we have the estimate   |||(u−uN,p−pN)|||≤C(1+1B(N,ε,Jp−pN))(N+(ln⁡N)1/2)E(N,ε)k, i.e., the convergence order reduces by 1 compared to the stabilized case. Proof The only difference to the estimation before is the bounding of the term $$|(J p-p,{\text{div}} {\boldsymbol{w}}_N)|$$. This term cannot be estimated against $$\|\gamma{\text{div}} {\boldsymbol{w}}_N\|$$ due to $$\gamma=0$$. Alternatively, we get   |(Jp−p,divwN)| =|(Jp−p,divwN)| ≤C(‖Jp−p‖Ω11‖∇wN‖Ω11+‖Jp−p‖Ω∖Ω11‖∇wN‖Ω∖Ω11) ≤C(N‖Jp−p‖Ω11‖wN‖Ω11+(εln⁡N)1/2‖Jp−p‖∞,Ω∖Ω11‖∇wN‖Ω∖Ω11) ≤C(N‖Jp−p‖+(ln⁡N)1/2‖Jp−p‖∞)|||(wN,rN)||| ≤C(N+(ln⁡N)1/2)E(N,ε)k|||(wN,rN)||| and the reduced order can be seen. □ Corollary 3.13 If we choose $$\gamma=0$$ only in the layer region, we can recover the convergence rate $$k$$ compared with $$\gamma=0$$ on the whole domain. Proof See the proof before. In the layer region, we obtain an additional factor of only $$(\ln N)^{1/2}$$ and not $$N$$. □ 4 Examples for suitable spaces We will consider two families of finite element pairs on the introduced meshes. The first one is the Taylor–Hood family $$Q_k\times Q_{k-1}$$, $$k\ge 2$$, and the second one is $$Q_k\times P_{k-1}^{\text{disc}}$$, $$k\ge 2$$. On isotropic meshes, the pairs $$Q_k\times Q_{k-1}$$, $$k\ge 2$$, and $$Q_k\times P_{k-1}^{\text{disc}}$$, $$k\ge 2$$, are known to fulfil Assumption 3.1 even with a constant $$\beta_0>0$$ that is independent of $$N$$; see Girault & Raviart (1986, Ch. II), and Assumption 3.1 becomes just the well-known inf–sup or Babuška–Brezzi condition. Since the $${\varepsilon}$$-dependence of $$\beta_0(N,{\varepsilon})$$ comes from the underlying mesh only, there is also no $${\varepsilon}$$-dependence on general isotropic meshes. On anisotropic meshes, the situation is much more complicated. The existence of an inf–sup constant that is independent of the maximal aspect ratio of the underlying mesh depends strongly on the precise mesh structure and the chosen spaces for approximating velocity and pressure. Proofs for the existence of inf–sup constants are mainly available for situations where the difference between the polynomial orders of velocity and pressure is $$2$$; see Sch¨otzau & Schwab (1998), Sch¨otzau et al. (1999) for the cases $$Q_k\times Q_{k-2}$$ on two-dimensional domains. A review of results on inf–sup constants for the Stokes problem on several types of anisotropic meshes can be found in Apel & Maharavo Randrianarivony (2003). Nonconforming discretizations of the Stokes problem on anisotropic rectangular meshes were investigated in Apel & Matthies (2008). Although the four nonconforming velocity approximation spaces in Apel & Matthies (2008) differ only slightly, just two of them provide an inf–sup constant independent of the aspect ratio. Another approach for generating inf–sup-stable pairs by eliminating certain pressure modes is presented in Ainsworth et al. (2015). In Section 5.1.1, we investigate the values of $$\beta_0=\beta_0(N,{\varepsilon})$$ and $$B_0=B_0(N,{\varepsilon},q_N)$$ for both element pairs. It can be observed that for the continuous pressure space, $$\beta_0$$ is bounded independent of $${\varepsilon}$$ but not independent of $$N$$, while for the discontinuous pressure space it goes to zero for $${\varepsilon}\to 0$$. In contrast to this negative observation, the values of $$B_0$$ stay bounded away from zero, independent of $$N$$ and $${\varepsilon}$$, for both cases if $$q_N$$ is the difference between the interpolation of $$p$$ and the discrete solution $$p_N$$. This shows that $$B_0$$ might be more appropriate to estimate the errors although we do not have a theoretical insight into the boundedness of $$B_0$$. In order to define the interpolation operators, let $$I_k$$ be the scalar nodal interpolator of order $$k$$. Then we set $$I=(I_k,I_k)$$. In Apel (1999, Theorem 2.7) the anisotropic interpolation error bounds for Assumption 3.5 are given. The stability estimates are standard and can be shown directly. For the interpolation operator into the discrete pressure space $$Q_N$$, we have to distinguish two cases. In the case $$Q_N=P_{k-1}^{\text{disc}}$$, we can use the $$L^2$$-projection into $$Q_N$$ as $$J$$. The global $$L^2$$-projection furthermore localizes to $$L^2$$-projections on each mesh cell. Hence, the proof of Assumption 3.6 can exploit the anisotropic error estimates for the $$L^2$$-projection that follow the same reasoning as in Apel (1999). In the case $$Q_N=Q_{k-1}$$, the construction of $$J$$ is two stage. Let $$I_{k-1}$$ denote the scalar nodal interpolation into the space of continuous, piecewise $$Q_{k-1}$$ functions. Then we set   Jq=Ik−1q−c¯,wherec¯=1|Ω|∫ΩIk−1q to ensure that $$J$$ maps into $$L^2_0({\it{\Omega}})$$. Using again the anisotropic interpolation error estimates by Apel (1999), we have   ‖p−Ik−1p‖Ω11≤CN−kand‖p−Ik−1p‖Ω∖Ω11≤CE(N,ε)k. From the definition of $$J$$ and the fact that constants are $$L^2({\it{\Omega}})$$-orthogonal to functions from $$L^2_0({\it{\Omega}})$$, we get   ‖p−Ik−1p‖2=‖p−Jp−c¯‖2=‖p−Jp‖2+‖c¯‖2 and therefore   ‖p−Jp‖≤‖p−Ik−1p‖≤CE(N,ε)k. Furthermore, we obtain in $$L^\infty$$,   ‖p−Jp‖∞ ≤‖p−Ik−1p‖∞+|c¯|≤‖p−Ik−1p‖∞+|∫ΩIk−1p−p|≤CE(N,ε)k which provides Assumption 3.6. Corollary 4.1 For both choices of the discrete pressure space, Corollary 3.12 can be improved slightly. If $$Q_N=P_{k-1}^{\rm disc}$$ then $$\|J p-p\|_{{\it{\Omega}}_{11}}\leq CN^{-k}$$ and the final lines in the proof of Corollary 3.12 become   |(Jp−p,divwN)| ≤C(N‖Jp−p‖Ω11+(ln⁡N)1/2‖Jp−p‖∞,Ω∖Ω11)|||(wN,rN)||| ≤C(N−k+1+(ln⁡N)1/2E(N,ε)k)|||(wN,rN)||| ≤CE(N,ε)k−1|||(wN,rN)|||. If $$Q_N=Q_{k-1}$$ then $$(J p-p,{\text{div}} {\boldsymbol{w}}_N)=(I_{k-1}p-p,{\text{div}}{\boldsymbol{w}}_N)$$ and the proof of Corollary 3.12 can be done with $$I_{k-1}$$ instead of $$J$$. Using $$\|I_{k-1} p-p\|_{{\it{\Omega}}_{11}}\leq CN^{-k}$$ leads to   |(Jp−p,divwN)| =|(Ik−1p−p,divwN)| ≤C(N−k+1+(ln⁡N)1/2E(N,ε)k)|||(wN,rN)||| ≤CE(N,ε)k−1|||(wN,rN)|||. 5 Numerical results In the numerical section, we show the results of computations supporting our theoretical findings. We start with a problem where the exact solution is known and end with some considerations on problems with unknown exact solutions where our theory does not hold to its full extent. 5.1 Problem with known solution We consider the singularly perturbed Oseen equations   −εΔu−((2+3x3+2y2)⋅∇)u+(1+2xy)u+∇p =f in Ω=(0,1)2,divu =0 in Ω,u =0 on ∂Ω, (5.1) where the right-hand side $${\boldsymbol{f}}$$ of (5.1) was chosen such that   u1(x,y)=∂Ψ(x,y)∂y,u2(x,y)=−∂Ψ(x,y)∂x,p(x,y)=2cos⁡(x)sin⁡(y)−2sin⁡(1)(1−cos⁡(1)) with $$\Psi(x,y):=F(x)G(y)$$ and   F(x) :=c3(1−x)3+c2(1−x)2+c1(1−x)−ε(e−2x/ε−e−2/ε)−sin⁡(1−x),G(y) :=d3(1−y)3+d2(1−y)2+d1(1−y)−1+ε(e−3y/ε−e−3/ε)+cos⁡(1−y) and constants   c1 :=(1+2e−2/ε),c2:=−((4+3ε)e−2/ε+4+cos⁡(1)−3sin⁡(1)−3ε),c3 :=((2+2ε)e−2/ε+3+cos⁡(1)−2ε−2sin⁡(1)),d1 :=−3e−3/ε,d2:=((6+3ε)e−3/ε+6−3cos⁡(1)−sin⁡(1)−3ε),d3 :=−((3+2ε)e−3/ε+5−2cos⁡(1)−sin⁡(1)−2ε) is the solution of (5.1). The function $$u_1$$ shows a strong exponential boundary layer at $$y=0$$ and a weak exponential boundary layer at $$x=0$$, and vice versa for $$u_2$$. The pressure shows no layer behaviour. Moreover, Assumption 2.1 is satisfied. For a visualization of the solution, see Fig. 1. All calculations were done in MATLAB using the finite element suite $$\mathbb{SOFE}$$ developed by Lars Ludwig; see Ludwig (2016). In the following, ‘order’ will always denote the exponent $$r$$ in a convergence order of form $$\mathcal{O}(N^{-r})$$ while ‘$$\ln$$-order’ corresponds to the exponent $$r$$ in a convergence order of form $$\mathcal{O}((N^{-1}\ln N)^r)$$. If not stated otherwise then the results on the two finest meshes were used to calculate the convergence order. We have chosen $$\gamma=1$$ on the whole domain, unless otherwise stated. The parameter $$\sigma$$ was chosen to be $$\sigma=k+2$$, where $$k$$ is the order of the elements for the velocity $${\boldsymbol{u}}$$. The errors will be measured in the norm   |||(u,p)|||ε2:=ε‖∇u‖2+‖u‖2+‖p‖2+‖γdivu‖2; (5.2) this means that $$\alpha$$ is set to $$1$$. 5.1.1 Dependence of $${\beta_0}$$ and $${B_0}$$ on $${N}$$ and $${{\varepsilon}}$$ First we will consider the dependence of $$\beta_0=\beta_0(N,{\varepsilon})$$ and $$B_0=B_0(N,{\varepsilon},Jp-p_N)$$ on $$N$$ and $${\varepsilon}$$. The pictures in Fig. 3 show the dependence of $$\beta_0$$ on $$N$$ and $${\varepsilon}$$ for the discretizations $$Q_2\times Q_1$$ and $$Q_3\times Q_2$$ on Shishkin meshes. Clearly, there is an $${\varepsilon}$$-uniform bound for any fixed $$N$$ while the value of $$\beta_0$$ will decrease with increasing $$N$$ if $${\varepsilon}$$ is fixed. If $$N$$ and $${\varepsilon}$$ are fixed, the constant $$\beta_0$$ decreases with increasing approximation order $$k$$. This effect appears also on isotropic meshes and spectral methods; see Bernardi & Maday (1997). The behaviour of $$\beta_0$$ on Bakhvalov–Shishkin meshes is presented in Fig. 4. We observe that $$\beta_0$$ takes much larger values and the dependence on $$N$$ is much weaker compared to Shishkin meshes. Numerical experiments for discretizations $$Q_k\times P_{k-1}^{\text{disc}}$$ show that $$\beta_0$$ approaches zero if $${\varepsilon}$$ tends to zero. Fig. 3. View largeDownload slide Problem (5.1): The inf–sup constants for different values of $$N$$ as a function of $$\varepsilon$$ on Shishkin meshes: $$Q_2\times Q_1$$ (left) and $$Q_3\times Q_2$$ (right). Fig. 3. View largeDownload slide Problem (5.1): The inf–sup constants for different values of $$N$$ as a function of $$\varepsilon$$ on Shishkin meshes: $$Q_2\times Q_1$$ (left) and $$Q_3\times Q_2$$ (right). Fig. 4. View largeDownload slide Problem (5.1): The inf–sup constants for different values of $$N$$ as a function of $$\varepsilon$$ on Bakhvalov–Shishkin meshes: $$Q_2\times Q_1$$ (left) and $$Q_3\times Q_2$$ (right). Fig. 4. View largeDownload slide Problem (5.1): The inf–sup constants for different values of $$N$$ as a function of $$\varepsilon$$ on Bakhvalov–Shishkin meshes: $$Q_2\times Q_1$$ (left) and $$Q_3\times Q_2$$ (right). Figure 5 shows the dependence of $$\beta_0$$ and $$B_0$$ on $$N$$ for $${\varepsilon}=10^{-8}$$. We observe that the values of $$B_0$$ are much larger than the corresponding values of $$\beta_0$$. Moreover, the value of $$B_0$$ is bounded away from zero for both continuous and discontinuous pressure spaces and is nondecreasing in $$N$$. Additional studies showed that there is also an $${\varepsilon}$$-uniform lower bound for $$B_0$$. Hence, Theorem 3.11 provides an $${\varepsilon}$$-uniform error bound. Fig. 5. View largeDownload slide Problem (5.1): The behaviour of $$\beta_0$$ and $$B_0$$ with $${\varepsilon}=10^{-8}$$ as a function of $$N$$ for $$Q_3\times Q_2$$ (left) and $$Q_3\times\,P_2^{\text{disc}}$$ (right), respectively. Fig. 5. View largeDownload slide Problem (5.1): The behaviour of $$\beta_0$$ and $$B_0$$ with $${\varepsilon}=10^{-8}$$ as a function of $$N$$ for $$Q_3\times Q_2$$ (left) and $$Q_3\times\,P_2^{\text{disc}}$$ (right), respectively. 5.1.2 Uniformity in $${{\varepsilon}}$$ For this study, we have chosen $$\gamma=1$$ on the whole domain in the computations. Table 2 presents the numerical results on Shishkin meshes (S-mesh) and Bakhvalov–Shishkin meshes (B–S-mesh) with $$N=32$$ and $${\varepsilon}$$ varying from $$10^{-1}$$ to $$10^{-10}$$. We show second-order approximations by $$Q_2\times Q_1$$ and $$Q_2\times P_1^{\text{disc}}$$ as well as third-order approximations by $$Q_3\times Q_2$$ and $$Q_3\times P_2^{\text{disc}}$$. Note that for large values of $${\varepsilon}$$ the mesh is equidistant while the characteristic meshes are obtained for smaller values of $${\varepsilon}$$. We clearly see from Table 2 that the error for small $${\varepsilon}$$ is independent of $${\varepsilon}$$, as predicted by our theory and the $${\varepsilon}$$-uniform bound for $$B_0$$. Furthermore, the errors for discretizations of the same order and on the same mesh family differ only slightly. Moreover, the errors on Bakhvalov–Shishkin meshes are smaller than those on Shishkin meshes. Table 2 Problem (5.1): Error in energy norm $${\left|\!\!\;\left|\!\!\;\left| {(\cdot,\cdot)} \right|\!\!\;\right|\!\!\;\right|}_{\varepsilon}$$ for $$N=32$$    $$Q_2\times Q_1$$  $$Q_2\times P^{\text{disc}}_1$$  $$Q_3\times Q_2$$  $$Q_3\times P^{\text{disc}}_2$$  $${\varepsilon}$$  S-mesh  B–S-mesh  S-mesh  B–S-mesh  S-mesh  B–S-mesh  S-mesh  B–S-mesh  $$10^{-1}$$  1.7439-02  8.6659-03  1.7439-02  8.6661-03  1.4779-03  1.4779-03  1.4779-03  1.4779-03  $$10^{-2}$$  2.4899-02  3.5717-03  2.4901-02  3.5810-03  3.4299-03  1.2697-04  3.4299-03  1.2702-04  $$10^{-3}$$  2.5197-02  3.7450-03  2.5198-02  3.7565-03  3.4685-03  1.2840-04  3.4685-03  1.2850-04  $$10^{-4}$$  2.5230-02  3.8047-03  2.5233-02  3.8166-03  3.4723-03  1.2854-04  3.4723-03  1.2865-04  $$10^{-5}$$  2.5234-02  3.8146-03  2.5237-02  3.8273-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04  $$10^{-6}$$  2.5235-02  3.8156-03  2.5237-02  3.8292-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04  $$10^{-7}$$  2.5235-02  3.8157-03  2.5237-02  3.8282-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04  $$10^{-8}$$  2.5235-02  3.8157-03  2.5237-02  3.8272-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04  $$10^{-9}$$  2.5235-02  3.8157-03  2.5237-02  3.8271-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04  $$10^{-10}$$  2.5235-02  3.8157-03  2.5237-02  3.8271-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04     $$Q_2\times Q_1$$  $$Q_2\times P^{\text{disc}}_1$$  $$Q_3\times Q_2$$  $$Q_3\times P^{\text{disc}}_2$$  $${\varepsilon}$$  S-mesh  B–S-mesh  S-mesh  B–S-mesh  S-mesh  B–S-mesh  S-mesh  B–S-mesh  $$10^{-1}$$  1.7439-02  8.6659-03  1.7439-02  8.6661-03  1.4779-03  1.4779-03  1.4779-03  1.4779-03  $$10^{-2}$$  2.4899-02  3.5717-03  2.4901-02  3.5810-03  3.4299-03  1.2697-04  3.4299-03  1.2702-04  $$10^{-3}$$  2.5197-02  3.7450-03  2.5198-02  3.7565-03  3.4685-03  1.2840-04  3.4685-03  1.2850-04  $$10^{-4}$$  2.5230-02  3.8047-03  2.5233-02  3.8166-03  3.4723-03  1.2854-04  3.4723-03  1.2865-04  $$10^{-5}$$  2.5234-02  3.8146-03  2.5237-02  3.8273-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04  $$10^{-6}$$  2.5235-02  3.8156-03  2.5237-02  3.8292-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04  $$10^{-7}$$  2.5235-02  3.8157-03  2.5237-02  3.8282-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04  $$10^{-8}$$  2.5235-02  3.8157-03  2.5237-02  3.8272-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04  $$10^{-9}$$  2.5235-02  3.8157-03  2.5237-02  3.8271-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04  $$10^{-10}$$  2.5235-02  3.8157-03  2.5237-02  3.8271-03  3.4727-03  1.2856-04  3.4727-03  1.2866-04  5.1.3 Convergence studies Table 3 shows for Shishkin meshes and $${\varepsilon}=10^{-8}$$, which is sufficiently small to have a solution with layers, the errors in the energy norm for pairs of second-, third-, and fourth-order with continuous and discontinuous pressure approximations. It can be seen that the theoretically predicted convergence orders are achieved. The differences in the energy norm between the two different pressure approximation spaces are again very small. Here the different pressure errors are similar and in both cases dominated by the velocity errors. Table 3 Problem (5.1): Error in energy norm $${\left|\!\!\;\left|\!\!\;\left| {(\cdot,\cdot)} \right|\!\!\;\right|\!\!\;\right|}_{\varepsilon}$$ on Shishkin meshes with $${\varepsilon}=10^{-8}$$ $$N$$  $$Q_2\times Q_1$$  $$Q_2\times P_1^{\text{disc}}$$  $$Q_3\times Q_2$$  $$Q_3\times P_2^{\text{disc}}$$  $$Q_4\times Q_3$$  $$Q_4\times P_3^{\text{disc}}$$  4  2.7153-01  2.5880-01  7.2169-02  7.2220-02  2.5802-02  2.5803-02  8  1.3134-01  1.3131-01  3.6858-02  3.6859-02  1.0660-02  1.0660-02  16  6.1849-02  6.1859-02  1.2918-02  1.2918-02  2.6828-03  2.6828-03  32  2.5235-02  2.5237-02  3.4727-03  3.4727-03  4.7297-04  4.7297-04  64  9.2647-03  9.2650-03  7.8324-04  7.8324-04  6.5527-05  6.5527-05  $$\ln$$-order  1.96  1.96  2.92  2.92  3.87  3.87  Theory  2  2  3  3  4  4  $$N$$  $$Q_2\times Q_1$$  $$Q_2\times P_1^{\text{disc}}$$  $$Q_3\times Q_2$$  $$Q_3\times P_2^{\text{disc}}$$  $$Q_4\times Q_3$$  $$Q_4\times P_3^{\text{disc}}$$  4  2.7153-01  2.5880-01  7.2169-02  7.2220-02  2.5802-02  2.5803-02  8  1.3134-01  1.3131-01  3.6858-02  3.6859-02  1.0660-02  1.0660-02  16  6.1849-02  6.1859-02  1.2918-02  1.2918-02  2.6828-03  2.6828-03  32  2.5235-02  2.5237-02  3.4727-03  3.4727-03  4.7297-04  4.7297-04  64  9.2647-03  9.2650-03  7.8324-04  7.8324-04  6.5527-05  6.5527-05  $$\ln$$-order  1.96  1.96  2.92  2.92  3.87  3.87  Theory  2  2  3  3  4  4  The results on Bakhvalov–Shishkin meshes for $${\varepsilon}=10^{-8}$$ are shown in Table 4. The same approximation spaces as on Shishkin meshes have been used. The convergence orders predicted by our theory are obtained. Moreover, one sees that the errors on Bakhvalov–Shishkin meshes are smaller than those on the corresponding Shishkin meshes. This is due to the logarithmic factor which is present only for Shishkin meshes. The difference between Shishkin and Bakhvalov–Shishkin meshes becomes larger with increasing approximation order and already reaches two orders of magnitude for the fourth-order pairs $$Q_4\times Q_3$$ and $$Q_4\times P_3^{\text{disc}}$$. Table 4 Problem (5.1): Error in energy norm $${\left|\!\!\;\left|\!\!\;\left| {(\cdot,\cdot)} \right|\!\!\;\right|\!\!\;\right|}_{\varepsilon}$$ on Bakhvalov–Shishkin meshes with $${\varepsilon}=10^{-8}$$ $$N$$  $$Q_2\times Q_1$$  $$Q_2\times P_1^{\text{disc}}$$  $$Q_3\times Q_2$$  $$Q_3\times P_2^{\text{disc}}$$  $$Q_4\times Q_3$$  $$Q_4\times P_3^{\text{disc}}$$  4  2.2054-01  2.0379-01  3.1063-02  3.1183-02  8.3515-03  8.3526-03  8  5.6584-02  5.6509-02  6.1054-03  6.1147-03  9.5490-04  9.5494-04  16  1.4800-02  1.4845-02  9.3378-04  9.3470-04  7.8247-05  7.8249-05  32  3.8157-03  3.8272-03  1.2856-04  1.2866-04  5.5671-06  5.5672-06  64  9.7042-04  9.7364-04  1.6851-05  1.6864-05  3.7078-07  3.7079-07  Order  1.98  1.97  2.93  2.93  3.91  3.91  Theory  2  2  3  3  4  4  $$N$$  $$Q_2\times Q_1$$  $$Q_2\times P_1^{\text{disc}}$$  $$Q_3\times Q_2$$  $$Q_3\times P_2^{\text{disc}}$$  $$Q_4\times Q_3$$  $$Q_4\times P_3^{\text{disc}}$$  4  2.2054-01  2.0379-01  3.1063-02  3.1183-02  8.3515-03  8.3526-03  8  5.6584-02  5.6509-02  6.1054-03  6.1147-03  9.5490-04  9.5494-04  16  1.4800-02  1.4845-02  9.3378-04  9.3470-04  7.8247-05  7.8249-05  32  3.8157-03  3.8272-03  1.2856-04  1.2866-04  5.5671-06  5.5672-06  64  9.7042-04  9.7364-04  1.6851-05  1.6864-05  3.7078-07  3.7079-07  Order  1.98  1.97  2.93  2.93  3.91  3.91  Theory  2  2  3  3  4  4  5.1.4 Necessity of stabilization The size of the stabilization parameter $$\gamma$$ is arbitrary and our analysis gives only the restriction that it should be constant with respect to $$N$$ and $${\varepsilon}$$. Figure 6 shows for the energy norm of the solution and the $$L^2$$-norm of the pressure alone the dependency on $$\gamma$$ for a Bakhvalov–Shishkin mesh with fixed $$N=64$$, $${\varepsilon}=10^{-8}$$, $$Q_3\times Q_2$$ and stabilization in the full domain. It can clearly be seen that a moderate value of $$\gamma$$ gives best results. Furthermore, we observe almost constant error values on a wide range of values for $$\gamma$$ that include $$\gamma=1$$. This is the reason for choosing this particular value in the numerical studies. Fig. 6. View largeDownload slide Problem (5.1): Dependency of errors on the stabilization parameter $$\gamma$$. Fig. 6. View largeDownload slide Problem (5.1): Dependency of errors on the stabilization parameter $$\gamma$$. Table 5 shows the influence of the grad-div stabilization for $${\varepsilon}=10^{-8}$$ and $$Q_3\times Q_2$$. We consider three choices for $$\gamma$$: first $$\gamma=1$$ in $${\it{\Omega}}$$, second $$\gamma=1$$ in $${\it{\Omega}}_{11}$$ only, and third $$\gamma=0$$ in $${\it{\Omega}}$$. The last choice gives the standard Galerkin method, thus no stabilization. The errors are measured in the corresponding energy norm $$\left|\!\!\;\left|\!\!\;\left| {(\cdot,\cdot)} \right|\!\!\;\right|\!\!\;\right|_{\varepsilon}$$. We observe that the errors for discrete problems without grad-div stabilization are larger than those for both stabilized discrete problems. Moreover, the convergence order for problems without stabilization is less than for problems with stabilization. This is in agreement with Corollary 3.12. Furthermore, we observe that the errors of both stabilized methods are almost indistinguishable, confirming Corollary 3.13. Table 5 Problem (5.1): Error in energy norm $${\left|\!\!\;\left|\!\!\;\left| {(\cdot,\cdot)} \right|\!\!\;\right|\!\!\;\right|}_{\varepsilon}$$ for $$Q_3\times Q_2$$ and several choices of $$\gamma$$ on Bakhvalov–Shishkin meshes with $${\varepsilon}=10^{-8}$$    $$\gamma=1$$ in $${\it{\Omega}}$$  $$\gamma=1$$ in $${\it{\Omega}}_{11}$$  $$\gamma=0$$ in $${\it{\Omega}}$$  $$N$$  Error  Order  Error  Order  Error  Order  4  3.1063-02     3.1075-02     3.1302-02     8  6.1054-03  2.35  6.1061-02  2.35  6.1727-03  2.34  16  9.3378-04  2.71  9.3381-04  2.71  9.6048-04  2.68  32  1.2856-04  2.86  1.2856-04  2.86  1.4026-04  2.78  64  1.6851-05  2.93  1.6851-05  2.93  2.1901-05  2.68  128  2.1567-06  2.97  2.1567-06  2.97  4.1048-06  2.42  196  6.0560-07  2.98  6.0560-07  2.98  1.6071-06  2.20     $$\gamma=1$$ in $${\it{\Omega}}$$  $$\gamma=1$$ in $${\it{\Omega}}_{11}$$  $$\gamma=0$$ in $${\it{\Omega}}$$  $$N$$  Error  Order  Error  Order  Error  Order  4  3.1063-02     3.1075-02     3.1302-02     8  6.1054-03  2.35  6.1061-02  2.35  6.1727-03  2.34  16  9.3378-04  2.71  9.3381-04  2.71  9.6048-04  2.68  32  1.2856-04  2.86  1.2856-04  2.86  1.4026-04  2.78  64  1.6851-05  2.93  1.6851-05  2.93  2.1901-05  2.68  128  2.1567-06  2.97  2.1567-06  2.97  4.1048-06  2.42  196  6.0560-07  2.98  6.0560-07  2.98  1.6071-06  2.20  5.2 Problems with unknown solution In this part of our numerical study, we do not prescribe a solution but give only the data for the problem. Therefore it is not clear whether Assumption 2.1 holds. In order to calculate the errors and the value of $$B_0$$, we replace the exact solution with a reference solution computed on a Bakhvalov–Shishkin mesh with $$N=128$$ and $$Q_4\times Q_3$$. This approach is a variant of the double-mesh principle; see Farrell et al. (2000). Let us start with   −εΔu−((11)⋅∇)u+u+∇p =f in Ω=(0,1)2,divu =0 in Ω,u =0 on ∂Ω, (5.3) where the right-hand side $${\boldsymbol{f}}=(1,x)^{\rm T}$$ was chosen to be simple but not a gradient. First we look at the convergence of our method with full stabilization for $${\varepsilon}=10^{-8}$$ on a Bakhvalov–Shishkin mesh. We present the results for $$Q_2\times Q_1$$ and $$Q_3\times Q_2$$ only as the counterparts with discontinuous pressure approximations behave similarly. The left part of Table 6 shows the error behaviour measured in the energy norm. Clearly the convergence order is far from the theoretically predicted values and does not improve upon increasing the polynomial order. A possible explanation is that Assumption 2.1 is not fulfilled and some corner singularities arise. This is similar to convection–diffusion problems where additional compatibility conditions of the right-hand side have to be fulfilled; see, e.g., Linß & Stynes (2001), Kellogg & Stynes (2005, 2007). Thus, we modify our problem by changing the right-hand side to   f2:=(1−1)(1−x−y)(x−y). Table 6 Problem (5.4): Error in energy norm $${\left|\!\!\;\left|\!\!\;\left| {(\cdot,\cdot)} \right|\!\!\;\right|\!\!\;\right|}_{\varepsilon}$$ on Bakhvalov–Shishkin meshes with $${\varepsilon}=10^{-8}$$ $${\boldsymbol{f}}=(1,x)^{\rm T}$$  $$N$$  $$Q_2\times Q_1$$  $$Q_3\times Q_2$$  8  3.178-02  2.244-02  16  1.787-02  1.258-02  32  1.121-02  8.035-03  64  6.780-03  4.891-03  Order  0.73  0.72  $${\boldsymbol{f}}=(1,x)^{\rm T}$$  $$N$$  $$Q_2\times Q_1$$  $$Q_3\times Q_2$$  8  3.178-02  2.244-02  16  1.787-02  1.258-02  32  1.121-02  8.035-03  64  6.780-03  4.891-03  Order  0.73  0.72  $${\boldsymbol{f}}=(1-x-y)(x-y)(1,-1)^{\rm T}$$  $$N$$  $$Q_2\times Q_1$$  $$Q_3\times Q_2$$  8  9.403-03  1.300-03  16  2.488-03  2.495-04  32  6.322-04  4.794-05  64  1.571-04  9.191-06  Order  2.00  2.38  $${\boldsymbol{f}}=(1-x-y)(x-y)(1,-1)^{\rm T}$$  $$N$$  $$Q_2\times Q_1$$  $$Q_3\times Q_2$$  8  9.403-03  1.300-03  16  2.488-03  2.495-04  32  6.322-04  4.794-05  64  1.571-04  9.191-06  Order  2.00  2.38  This function $${\boldsymbol{f}}_2$$ is zero in the corners of $${\it{\Omega}}$$ which is a minimal compatibility condition in convection–diffusion problems. The convergence results are presented in the right part of Table 6. We do observe the optimal convergence order for $$Q_2\times Q_1$$ but still a reduced order for the higher-order pair. Thus again, not enough compatibility conditions are fulfilled to ensure that Assumption 2.1 holds for the higher derivatives. It remains an open question which conditions are needed. Using $${\boldsymbol{f}}_2$$ we can also look into the constant $$B_0$$. The obtained results are shown in Fig. 7. Again we observe for all values of $$N$$ and $${\varepsilon}$$ that $$B_0$$ is bounded away from zero. In order to get more insight into the behaviour of $$B_0$$, we now change the vector $${\boldsymbol{b}}$$ in (5.3) to   b=(cos⁡(φ)sin⁡(φ)) for $$\varphi\in(0,\pi/4]$$. For fixed values of $$N=32$$ and $${\varepsilon}=10^{-8}$$, we vary in Fig. 8 the angle $$\varphi$$. Note that $$\varphi=0$$ would give characteristic layers that we will investigate in a moment. Although the value of $$B_0$$ depends on $$\varphi$$, it is bounded away from zero over the full range of $$\varphi$$. Fig. 7. View largeDownload slide Problem (5.3): The value of $$B_0$$ for different values of $$N$$ as a function of $$\varepsilon$$ on Bakhvalov–Shishkin meshes: $$Q_2\times Q_1$$ (left) and $$Q_3\times Q_2$$ (right). Fig. 7. View largeDownload slide Problem (5.3): The value of $$B_0$$ for different values of $$N$$ as a function of $$\varepsilon$$ on Bakhvalov–Shishkin meshes: $$Q_2\times Q_1$$ (left) and $$Q_3\times Q_2$$ (right). Fig. 8. View largeDownload slide Problem (5.3): The value of $$B_0$$ as function of $$\varphi$$ on Bakhvalov–Shishkin meshes. Fig. 8. View largeDownload slide Problem (5.3): The value of $$B_0$$ as function of $$\varphi$$ on Bakhvalov–Shishkin meshes. Let us now consider the case of characteristic layers, e.g.,   −εΔu−((10)⋅∇)u+u+∇p =(1−1)(1−x−y)(x−y) in Ω=(0,1)2,divu =0 in Ω,u =0 on ∂Ω. (5.4) Here our theory about the decomposition of $${\boldsymbol{u}}$$ does not hold at all. Nevertheless, we can look at the corresponding convection–diffusion case where a decomposition is known; see Kellogg & Stynes (2005, 2007), and use the analogy of the layer structure. Thus we expect $$u_1$$ to have characteristic layers at $$y=0$$ and $$y=1$$ and $$u_2$$ to have an exponential layer at $$x=0$$. An approximation to the solution is shown in Fig. 9. With this assumption on the layer structure, a suitable layer-adapted mesh can be constructed in the usual way and the analysis presented could also be done. We present only the convergence results in Table 7 that are similar to those of only exponential layers with reduced regularity. Also the constant $$B_0$$ is bounded away from zero for all $${\varepsilon}$$ and $$N$$. Table 7 Problem (5.4): Error in energy norm $${\left|\!\!\;\left|\!\!\;\left| {{(\cdot,\cdot)}} \right|\!\!\;\right|\!\!\;\right|}_{\varepsilon}$$ on Bakhvalov–Shishkin meshes with $${\varepsilon}=10^{-8}$$ $$N$$  $$Q_2\times Q_1$$  $$Q_3\times Q_2$$  8  9.403-03  1.300-03  16  2.488-03  2.495-04  32  6.322-04  4.794-05  64  1.571-04  9.191-06  Order  2.00  2.38  $$N$$  $$Q_2\times Q_1$$  $$Q_3\times Q_2$$  8  9.403-03  1.300-03  16  2.488-03  2.495-04  32  6.322-04  4.794-05  64  1.571-04  9.191-06  Order  2.00  2.38  Fig. 9. View largeDownload slide Problem (5.4): Solution of (5.4) with exponential and characteristic layers in the velocity. Fig. 9. View largeDownload slide Problem (5.4): Solution of (5.4) with exponential and characteristic layers in the velocity. 6. Outlook This article is a starting point for the analysis of more complex domains and problems. The analysis relies on a solution decomposition and therefore a precise knowledge of the layer structure of $${\boldsymbol{u}}$$. Having such a structure, it is relatively easy to extend the techniques of this article to the three-dimensional case or more general domains. Numerical experiments do actually indicate that the solution structure in three dimensions is similar to the two-dimensional situation. In the case of second-order reaction–diffusion problems, more general but smooth domains were investigated in, e.g., Xenophontos & Fulton (2003), Kopteva (2009). Their techniques could be incorporated into the system case. For more general domains and problems, the type of layers can change from exponential to characteristic (see Roos et al., 2008 for these types of layers), but the exact transition from one type to another is not obvious. A technique to handle characteristic layers is the use of problem-specific enriched function spaces where the enrichment could be based on the law of wall; see Krank & Wall (2016) for a recent development. For flows in channels and around obstacles in unbounded domains, one tries to prevent strong boundary layers at the outflow. The idea is to apply suitable boundary conditions that allow us to restrict the domains without significant changes of the flow solutions. Such methods can be found, e.g., in Griffiths (1997), Braack & Mucha (2014). Another topic for further research is to show that $$B_0$$ is actually bounded away from zero, independent of $$N$$ and $${\varepsilon}$$. Having this would provide a tool for proving the uniform error estimation. As indicated in Fig. 8, the value of $$B_0$$ stays away from zero even for the transition from exponential to characteristic layers. Assuming a suitable layer structure for the solution of the Oseen problem with characteristic layers, a convergence behaviour similar to the case with only exponential layers is obtained; compare Tables 6 (right) and 7. Funding K.H. was funded by the grant FR 3052/2-1 of the German Research Foundation (DFG). References Ainsworth M., Barrenechea G. R. & Wachtel A. ( 2015) Stabilization of high aspect ratio mixed finite elements for incompressible flow. SIAM J. Numer. Anal. , 53, 1107– 1120. Google Scholar CrossRef Search ADS   Apel T. ( 1999) Anisotropic Finite Elements: Local Estimates and Applications . Advances in Numerical Mathematics. Stuttgart: B. G. Teubner. Apel T., Knopp T. & Lube G. ( 2008) Stabilized finite element methods with anisotropic mesh refinement for the Oseen problem. Appl. Numer. Math. , 58, 1830– 1843. Google Scholar CrossRef Search ADS   Apel T. & Maharavo Randrianarivony H. ( 2003) Stability of discretizations of the Stokes problem on anisotropic meshes. Math. Comput. Simul. , 61, 437– 447. Google Scholar CrossRef Search ADS   Apel T. & Matthies G. ( 2008) Nonconforming, anisotropic, rectangular finite elements of arbitrary order for the Stokes problem. SIAM J. Numer. Anal. , 46, 1867– 1891. Google Scholar CrossRef Search ADS   Bakhvalov N. S. ( 1969) The optimization of methods of solving boundary value problems with a boundary layer. U.S.S.R. Comput. Math. Math. Phys. , 9, 139– 166. Google Scholar CrossRef Search ADS   Bernardi C. & Maday Y. ( 1997) Spectral methods. Handbook of Numerical Analysis , vol. V, pp. 209– 485. Amsterdam: North-Holland. Google Scholar CrossRef Search ADS   Braack M. & Mucha P. B. ( 2014) Directional do-nothing condition for the Navier-Stokes equations. J. Comput. Math. , 32, 507– 521. Google Scholar CrossRef Search ADS   Brooks A. N. & Hughes T. J. R. ( 1982) Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations. Comput. Methods Appl. Mech. Eng. , 32, 199– 259. Google Scholar CrossRef Search ADS   Burman E. & Linke A. ( 2008) Stabilized finite element schemes for incompressible flow using Scott–Vogelius elements. Appl. Numer. Math. , 58, 1704– 1719. Google Scholar CrossRef Search ADS   Dallmann H., Arndt D. & Lube G. ( 2016) Local projection stabilization for the Oseen problem. IMA J. Numer. Anal. , 36, 796– 823. Google Scholar CrossRef Search ADS   Farrell P. A., Hegarty A. F., Miller J. J. H., O’Riordan E. & Shishkin G. I. ( 2000) Robust Computational Techniques for Boundary Layers . Applied Mathematics (Boca Raton), vol. 16. Boca Raton, FL: Chapman & Hall/CRC. Franca L. P. & Frey S. L. ( 1992) Stabilized finite element methods: II. The incompressible Navier–Stokes equations. Comput. Methods Appl. Mech. Eng. , 99, 209– 233. Google Scholar CrossRef Search ADS   Franz S. ( 2008) Continuous interior penalty method on a Shishkin mesh for convection-diffusion problems with characteristic boundary layers. Comput. Meth. Appl. Mech. Eng. , 197, 3679– 3686. Google Scholar CrossRef Search ADS   Franz S., Höhne K. & Waurick M. ( 2016) A solution decomposition for a singularly perturbed 4th order problem (submitted). Franz S., Linß T. & Roos H.-G. ( 2008) Superconvergence analysis of the SDFEM for elliptic problems with characteristic layers. Appl. Numer. Math. , 58, 1818– 1829. Google Scholar CrossRef Search ADS   Franz S. & Matthies G. ( 2010) Local projection stabilisation on S-type meshes for convection-diffusion problems with characteristic layers. Computing , 87, 135– 167. Google Scholar CrossRef Search ADS   Gelhard T., Lube G., Olshanskii M. A. & Starcke J.-H. ( 2005) Stabilized finite element schemes with LBB-stable elements for incompressible flows. J. Comput. Appl. Math. , 177, 243– 267. Google Scholar CrossRef Search ADS   Girault V. & Raviart P.-A. ( 1986) Finite Element Methods for Navier–Stokes Equations . Berlin: Springer. Google Scholar CrossRef Search ADS   Griffiths D. F. ( 1997) The “no boundary condition” outflow boundary condition. Int. J. Numer. Methods Fluids , 24, 393– 411. Google Scholar CrossRef Search ADS   Hansbo P. & Szepessy A. ( 1990) A velocity-pressure streamline diffusion method for the incompressible Navier–Stokes equations. Comput. Methods Appl. Mech. Eng. , 84, 175– 192. Google Scholar CrossRef Search ADS   Hughes T. J. R., Franca L. P. & Balestra M. ( 1986) A new finite element formulation for computational fluid dynamics. V: circumventing the Babuska–Brezzi condition: a stable Petrov–Galerkin formulation of the Stokes problem accommodating equal-order interpolations. Comput. Methods Appl. Mech. Eng. , 59, 85– 99. Google Scholar CrossRef Search ADS   Jenkins E. W., John V., Linke A. & Rebholz L. G. ( 2014) On the parameter choice in grad-div stabilization for the Stokes equations. Adv. Comput. Math. , 40, 491– 516. Google Scholar CrossRef Search ADS   Johnson C. & Saranen J. ( 1986) Streamline diffusion methods for the incompressible Euler and Navier–Stokes equations. Math. Comput. , 47, 1– 18. Google Scholar CrossRef Search ADS   Kellogg R. B. & Stynes M. ( 2005) Sharpened and corrected version of: corner singularities and boundary layers in a simple convection-diffusion problem. J. Differ. Equ. , 213, 81– 120. Google Scholar CrossRef Search ADS   Kellogg R. B. & Stynes M. ( 2007) Sharpened bounds for corner singularities and boundary layers in a simple convection-diffusion problem. Appl. Math. Lett. , 20, 539– 544. Google Scholar CrossRef Search ADS   Kopteva N. ( 2009) Numerical analysis of a 2d singularly perturbed semilinear reaction-diffusion problem. Numerical Analysis and Its Applications  ( Margenov S. Vulkov L. G. & Waśniewski J. eds). Lecture Notes in Computer Science, vol. 5434. Berlin: Springer, pp. 80– 91. Google Scholar CrossRef Search ADS   Krank B. & Wall W. A. ( 2016) A new approach to wall modeling in LES of incompressible flow via function enrichment. J. Comput. Phys. , 316, 94– 116. Google Scholar CrossRef Search ADS   Linke A. ( 2014) On the role of the Helmholtz decomposition in mixed methods for incompressible flows and a new variational crime. Comput. Methods Appl. Mech. Eng. , 268, 782– 800. Google Scholar CrossRef Search ADS   Linke A. & Rebholz L. G. ( 2013) On a reduced sparsity stabilization of grad-div type for incompressible flow problems. Comput. Methods Appl. Mech. Eng. , 261/262, 142– 153. Google Scholar CrossRef Search ADS   Linß T. ( 1999) An upwind difference scheme on a novel Shishkin-type mesh for a linear convection-diffusion problem. J. Comput. Appl. Math. , 110, 93– 104. Google Scholar CrossRef Search ADS   Linß T. ( 2000) Analysis of a Galerkin finite element method on a Bakhvalov-Shishkin mesh for a linear convection-diffusion problem. IMA J. Numer. Anal. , 20, 621– 632. Google Scholar CrossRef Search ADS   Linß T. ( 2010) Layer-Adapted Meshes for Reaction-Convection-Diffusion Problems . Lecture Notes in Mathematics, vol. 1985. Berlin: Springer. Google Scholar CrossRef Search ADS   Linß T. & Stynes M. ( 2001) Asymptotic analysis and Shishkin-type decomposition for an elliptic convection-diffusion problem. J. Math. Anal. Appl. , 261, 604– 632. Google Scholar CrossRef Search ADS   Ludwig L. ( 2017) https://github.com/SOFE-Developers/SOFE, last accessed 17 March 2017. Matthies G., Lube G. & Röhe L. ( 2009) Some remarks on residual-based stabilisation of inf-sup stable discretisations of the generalised Oseen problem. Comput. Methods Appl. Math. , 9, 368– 390. Google Scholar CrossRef Search ADS   Matthies G. & Tobiska L. ( 2015) Local projection type stabilization applied to inf–sup stable discretizations of the Oseen problem. IMA J. Numer. Anal. , 35, 239– 269. Google Scholar CrossRef Search ADS   Miller J. J. H., O’Riordan E. & Shishkin G. I. ( 1996) Fitted Numerical Methods for Singular Perturbation Problems: Error Estimates in the Maximum Norm for Linear Problems in One and Two Dimensions.  River Edge, NJ: World Scientific Publishing Co. Google Scholar CrossRef Search ADS   Olshanskii M. A. ( 2002) A low order Galerkin finite element method for the Navier-Stokes equations of steady incompressible flow: a stabilization issue and iterative methods. Comput. Methods Appl. Mech. Eng. , 191, 5515– 5536. Google Scholar CrossRef Search ADS   Olshanskii M. A., Lube G., Heister T. & Löwe J. ( 2009) Grad-div stabilization and subgrid pressure models for the incompressible Navier-Stokes equations. Comput. Methods Appl. Mech. Eng. , 198, 3975– 3988. Google Scholar CrossRef Search ADS   Olshanskii M. A. & Reusken A. ( 2004) Grad-div stabilization for the Stokes equations. Math. Comput. , 73, 1699– 1718. Google Scholar CrossRef Search ADS   Roos H.-G. ( 2012) Robust numerical methods for singularly perturbed differential equations: a survey covering 2008–2012. ISRN Appl. Math. , 2012, Article ID 379547. Roos H.-G. & Linß T. ( 1999) Sufficient conditions for uniform convergence on layer-adapted grids. Computing , 63, 27– 45. Google Scholar CrossRef Search ADS   Roos H.-G., Stynes M. & Tobiska L. ( 2008) Robust Numerical Methods for Singularly Perturbed Differential Equations . Springer Series in Computational Mathematics, vol. 24, 2nd edn. Berlin: Springer. Schötzau D. & Schwab C. ( 1998) Mixed $$hp$$-FEM on anisotropic meshes. Math. Models Methods Appl. Sci. , 8, 787– 820. Google Scholar CrossRef Search ADS   Schötzau D., Schwab C. & Stenberg R. ( 1999) Mixed $$hp$$-FEM on anisotropic meshes. II. Hanging nodes and tensor products of boundary layer meshes. Numer. Math. , 83, 667– 697. Google Scholar CrossRef Search ADS   Stynes M. & O’Riordan E. ( 1997) A uniformly convergent Galerkin method on a Shishkin mesh for a convection-diffusion problem. J. Math. Anal. Appl. , 214, 36– 54. Google Scholar CrossRef Search ADS   Stynes M. & Tobiska L. ( 2003) The SDFEM for a convection-diffusion problem with a boundary layer: optimal error analysis and enhancement of accuracy. SIAM J. Numer. Anal. , 41, 1620– 1642. Google Scholar CrossRef Search ADS   Stynes M. & Tobiska L. ( 2008) Using rectangular $$Q_p$$ elements in the SDFEM for a convection-diffusion problem with a boundary layer. Appl. Numer. Math. , 58, 1709– 1802. Google Scholar CrossRef Search ADS   Tezduyar T. E., Mittal S., Ray S. E. & Shih R. ( 1992) Incompressible flow computations with stabilized bilinear and linear equal order interpolation velocity pressure elements. Comput. Methods Appl. Mech. Eng. , 95, 221– 242. Google Scholar CrossRef Search ADS   Tobiska L. & Lube G. ( 1991) A modified streamline diffusion method for solving the stationary Navier–Stokes equations. Numer. Math. , 59, 13– 29. Google Scholar CrossRef Search ADS   Xenophontos C. & Fulton S. R. ( 2003) Uniform approximation of singularly perturbed reaction-diffusion problems by the finite element method on a Shishkin mesh. Numer. Methods Partial Differ. Equ. , 19, 89– 111. Google Scholar CrossRef Search ADS   © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.

### Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: Jan 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations