An extension of assumed stress finite elements to a general hyperelastic framework

Assumed stress finite elements are known for their extraordinary good performance in the framework of linear elasticity. In this contribution we propose a mixed variational formulation of the Hellinger–Reissner type for hyperelasticity. A family of hexahedral shaped elements is considered with a classical trilinear interpolation of the displacements and different piecewise discontinuous interpolation schemes for the stresses. The performance and stability of the new elements are investigated and demonstrated by the analysis of several benchmark problems. In addition the results are compared to well known enhanced assumed strain elements.


Introduction
An enormous effort was invested in the development of finite element methods based on the variational approach going back to Galerkin [28], which are in general considering the approximation of one basic variable. It is well known that low order elements, based on this variational approach yield poor results in the framework of nearly incompressibility and bending dominated problems. The reason for this are locking effects, see [12]. Developments considering the approximation of additional fields in the variational setup are for example given in Reissner [44] (compare also [29] and [43]), here an independent stress approximation is applied in addition to the displacements, which acts as a Lagrange multiplier. We refer to this type of formulation, which is based on a complementary stored energy function, as Hellinger-Reissner (HR) formulation. A few years later, [30] and [57] proposed independently a variational principle related to displacements, stresses and strains, based on the so-called Hu-Washizu functional. In the framework of finite element analysis, mixed formulations lead to saddle-point problems and therefore a major constraint of these approaches is the restriction to the so-called LBB-conditions, see [11,19,20,35] and [10]. Mathematical aspects concerning the mixed finite element formulation for elasticity based on the Hellinger-Reissner (HR) principle are discussed in Arnold and Winther [4], Auricchio et al. [10], Lonsing and Verfürth [37], Arnold et al. [3], Boffi et al. [18] and Cockburn et al. [25]. In contrast to the Stokes problem, where the enrichment of the discrete space of the primal variable leads always to stable elements in the linear range, the situation is more difficile in the framework of the HR formulation. Here, the discrete spaces have to be carefully balanced.
Within the principle of HR two methods can be distinguished. A shift of the derivatives from the displacements to the stresses, using integration by parts, leads to a formulation with the displacements in the L 2 and the stresses in the H(Div) Sobolev spaces. In the literature, this is often denoted as the dual HR formulation. Stable elements in the linear range based on the dual HR formulation with a strong enforcement of the symmetry can only be achieved at high polynomial order, see e.g. [6,31] in 2D and [1,7] in 3D. A reduction of the symmetry constraint leads to some flexibility in the construction of the finite elements, see [5,18,53,54] and [32].
In the primal HR formulation the derivatives are correlated to the displacements, such that the displacements are in the H 1 and the stresses in L 2 Sobolev spaces. Elements, based on this formulation, are called assumed stress elements going back to the pioneering work of Pian [39]. In 2D, often a discontinuous stress approximation using a 5-parameter ansatz proposed by Pian and Sumihara [40] is used. The advantages of this approach are characterized by a remarkable insensitivity to mesh distortion, locking free behavior for plane strain quasi-incompressible elasticity and superconvergent results for bending dominated problems, see e.g. [40,52] and [23]. Its stability with regard to the LBB-conditions and an a posteriori error estimation has been shown by Yu et al. [61] and Li et al. [36]. A well known extension to the 3D case is the element by Pian and Tong [41]. The family of assumed stress elements is closely related to the family for enhanced assumed strains elements, see e.g. [2,60] and [16].
The "direct" extension of the HR principle, which requires a complementary energy function in terms of stress measures, is in general not possible. Some effort has been pursued in the extension of the dual version of the HR principle considering small-strain elasto-plasticity, [46,48]. Considering a large-strain setup the workgroup of Atluri (see [8,9,49]), has been proposed an incremental variational formulation, involving the discretization of the displacement, rotations and the hydrostatic pressure considering the primal HR formulation. In addition [42] extended the underlying variational formulation to a Hu-Whashizu type form, closely related to the family of enhanced assumed strains. Another approach has been discussed in Wriggers [58], where the complementary constitutive relation has been derived in an explicit form for a special Neo-Hookean type.
In the proposed work, a general framework is introduced, which extends the family of assumed stress elements to hyperelasticity. It is based on an iterative solution procedure of the constitutive law at the level of the element integration points. In this work, the framework is adopted to a family of 3D hexahedral elements, with a varying interpolation scheme of the stresses. It should be noted that in the framework of hyperelasticity, the corresponding solution spaces are shifted from the classical Lebesgue spaces to the more general framework of the Hilbert spaces, see [24] for a detailed discussion. Thus, the related solution spaces for the proposed formulations are represented for the displacements and stresses by W 1,p and W 0,p , whereas p ≥ 2 depends on the constitutive relation. The stress interpolation is discussed with respect to volumetric locking, shear locking and hourglassing. In addition the related EAS elements are mentioned and the properties and characteristics of the various elements are discussed on the example of different numerical benchmarks.

Kinematics and variational formulation
Let B ⊂ IR 3 be the body of interest in the reference placement, parametrized in X. Its boundary ∂B is decomposed into a nontrivial Dirichlet part ∂B u and a Neumann part ∂B t with ∂B u ∪ ∂B t = ∂B and ∂B u ∩ ∂B t = ∅. The nonlinear deformation map ϕ, which maps points X ∈ B onto points x of the actual placement, is given by x = ϕ(X). As basic kinematical quantities we define the deformation gradient, the right Cauchy-Green tensor, and the Green-Lagrange strain tensor respectively. Here, I denotes the second-order identity tensor. The Jacobian of the deformation gradient has to satisfy J := det F > 0. In addition we introduce the symmetric second Piola-Kirchhoff stresses S as an additional independent variable. The relevant continuum-mechanical quantities are listed in Table 1. Before we extend the discrete formulation to a general hyperelastic framework, we first assume for simplicity here the existence a complementary stored energy function χ(S) which describes the constitutive equation by E.g. in the case of St. Venant type nonlinear elasticity E = C −1 : S, we simply obtain the explicit expression where the compliance tensor C −1 is defined as the inverse of the fourth-order elasticity tensor. In case of isotropy the compliance tensor is given by with the dyadic product defined as (I ⊗ I) ijkl = δ ij δ kl , the fourth-order identity tensor II ijkl = 1 2 (δ il δ jk + δ ik δ jl ) and the Lamé constants and μ. It should be mentioned that in general an explicit complementary stored energy does not exist for arbitrary hyperelastic constitutive laws. This issue will be discussed in the next chapter. The balance of momentum closes, under the assumption of suitable boundary conditions, the set of equations for the boundary value problem in hyperelasticity where f denotes the body force vector and Div the divergence operator with respect to X. For a more detailed overview on the underlying continuum mechanics, the reader is referred e.g. to Ciarlet [24]. The solution of Eqs. (2) and (5) for the displacements u and the second Piola-Kirchhoff stresses S are equivalent to the stationary point of the HR functional with the external potential ext (u) given by where t denotes the prescribed traction vector on the Neumann boundary and u satisfies a priori the Dirichlet boundary conditions. In order to find a stationary point of the functional, the roots of the first variations with respect to the unknown fields u and S have to be calculated. In detail, we obtain with the virtual deformation δu and the virtual stress field δS. Furthermore, the virtual strains are defined by δE = 1 2 (δF T F + F T δF ) with δF = ∇δu.

Discrete formulation and interpolation
We consider a standard decomposition of the reference body B into an assemblage of hexahedral elements B e such that B ≈ B h = n ele e=1 B e , where n ele is the number of elements. In the following, underlined quantities denote the application of a matrixnotation, which reduces second-and fourth-order tensors to suitable matrices using the Voigt notation. The displacements and its variations are accordingly approximated by where IN denotes a matrix containing the Lagrangian trilinear shape functions and d the vector of elementwise nodal displacements. This leads to a displacement interpolation which is C 0 -continuous on B h . The interpolation of the assumed stress field on the isoparametric reference element is given in vector notation by where β is the vector of element-wise unknowns and L ξ the matrix with the corresponding interpolation functions, with the general structure A variety of suitable interpolation matrices are discussed in the following section. The transformation from the isoparametric domain to the reference configuration for the second Piola-Kirchhoff stresses is described by where the Jacobian matrix, mapping between the isoparametric coordinates ξ and the reference coordinates X follows as In order to pass the patch test, see Figs. 2 and 3, it is necessary to use the values of the Jacobian at the origin {ξ , η, ζ } = {0, 0, 0} as it is discussed in Pian and Sumihara [40] and Pian and Tong [41]. Therefore, the second Piola-Kirchhoff stress in the physical space is given by S = L β , with L = T L ξ where the transformation matrix T is given by The discretized weak forms appear for a typical element as where we have utilized δE = B δd with B being a suitable matrix containing the derivatives of the shape functions. In general cases, where complementary stored energy function is not known, the partial derivative ∂ S χ(S) has to be computed iteratively. The Green-Lagrange strain tensor is given by the approximation of the geometry whereas ∂ S χ(S) =: E cons is implicitly evaluated from the constitutive equation. Therefore, we compute at each integration point E cons from the residual at fixed S. Using a Newton scheme we have to update until r(E cons ; S) ≈ 0. Table 2 sketches the nested algorithmic treatment for a typical element for the case that the complementary stored energy cannot be computed analytically. The linearization, LinG e = G e (d, β) + G e ( d, β), yields the increments (2) Compute stresses S and Green-Lagrange strain tensor E at each Gauss Point: where is defined as B = d. We introduce for convenience the element matrices and right hand side vectors This leads to the system of equations Assembling over the number of elements n ele leads to the global system of equations and therefore the nodal unknowns are computed via Due to the elementwise discontinuous interpolation of the stresses, the unknowns β in (20) can already be eliminated at element level. This leads to a global system of equations with the same number of unknowns, and almost the same computational cost, as a displacement based trilinear element.

Stress interpolation
Following the limitation principle by Fraeijs de Veubeke [27], it can be shown that the following 39 parameter based interpolation is in the linear elastic framework equivalent to a primal displacement formulation with a trilinear interpolation for the displacements because the resulting stress space are equivalent in both elements. For this case the interpolation vectors for the stresses are given by It will be shown in the numerical examples, that this formulation suffers due to volumetric locking, which is caused by the choice of parameter corresponding to L ξξ , L ηη and L ζ ζ .
Here, a meaningful reduction will suppress the artificial stiffness. In addition, it will be shown that this interpolation also leads to shear locking, which can be attributed to the choice of L ξη , L ηζ and L ξζ . An expedient reduction of introduced stress-modes, lead to a softening of the formulation. However, care must be taken not to relax the formulation too much, since this could lead to artificial deformation states such as hourglassing. A well known and very efficient stress discretization for the linear elastic counterpart of the proposed HR formulation is the 18 parameter based interpolation scheme proposed by Pian and Tong [41], which is a 3D extension of the element by Pian and Sumihara [40]. It has been shown by Andelfinger and Ramm [2] that this interpolation leads in the small deformation framework to an equivalent EAS formulation with 21 additional modes. Here the individual interpolation vectors are given by Note, that this stress interpolation considers a minimal number of stress unknowns, since any further reduction would lead to singular system matrices due to a violation of the count condition proposed by Zienkiewicz et al. [62]. In the framework of the primal HR formulation this count condition requires that the number of unknowns related to the stresses is always larger or equal to the number of unknowns related to the displacements. Note that in this case a single element constitutes the most critical patch, due to the assembling procedure and the related tying of the displacement related degrees of freedom.
The numerical examples will show that this formulation does not suffer due to volumetric and shear locking. But in case of large uniaxial homogeneous stress states hourglassing modes can be detected. It is known (see e.g. [34] and the references therein) that the hourglass modes are caused by a violation of the kinematics related to the shear deformation terms. The EAS formulation with 21 enhanced modes depict the same characteristics in the hyperelastic framework. In a straight forward manner, further discretization schemes can be adopted from EAS formulations. A promising approach has been published recently by Krischok and Linder [34] where a 9 parameter based EAS element was proposed, which does not suffer to hourglassing modes and is free from volumetric locking. The related stress interpolation in the framework of the HR formulation needs the introduction of a 30-parameter based discretization It can be recognized that the only difference between the interpolation schemes of Eqs. (24) and (25) are the additional terms for the interpolation of the shear stresses. These additional stress modes lead to a stabilization such that hourglassing modes have not been detected in the numerical examples. Unfortunately they also lead to shear locking phenomena as it can be expected by the comparison L ξη − L ξζ of Eqs. (25) and (23).
With this knowledge, it makes sense to investigate one further interpolation scheme. Since the volumetric locking problem seems to be solved by a suitable interpolation of S 11 , S 22 and S 33 , only the interpolation matrices corresponding to the shear stresses are modified. In particular we consider the interpolation scheme which can be nested in between the interpolations from Eqs. (24) and (25). The corresponding matrices follow by

Numerical simulations
In the following numerical examples the proposed family of finite elements will be investigated with respect to locking phenomena, stability, robustness and efficiency. Each considered benchmark focuses on one of these aspects in order to be able to classify each element.
In addition to the proposed assumed stress elements the numerical studies will also be discussed for a non-mixed lowest order element, the well known H1P0 (see [50]) and a couple of enhanced assumed strain element formulations. In particular we consider exactly those elements which are, in the framework of linear elasticity and a constant Jacobian, equivalent to the considered assumed stress elements. A list of the considered elements is given Enhanced assumed strain element with 15 modes, see [38] H1P0 Displacement-pressure approach with piecewise constant pressure, see [50] Boundary Conditions:  Table 3. Please note that the *EAS-9 element is not the popular enhanced assumed strain element by Simo and Rifai [51], which also comprises 9 enhanced modes.

Equivalence of the formulations in linear elasticity
All considered assumed stress element formulations have an equivalent counterpart in the linear elastic framework if the Jacobian mapping from the isoparametric to the physical domain is piecewise constant. This equivalency has been already discussed and proven in a variety of publications, see e.g. [16]. Here, a numerical test is used in order to verify the implementation of the finite element code. Therefore we consider a Neo Hookean type free energy function of the form which fulfills = C, where C is the elasticity tensor related to Hooke's law. In order to obtain the result referring to the linear elastic framework we consider the result after the first Newton-iteration. The considered boundary value problem is a block with the dimension 100×100×50 which is subjected to a constant surface load at the top central quarter t = (0, 0, −3) T . Due to the axial symmetry of the problem, only a quarter of the block is discretized. The geometry, a representative mesh and the boundary conditions are given in Fig. 1. The material parameters are assumed to be compressible, with a Young's modulus E = 5 and a Poisson's ratio of ν = 0.3. In terms of the Lamé parameter this is related to μ = 1.92308 and = 2.88462. A regular hexahedral mesh with uniform refinement is considered, such that the Jacobian is always constant in each element. Table 4 depicts the convergence of the nodal u 3 displacements at point x = (0, 0, 50). The

Patch test
The Patch test is a necessary condition for the convergence of finite elements. It demands that an arbitrary patch of assembled elements is able to reproduce a constant state of stress and strain if subjected to boundary displacements consistent with constant straining. This condition is necessary since with respect to mesh refinement, where h → 0, all boundary value problems tend to constant stress and strains in each element. This test is mainly attributed to the work of Bruce Iron, first presented in Bazeley et al. [15]. A summary on its theory, practice and possible conclusions on its satisfaction can be found in Taylor et al. [55]. Following Korelc et al. [33], two different load scenarios are considered, described in Figs. 2 and 3. Load case (A) prescribes a pure rigid body motion by a rotation around the zaxis. All proposed elements are free of resulting stresses and strains and thus fulfill the first patch test. Load case (B) prescribes a combined deformation of shear and uniaxial strain, which analytical solution leads to a constant strain and stress over the whole domain. All proposed elements result in the expected constant stress and strain field and therefore fulfill the second patch test. Note that the patch tests verify only the consistency of the finite element and thus represents only a necessary but not a sufficient condition for the stability of the formulation.

Inf-sup test
The theorems of Babuška [11] and Brezzi [20] ensure stability, existence and uniqueness of the solution of mixed finite elements in the case of linear elasticity. The crucial aspect in the validation of these theorems for particular formulations is mainly encountered in the proof of the discrete inf-sup condition. In many cases its analytical proof is cumbersome and a direct numerical evaluation is impossible because an infinite number of problems must be taken into account. For the engineering praxis, a numerical inf-sup test has been proposed in Chapelle and Bathe [22] and Bathe [13,14]. The objective of this test is to estimate the inf-sup constant numerically on a set of refined meshes. In the following this inf-sup test is adopted for the proposed Hellinger-Reissner formulations considering a stress free and undeformed (F = I) configuration which represents the framework of linear elasticity. A simply supported unit cube with the dimension 1 × 1 × 1 is considered with a consecutive regular mesh refinement. It has been shown by Brezzi and Fortin [21] that the square root of the smallest eigenvalue, denoted by λ p , of the problem is equivalent to the inf-sup constant, where the matrices T and A are defined as In addition to the evaluation of the inf-sup values for the proposed assumed stress elements the corresponding test is verified by means of the H1P0. This element is a well known textbook example which does not satisfy the inf-sup condition. A detailed instruction on the accomplishment of the inf-sup test for displacement-pressure based elements is given in Chapelle and Bathe [22]. Figure 4 shows the development of the inf-sup value for the chosen boundary value problem over the number of elements considering a regular structured mesh refinement. The value seems to be bounded from below for the AS-39, AS-30 and AS-24. In case of the AS-18 the inf-sup value shows asymptotical convergence, which also indicates a distinct lower bound and thus passes the inf-sup test. Both described behavior are in sharp contrast to the results of the H1P0. Here, the inf-sup value decreases with respect to mesh refinement by an almost constant rate, clearly indicating the failure of the inf-sup test. It should be mentioned, that such a numerical verification does not replace the need of an analytical investigation in order to ensure the statements on the stability. Nonetheless, to the best knowledge of the authors, not a single example has been found, where the inf-sup test does not give similar results as the analytical proofs.

Eigenvalue analysis of initial element matrix
The eigenvalue spectrum of a single element matrix in the initial state is inspected for the case of nearly incompressibility and in addition for a slender domain. Note that in the proposed setup the initial state represents the case of linear elasticity. For the nearly incompressible example, we follow the example of Andelfinger and Ramm [2] and consider a single square element with a side length of l = 1, a Young's modulus of E = 1 and a Poisson's ratio of ν = 0.49999. Table 5 shows the eigenvalue spectrum of the different elements, whereas eigenvalues larger than 10 3 are denoted by ∞ and the six zero eigenvalues are omitted. The incompressibility should affect only a single eigenvalue to tend to infinity which is related to the mode of volumetric dilation. It can be recognized, that in case of the H1 and the AS-39 a total number of seven eigenvalues tend to infinity which consequences in volumetric locking. All other proposed elements depict the correct eigenvalue spectrum. A second study on the eigenvalue spectrum is related to an increasing slenderness of the domain of interest. In Fig. 5 the development of the smallest eigenvalue of a single element cube with the dimensions 2 × 2 × t is investigated. In addition the domain is assumed to be clamped at x = 0. The Domain on the left of Fig. 5 depicts the related eigenmode, which is clearly the mode related to bending deformation. The evolution of the smallest eigenvalue of the element matrix shows two distinct characteristics. In case of the EAS-15, EAS-21, AS-24 and AS-18 the eigenvalue is decreasing rapidly for a shrinking thickness t. In contrast this progress is significant slower for the elements H1, AS-39, *EAS-9 and H1P0. Based on this observation, it can be recognized that bending deformation states are energetically preferable for the first group of elements, yielding to a locking behavior for the latter group. The same analysis is carried for a skew-shaped domain depicted in Fig. 6. It is interesting to mention, that for this domain the elements EAS-15 and AS-24 show a reduced ratio of decrease after a critical value of t. In contrast, the behavior of the remaining elements is almost unaffected.

Hyperelastic nearly incompressible Cook's membrane
A tapered cantilever beam known as the Cook's membrane problem is considered, representing a bulk related boundary value problem with nearly incompressible material. The geometry, boundary conditions and material parameter are summarized in Fig. 7. The cantilever, with a thickness of t = 10 is clamped on the left and a constant shear stress is applied on the right face. In addition the displacements in out of place direction are restricted, imitating a plane strain setting. The material parameters are chosen to be nearly incompressible with a Young's modulus of E = 200 and a Poisson's ratio of ν = 0.499. This corresponds to the Lamé parameters given by = 33288.9 and μ = 66.711. The convergence of the tip displacement for a regular mesh refinement in x and y direction are shown in Fig. 8a. Note that only a single element in z direction is considered. The suffering due to volumetric locking can be recognized for the displacement based elements H1. Unfortunately, the AS-39 does not converge at all for this numerical example. In contrast, all other mixed finite elements achieve a comparably good convergence behavior of the tip displacements, since they do not suffer due to volumetric locking, as expected from the  Table 5. Note that the elements which are equivalent in the linear elastic framework do not yield exactly the identical solution for the displacements in the finite deformation case. However their results are still close to each other. Consideration of the necessary load steps, depicted in Fig. 8b, indicate that the assumed stress elements are able to deal with large load steps in case of nearly incompressibility. For this boundary value problem the considered elements require only a single load step, independent of the mesh size. In contrast, the enhanced assumed strain elements and the H1P0 element require a significant higher number of load increments, leading to a substantially larger computation time. The level of necessary load steps in case of the H1 element is moderate but it should kept in mind that their performance by means of displacement accuracy is insufficient.

Hyperelastic bending of a compressible clamped plate
We consider a thin rectangular plate with the dimension 10 × 10 × 0.1 which is clamped at one end and loaded by a bending force at the opposing end. The geometry, boundary conditions and a representative mesh are depicted in Fig. 9. Two different meshing strategies are considered in the following, note that both are restricted to a single element in thickness direction. On the one hand a regular and structured meshing procedure is adopted, where the number of elements in x and y direction coincide. A detailed visualization of the appropriated meshes is neglected for the sake of brevity. In addition an unstructured inplane mesh is adopted, whereas the considered refinement steps are explicitly depicted in Fig. 10. In a first step a perfect compressible material, characterized by the Neo Hookean energy, see (27), and a Young's modulus of E = 200 and a Poisson's ratio of ν = 0 is taken into account. In terms of the Lamé parameter this is related to μ = 100 and = 0. Due to the material parameter and the boundary conditions, this boundary value problem considers a pure bending mode. The convergence of the displacements with respect to both mesh refinement strategies are depicted in Fig. 11. The elements can be mainly distinguished into three groups, whereas the close relationship between the assumed stress and enhanced assumed strain elements can be recognized again. In case of the structured meshing the AS-18, EAS-21, AS-24 yield EAS-15 yield optimal results for the tip-displacement already for the coarsest meshing. However, taking into account the unstructured meshes the quality of the results is weakened, especially for the coarsest level. However, the AS-18 and EAS-21 depict the best mesh convergence of all considered elements, independent of the discretization strategy. The AS-24 and EAS-15 are slightly weaker in case of unstructured coarse meshes. In contrast to these four elements the H1P0, H1, AS-39, AS-30 and *EAS-9 demonstrate a distinct locking behavior in this numerical example. Note that these observations are consent with the results of a b Fig. 11 Bending of a compressible clamped plate; convergence of tip displacement and number of necessary load steps over the number of elements for a structured (a) and an unstructured (b) mesh refinement strategy the eigenvalue study related to Fig. 6, where locking behavior has been predicted for the latter elements in case of boundary value problems with slender domains. Interestingly, their tip-displacement convergence is slightly better for the unstructured case, which is in contrast to the behavior of the remaining elements.
In addition the number of necessary load steps are depicted on the right of Fig. 11. It can be recognized that also in this bending dominated problem, the number of necessary load steps is significantly smaller for the proposed family of assumed stress elements. Especially the enhanced assumed strain elements, which do not show locking effects (EAS-15, EAS-21) suffer due to the need of smaller load increments compared to their AS counterparts (AS-18, EAS-24).

Hyperelastic bending of a clamped plate-nearly incompressible
We investigate the boundary value problem with the same geometry, boundary conditions and constitutive relation as in the prior example. Only, the material is assumed to be nearly incompressible increasing the Poisson's ratio to ν = 0.495. In terms of the Lamé parameter this is related to μ = 66.8896 and = 6622.07.
Considering the displacement convergence, shown on the left of Fig. 12, it can be noted that the qualitatively response of the elements is equivalent to the compressible case, except for the AS-39 and H1. These elements clearly depict a conspicuous additional volumetric locking. A similar picture as for the compressible case is obtained considering the number of necessary load steps, shown in Fig. 12b. Even if the displacement results a b Fig. 12 Bending of a nearly incompressible plate; convergence of tip displacement and number of necessary load steps over the number of elements for a structured (a) and an unstructured (b) mesh refinement strategy for the EAS-21 and EAS-15 seem to be satisfactory, they suffer due to the need of smaller load steps (a factor of 20), compared to their assumed stress counterparts.

Hyperelastic hourglassing test
The deficiency of a couple of enhanced assumed strain elements to artificial modes, often denoted as hourglassing, have first been detected in case of homogeneous stress states by Wriggers and Reese [59]. The investigation of these unphysical free energy modes have been subject of many following publications, see e.g. [26] and [56]. In order to investigate the assumed stress elements with regard to hourglassing modes we consider the boundary value problem depicted in Fig. 13.
The numerical test considers a cube with a side length of 50 where a compression is applied by displacement boundary conditions. In addition displacement boundary conditions are considered such that the edges are fixed to be straight. It is an easy task to see, that for this problem surface buckling is excluded due to the choice of boundary conditions and therefore the solution for the displacements has to be unique. Therefore each point of instability is introduced by the numerical discretization scheme and can be considered as artificial. The same constitutive law as in the prior examples is used, whereas the material parameter are given by E = 4.337 and ν = 0.355. In terms of the Lamé parameter this is related to μ = 1.60037 and = 3.91814. Figure 14 show  formulations depict stable behavior. Note that the unstable elements are exactly the shear locking free elements. Interestingly, the load levels where the instability occurs is not close to each other for the formulations which are equivalent in the linear elastic framework. In case of the assumed stress discretization the point of hourglassing is on a higher state of stress, compared to its enhanced assumed strain counterpart. The relationship between these elements occurs again in the related eigenforms see Fig. 15. It can be recognized, that the hourglassing modes itself are equivalent for the AS-18 and EAS-21 and also for the AS-24 and EAS-15.

Hyperelastic fiber reinforced Cook's membrane problem
In order to show that the proposed algorithm is also applicable to more complex constitutive equations the Cook's membrane with a fiber reinforced material is considered. Therefore, the underlying strain energy function is split into an isotropic and an anisotropic part

Boundary
Conditions: where α, β, γ , 1 and 2 are material parameter and the Cofactor of a second order tensor is defined by CofA = det[A] A −T . For the formulation of anisotropic free energies as isotropic tensor functions we apply the concept of structural tensors, see e.g. [17]. Considering here the case of transverse isotropy we introduce a preferred direction vector a of unit length and the structural tensor M = a ⊗ a. The anisotropic part of the strain energy, originally introduced in Schröder et al. [47], is given by with the material parameter g 0 , g c , g h , g θ , see in this context also [45]. The geometry, boundary conditions and the material parameter are depicted in Fig. 16. The convergence of the displacements at node x = (48, 60, 0) T and the corresponding necessary load steps are depicted in Fig. 17. First, it should be mentioned that for the AS-39 we have not been able to obtain a solution for the considered boundary value problem. In contrast to that, the previously mentioned statements are confirmed. The AS-30, *EAS-9 and the H1 suffer due to the bending of the boundary value problem. In addition, the H1 suffers due to the constraint on the volumetric deformation, corresponding to the material parameter

Concluding remarks
We proposed an extension of a mixed finite element formulation based on the Hellinger-Reissner principle to the framework of hyperelasticity and investigated it numerically. This principle requires an independent interpolation of the stresses and displacements. The displacements are interpolated by classical trilinear shape functions, whereas for the stresses four different interpolation strategies are discussed. The numerical results for each interpolation can be summarized by the following. The AS-39 correlates to the H1 element also in the nonlinear regime and shows volumetric as well as shear locking. The AS-18 represents the hyperelastic extension of the element proposed by Pian and Tong [41]. A close relationship can be drawn to the EAS-21. The AS-18 performs very well in bending dominated and nearly incompressible problems. Unfortunately, it suffers due to hourglassing modes which could lead to instabilities also in more complex boundary value problems. The AS-30 element is inspired by the *EAS-9 enhanced formulation, proposed in Krischok and Linder [34]. It is free of volumetric locking, does not show hourglass instabilities in the discussed numerical example. Unfortunately, it is not free of shear locking and therefore behaves poor in the bending dominated problems. The AS-24 element is closely related to the EAS-15 proposed by Pantuso and Bathe [38]. It is free of volumetric and shear locking. However, it suffers due to hourglass instabilities. Even if non of the investigated interpolation schemes are free of all drawbacks, the proposed procedure of Hellinger-Reissner formulations for large deformations has emerged as a promising approach. The proposed elements seem to be significantly improved in terms of robustness and large load increments. This leads to an enormous gain in terms of computational cost comparing to the widely used EAS formulations.