Pressure-robust error estimate of optimal order for the Stokes equations: domains with re-entrant edges and anisotropic mesh grading

The velocity solution of the incompressible Stokes equations is not affected by changes of the right hand side data in form of gradient fields. Most mixed methods do not replicate this property in the discrete formulation due to a relaxation of the divergence constraint which means that they are not pressure-robust. A recent reconstruction approach for classical methods recovers this invariance property for the discrete solution, by mapping discretely divergence-free test functions to exactly divergence-free functions in the sense of H(div)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{H}}({\text {div}})$$\end{document}. Moreover, the Stokes solution has locally singular behavior in three-dimensional domains near concave edges, which degrades the convergence rates on quasi-uniform meshes and makes anisotropic mesh grading reasonable in order to regain optimal convergence characteristics. Finite element error estimates of optimal order on meshes of tensor-product type with appropriate anisotropic grading are shown for the pressure-robust modified Crouzeix–Raviart method using the reconstruction approach. Numerical examples support the theoretical results.


Introduction
When considering polyhedral domains, the solution of the Stokes equations shows in general singular behavior near corners and edges. On quasi-uniform meshes, this leads to sub-optimal performance of standard numerical methods, which can be remedied by local mesh refinement near the singular sections of the boundary. Isotropic refinement can compensate the negative effect of the singular solution, but also leads to over-refinement near edges and thus a waste of computational resources. Anisotropic refinement on the other hand can recover the optimal convergence rate [6,11,12], while the number of elements N in the mesh still satisfies N ∼ h −3 , where h is the mesh size parameter.
Unfortunately, many classical mixed methods do not fulfill the discrete inf-sup stability condition independently of the aspect ratio of the triangulation, which may be unbounded in the case of anisotropic grading. For instance, the standard proof of the inf-sup condition for the Taylor-Hood and Mini-element leads to a constant that depends on the aspect ratio. While for the lowest order Taylor-Hood pair a new proof has been found recently that shows inf-sup stability on a class of anisotropic meshes [13], the Mini-element is reported to become unstable with decreasing minimum angle in the triangulation [2].
However, several inf-sup stable methods are known for anisotropic triangulations in two dimensions, e.g. the Bernardi-Raugel and related elements [10] and the stabilized Q 1 ∕Q 1 , Q 1 ∕Q 0 and rotated Q 1 ∕P 0 elements for quadrilaterals [14,15]. Additionally, results are available for the hp-version finite element method, see e.g. [4,5,36]. The Crouzeix-Raviart element [18], which we will focus on in this contribution, is inf-sup stable on simplicial triangulations in two and three dimensions, without any condition on the mesh [12,Lemma 3.1].
In addition to its low regularity near concave edges, the velocity solution of the Stokes problem is not affected by changes in form of gradient fields on the right hand side. This property leads to the notion velocity-equivalence of forces, i.e. f 1 , f 2 ∈ L 2 (Ω) are velocity-equivalent, f 1 ≃ f 2 , if they lead to the same velocity solution of (1), see [23]. That is the case if and only if they differ by a gradient field, see e.g. [3,27,31]. Reproducing this continuous property on the discrete level poses an additional difficulty for discretization schemes, and most classical methods do not overcome it. In fact, error estimates for classical H 1 -conforming methods are in general of the form, see e.g. [24,27], where ∥ ⋅ ∥ k is the standard H k (Ω)-Sobolev norm and C F is the stability constant of the Fortin operator of the mixed method. Some divergence-free H 1 -conforming methods exist, e.g. the Scott-Vogelius element, see [39], and the rational bubble enriched methods from [25,26], where the error estimate does not contain the second, pressure-dependent term, for more references see also [27]. For the non-divergence-free methods however, this type of estimate implies that in settings where the (1a) − Δu + ∇p = f , (1b) ∇ ⋅ u = 0, continuous pressure is more difficult to approximate compared to the velocity, the velocity approximation can be highly inaccurate. Consider for example a hydrostatic case where the exact velocity is given as u ≡ 0 and the continuous pressure is a polynomial of order k. Then for classical methods using piecewise polynomials of order less than k for the pressure approximation, in general inaccurate discrete velocity solutions u ≢ 0 and a locking effect for → 0 can be observed, where the errors may become arbitrarily large. In contrast, so called pressure-robust methods, i.e. methods that see the velocity-equivalence of forces, yield the exact velocity solution, even for lowest-order mixed methods with piecewise constant pressure approximation, see [23,Section 2.5].
Additionally to the naturally pressure-robust class of H(div)-conforming finite element methods, see e.g. [16,37,38], and the already mentioned divergence-free H 1 -conforming methods from [25,26,39], a recent approach using a reconstruction operator on the velocity test functions in the linear form showed that most classical mixed methods can be made pressure-robust at the cost of an additional consistency error, see e.g. [27,28,30]. In [9] the pressure-robust modified Crouzeix-Raviart element was analyzed on anisotropic triangulations, using the assumption of a regular solution, i.e. (u, p) ∈ H 2 (Ω) × H 1 (Ω) . In this article we extend those results to the case of domains with concave edges and low regularity of the exact solutions. Recently in [31], the reconstruction approach was investigated with minimal assumptions on the regularity of the solution, but without admitting the use of anisotropic triangulations.
The main contribution of this paper is a pressure-robust estimate for the velocity solution of the modified Crouzeix-Raviart method in low-regularity settings due to non-smooth domains. The estimate shows that when appropriate anisotropic mesh grading towards a non-convex edge is used, an optimal convergence rate can be achieved by the pressure robust method. We provide an estimate on the pressure error for the modified method in the anisotropic setting as well. Numerical examples support the theoretical results, and show that in cases with low regularity data the pressure-robust variants can be superior to the standard method.
The article is structured as follows. Section 2 introduces the problem and basic notation. The type of mesh and the modified Crouzeix-Raviart method is described in Sect. 3, and Sect. 4 shows some aspects of the Helmholtz-Hodge decomposition of vector fields which are important to the analysis. Section 5 contains the a-priori error analysis, Sect. 6 shows the performance of the method with the help of two numerical examples.

Continuous Stokes problem
Consider a prismatic domain Ω = G × Z , where G is a polygonal shape with one concave corner at which the interior angle is ∈ ( , 2 ) , and Z is a bounded interval. To facilitate notation we assume that the non-convex corner of G is placed at the origin, i.e. the relevant edge of Ω is located on the z-axis. On the domain Ω , consider the incompressible Stokes equations (1) with homogeneous Dirichlet boundary condition where is the kinematic viscosity and vector valued quantities are denoted in bold symbols. For f ∈ L 2 (Ω) , the corresponding weak formulation given by has a unique solution (u, p) ∈ X × Q , see e.g. [24,Section I.5.1], where and (⋅, ⋅) denotes the L 2 (Ω) scalar product. With the space of divergence-free functions we can reformulate the problem, see [24, Section I.5.1]: find u ∈ V 0 , so that Additionally to the well known Stokes theory, see e.g. [24], which states the regularity of the solution in the Hilbert space case as above, Theorem 2.1 in [22] classifies the solution in a more general setting: for f ∈ W −1,q (Ω) and appropriate regularity of the boundary condition we have (u, p) ∈ W 1,q (Ω) × L q 0 (Ω) , with 1 < q < ∞. For the special case of convex prismatic polyhedral domains, we can assume that the solution of problem (3) satisfies (u, p) ∈ H 2 (Ω) × H 1 (Ω) , see [19]. This is in general not the case for non-convex geometries like the ones we are considering, but Theorem 2.1 in [12] gives a regularity result for our case in weighted Sobolev spaces. In particular, the derivatives of the solution in the direction parallel to the concave edge have the standard regularity, i.e. z u ∈ H 1 (Ω) and z p ∈ L 2 (Ω) . The global regularity of u is however characterized by r , where r is the distance to the singular edge and is the smallest positive solution of for which 1 ∕2 < < ∕ holds, see [19]. Figure 1 shows the type of anisotropically graded tensor-product mesh used for the discretization of the problem, and we briefly describe the process of mesh generation in the following paragraph. This type of mesh was introduced in [7] for the treatment of edge singularities that occur in the Poisson problem, and was used in subsequent works also for the Stokes and Maxwell equations [11,12,21,35].

Discretization
Let D h be a conforming, shape regular triangulation of the two-dimensional domain G, which has a mesh size parameter h = max D∈D h h D , where h D = diam (D) . This mesh is graded towards the non-convex corner, so that the size of every element satisfies where r D = inf x∈D {dist(x, 0)} is the distance of an element D ∈ D h to the concave corner, ∈ (0, 1] is a grading parameter and R > 0 is the radius of the refinement zone. The graded two dimensional mesh is extended into the z-direction with uniform mesh size h 3 ∼ h . The resulting prisms are subdivided into tetrahedra, which form the simplicial mesh T h . With r T being the distance of an element T ∈ T h to the z-axis and h 1,T , h 2,T , and h 3,T the length of the projection on the x-, y-, and z-axis, respectively, the procedure yields a mesh where and the number of elements satisfies N ∼ h −3 . By F(T h ) we denote the set of facets of the mesh T h .

Remark 1
By construction, this type of tensor-product mesh satisfies a maximum angle condition, i.e. all angles between edges and faces of the triangulation are bounded by a constant ̄< . The subsequent analysis depends on this regularity assumption on the tetrahedra, which means that meshes like the ones used in [29], where the maximum angle condition is violated, can not easily be included in the theory. Our discretization is nonconforming, thus we need tools to handle potential discontinuities at the interfaces. Let [[v]] F be the jump of a function v over a facet F, which is defined for an interior facet belonging to two elements T 1  where P k (T) denotes the space of all polynomials with maximal degree k on the element T. We also need the broken gradient , which define the derivatives elementwise for all T ∈ T h by and which are on X equivalent to the standard operators, see e.g. [20, Sections 1.2.5, 1.2.6]. The discrete gradient norm for the space X ⊕ X h is defined by For the next part we need the function spaces where n denotes the unit outward normal vector on Ω . Our discretization uses a reconstruction operator on the velocity test functions in the linear form, and in order to yield a pressure-robust method the operator needs to satisfy some properties, see e.g. [17,30,31], which we summarize in the following assumption.

Assumption 1 Assume there is a reconstruction operator
When using the approximation spaces X h and Q h , the lowest-order Raviart-Thomas and Brezzi-Douglas-Marini interpolation operators satisfy this assumption, see [17,30] for the isotropic and [9] for the anisotropic case. For our intended use of the method in an anisotropic setting, the constant in estimate (5) must be independent of the aspect ratio of the mesh. Under the mild assumption of the maximum angle condition which is satisfied for the type of mesh described above, see Remark 1, this is the case for both Raviart-Thomas interpolation, see [1], and Brezzi-Douglas-Marini interpolation, see [8].

Remark 2
In contrast to our reconstruction operator, that maps into H 0 (div, Ω) , there are approaches to use divergence preserving operators which map into H 1 0 (Ω) , see e.g. [32], where such an operator is used for theoretical purposes or [40], where the operator is used as reconstruction operator like in our case. However, these operators do not seem to work for highly anisotropic triangulations. It is an open problem to find an operator mapping discretely divergence-free functions to exactly divergencefree functions in H 1 0 (Ω) on anisotropic meshes in a stable way.
Using the discrete bilinear forms we get the discrete weak formulation where f ∈ L 2 (Ω) and I h must satisfy Assumption 1. As in the continuous case, using the space of discretely divergence constrained functions we can rewrite the problem, see [17,24,33]. Thus u h ∈ V 0 h is uniquely defined by To conclude this section, we state the well-known discrete inf-sup stability for the Crouzeix-Raviart element, see e.g. [12, Lemma 3.1].

Lemma 1 The pair of function spaces X h × Q h satisfies the discrete inf-sup condition
where the discrete inf-sup constant ̃ does not depend on the mesh size parameter h or the regularity of the mesh.

Helmholtz-Hodge decomposition
This section introduces some aspects of the Helmholtz-Hodge decomposition of vector fields, which is needed for overall context and explanation. The main idea of this section is from [31, Section 3]. Every vector field g ∈ L 2 (Ω) can be uniquely decomposed into g = ℙ(g) + ∇ , where ∈ H 1 (Ω)∕ℝ and The function ℙ(g) is called Helmholtz-Hodge projection of g , see e.g. [24, We can extend the domain of the Helmholtz-Hodge projection operator from L 2 (Ω) to H −1 (Ω) with range in (V 0 ) � , the dual space of V 0 , by defining the projection for every g ∈ H −1 (Ω) as the restriction to V 0 , i.e.
Note that it holds V 0 ⊂ L 2 (Ω) . For a more detailed and technical introduction of this extension we refer to [34,Section 2] Proof The Stokes equations satisfy a fundamental invariance property, i.e. adding a gradient field to the data only changes the pressure solution, see [27]. Thus (u, p − ) is the solution with data function ℙ(f ) = f − ∇ . Dividing the momentum equation by , we get the statement of the lemma.

A-priori error estimates
For an estimate on the finite element error, the consistency error of the method has to be estimated. For self-containedness we restate [12, Lemma 3.3] which estimates the consistency error for the standard method in the case = 1. (u, p) be the solution of the Stokes problem with = 1 and data f ∈ L 2 (Ω) . Then if the mesh is refined according to < , with from (4), the estimate We can now state the error estimate of the velocity solution of our method. It shows that for appropriately refined meshes the method has an optimal order of convergence and is pressure-robust, i.e. the estimate does not depend on the viscosity or the pressure approximability.

Lemma 3 Let
Theorem 1 Let (u, p) be the solution of (3), (u h , p h ) the solution of (6), and let the mesh be refined according to < , with from (4). In addition, let the reconstruction operator satisfy Assumption 1. Then we have the estimate Due to the Pythagoras theorem we get Using (7) we get where in the last step we used a h (u − v h , w h ) = 0 . Combined with (11) we get Denote the Helmholtz-Hodge decomposition of the data by f = ℙ(f ) + ∇ and note that ∇ ⋅ I h w h = 0 due to Assumption 1 and w h ∈ V 0 h . Using (∇ , I h w h ) = 0 , we get Factoring out in the first term of (13) we get and, due to Lemma 2, see that u is also the solution to a Stokes problem with unit viscosity and data −1 ℙ(f ) , which means we can now use Lemma 3 and estimate Using the Cauchy-Schwarz inequality and the interpolation error estimate for the operator I h from Assumption 1 we estimate for the second term of (13) Combining estimates (14), (15) with (13), inserting the result in (12), using (9) and by seeing that v h was chosen as the best-approximation of u in V 0 h , we get the final estimate ◻ The term for the approximation error can be easily bounded using known results:

Corollary 1 With the assumptions from Theorem 1 the estimate holds.
Proof Using Lemma 3.5 from [9] we get By Lemma 2 u is also the velocity solution of the Stokes problem with unit viscosity and right hand side data −1 ℙ(f ) , and thus using Lemma 3.2 from [12] and (9) we get which combined with (10) proves the statement. ◻

Remark 3
Considering Lemma 2, the relationship between the data f and the velocity solution u with regard to the viscosity parameter can be looked at from different points of view. On the one hand, in Theorem 1 we establish a velocity error estimate in terms of the divergence-free part of the Laplacian of the exact velocity ℙ(−Δu) . In this form, the estimate is pressure-robust, i.e. it does not depend on the irrotational part of the data, and it does not have an apparent dependence on the viscosity. If on the other hand by using (9) we would put the estimate in terms of the Helmholtz-Hodge projection ℙ(f ) of the data, it would still be pressure-robust, but we would see a dependence on −1 . The difference is of interest for numerical examples and the information we want to extract from them. Consider e.g. the examples from [30,Section 5]. Here the exact velocity and pressure solutions are fixed, and the data function f changes when the viscosity parameter is adjusted due to the factor in front of the Laplacian. This can nicely show the effect of pressure-robustness, since non-pressure-robust methods show a viscosity induced locking effect, i.e. the velocity error scales with −1 , while pressure-robust methods do not, as the discrete velocity solution is the same for all values of . If however the data function f is fixed, we also see a dependence on the viscosity in the error for pressure-robust methods, since the velocity solution now scales with −1 . When altering the viscosity parameter while using fixed data, pressure-robustness can still be observed by changing the irrotational part of f , i.e. adding a gradient field, which has no effect on the numerical velocity solution of pressure-robust methods.
For the pressure error we get the following estimate.

Proposition 1 With the assumptions of Theorem 1 we have the estimate
Proof Let h ∶ L 2 0 (Ω) → Q h be the L 2 -orthogonal projection onto the discrete pressure space. We start with a triangle inequality, which gives where we see that for the first term it holds Because of h p − p h ∈ Q h , we can use the inf-sup condition (8) and estimate h . Since p h is the discrete pressure solution of (6) we can further calculate where in the last step we used Lemma 3, the Cauchy-Schwarz inequality and the interpolation error estimate (5) for the reconstruction operator. Now combining the estimates, using Theorem 1 and (9) yields the desired inequality. ◻

Corollary 2 With the assumptions of Theorem 1 we have the estimate
Proof The estimate is obtained from (17) by using [12,Lemma 3.2] for the first term and (16) in combination with (9) for the second term. ◻

Remark 4
We consider only the three dimensional case, since the focus of this paper is on anisotropic elements. The main results are nevertheless valid for the corresponding two-dimensional problem in a domain with a re-entrant corner, as long as adequate local mesh grading near the corner, as described in the first part of Section 3, is applied. The proofs for the intermediate results from [11] can be adapted to fit the twodimensional setting. With them, the consistency error for the standard method, see Lemma 3, can be proved analogously to the first part of the proof of Lemma 3.2 in [12], without the additional difficulty for the third component. From there, our proofs in this section apply analogously.

Numerical examples
With the following two examples we show the performance of the pressurerobust modified Crouzeix-Raviart method with Raviart-Thomas (CR-RT) and Brezzi-Douglas-Marini (CR-BDM) reconstruction compared to the standard Crouzeix-Raviart (CR) method on anisotropically graded meshes. Considering Remark 3, we first choose the approach of fixing an exact solution, where the data changes when altering the viscosity. However, for our specified exact solution we get f ∉ L 2 (Ω) for ≠ 1 , which does not comply with the assumptions of our theory. Thus for the second example we use the other approach, where the divergence-free part ℙ(f ) of the data is fixed and only the irrotational part of f is changed in order to show pressure-robustness.

Example with fixed exact solution
Consider the inhomogeneous Stokes problem, i.e. problem (1) with boundary condition u = g on Ω , on the domain where = 3 2 . The results below show that the change to inhomogeneous boundary conditions does not impact the performance of the numerical method. We use the exact velocity and exact pressure solutions defined by where we use From (4) we get ≈ 0.54448 . The velocity solution and the singular nature of the exact pressure along the edge at r = 0 are illustrated in Figure 2.
The data function f for the numerical calculations is obtained by evaluating (1a), from which we get where f 1 = f 2 = 0 for = 1.
In [9] this example was used to show that the modified Crouzeix-Raviart method can be used for anisotropic meshes. However, no theoretical foundation for the numerical results was given, since due to the low regularity of the solution in this example, i.e. (u, p) ∉ H 2 (Ω) × H 1 (Ω) , Δu ∉ L 2 (Ω) , the results from [9] are not directly applicable. This gap in the theory is closed by this contribution, at least for the case = 1 where f ∈ L 2 (Ω).
As mentioned in Sect. 4 we know that ℙ(−Δu) ∈ L 2 (Ω) , since for = 1 by [12, Theorem 2.1] it holds z p ∈ L 2 (Ω) and thus the data function f = (0, 0, z p) is in L 2 (Ω) . As u is fixed, this does not change for other values of , even though f is no longer in L 2 (Ω) for ≠ 1.
The calculations were performed with parameter values ∈ {10 −1 , 1} and ∈ {0.4, 1} . Tables 1 and 2 contain the computed errors. Comparing the estimated order of convergence (eoc) for meshes without grading, = 1 , and with grading towards the edge, = 0.4 , shows that anisotropic grading recovers the optimal convergence rate for all methods. The results with viscosity = 10 −1 show the pressure-robustness of the modified method, as the absolute value of the velocity error does not depend on , contrary to the standard method. The modified method seems to perform optimally in the anisotropic setting even with the low regularity data in the case ≠ 1 , where the optimal convergence rate could not be observed with the standard method.

Remark 5
The data function (19) is not in L 2 (Ω) for ≠ 1 , but the right hand side integrals for our methods are still finite. However, in order to produce the shown results in Table 2 the numerical quadrature for the linear form had to be highly accurate. For our CR-BDM calculations, additionally to choosing a high quadrature degree as for the other methods, we used local mesh refinement near the singular axis.
Neither the quadrature procedure described in Remark 5 nor changing the grading parameter improved the convergence results of the standard method on graded meshes with ≠ 1 , where the optimal rate could not be observed. We do not have a proof, and irregular data do not fit our theory, as for Lemma 2 and Theorem 1 we assume f ∈ L 2 (Ω) , but the differing behavior of the methods seems to be a result of f ∉ L 2 (Ω): The velocity error estimate from [12] for the standard method, which is shown for f ∈ L 2 (Ω) , comprises the consistency error and the best approximation error, the latter being bounded in terms of the interpolation error of the Crouzeix-Raviart interpolation. While we could see the interpolation error in this test converging optimally on the graded meshes, the consistency error does not seem to converge for irregular data. In contrast to the standard Crouzeix-Raviart method, the proof of our pressure robust estimate from Sect. 5 only needs to bound the consistency error for the Helmholtz-Hodge projection of the data, which, for this example, are in L 2 (Ω) . This is the reason why the modified methods work for this example.
Since the consistency error estimate from [12] was prepared in [11] with a similar estimate for the Poisson equation, we did a further test computation for the Poisson problem with exact solution u = r 1 ∕2 sin( 2 ∕3 ) on the same meshes. For this exact solution the data are not in L 2 (Ω) as in the Stokes case, and the results showed a similar convergence behavior as the Stokes example. This is another indication that the consistency error of the Crouzeix-Raviart method causes the bad numerical performance.
The discussion in the previous paragraphs shows that the modified, pressurerobust Crouzeix-Raviart method, in addition to the advantages of the robust behavior concerning small viscosities and irrotational forces, also seems to perform better in certain cases with low regularity data and anisotropic mesh grading. Table 1 Errors and experimental convergence orders of the standard and modified Crouzeix-Raviart methods on uniform and graded meshes, = 1

Example with fixed data
Consider the same general setting as in the previous example. We now use the data where f 0 and i are chosen as with Φ( ) from (18). The function f 0 is aquired by setting = 1 in (19) and the functions i are used to show the pressure-robustness in the case of a scaled exact velocity solution: the errors for the CR-RT and CR-BDM methods do not change when adding gradient fields to the data. The exact solutions for the convergence analysis can be deduced from the first example using the considerations from Lemma 2 and Remark 3.  As before we have −Δu ∉ L 2 (Ω) , but due to our choice of the functions f 0 and i we now get f ∈ L 2 (Ω) for all calculations. The calculations were performed with viscosity parameter ∈ {10 −3 , 1} and, since the difference in convergence orders was already demonstrated in the previous example, we only use anisotropic meshes with grading parameter = 0.4. Tables 3 and 4 show the errors for both choices of i . We see that while the asymptotic convergence rates are optimal for all methods, the additional gradient part ∇ 2 in the data has a significant influence on the value of the velocity error of the standard method. In contrast, the modified methods show their pressure-robustness by yielding the same velocity solution, and thus unchanged velocity errors. The scaling of the velocity solution with −1 for fixed f is clearly visible when comparing the two tables.