Versatile stabilized finite element formulations for nearly and fully incompressible solid mechanics

Computational formulations for large strain, polyconvex, nearly incompressible elasticity have been extensively studied, but research on enhancing solution schemes that offer better tradeoffs between accuracy, robustness, and computational efficiency remains to be highly relevant. In this paper, we present two methods to overcome locking phenomena, one based on a displacement-pressure formulation using a stable finite element pairing with bubble functions, and another one using a simple pressure-projection stabilized \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbb {P}_1 - \mathbb {P}_1$$\end{document}P1-P1 finite element pair. A key advantage is the versatility of the proposed methods: with minor adjustments they are applicable to all kinds of finite elements and generalize easily to transient dynamics. The proposed methods are compared to and verified with standard benchmarks previously reported in the literature. Benchmark results demonstrate that both approaches provide a robust and computationally efficient way of simulating nearly and fully incompressible materials.


Introduction
Locking phenomena, caused by ill-conditioned global stiffness matrices in finite element analyses, are an often observed and extensively studied issue when modeling nearly incompressible, hyperelastic materials [10,18,46,84,87]. Typically, methods based on Lagrange multipliers are applied to enforce incompressibility. A common approach is the split of the deformation gradient into a volumetric and an isochoric part [38]. Here, locking commonly arises when unstable standard displacement formulations are used that rely on lin-B Christoph M. Augustin christoph.augustin@medunigraz.at 1 ear shape functions to approximate the displacement field u and piecewise-constant finite elements combined with static condensation of the hydrostatic pressure p, e.g., P 1 − P 0 elements. It is well known that in such cases solution algorithms may exhibit very low convergence rates and that variables of interest such as stresses can be inaccurate [41].
From mathematical theory it is well known that approximation spaces for the primal variable u and p have to be well chosen to fulfill the Ladyzhenskaya-Babuŝka-Brezzi (LBB) or inf-sup condition [9,19,26] to guarantee stability. A classical stable approximation pair is the Taylor-Hood element [78], however, this requires quadratic ansatz functions for the displacement part. For certain types of problems higher order interpolations can improve efficiency as higher accuracy is already reached with coarser discretizations [25,57]. In many applications though, where geometries are fitted to, e.g., capture fine structural features, this is not beneficial due to a possible increase in degrees of freedom and consequently a higher computational burden. Also for coupled problems such as electromechanical or fluid-structure-interaction models high-resolution grids for mechanical problems are sometimes required when interpolations between grids are not desired [5,51]. As a remedy for these kind of applications quasi Taylor-Hood elements with an order of 3 2 have been considered, see [62], as well as equal order linear pairs of ansatz functions which has been a field of intensive research in the last decades, see [6,48] and references therein. Unfortunately, equal order pairings do not fulfill the LBB conditions and hence a stabilization of the element is of crucial importance. There is a significant body of literature devoted to stabilized finite elements for the Stokes and Navier-Stokes equations. Many of those methods were extended to incompressible elasticity, amongst other approaches by Hughes, Franca, Balestra, and collaborators [39,47]. Masud and coauthors followed an idea by means of variational multiscale (VMS) methods [58][59][60]85], a technique that was recently extended to dynamic problems (D-VMS) [66,71]. Further stabilizations of equal order finite elements include orthogonal sub-scale methods [24,27,30,54] and methods based on pressure projections [33,86]. Different classes of methods to avoid locking for nearly incompressible elasticity were conceived by introducing nonconforming finite elements such as the Crouzeix-Raviart element [32,37] and Discontinuous Galerkin methods [49,80]. Enhanced strain formulations [64,79] have been considered as well as formulations based on multi-field variational principles [17,68,69].
In this study we introduce a novel variant of the MINI element for accurately solving nearly and fully incompressible elasticity problems. The MINI element was originally established for computational fluid dynamics problems [3] and pure tetrahedral meshes and previously used in the large strain regime, e.g. in [25,56]. We extend the MINI element definition for hexahedral meshes by introducing two bubble functions in the element and provide a novel proof of stability and well-posedness in the case of linear elasticity. The support of the bubble functions is restricted to the element and can thus be eliminated from the system using static condensation. This also allows for a straightforward inclusion in combination with existing finite element codes since all required implementations are purely on the element level. Additionally, we introduce a pressure-projection stabilization method originally published for the Stokes equations [14,33] and previously used for large strain nearly incompressible elasticity in the field of particle finite element methods and plasticity [22,65]. Due to its simplicity, this type of stabilization is especially attractive from an implementation point of view.
Robustness and performance of both the MINI element and the pressure-projection approach are verified and compared to standard benchmarks reported previously in literature. A key advantage of the proposed methods is their high versatility: first, they are readily applicable to nearly and fully incompressible solid mechanics; second, with little adjustments the stabilization techniques can be applied to all kinds of finite elements, in this study we investigate the performance for hexahedral and tetrahedral meshes; and third, the methods generalize easily to transient dynamics.
Real world applications often require highly-resolved meshes and thus efficient and massively parallel solution algorithms for the linearized system of equations become an important factor to deal with the resulting computational load. We solve the arising saddle-point systems by using a GMRES method with a block preconditioner based on an algebraic multigrid (AMG) approach. Extending our previous implementations [5] we performed the numerical simulations with the software Cardiac Arrhythmia Research Package (CARP) [82] which relies on the MPI based library PETSc [12] and the incorporated solver suite hypre/BoomerAMG [43]. The combination of these advanced solving algorithms with the proposed stable elements which only rely on linear shape functions proves to be very efficient and renders feasible simulations on grids with high structural detail.
The paper is outlined as follows: Sect. 2 summarizes in brief the background on the methods. In Sect. 3, we introduce the finite element discretization and discuss stability. Subsequently, Sect. 4 documents benchmark problems where our proposed elements are applied and compared to results published in the literature. Finally, Sect. 5 concludes the paper with a discussion of the results and a brief summary.

Nearly incompressible nonlinear elasticity
Let Ω 0 ⊂ R 3 denote the reference configuration and let Ω t ⊂ R 3 denote the current configuration of the domain of interest. Assume that the boundary of Ω 0 is decomposed into ∂Ω 0 = Γ D,0 ∪ Γ N,0 with |Γ D,0 | > 0. Here, Γ D,0 describes the Dirichlet part of the boundary and Γ N,0 describes the Neumann part of the boundary, respectively. Further, let n 0 be the unit outward normal on ∂Ω 0 . The nonlinear mapping φ : X ∈ Ω 0 → x ∈ Ω t , defined by φ := X + u(X, t), with displacement u, maps points in the reference configuration to points in the current configuration. Following standard notation we introduce the deformation gradient F and the Jacobian J as F := Grad φ = I + Grad u, J := det(F), and the left Cauchy-Green tensor as C := F F. Here, Grad(•) denotes the gradient with respect to the reference coordinates X ∈ Ω 0 . The displacement field u is sought as infimizer of the functional over all admissible fields u with u = g D on Γ D,0 , where, Ψ denotes the strain energy function; ρ 0 denotes the material density in reference configuration; f denotes a volumetric body force; g D denotes a given boundary displacement; and h denotes a given surface traction. For ease of presentation it is assumed that ρ 0 is constant and f, g D , and h do not depend on u. Existence of infimizers is, under suitable assumptions, guaranteed by the pioneering works of Ball, see [13].
In this study we consider nearly incompressible materials, meaning that J ≈ 1. A possibility to model this behavior was originally proposed by Flory [38] using a split of the deformation gradient F such that Here, F vol describes the volumetric change while F describes the isochoric change. By setting F vol := J 1 3 I and F := J − 1 3 F we get det(F) = 1 and det(F vol ) = J . Analogously, by setting C vol := J 2 3 I and C := J − 2 3 C, we have C = C vol C. Assuming a hyperelastic material, the Flory split also postulates an additive decomposition of the strain energy function where κ is the bulk modulus. The function U (J ) acts as a penalization of incompressibility and we require that it is strictly convex and twice continuously differentiable. Additionally, a constitutive model for U (J ) should fulfill that (i) it vanishes in the reference configuration and that (ii) an infinite amount of energy is required to shrink the body to a point or to expand it indefinitely, i.e., For the remainder of this work we will focus on functions U (J ) that can be written as In the literature many different choices for the function Θ(J ) are proposed, see, e.g., [34,42,66] for examples and related discussion.
As we also want to study the case of full incompressibility, meaning κ → ∞, we need a reformulation of the system. In this work we will use a perturbed Lagrange-multiplier functional, see [4,21,77] for details, and we introduce We will now seek infimizers (u, p) ∈ V g D ×Q of the modified functional To guarantee that the discretization of (4) is well defined, we assume that with H 1 (Ω 0 ) being the standard Sobolev space consisting of all square integrable functions with square integrable gradient, and Q = L 2 (Ω 0 ). Existence of infimizers in V g D cannot be guaranteed in general. However, assuming suitable growth conditions on the strain energy function Ψ , and assuming that the initial data keeps the material in the hyperelastic range, one can conclude that the space V for the infimizer u contains V g D as a subset, see [13] for details.
To solve for the infimizers of (4) we require to compute the variations of (4) with respect to Δu and Δp with abbreviations as, e.g., in [45] S isc := J − 2 3 Dev(S), where S := ∂Ψ (C) ∂ C (7) Next, with notations we formulate the mixed boundary value problem of nearly incompressible nonlinear elasticity via a nonlinear system of equations. This yields a nonlinear saddle-point problem, find for all (Δv, Δq) ∈ V 0 × Q.

Consistent linearization
To solve the nonlinear variational Eqs. (18) and (19), with a finite element approach we first apply a Newton-Raphson scheme, for details we refer to [31]. Given a nonlinear and continuously differentiable operator F : X → Y a solution to F(x) = 0 can be approximated by which is looped until convergence. In our case, we have . We obtain the following linear saddle-point problem for each (u k , where a k (Δu, Δv) := Ω 0 Grad ΔvS tot,k : Grad Δu dX where C isc is given in (57). The derivation of the consistent linearization is lengthy but standard, we refer to [45,Chapter 8] for details. The definition of the higher order tensor and other abbreviations are given in "Appendix".

Review on solvability of the linearized problem
Since (20) and (21) is a linear saddle-point problem for each given (u k , p k ) we can rely on the well-established Babuška-Brezzi theory, see [15,36,67,70]. The crucial properties to guarantee that the problem (20) and (21) is well-posed are continuity of all involved bilinear forms and the following three conditions: (ii) The coercivity on the kernel condition: there exists c 2 > 0 such that where (iii) Positivity of c: it holds Upon observing that F − : Grad v = div v, see [45], we rewrite the bilinear form b k (q, v) as Assuming that Θ (J ) ≥ 1, we can conclude the inf-sup condition from standard arguments, see [83,Section 5.2]. The positivity of the bilinear form c is always fulfilled. However, it is not possible to show the coercivity condition (25) for a general hyperelastic material or load configuration. Nevertheless, for some special cases it is possible to establish a result. We refer to [7,8,83] for a more detailed discussion. Henceforth, we will assume that our given input data is such that we stay in the range of stability of the problem. Examples for cases in which bilinear form a k lacks coercivity can be found in [83,Chapter 9] and [7, Section 4].

Finite element approximation and stabilization
Let T h be a finite element partitioning of Ω into subdomains, in our case either tetrahedral or convex hexahedral elements. The partitioning is assumed to fulfill standard regularity conditions, see [29]. LetK be the reference element, and for K ∈ T h denote by F K the affine, or trilinear mapping from K onto K . We assume that F K is a bijection. For a tetrahedral element K this can be assured whenever K is non-degenerate, however, for hexahedral elements this may not necessary be the case, see [53] for details. Further, letV andŶ denote two polynomial spaces defined overK . We denote by the spaces needed for further analysis in the following sections.
As a model problem we study the well-known equations for nearly incompressible linear elasticity. In this case it is assumed that Ω := Ω 0 ≈ Ω t . Then, the linear elasticity problem reads: Here, μ > 0 and λ denote the Lamé parameters, and ε(v) := sym(grad v).
The regularity of (31) and (32) is a classical result [75] and follows with the same arguments as for the Stokes equations. The discretized analogue of (31) and (32) Coercivity on the kernel condition (25) is a standard result for the case of nearly incompressible linear elasticity posed in the form (31)- (32) and (33)- (34). In the nonlinear case this is not true in general and will be addressed in Sect. 3.4. The crucial point for checking well-posedness of the discrete Eqs. (33) and (34) is the fulfillment of the discrete inf-sup condition, reading The discrete inf-sup condition puts constraints on the choice of the spaces V h,0 and Q h . A finite element pairing fulfilling (35) is called a stable pair. A classic example for tetrahedral meshes would be the Taylor-Hood element. In this paper, we will focus on two different finite element pairings, the MINI element and a stabilized equal order element. The stabilized equal order pairing has been used in this context for pure tetrahedral meshes, see [22,65]. To the best of the authors knowledge those elements have not been used in the present context for general tesselations.

The pressure-projection stabilized equal order pair
In the following, we present a stabilized lowest equal order finite element pairing, adapted to nonlinear elasticity from the pairing originally introduced by Dohrmann and Bochev [14,33] for the Stokes equations. We chooseV andŶ in (28) and (29) as the space of linear (or trilinear) functions overK . This choice of spaces is a textbook example of an unstable element, however, following [33], we can introduce a stabilized formulation of (33) and (34) and μ * > 0 a suitable parameter. We note that the integral in (38) has to be understood as sum over integrals of elements of the tessellation. The projection operator We can state the following results for this discrete problem: Theorem 1 There exists a unique bounded solution to the discrete problem (36). (31) and (32). Further, assume that (u h , p h ) are the solutions to the stabilized problem (36). Then there exists a constant c 3 independent of the mesh size h and it holds:

Tetrahedral elements
One of the earliest strategies in constructing a stable finite element pairing for discrete saddle-point problems arising from Stokes Equations is the MINI-Element, dating back to the works of Brezzi et al., see for example [3,20]. In the case of Stokes the velocity ansatz space is enriched by suitable polynomial bubble functions. More precisely, if we denote byP 1 the space of polynomials with degree ≤ 1 over the reference tetrahedronK , we will choosê where (ξ 0 , ξ 1 , ξ 2 ) ∈K , see also [15]. Classical results [15] guarantee the stability of the MINI-Element for tetrahedral meshes. Due to compact support of the bubble functions, static condensation can be applied to remove the interior degrees of freedom during assembly. A short review on the static condensation process is given in "Appendix". Hence, these degrees of freedom are not needed to be considered in the full global stiffness matrix assembly which is a key advantage of the MINI element.

Hexahedral meshes
In the literature mostly two dimensional quadrilateral tessellations, see for example [11,15,55], were considered for MINI element discretizations. In this case, the proof of stability relies on the so-called macro-element technique proposed by Stenberg [76].
To motivate our novel ansatz for hexahedral bubble functions, we will first give an overview of Stenbergs main results. A macro-element M is a connected set of elements in T h . Moreover, two macro-elements M 1 and M 2 are said to be equivalent if and only if they can be mapped continuously onto each other. Additionally, for a macro element M we define the spaces Denote by Γ h the set of all edges in T h interior to Ω. The macro-element partition M h of Ω then consists of a (not necessarily disjoint) partitioning into macro-elements The macro element technique is then described by the following theorem, see [76].  We will next show, that Assumption (M1) depends on the choice of the bubble functions inside every K ∈ M i . For ease of presentation and with no loss of generality we will assume that M i is a parallelepiped. This means that the mapping F M i fromK onto M i is affine, so there exists an invertible matrix where ξ ∈K = [−1, 1] 3 and x 0 is a given node of M i . The case of M i not being the image of an affine mapping ofK can be handled analogously, however, there are constraints on the invertibility of F M i , see [53]. Let {ψ j } 8 j=1 denote the standard trilinear basis functions on the unit hexahedron. These functions will serve as a basis for P M i . For the space V 0,M i we will chose one piecewise continuous trilinear ansatz function defined in x i and for each sub-hexahedron we will add two bubble functions as degrees of freedom. The distribution  Fig. 2. OnK we will define the following two bubble functionŝ where the indices α, β are chosen such thatψ α andψ β are two ansatz functions belonging to two diagonally opposite nodes. Having this, we will form a basis for V 0,M i by gluing together the images of the basis functions of each sub-hexahedron. So we can write a basis for V 0,M i as Here, ψ x i corresponds to a piecewise trilinear ansatz function that has unit value in x i and zero in all other nodes of M i . Thus, we can calculate that dim(P M i ) = 27 and dim(V 0,M i ) = 51. For ease of presentation we will rename the elements of (43) as Next, we use the chain rule to get grad x φ k = J − iĝ rad ξ and a change of variables to obtain This means we can find a matrix D ∈ R 27×51 such that where q and v encode the nodal values of q h and v h . The following ordering will be employed for v To proof (M1) we need to show that the rank of the matrix D is 26. Due to the invertibility of J i the rank of the matrix D will remain unchanged by replacing M i byK . Thus, it suffices to compute the rank of the matrix D whose j th row is defined by By this formula the matrix D can be explicitly calculated, e.g., by using software packages like Mathematica TM and further analyzed. We can conclude that the rank of D is 26 and thus (M1) holds and we can apply Theorem 3. A Mathematica TM notebook containing computations discussed in this section is available upon request.

Remark 1
Contrary to the two-dimensional case studied in [11,55] it is not sufficient to enrich the standard isoparametric finite element space for hexahedrons with only one bubble function. In this case both the spaces V 0,M i and P M i have a dimension of 27, however, matrix D has only rank 24.

Remark 2
Although not mentioned explicity, the stability of the MINI element holds also for mixed discretizations.

Changes and limitations in the nonlinear case
One of the main differences between the linear and nonlinear case stems from the definition of the pressure p as remarked in [16]. Consider, as an example, the strain energy function for a nearly incompressible neo-Hookean material where with μ > 0 a material parameter. Then, S tot and C tot , evaluated at (u k , p k ) = (0, 0), are given by independent of the choice of Θ(J ). Assuming that Ω := Ω 0 ≈ Ω t we obtain from Eqs. (20) and (21) the following linear system where ε d (u) := ε(u) − 1 3 div(u)I. While the pressure in formulation (31) and (32) is usually denoted as Herrmann pressure [44], above formulation (44) and (45) uses the socalled hydrostatic pressure.
The arguments to prove the inf-suf condition for this linear problem remains the same as for (31) and (32). For the extension of the inf-suf condition to the nonlinear case we already stated earlier in Eq. (27) Here, Ω t,h is the approximation of the real current configuration Ω t . Our conjecture is that stability of the chosen elements is given provided sufficient fine discretizations and volumetric functions Θ(J ) fulfilling Θ (J ) ≥ 1. However, we can not offer a rigorous proof of this, and rely on our numerical results which showed no sign of numerical instabilities.
Concerning well-posedness of (44)- (45), it was noted in [16], that the coercivity on the kernel condition (25) does not hold in general, which makes the formulation with hydrostatic pressure not well-posed in general. However, it remains well-posed for strictly divergence-free finite elements or pure Dirichlet boundary conditions. This has also been observed by other authors, see [52,81]. Even if the coercivity on the kernel condition can be shown for the hydrostatic, nearly incompressible linear elastic case this result may not transfer to the nonlinear case. Here, this condition is highly dependent on the chosen nonlinear material law and for the presented benchmark examples (Sect. 4) we did not observe any numerical instabilities.
For an in-depth discussion we refer the interested reader to [7,8]. A detailed discussion on Herrmann-type pressure in the nonlinear case is presented in [72,73].
To show well-posedness for the special case of the presented MINI element discretizations we rely on results given in [16,Section 4]. There it is shown, that discrete coercivity on the kernel holds, provided that a rigid body mode is the only function that renders (44) to (45) zero. We could obtain this result following the same procedure outlined in [16] for both hexahedral and tetrahedral MINI elements. A Mathematica TM notebook containing the computations discussed is available upon request.
In the case of the pressure-projection stabilization we will modify Eq. (17) using the stabilization term (38) Here, the stabilization parameter μ * > 0 is supposed to be large enough and will be specifically defined for each nonlinear material considered. Note, that by modifying the definition of the lower residual, we introduced a mesh dependent perturbation of the original residual. An estimate of the consistency error caused by this is not readily available and will be the topic of future research. However, results and comparisons to benchmarks in Sect. 4 suggest that this error is negligible for the considered problems as long as μ * is well-chosen. If not specified otherwise we chose μ * = μ for neo-Hookean materials and μ * = c 1 for Mooney-Rivlin materials in the results section. For the pressure-projection stabilized equal order pair we can not transfer the results from the linear elastic case to the non-linear case, as the proof of wellposedness relies on the coercivity of a k (u, v) which can not be concluded for this formulation. However, no convergence issues occured in the numerical examples given in Sect. 4. The considerable advantage of the MINI element is that there are no modifications needed and that no additional stabilization parameters are introduced into the system.

Changes and limitations in the transient case
The equations presented in Sect. 2 are not yet suitable for transient simulations. To include this feature we modify the nonlinear variational problem (18) in the following way: R trans lower (u, p; Δq) := R lower (u, p; Δq).
For time discretization we considered a generalized-α method, see [28] and also the "Appendix" for a short summary. Due to the selected formulation, the resulting ODE system turns out to be of degenerate hyperbolic type. Hence, we implemented a variant of the generalized-α method as proposed in [50] and using that we did not observe any numerical issues in our simulations. Note, that other groups have proposed a different treatment of the incompressibility constraints in the case of transient problems, see [66,71] for details.

Numerical examples
While benchmark cases presented in this section are fairly simple, mechanical applications often require highly resolved meshes. Thus, efficient and massively parallel solution algorithms for the linearized system of equations become an important factor to deal with the resulting computational load. After discretization, at each Newton-Raphson step a block system of the form has to be solved. In that regard, we used a generalized minimal residual method (GMRES) and efficient preconditioning based on the PCFIELDSPLIT 1 package from the library PETSc [12] and the incorporated solver suite hypre/BoomerAMG [43]. By extending our previous work [5] we implemented the methods in the finite element code Cardiac Arrhythmia Research Package (CARP) [82].

Analytic solution
To verify our implementation we consider a very simple uniaxial tension test, see also [83,Sec. 10.1]. The computational domain is described by one eighth part of a cylinder with

Pressure [MPa]
Hexahedra -MINI Hexahedra -Projection Tetrahedra -MINI Tetrahedra -Projection   Fig. 3. This cylinder is stretched to a length of L + ΔL, with ΔL = 2 mm. We chose a neo-Hookean material with μ = 7.14 MPa and impose full incompressiblity, i.e., 1/κ = 0. For this special case, an analytic solution can be computed by where t ∈ [0, 1] corresponds to the load increment. Two meshes consisting of 5420 points and 4617 hexahedral or 27,702 tetrahedral elements were used. We performed 20 incremental load steps with respect to ΔL. In Fig. 4 it is shown that the results of the numerical simulations render  hexahedral (a, b) and tetrahedral (c, d) elements for the = 2 mesh in Table 1

Block under compression
The computational domain, studied by multiple authors, see, e.g., [23,58,63], consists of a cube loaded by an applied pres-  The same neo-Hookean material model as in [58] is used: with material parameters λ = 400,889.806 MPa, μ = 80.194 MPa. To test mesh convergence the simulations were computed on a series of uniformly refined tetrahedral and hexahedral meshes, see Table 1. Figure 6 shows the deformed meshes for the level = 2 with loads p 0 = 320 MPa and p 0 = 640 MPa, respectively. In all cases discussed in this section we used 10 loading steps to arrive at the target pressure p 0 . As a measure of the compression level the vertical displacement of the node at the center of the top surface, i.e. the edge point A of the quarter of the cube, is plotted in Fig. 7. Small discrepancies can be attributed to differences in the meshes for tetrahedral and hexahedral grids, however, the different stabilization techniques yield almost the same results for finer grids. Note, that the dis- placements at the edge point A obtained using the simple Q 1 − P 0 hexahedral and P 1 − P 0 tetrahedral elements seem to be in a similar range compared to the other approaches. The overall displacement field, however, was totally inaccurate rendering Q 1 − P 0 and P 1 − P 0 elements an inadequate choice for this benchmark problem. The solution for Taylor-Hood (P 2 − P 1 ) tetrahedral elements was obtained using the FEniCS project [2]. Here, as a linear solver, we used a GMRES solver with preconditioning similar to the MINI and projection-based approach, see first paragraph of Sect. 4. The PCFIELDSPLIT and hypre/BoomerAMG settings were slightly adapted to optimize computational performance for quadratic ansatz functions. We comparing simulations with about the same number of degrees of freedom, not accuracy as, e.g., in [25]. For coarser grids computational times were in the same time range for all approaches; see, e.g., the cases with approximately 10 6 degrees of freedom and target pressure of p 0 = 320 mmHg in Table 2(a). For the simulations with the finest grids with approximately 10 7 degrees of freedom, however, we could not find a setting for the Taylor-Hood elements that was competitive to MINI and pressure-projection stabilizations. The computational times to arrive at the target pressure of p 0 = 320 mmHg using 192 cores on ARCHER, UK were about 10 times higher for Taylor-Hood elements using FEniCS, see Table 2(b). We attribute that to a higher communication load and higher memory requirements of the Taylor-Hood elements: memory to store the block stiffness matrices was approximately 2.5 times higher for Taylor-Hood elements compared to MINI and projection-stabilization approaches (measured using the MatGetInfo 2 function provided by PETSc). Note, that although we used the same linear solvers, the time comparisons are not totally just as results were obtained using two different finite element solvers, CARP and FEniCS. Note also, that timings are usually very problem dependent and for this block under compression benchmark high accuracy was already achieved with coarse grids for hexahedral and Taylor-Hood discretizations. For a further analysis regarding computational costs of the MINI element and the pressure-projection stabilization, see Sect. 4.4.
In Fig. 8 the hydrostatic pressure is plotted for the MINI element and the projection-based stabilization. These results are very smooth in all cases and agree well with those published in [23,35,58,63].

Cook-type cantilever problem
In this section, we analyze the same Cook-type cantilever beam problem presented in [17,69], see also Fig. 9. Displacements at the plane x = 0 mm are fixed. At the plane x = 48 mm a parabolic load, which takes its maximum at t 0 = 300 kPa, is applied. Note, that this in-plane shear force Fig. 9 Cook-type cantilever problem: geometry and boundary conditions in y-direction is considered as a dead load in the deformation process. To compare to results in [69] the same isotropic strain energy function was chosen with material properties c 1 = 21 kPa, c 2 = 42 kPa, and γ = 12c 1 + 24c 2 to satisfy the condition of a stress-free reference geometry.
We chose a fully incompressible material, hence, with 1/κ = 0. First, mesh convergence with respect to resulting displacements is analyzed for the tetrahedral and hexahedral meshes with discretization details given in Table 3.

z-displacement at point C [mm]
Hexahedra -MINI Hexahedra -Projection Tetrahedra -MINI Tetrahedra -Projection Displacements u x , u y , and u z at point C are shown in Fig. 10. The proposed stabilization techniques give comparable displacements in all three directions and also match results published in [17,69]. Mesh convergence can also be observed for the stresses σ x x at point A and B and σ yy at point B, see Fig. 11. Again, results match well those presented in [17,69]. Small discrepancies can be attributed to the fully incompressible formulation used in our work and differences in grid construction. Additionally, the minimal and maximal value, as well as the mean (μ) and the standard deviation (σ ) is given for each setting cantilever remained unchanged at 14,400 mm 3 , rendering the material fully incompressible on the domain level. Figure 13 gives a comparison of several computed values in the deformed configuration of Cook's cantilever for the finest grids ( = 4). Slight pressure oscillations in Fig. 13b on the domain boundary for the MINI element are to be expected, see [74]; this also affects the distribution of J in Fig. 13a. A similar checkerboard pattern is present for the projection based stabilization.
In the third row of Fig. 13 we compare the stresses σ x x for the different stabilization techniques. We can observe slight oscillations for the the projection-based approach, whereas  Figure 10] the σ x x stresses have a similar contour but are slightly higher. As before, we attribute that to the fully incompressible formulation in our paper compared to the quasi-incompressible formulation in [69].

Twisting column test
Finally, we show the applicability of our stabilization techniques for the transient problem of a twisting column [1,40,71]. The initial configuration of the geometry is depicted in Fig. 14. There is no load prescribed and the column is restrained against motion at its base. A twisting motion is applied to the domain by means of the following initial condition on the velocity v(x, 0) = v(x, y, z, 0) = 100 sin π y 12 (z, 0, −x) m/s, for y ∈ [0, 6] m. To avoid symmetries in the problem the column is rotated about the z-axes by an angle of θ = 5.2 • . We chose the neo-Hookean strain-energy with parameters μ = 5704.7 kPa and κ = 283,333 kPa for the nearly incompressible and 1/κ = 0 for the truly incompressible case. For the results presented, we considered hexahedral and tetrahedral meshes with five levels of refinement, respectively; for discretization details of the column meshes see Table 4. In Fig. 15, mesh convergence with respect to tip displacement (u x , u y , u z ) at point D is analyzed. While differences at lower levels of refinement = 1, 2 are severe, the displacements converge for higher levels of refinement = 3, 4, 5. For finer grids the curves for tetrahedral and hexahedral elements are almost indistinguishable, see also Fig. 16, and the results are in good agreement with those presented in [71]. While this figure was produced using MINI elements we also observed a similar behavior of mesh convergence for the projection-based stabilization. In fact, for the finest grid, all the proposed stabilization techniques and elements gave virtually identical results, see Fig. 16. Further, as already observed by Scovazzi et al. [71], the fully and nearly incompressible formulations gave almost identical deformations, see Fig. 17.
In Fig. 18 stress σ yy and pressure p contours are plotted on the deformed configuration for the incompressible case at time instant t = 0.3 s. Minor pressure oscillations can be observed for tetrahedral elements. Again, results match well those presented in [71, Figure 22].
Finally, in Fig. 19, we compare the magnitude of velocity and acceleration at time instant t = 0.3 s. Results for these variables are very smooth and hardly distinguishable for all the different approaches.
The computational costs for this nonlinear elasticity problem were significant due to the required solution of a saddle-point problem in each Newton step and a large number of time steps. However, this challenge can be addressed by using a massively parallel iterative solving method and exploiting potential of modern HPC hardware. The most expensive simulations were the fully incompressible cases for the finest grids with a total of 840,708 degrees of freedom and 400 time steps. These computations were executed at the national HPC computing facility ARCHER in the United Kingdom using 96 cores. Computational times were as follows: 239 min for tetrahedral meshes and projection-based stabilization; 283 min for tetrahedral meshes and MINI elements; 449 min for hexahedral meshes and projection-based stabilization; and 752.5 min for hexahedral meshes and MINI elements. Simulation times for nearly incompressible problems were lower, ranging from 177 to 492 min. This is due to the additional matrix on the lower-right side of the block stiffness matrix which led to a smaller number of linear iterations. Simulations with hexahedral meshes were, in general, computationally more expensive compared to simulations with tetrahedral grids; the reason being mainly a higher number of linear iterations. Computational burden for MINI elements was larger due to higher matrix assembly times. However, this assembly time is highly scalable as there is almost no communication cost involved in this process.

Conclusion
In this study we described methodology for modeling nearly and fully incompressible solid mechanics for a large variety of different scenarios. A stable MINI element was presented which can serve as an excellent choice for applied problems where the use of higher order element types is not desired, e.g., due to fitting accuracy of the problem domain. We also proposed an easily implementable and computationally cheap technique based on a local pressure projection. Both approaches can be applied to stationary as well as transient problems without modifications and perform excellent with both hexahedral and tetrahedral grids. Both approaches allow a straightforward inclusion in combination with existing finite element codes since all required implementations are purely on the element level and are well-suited for simple single-core simulations as well as HPC computing. Numerical results demonstrate the robustness of the formulations, exhibiting a great accuracy for selected benchmark problems from the literature.
While the proposed projection method works well for relatively stiff materials as considered in this paper, the setting of the parameter μ * has to be adjusted for soft materials such as biological tissues. A further limitation is that both formulations render the need of solving a block system, which is computationally more demanding and suitable preconditioning is not trivial. However, the MINI element approach can be used without further tweaking of artificial stabilization coefficients and preliminary results suggested robustness, even for very soft materials. Consistent linearization as presented ensures that quadratic convergence of the Newton-Raphson algorithm was achieved for all the problems considered. Note that all computations for forming the tangent matrices and also the right hand side residual vectors are kept local to each element. This benefits scaling properties of parallel codes and also enables seamless implementation in standard finite element software.
The excellent performance of the methods along with their high versatility ensure that this framework serves as a solid platform for simulating nearly and fully incompressible phenomena in stationary and transient solid mechanics. In future studies, we plan to extend the formulation to anisotropic materials with stiff fibers as they appear for example in the simulation of cardiac tissue and arterial walls.
In all our simulations we used a value of ρ ∞ = 0.5.

Remark on the implementation of the pressure-projection stabilized equal order pair
Considering the bilinear form s h ( p h , q h ) defined in (38) we can rewrite this with a simple calculation into Denoting by {φ i } n i=1 the chosen ansatz functions the element contribution for an arbitrary element K to the matrix C h is given by This corresponds to an element mass matrix minus a rank-one correction.

Static condensation
For completeness we provide a summary for the static condensation used for the MINI element. Consider a finite element K ∈ T h with a local ordering of the unknowns u and p as p = p 1 , p 2 , . . . , p ndofs N .
Here, ndofs N corresponds to the nodal degrees of freedom per element and ndofs B to the bubble degrees of freedom (one for tetrahedral elements and two for hexahedral elements). Then the element contribution to the global saddle-point system can be written as The bubble part of the stiffness matrix, K BB is local to the element and can be directly inverted. This gives the condensed system where the effective matrices and vectors are given as The effective matrices and vectors can then be assembled in a standard way into the global system. The bubble update contributions can be calculated once Δu N and Δp N are know as

Tensor calculus
We use the following results from tensor calculus, for more details we refer to, e.g., [45,84]. The isochoric part of the second Piola-Kirchhoff stress tensor as well as the isochoric part of the fourth order elasticity tensor are given as