Gradient Robust Mixed Methods for Nearly Incompressible Elasticity

Within the last years pressure robust methods for the discretization of incompressible fluids have been developed. These methods allow the use of standard finite elements for the solution of the problem while simultaneously removing a spurious pressure influence in the approximation error of the velocity of the fluid, or the displacement of an incompressible solid. To this end, reconstruction operators are utilized mapping discretely divergence free functions to divergence free functions. This work shows that the modifications proposed for Stokes equation by Linke (Comput Methods Appl Mech Eng 268:782–800, 2014) also yield gradient robust methods for nearly incompressible elastic materials without the need to resort to discontinuous finite elements methods as proposed in Fu et al. (J Sci Comput 86(3):39–30, 2021).


Introduction
The Stokes equation for steady flow of an incompressible fluid is given as in a, polygonal, domain ⊂ R d ; d = 2, 3 for given data f ∈ L 2 ( ) and ν > 0, where u denotes the fluid velocity and p denotes the pressure. Under the famous inf-sup condition for the finite element spaces V h and Q h , the use of mixed finite elements allows to obtain discrete approximations u h ∈ V h and p h ∈ Q h satisfying an error estimate of the form see, e.g., [3] Here β is the inf-sup constant associated to the choice of V h and Q h , · 1 and · 0 denote the H 1 and L 2 norm on , respectively. Further, here and throughout the paper c denotes a generic constant which is independent of all relevant quantities of the estimate but may take a different value at each appearance.
While the estimate can yield asymptotically optimal orders, provided suitable finite element spaces are taken, without the need to utilize exactly divergence free finite element functions for the approximation of u h the right hand side of the estimate hints towards an undesirable influence of the pressure on the approximation error of the velocity. In fact, it has been observed, e.g., in [1] that indeed complicated pressures can give rise to a large error in the velocity approximation, even in situations where the true velocity can be represented in the discrete space V h .
A potential remedy, allowing for arbitrary inf-sup stable element pairs while providing pressure independent velocity has been proposed by [1]. He proposed the use of reconstruction operators on the right hand side of the equation to map discretely divergence free functions to divergence free functions. This proposed method has been implemented to a range of problems and a variety of finite element pairs for the discretization of Stokes equation, such as non-conforming Crouzeix-Raviart element [4], Taylor-Hood and MINI elements with continuous pressure spaces [5], on rectangular elements [6], for embedded discontinuous Galerkin methods (EDG) [7]. For 3-d polyhedral domains with concave edges a pressure robust reconstruction is given in [8]. While the obtained convergence orders are optimal, the price to pay, for these methods is a loss of quasi optimality of the method due to Strang's first lemma. Recently, [9] showed that a more involved construction of the reconstruction operator allows for a quasi-optimal discretization.
In this paper, we consider the extension of these results to nearly incompressible linear elasticity, e.g., where ε(u) denotes the symmetric gradient, and μ, λ > 0 are the Lamé parameters. To avoid the locking phenomenon, e.g., [10,Chapter VI.3], typically a mixed form is considered. Here the incompressible case, i.e., λ = ∞, can easily be included by dropping the term − 1 λ p in the second line. It is clear conceptually that the same difficulties as for the Stokes problem will occur in the incompressible limit. However, the treatment of the nearly incompressible case requires additional care. To this end, [2] defined a discretization to be "gradient robust", if the influence of gradient forces f = ∇φ in the discrete solution vanishes sufficiently fast as λ → ∞. [2] showed that a standard mixed discretization of (2) is not gradient robust and provided a gradient robust hybrid discontinuous Galerkin (HDG) scheme. Within this article, we will show that mixed methods can be made gradient robust using the approach proposed by [1] for the mixed discretization of (2).
The rest of the paper is structured as follows. In Sect. 2, we introduce the notion of gradient robustness and discuss the discretization of (2). Next, in Sect. 3, we show that the proposed discretization is indeed gradient robust and provide error estimates. We conclude the paper with a series of examples highlighting the derived results in Sect. 4.

Gradient Robustness
We define the spaces V 0 of divergence free function and its orthogonal complement V ⊥ as with the L 2 ( )-scalar product ( ·, · ). Now, any function u ∈ V can be uniquely written as Using Helmholtz decomposition, f ∈ L 2 ( ; R d ) can be uniquely decomposed as where φ ∈ H 1 ( )/R is irrotational, w is divergence free and both are orthogonal with respect to the L 2 ( )-scalar product, i.e., With these definitions, the decay of the influence of gradient forces, i.e., w = 0, onto the solutions u of (2) can be quantified as the following result from [2, Theorem 1] shows: Then for the solution u = u 0 + u ⊥ of (2) it holds u 0 = 0 and In particular, Since this bound need not hold for arbitrary, inf-sup stable, discretizations, [2, Definition 2] introduced the following notion: Definition 1 A discretization of (2) is called gradient robust, if for any fixed f = ∇φ with φ ∈ L 2 ( ), μ > 0 and any discretization parameter h there is a constant c h such that the approximate solution u h ∈ V ⊥ h and satisfies

Abstract Discretization
In order to discretize (2), we define a second bilinear form b : Now we select subspaces V h ⊂ V and Q h ⊂ Q such that there is a positive constant β satisfying the inf-sup condition Now, the standard, in general not gradient robust, weak formulation is given as follows: Under the well known inf-sup condition (7) on V h and Q h , the system (8) Following [1], we assume that there exists a reconstruction operator to be specified later in Sect. 2.3, mapping discretely divergence free functions to divergence free functions. Then the modified problem is given as: Clearly, by construction, the modified problem (10) admits a solution under the same conditions as (8), since only the right hand side has been modified. In Theorem 4, we will see that the discretization (10) is gradient robust, under appropriate assumptions on π div . Further, in Theorem 5, we show the gradient robust displacement error estimate where · k denotes the norm on H k ( ) or H k ( ; R d ); of course assuming sufficient regularity of u and p and approximation order of V h and Q h . While the introduction of a variational crime in (10) means that instead of a quasi-best approximation error we only provide an estimate of optimal convergence order the estimate (11) is clearly better than (9) in view of the asymptotics as λ → ∞ and μ → 0.

Reconstruction Operator and Assumptions
The construction of the reconstruction operator π div proposed by [1] is based on the choice of a suitable subspace M h ⊂ H div ( ;R d ) satisfying the commuting diagram in Fig. 1 where The commuting diagram is equivalently expressed by the equation Then clearly, by (12) we have that the restriction of π div to discretely divergence free functions maps into divergence free functions, i.e., and further for any v h ∈ V h it holds where, n is the unit outward normal vector. Analogously to the continuous setting, we can define the orthogonal complement V ⊥ h by Before we continue, let us make some, generic assumptions on the considered spaces V h and Q h defined on a shape regular family T h of decompositions of . Assumption 1 Following [6, Assumptions A1,A2, and A3], we assume, that for some k ≥ 2 and i = 0, 1 the finite element space V h is equipped with an interpolation operator where · i,T denotes the respective norm on the element T , and h T is the element diameter. For the space Q h , we assume that the L 2 -projection π L 2 : Further, it is assumed that V h and Q h satisfy the inf-sup inequality (7). Finally, we assume that there exists a subspace Q h ⊂ L 2 ( ; R d ) such that the respective L 2 -projection π L 2 satisfies Further requirements on Q h will be made in Assumption 2.
With these preparations, we can now state the additional assumptions on the reconstruction operator.
Assumption 2 Following [6, Assumption A4], we first assume, that the reconstruction operator satisfies the following orthogonality relation where Q h ⊂ L 2 ( ; R d ) is given in Assumption 1. Second, we assume the following local approximation property to hold Before concluding the assumption, let us note that the assumptions can indeed be satisfied.
To this end, we give an example which we will also use for the numerical results in Sect. 4.

Example 1
Let us assume that the domain can be decomposed into a family T h of shape regular we consider, parametric, piecewise Q k and globally continuous finite elements with k ≥ 2. For the discretization of Q h = Q k−1 h , we select the space of discontinuous piecewise P k−1 functions. Indeed theses pairs satisfy the inf-sup condition (7), see, e.g., [ (12) together with the canonical interpolation π div . Further, they showed [6, Lemma 2.1], that the restriction of π div to discretely divergence free functions maps into divergence free functions, i.e., Remark 1 Infact, [6] showed that (12) follow from a set of assumed orthogonality properties and surjectivity of divergence and normal traces from which suitable choices of M h and constructions of π div can be obtained.

Error Analysis
In this section, we proceed with error analysis of the modified weak form (10). We split the analysis in two parts for incompressible materials (λ = ∞) and nearly incompressible materials (λ = ∞).

Incompressible Materials
We proceed to the error analysis of incompressible materials, where λ = ∞ and the term involving 1 λ is dropped in (10). The analysis follows, at large, the arguments in [4] with some minor adjustments to the elasticity case.
where | · | k denotes the H k -semi-norm, where k ≥ 2 is given by Assumption 1.
Before proving the above theorem, we would like to prove an important lemma which is need to prove the theorem.

Lemma 3 Let Assumptions 1 and 2 be satisfied and λ = ∞. Then for any functions
where | · | k,T denotes the H k -semi-norm on T , where k ≥ 2 is given by Assumption 1.
Proof We add and subtract (∇ · ε(u), w h ) on the left to obtain Since ∇ · ε(u) ∈ L 2 ( ; R d ), we can apply the projection π L 2 , from Assumption 1, to get π L 2 ∇ · ε(u) ∈ Q h . By the assumed orthogonality in (17), we have Using Assumption 1 and (18), we obtain, for the first summand on the right of (20), For the last two summands of (20), we apply Gauss divergence theorem to get since w h = 0 on ∂ . Combining (20) with the bounds (21) and (22) the assertion is shown. Now, we continue to prove Theorem 2 proof of Theorem 2 Let u h be the solution of (10), with λ = ∞, and let v h ∈ V 0 h be arbitrary.
h and applying the triangle inequality gives In view of the interpolation estimate in Assumption 1, we are left to estimate w h 1 . From Korn's inequality, we have For the first summand on the right of (24) we use Cauchy-Schwartz inequality to get Before we come to the bound of the second summand in (24), we make some preliminary calculations. Since u h is the solution of (10), Further, since u is the solution to the equation (2) multiplication with π div w h and integration yields −2μ ∇ · ε(u)π div w h dx − ∇ p π div w h dx = fπ div w h dx by the compatibility of the reconstruction with the kernel of the divergence, i.e., (15) and (16), this gives Combining this with (26), we get Now, we can bound the second summand on the right of (24), using (27) we get By the previously shown lemma, i.e., (19), we can bound the right hand side to get Now combining (24) with the two bounds (25) and (28), we get Substituting this in (23) yields To bound the best approximation error on V 0 h in this inequality, we proceed using inf-sup condition as in [3, Chapter 2, (1.16)] and the assumed interpolation estimate on V h in Assumption 1, to get the estimate Using this in (29) gives the desired estimate.

Nearly Incompressible Materials
For the nearly incompressible case, i.e., (λ = ∞), we start by assuming a gradient force f = ∇φ, for some φ ∈ L 2 ( ). From Lemma 1, we have that the solution of (2) for such an f is u = u ⊥ . The following result shows, that our mixed discretization (10) is gradient robust in the sense of Definition 1. (10) is a gradient field, i.e., f = ∇φ, for some φ ∈ L 2 ( ), then the solution (u h , p h ) ∈ V h × Q h of (10) with λ = ∞ satisfies u h ∈ V ⊥ h and the gradient robust bound

Theorem 4 Let Assumptions 1 and 2 be satisfied. If the right hand side
with a constant c independent of h.

Proof
Consider v h = u h in equation (10) with f = ∇φ. Then integration by parts for the right hand side, using the zero trace from (16), we get Since ∇ · π div u h ∈ Q h we can rewrite the right hand side as Since π L 2 ∇ · u h ∈ Q h , we can use it to test the second line in (10) giving Substituting (32) and (33) in (31), we get Now π L 2 φ ∈ Q h and u h ∈ V h hence, by (12), it holds Filling this into (34) gives Using Cauchy-Schwartz inequality, we get Now, to estimate the H 1 -norm of u h , we notice that by the choice of f and (15), testing the first equation in (10) and thus u h ∈ V ⊥ h . Hence by, e.g., [13,Lemma 3.58] it holds with a constant c depending on the inf-sup constant β from (7), since π L 2 ∇ · u h ∈ Q h . Using Korn's inequality, (36), and (37), we get and thus the assertion is shown.
Theorem 5 Let Assumptions 1 and 2 be satisfied. Then the solutions (u, p) ∈ V × Q, of the problem (2) and (u h , p h ) ∈ V h × Q h of (10) satisfy the error estimate provided the regularity (u, p) ∈ H k+1 ( ; R d ) × H k ( ) is given, where k ≥ 2 is given by Assumption 1.
Proof As in the proof of Theorem 2, we could split the error with arbitrary v h ∈ V h and q h ∈ Q h . However, as it will turn out to be useful, we will select q h = π L 2 p and v h as a particular Fortin operator applied to u, i.e., satisfying the following equation Clearly, the solution to the continuous counterpart is (v,p) = (u, 0). Since the above equation is uniquely solvable, see, e.g. [11,Theorem 4.2.3], we have the orthogonality b(π L 2 p − p h , u − v h ) = 0 and the approximation error satisfies, e.g., [11,Theorem 5.2.2].
which gives Due to the interpolation estimates in Assumption 1, we are left with bounding By definition of the bilinear forms a and b, i.e., (3) and (6), and the first line in (10) and (2), the remainder w h and r h satisfy, for any discrete function Analogously, from the second line in (10) and (2), we get for arbitrary Testing (43) and (44) with ' h = w h and s h = r h we get by our choice of v h . This provides the bound Of course (47) provides a bound on w h but as it is suboptimal, in view of the λ dependence, we continue by splitting w h = w 0 h + w ⊥ h . We first bound w 0 h 1 . Consider cμ w 0 h 1 and using that a(w ⊥ h , w 0 h ) = 0, we have, using (13), (43), and the choice of v h by (40) that Lemma 3.58], we get with a constant c depending on the inf-sup constant

Thus, by Lemma 3, we conclude
from the definition of v h in (40). Hence, noting that ∇ · u = 1 λ p, we obtain We conclude from (47) and thus we obtain the final bound on w ⊥ Now, we can bound w h 1 using (48) and (49) Finally, we arrive at the desired bound by definition of q h and Assumption 1.

Numerical Results
For our computation, we use DOpElib [14] based on the deal.II [15] finite element library with rectangular meshes. All examples are posed on square domains and the meshes are obtained by bisection. For the computation we considered the inf-sup stable Taylor-Hood element (Q 2 × Q 1 ), for comparison of our results with [2]. Further, we utilized the inf-sup stable discretization Q 2 × DGP 1 (discontinuous P 1 pressure) and its gradient robust modification by interpolation into BDM 2 as discussed in Example 1.
First, we present an example for incompressible materials.
Example 2 For the first numerical example, we consider a small variation of Example 5.1 in [1], where the displacement and pressure on the domain = (0, 1) 2 is given as for the incompressible linear elasticity equation with homogeneous boundary conditions on u and thus define f, of course the pressure is defined up to a constant only. For future examples, we consider nearly incompressible materials given by equation (2). From Lemma 1, the solution for Example 3 in the limiting case (λ = ∞) is given as u ∞ = 0 and p ∞ = x 6 + y 6 . From equation (30), we have the bound for the solution u λ h as on the discrete function for a gradient robust discretization. For μ = 10 −5 , we have λ + μ ≈ λ, ∀λ ≥ 1. Hence, we see a green line with positive slope in Fig. 3a for the gradient robust method, while the non robust method shows an almost constant u λ h 1 = 0. However, for λ = 10 5 we have 1 λ+μ ≈ c(constant) ∀0 < μ ≤ 1, which is seen in the flat green line in Fig. 3b.
For non-gradient robust methods, we have from equation (9). For μ = 10 −5 , the term 1 λ + 1 → 1 as λ → ∞. The same is shown by the flat red line in Fig. 3a. However, for λ = 10 −5 , we have u λ h 1 ≤ c μ φ 0 . Which is shown by the red line with negative slope in Fig. 3b.
It should be noted in this example, that the (blue with triangles) line for the non-gradient robust Q 2 × DGP 1 method coincides with the (green with dots) line for the gradient robust However, this appears to be due to the particular problem data hiding the nongradient robustness of the Q 2 × DGP 1 discretization. That indeed, the standard Q 2 × DGP 1 method is not gradient robust is shown in the following example.
Example 4 For the third numerical example, we consider the right hand side f = ∇φ; φ = − 10(x − 0.5) 3 y 2 + (1 − x) 3 (y − 0.5) 3 + 1/8 in Example 3 while keeping the homogeneous boundary values, and the equation, for u. Figure 4 shows our previous statement, that Example 3 failed to show the missing gradient robustness of the standard Q 2 × DGP 1 discretization. Indeed, in this example, both Q 2 × Q 1 and Q 2 × DGP 1 discretization show the undesirable blowup for μ → 0 and the constant value as λ → ∞, while the gradient robust modification shows the desired convergence.
Example 5 For the fourth numerical example, we consider, again, the nearly incompressible case with homogeneous boundary conditions on u, the values of u ∞ and thus f are given as in Example 2 and the domain = (0, 1) 2 .
In this example, for λ = ∞, the solution u ∞ is known, i.e., it is given in (52). We denote the solution, for λ = ∞, as u λ , p λ . We compute the error u ∞ − u λ h 1 in our numerical results, where u λ h is the discrete approximated solution for a given value of λ. Since, Theorem 5 provides an estimate, for u λ − u λ h 1 , we use the triangle inequality to get (55) Figure 5b follows the same pattern as Fig. 4b. However, there is a slight difference between Figs. 4a and 5a, which can be explained by (55). In the limit λ → ∞ and fixed μ = 10 −5 , we have 1 + μ λ → 1 and u ∞ − u λ 1 → 0, and we observe  for fixed refinement as it is shown in Fig. 5a. Figure 5a further confirms (56) as we can see the order O(h 2 ) for u ∞ − u λ h 1 for large values of λ. In Fig. 5b, we observe the convergence u ∞ − u λ h 1 → u ∞ − u λ as h → 0 and the decay of the error u ∞ − u λ → 0 as λ → 0 Example 6 Finally, we would like to compare our results with the thermo-elastic solids example given in [2,Sect. 6]. The gradient force f is given by a temperature θ as The material used is a nearly incompressible hard rubber with Young's Modulus E = 5 × It is important to note that θ ∈ H 1 ( ) and thus f ∈ L 2 ( ; R 2 ). For numerical computation, we additionally solve the temperature equation by a standard H 1 -conforming finite element discretization. Hence, the finite element spaces now consist of three components, the first two denote the displacement and pressure discretization as before. The third element, always Q 2 , is used to solve the equation for the temperature θ .
In Fig. 7, we can see that we achieve a well represented solution for the displacement with only 64 elements using a gradient robust method, and the magnitude is already captured with only 16 elements. In comparison, the non gradient robust methods require 256 and 1024 elements, respectively, to get a solution of similar shape and magnitude, see Figs. 8 and 9.

Conclusion
In this paper, we have shown that a gradient robust modification of nearly incompressible elasticity is possible by the same techniques proposed for incompressible flows. For this gradient robust methods, we have shown convergence estimates of optimal order w.r.t the mesh size and optimal dependence on the Lamé-constants. Several numerical examples highlighted the proven convergence rates.