Rot‐free mixed finite elements for gradient elasticity at finite strains

Through enrichment of the elastic potential by the second‐order gradient of deformation, gradient elasticity formulations are capable of taking nonlocal effects into account. Moreover, geometry‐induced singularities, which may appear when using classical elasticity formulations, disappear due to the higher regularity of the solution. In this contribution, a mixed finite element discretization for finite strain gradient elasticity is investigated, in which instead of the displacements, the first‐order gradient of the displacements is the solution variable. Thus, the C1 continuity condition of displacement‐based finite elements for gradient elasticity is relaxed to C0. Contrary to existing mixed approaches, the proposed approach incorporates a rot‐free constraint, through which the displacements are decoupled from the problem. This has the advantage of a reduction of the number of solution variables. Furthermore, the fulfillment of mathematical stability conditions is shown for the corresponding small strain setting. Numerical examples verify convergence in two and three dimensions and reveal a reduced computing cost compared to competitive formulations. Additionally, the gradient elasticity features of avoiding singularities and modeling size effects are demonstrated.


INTRODUCTION
The elastic potential in classical solid mechanical modeling is usually a functional of the first order gradient of deformation. Both in small strain and in finite strain elasticity this approach is sufficient for many applications. However, they do not take into account the finite resolution of the microstructure of the material. If the modeled geometry has, for example, sharp corners this can lead to nonphysical singularities of the strains and stresses resulting from these local models. In consequence, corresponding finite element simulations yield results, which are dependent on the mesh size around these singularities and numerical problems may arise. Due to the local nature of the classical elasticity formulations, the influence of material heterogeneities on the global elastic behavior is not taken into account. This becomes relevant when modeling specialized materials (such as metamaterials), in which material heterogeneities approach the scale of the macroscopic mechanical fields. Nonlocal approaches such as gradient elasticity, where the elastic potential is enriched by a dependency on the second-order displacement gradient represent a remedy in these cases. Although nonlocal approaches have been in existence for many years, due to recent developments of specialized materials, the research of formulations has experienced an increased interest.
Various generalized continuum models are based on the theory of the Cosserat brothers, 1 in which rotational degrees of freedom derived from the balance of angular momentum at the material point are considered in order to take into account nonlocal effects. Recent contributions have been made in this field (such as Reference 2, among several others). In the frequently cited, more general theory of Reference 3, microdeformations are used to take into account microstructural effects on the material behavior. As a special case, in Reference 3 the second-order gradient elasticity approach is introduced. In this case, microstuctural effects are modeled through the enhancement of the elastic potential by a dependency on second-order gradients of displacements. For an overview of gradient elasticity approaches the reader is referred to Reference 4 (see also Reference 5 for a classification with respect to other nonlocal models). While the gradient elasticity approach itself is rather straightforward, a remaining research challenge is the identification and physical interpretation of the additional constitutive parameters that are appearing in the gradient elasticity formulation. Recent developments have been made by Reference 6 in which a determination of the parameters based on action principles has been proposed (see also Reference 7). Another field of study which has experienced recent increase of interest is the development of models, which incorporate flexoelectric effects. In these models, second-order displacement gradients are coupled with electric polarization fields, which are of increasing interest in the research of nanoelectronics and nanowire semiconductors. 8 The main challenge in the development of numerical solution procedures for gradient-enhanced models is the increased continuity requirement for the interpolation functions of the solution fields, namely, C 1 continuity requirement for displacement-based finite elements. Straightforward corresponding C 1 -continuous finite elements require high polynomial orders of the interpolation functions and additional element degrees of freedom corresponding to the gradients of the solution fields. An investigation of a C 1 -continuous finite element formulation for gradient elasticity can be found in Reference 9. Another common approach is the use of isogeometric analysis (IGA) as discretization scheme. See Reference 10 for an IGA model for finite strain gradient elasticity, whereas in References 11 and 12 flexoelectric material behavior is modeled. While the realization of the C 1 continuity condition is straightforward with IGA, the discretization of complex (practically relevant) structures remain a challenge. Alternatively, mixed finite elements can be used. Through the introduction of additional mixed solution variables C 0 -continuous interpolation functions can be used. In Reference 13 a mixed finite element formulation is investigated, which incorporates both displacements and Mindlins microdeformations as solution variables. In the rather recent contribution of Reference 14 various finite element discretizations incorporating the displacement and strain field as solution variables are investigated in a two-dimensional (2D) linear small strain setting. Another small strain finite element approach is introduced in Reference 15, where both first and second displacement gradients are introduced as mixed variables and corresponding numerical investigations focus on the incorporation of point and line forces. The small strain two-dimensional (2D) and three-dimensional (3D) mixed formulations of References 16 and 17, respectively, incorporate displacements, displacement gradients and Lagrange multipliers as solution variables. For a corresponding mathematical investigation with respect to inf-sup stability and an extension to finite strains see Reference 18. So far, the main challenges of these existing mixed formulations is the computational costliness due to a large number of degrees of freedom and the lack of discrete inf-sup stability in some cases.
In this contribution a decoupled approach, based on the polyharmonic operator splitting approach of Reference 19, is investigated for the hyperelastic finite-strain gradient elasticity problem. In Section 3.1 the continuous large deformation variational formulation is introduced. Specifically, the case is examined in which the constitutive parameter corresponding to the nonlocal part of the formulation is significantly smaller than for the local part. Through a stabilization term, a loss of stability in this case (as also apparent in the existing non-decoupled mixed approaches of Reference 18, following References 16 and 17) can be avoided. In the following Section 3.2 several finite element discretizations, both in 2D and in 3D, are investigated. Numerical experiments (Section 4) verify convergence of the proposed formulations. In the 2D case a cost reduction compared to the nondecoupled mixed approach of Reference 18 (which is based on References 16 and 17) is shown. Furthermore, the gradient elasticity features of avoiding singularities and modeling size effects are illustrated.

FUNDAMENTALS
Section 2.1 gives preliminary remarks on notation, finite strain kinematics, and the Helmholtz decomposition, which is used for the decomposed Lagrange multiplier method introduced in Section 3. An overview on the continuum mechanical framework of the gradient elasticity approach is given in Section 2.4.

Tensor operations
For the dimension d ∈ {2, 3} let T (n) be a tensor-valued function of order n ∈ N 0 over (R d ) n . We define the scalar product T ⋅ T ∶ (R d ) n × (R d ) n → R n . Let ∧ denote the cross product, which is applied row-wise to tensor-valued functions. Explicitly, the cross products of a second-order tensor T and a vector s take the forms for d = 3 and d = 2, respectively. We define the row-wise applied derivative operators Div where e i (•) denote Cartesian base vectors.

Finite strain kinematics
For each position vector X ∈  of a material point in the reference configuration there exists a corresponding position vector x ∈  in the deformed configuration defined by the deformation map ∶  → . Through (X, t) = X + u(X, t), the deformation map can be described in terms of the displacement function u. We denote the deformation gradient by Within the scope of this contribution all derivative operators are with respect to coordinates in the reference configuration.

Helmholtz decomposition
Let  be a bounded, simply connected domain. Then, a second-order tensor function can be decomposed (see Reference 20, corollary 2.31) into a solenoidal (divergence-free) part c and an irrotational (rot-free) part g : Consequently, the two parts can be expressed by gradient and rotation functions, respectively: c = Rot and g = −∇g, with g ∈ H 1 0 being a vector-valued function and ∈ , where  is the space of vector-valued L 2 functions with vanishing integral mean for d = 2, while  is the space for tensor-valued L 2 functions with vanishing divergence and normal boundary conditions for d = 3; see Section 3.1.6 for details.

Gradient elasticity theory
When considering a gradient-elastic boundary value problem, the boundary of a given body  can be decomposed as follows: Here, Γ D and Γ N denote the standard Dirichlet-and Neumann boundary, while Γ H and Γ M denote a Dirichlet-and Neumann boundary corresponding to higher-order quantities (cf. References 3,21). The gradient elasticity solution u is the minimizer of the elastic potential Here, int denotes the volume integral of the internal elastic energy density . ext denotes the external potential of volume load f and surface traction t. Due to the appearance of second order displacement gradients in (8), the minimizer is sought to be in the Sobolev space H 2 of twice differentiable functions. Moreover we denote the displacement boundary u = u and ∇u ∧ N = ∇u ∧ N on Γ D and ∇u N = ∇u N on Γ H . Here N denotes the unit normal surface vector of the body in reference configuration and u denotes a prescribed displacement function. In what follows, the higher order surface traction vector (cf. Reference 10, see also Reference 3) r = 0 is assumed to be homogeneous on Γ M . Applying the stationarity condition leads to the following variational problem: For given f, t, and the previously discussed essential boundary conditions, find u ∈ H 2 such that for all u. Here, we denote ∇ 2 (•) := ∇(∇ (•)) as the second-order gradient and introduce the first-and second-order stress tensors P ∶= F and G ∶= ∇F , respectively.

THE ROT-FREE MIXED FORMULATION
This section introduces a finite strain mixed finite element approach, in which, through a rot-free constraint, the displacement and displacement gradient solution variables appear in a decoupled set of variational equations promising increased computational efficiency. While Section 3.1 discusses the continuous variational setting, in Section 3.2 suitable finite element interpolations are discussed and the corresponding set of discrete matrix-vector notations are obtained on the element level. The proposed approaches are complemented by a mathematical analysis of the small strain counterparts in Sections 3.1.6 and 3.2.5 for the continuous and the discrete setting, respectively.

Continuous variational framework
In the following it is assumed that the body  has a finite size and is simply connected with a connected essential boundary Γ D .

Decomposed Lagrange multiplier method
In order to arrive at a formulation, which contains only first order gradients and thus, enables C 0 continuous finite elements, the total elastic potential (8) is reformulated to the following Lagrangian: Here, the internal elastic potential int is a functional of the displacement gradient variable H and compatibility with the displacement u is enforced by the constraint term lag (13). The Lagrange multiplier = −∇g + Rot is decomposed according to relations (6) and (7) into a gradient and a rotation. Thus, by making use of the divergence theorem, (13) can be written as Note the exclusion of ∇u in the second term of (14), since the scalar product of rot and gradient tensor functions vanishes. The second term of (14) can be interpreted as constraint term enforcing H to be rot-free, which is a necessary condition for gradient functions.
Remark 1. In the case of a vanishing nonlocal contribution of the strain energy density (G = 0), can be identified as the first Piola Kirchhoff stress tensor P and −∇g + Rot can be viewed as corresponding split of the latter. This can be verified by the Euler-Lagrange equations from variation of (12) according to the Hu-Washizu variational principle. 22

Stabilization for the limit case of vanishing nonlocal contribution
In order to take into account the case in which the nonlocal contribution becomes small, the stabilization term is added to (12). Through adding this augmentation term, the problem remains well posed even in the limit case of vanishing nonlocal contribution and the inf-sup stability condition is fulfilled. The problem (12) is not changed by the stabilization term since the solution satisfies H stab = 0. A corresponding mathematical analysis can be found in Section 3.1.6.

2D variational formulation
In the 2D case d = 2 with definition (4) Rot H simplifies to the vector Thus, in this case the Lagrange multiplier is also vector-valued and fulfills the inf-sup stability condition (cf. (B1) in Section 3.1.6) without further restrictions with respect to differential operations, namely as element of the Sobolev space L 2 with a fixed mean integral value ∫ B dV = 0. With (12), (14), and (15) the weak form corresponding to the Lagrangian = int + ext + lag + stab reads We seek u and g in the Sobolev space H 1 ∶= { u ∈ L 2 ∶ ∇ u ∈ L 2 } for all u, g ∈ H 1 . Moreover, we use Dirichlet boundary conditions u = u and g = 0 on Γ D with u being the prescribed displacement. Corresponding to the displacement boundary, the tangential direction of the displacement gradient variable is prescribed on Γ D with H ∧ N = ∇u ∧ N, where N denotes the surface normal vector in the reference configuration. The displacement gradient variable H ∈ H 1 (2) is sought for all H ∈ H 1 (2) . Furthermore, the normal direction of the displacement gradient corresponding to the higher-order essential boundary is prescribed by H N = ∇u N on Γ H .
Remark 2. The equation system (17) through (20) depicts a split of the second-order weak form (11) into a set of first-order equations. (cf. Reference 19, where this approach was proposed for polyharmonic problems). In the present system of equations, the displacement u is decoupled from the main problem (18) and constraint Equation (19). The displacements only appear in the simple Laplace-type Equations (17) and (20). In the following, (17) and (20) are referred to as pre-and postprocessing step, respectively. This set of equations follows first ideas in Reference 23 in the context of discretizing Kirchhoff's equations of thin plate bending by C0 finite elements. Note, that the higher-order tractions r, which were assumed to be zero could also be applied through extending (18) correspondingly.

3D variational formulation
In the 3D case d = 3, to fulfill the inf-sup condition, the Lagrange multiplier needs to be divergence-free, namely it needs to be an element of the Sobolev space, that is, ∈ H(Div 0 ) (2) with vanishing boundary trace N = 0 on Γ D (cf. Proposition 7). Thus, in this case, the second divergence-free constraint term is added to (12), where is a vector-valued, second Lagrange multiplier. For the space L 2 and a fixed mean integral value ∫ B dV = 0 is a suitable choice (cf. (39) and Proposition 3). Similar to the previous case, together with (21), the weak form corresponding to the Lagrangian = int + ext + lag + stab + div reads and u and g according to (17) and (20), and the spaces for u, g, and H unchanged.

Alternative 3D formulation
Starting point for an alternative approach is again the Lagrange multiplier decomposition of Section 3.1.1. However, the split takes the form = −∇g + c , where the variable c is sought to have the characteristics of a 3D rotational tensor field, namely to be divergence-free. Subsequently, the potentials (14) and (21) are modified as Note, that in the second term of (25) only the rotation part of H is controlled by the Lagrange multiplier c . This is due to the fact that with (26) c is specifically sought to be a rotation function and thus not having an influence on the divergence part of H. Proceeding analogously to Section 3.1.4, we obtain the modified weak forms and u and g according to (17) and (20), and the spaces for u, g, and H unchanged. For c a zero trace on the Neumann boundary c N = 0 on Γ N is prescribed.

Mathematical analysis of the continuous formulations
In this section we prove the stability of the proposed approach in the small strain framework. Proofs of the propositions, lemmas, and corollaries of this section are given in Appendix B1. For the mathematical analysis it will be assumed that the displacements are small and the internal energy in (9) is additively decomposed into a term quadratic in ∇u and a term quadratic in ∇ 2 u. Let L 2 (, (R d ) n ) denote the space of L 2 functions with values in (R d ) n and define the L 2 -inner product and the L 2 -norm for an n-order tensor valued function We furthermore define the space For the left-hand side of (11) we define the bilinear form where C denotes a constant fourth-order elasticity tensor and c 1 > 0 a constitutive parameter associated with the higher order stress response. For a homogeneous boundary  = Γ D = Γ H the small-strain gradient elasticity problem corresponding to (11) simplifies to finding u ∈  for a given volume force f ∈ L 2 so that For the mathematical investigation we define the spaces where H m 0 (; (R d ) n ) denotes the H m functions with values in (R d ) n satisfying zero Dirichlet boundary conditions on . Furthermore, the definition for the Lagrange multiplier involves the space Define the bilinear form with added stabilization term ã ∶  ×  → R by The small strain-analogon to (18) seeks (H, ) ∈  ×  for a given g ∈  such that for all ( H, ) ∈  × . We now show that (37) is stable and robust for c 1 → 0. The robustness is proved with respect to the following norm on  that depends on the parameter c 1 that describes the nonlocal contribution, Note that, if Rot H = 0, then H is a gradient of a function in H 1 0 . Therefore, Korn's inequality implies that for min{c 1 , } > 0 this is in fact a norm.
The following proposition proves the unique existence of solutions to (18) with stabilization in the linear strain case.
Remark 3. Note that Proposition 1 remains true, if the stabilization term is not added to the bilinear form, that is, if ã is replaced by a in (37). This comes from the fact that functions in the kernel of b are rot free and a is coercive on this kernel. However, this is no longer true for the discretization and therefore, we include the stabilization term.

Stability of the 3D variational formulation with two Lagrange multipliers
We now show that in the small-strain framework the solution of the 3D variational formulation (22) through (24) coincides with the solution of the original problem. Define the spaces of the Lagrange multipliers where The small strain analogon to (22), (23), and (24) reads: for all ( H, , ) ∈  ×Q × .

Stability of the alternative 3D variational formulation
In the following we discuss the stability of the small strain analogon to the alternative 3D formulation (27) through (29).
, that is, the space of all linear and continuous mappings from Further, define the norm We consider the problem: The small strain analogon to (27), (28), and (29) then reads: The following proposition states, that problem (44) is well posed, meaning it has a unique solution, which coincides with the solution of the original problem.

Finite element approximations
For the finite element discretization, a partition of  into a set of simplices  = ⋃ e T e is considered, where  is the set of corresponding element faces. By standard FE procedure, the solution variables and test functions of the continuous weak forms of Section 3.1 are replaced by piecewise polynomial functions, that are defined in the following subsections. Moreover, matrix vector representations corresponding to the discretized weak forms are given.

Discretization of the pre-and postprocessing step
In the following, vectors containing nodal degrees of freedom corresponding to the pre-and postprocessing problem (17) and (20) are denoted by d u and d g . We introduce the vector-matrix interpolation operators where N u and B u represent suitable finite element interpolation matrices containing Lagrange shape functions and corresponding derivatives (cf. Appendix A1). Inserting (45) and (46) and reformulating (17) and (20) leads to the following matrix equations where d H,ext contains the nodal degrees of freedom corresponding to H h computed by the main step. Note, that the preand postprocessing steps are linear and thus, even in the finite deformation regime, no further linearization of (47) and (48) is necessary.

Discretization of the main step in 2D
For the discretization of the weak form (18), the following matrix interpolation operators are introduced (cf. Appendix A1): Since the solution variable H is sought in the H 1(2) Sobolev space, Lagrange interpolation functions are used in matrices N H and B H . The rotation operator matrix R H is constructed from the suitable components of the gradient operator with relation (16) (cf. Appendix A1). In the 2D case the Lagrange multiplier ∈ L 2 is discretized with interpolation functions as As discussed in the mathematical stability analysis of Section 3.2.5 in 2D, any pairing of interpolation functions, which is stable for the Stokes problem (cf. Reference 24) is a suitable choice for the approximations H h and h . Therefore, the Mini interpolation scheme P1B H -P1 and the Taylor-Hood interpolation schemes P2 H -P1 and P3 H -P2 are used (cf. overview in Table 1). Note, that in the Mini interpolation scheme N H consists of linear Lagrange shape functions corresponding to the vertex nodes of the (linear) P 1 -triangle and the cubic Lagrange shape function corresponding to the interior node of the (cubic) P 3 -triangular element. In this context the latter is also referred to as volume bubble function (cf. Section 3.2.5, see also Reference 25). The discretization of (49) and (50) reads .
TA B L E 1 Overview of the two-dimensional finite element interpolation schemes used. For each scheme, the number of interpolation nodes is given in parantheses. The pre-and postprocessing elements are denoted by P2 g,u , P3 g,u , and P4 g,u . The order of this listing corresponds to the order of the elements appearing in the  (3) Quadratic (6) P2 H -P1 (Taylor-Hood) Quadratic (6) Linear (3) Cubic (10) P3 H -P2 (Taylor-Hood) Cubic (10) Quadratic (6) Quadratic (15) Here, k HΦ = k T ΦH represent the element submatrices corresponding to the rot-free constraint condition, whereas r ext g depicts the nodal load vector with d g,ext obtained from the solution of the preprocessing step (47). In general, the strain energy function is non-linear in H h and ∇H h and thus, nonlinear in d H . Therefore, the incremental Newton-Raphson loadstep solution procedure is used to obtain the solution. This requires the linearization of (51). The linearized problem then reads where Δd H and Δd Φ denote increments of the nodal solution vectors, and d H denotes the nodal solution from the previous step. The element tangent submatrix k H can be written as where the second term consists of the discretization of the rot-rot stabilization term. The integrals ∫ T • dV over the element are evaluated numerically via Gauss quadrature of the reference element.

Discretization of the main step in 3D
In order to arrive at a stable discrete set of equations, in the 3D case the Lagrange multiplier is discretized conforming to the Sobolev space H(Div) (cf. Section 3.1.6). Therefore, the lowest order Raviart-Thomas interpolation procedure is used for h . By using appropriate interpolation operator matrices, we denote the approximations as where S Φ and D Φ contain the Raviart-Thomas shape functions and their divergence, respectively. Details on the construction of S Φ and D Φ can be found in Appendix A1. Since the second Lagrange multiplier h appearing in (23) and (24) is sought in L 2 , a piecewise constant approximation is used. We define the interpolation matrix N (cf. Appendix A1) with Note, that since h is piecewise constant, the corresponding degrees of freedom d may be condensed at the element level. An overview of the used Lagrange-and Raviart-Thomas interpolation schemes that fulfill the mathematical stability analysis of Section 3.2.5 can be found in Table 2. Here, for the interpolation of H h , the matrix N H consists of the linear Lagrange shape functions corresponding to the four vertex nodes of the (linear) P 1 -tetrahedron and the cubic Lagrange shape functions corresponding to the four midface nodes of the (cubic) P 3 -tetrahedral element. In the mathematical analysis of Section 3.2.5, the latter are referred to as face bubble functions (see also Reference 19). Analogously to the  (56) we arrive at the following linearized system corresponding to the discretization of (22) through (24) Lin The element tangent submatrices k Φ = k T Φ and k HΦ = k T ΦH , which correspond to the rot-constraint and the div-constraint, respectively, are given by In the following, the discretization scheme (57) is denoted by P1FB H -RT0 -P0 .
Remark 4. In Reference 19, for a similar discretization the approach following Reference 26 is taken in which the first Lagrange multiplier is discretized with P 0 functions. In a second constraint condition, added to the discrete system of equations, through discrete Lagrange multipliers h continuity of in normal direction across the element faces is enforced. While this second constraint is not part of the continuous formulation, the global number of equations is similar to those of the proposed approach. For d = 3 the Lagrange multiplier and the additional constraint equation of one element contribute ndof + ndof h = 9 + 3 ⋅ 4 = 21 equations to the global system, where ndof may be condensed at the element level. Meanwhile, the P1FB H -RT0 -P0 -element counts ndof + ndof = 3 ⋅ 4 + 3 = 15 corresponding degrees of freedom, where ndof may be condensed at the element level.

3D discretization of the alternative formulation
For the discretization of the alternative approach (27) through (29), the same interpolation operators as in the previous section are used. Thus, the interpolation of c is given by Consequently, the structure of the linearized system of the discretized equations of (27) through (29) is analogous to the one in the previous section: However, here the tangent submatrix corresponding to the rot-constraint takes the form The submatrices k H and k Λ c = k Φ remain unchanged (cf. Section 3.2.3). In the following, the discretization scheme corresponding to (61) is denoted by P1FB H -RT0 -P0 (cf. Table 2).

Mathematical analysis of the finite element approximations
In this section we discuss the stability of the finite element approximation schemes discussed in the previous sections. Proofs of the propositions and lemmas of this section are given in appendix B1.

Stability of the 2D finite element approximations
Let A ≲ B abbreviate that there exists a generic constant C < ∞, that does not depend on the critical parameter like the mesh-size of a discretization or the parameter c 1 and defined below, such that A ≤ CB.
In the case d = 2, any finite element pairing that is stable for the Stokes equations is a suitable choice for the discretization of (18) and (19). This is proved in the following proposition.
For a simplicial triangulation  of , we introduce the following notation for the space of tensorial piecewise polynomial functions of degree ≤k: If T is a vector (n = 1), the corresponding superscript index will be omitted. We denote the discrete spaces corresponding to the MINI element by Here, B 3 ( , R 2 ) (2) is the space of cubic bubble functions defined on the triangle (cf. Reference 25). The discretization of (37) then seeks The inf-sup condition for the Stokes equations (Reference 24, theorem 8.8.1) together with Proposition 5 proves the stability of this discretization, and therefore, a unique solution exists and the error satisfies (63). The Taylor-Hood finite element subspaces (cf. Reference 25) are  24 together with Proposition 5 proves that this defines a stable discretization. For k = 1 and k = 2, respectively, the corresponding finite elements are referred to by P2 H -P1 and P3 H -P2 in the following (Table 1).

3.2.7
Stability of the 3D finite element approximations As discussed in Section 3.1.6, for d = 3 the Lagrange multiplier ∈  is required to be divergence free. Thus, the formulation with an additional Lagrange multiplier from Section 3.1.4 is employed with the following discrete subspaces: Here, B 3 ( , R 3 ) denotes the space of cubic face bubble functions vanishing on the element edges. Furthermore, RT 0 ( ; R 3 ) is the finite element space of Raviart and Thomas, which consists of polynomial functions, which are continuous in normal direction across the interelement boundarys. 27 The discretization of the modified problem (41) seeks This linear problem corresponds to the P1FB H -RT0 -P0 discretization of Section 3.2.3.
The following proposition proves the stability and an error estimate for this discretization.
The discretization of (41) with the above choice of spaces has a unique solution where (H, ) ∈  ×  is the solution to (37).

Stability of the alternative 3D finite element approximation
This section investigates a discretization of the alternative formulation (44), which is the linear analogon to (27) through (29), which showed appropriate performance in the numerical experiments. The space in (44) is discretized withQ ∩Q h withQ h ∶= RT 0 ( ; R 3 ) (2) ∩ H(Div) (2) . In particular, those functions have more smoothness than requested byQ, namely they are in H(Div) (2) . Therefore, the duality pairing can be replaced by the L 2 product. Moreover, since Div h ⊆ P 0 , the equation ( h , Div h c ) L 2 () = 0 for all h ∈ P 0 guarantees that a function h c ∈Q h is in fact inQ. Thus, the divergence free condition from can be incorporated via an additional Lagrange multiplier. In conclusion, the discrete problem for (44) seeks and  h from (68) from the previous section.The discretization of (71) corresponds to the P1FB H -RT0 -P0 formulation of Section 3.2.4. Note that h can be both identified as Lagrange multiplier enforcing h c to be divergence free and P 0 approximation of the displacement.
As mentioned above, since Div h ⊆ P 0 , the third equation of (71) shows that h c ∈ h ∩. Therefore, the discretization (71) is a conforming discretization of (44) in the space  h × ( h ∩) ⊆  ×. Proposition 4 proves the well-posedness of problem (44). Furthermore, the discretized problem is stable in the numerical experiments of section 4.

NUMERICAL TESTS
In this section, the proposed discretizations are numerically tested. An overview of the considered elements is given in Table 3. For the implementation, the AceGen/AceFEM software package has been used, which is based on automatic differentiation. Through construction of the elastic potential containing the approximated solution fields and differentiation with respect to the nodal degrees of freedom on the element level, both the element residual and tangent matrix are derived (cf. Reference 28). The integrals are evaluated numerically with Gauss integration over the corresponding reference coordinate space. In Table 3 the used Gauss integration orders are shown for each discretization. The standard Newton-Raphson load step solution procedure is used for the solution of the nonlinear global system. For the solution of the linearized system of equations, the PARDISO solver is used. For the local hyperelastic energy density the Neo Hookean ansatz of Reference 29 is used: and I 1 = tr C, J = det F. The coefficients c and d are set to c = ∕4 and d = ∕2 (see also Reference 30), where the parameters and are the Lamé constants which can be computed from the Young's modulus E and the Poisson's ratio by = E ∕((1 + )(1 − 2 )) and = E∕(2(1 + )). Based thereon, the Cauchy stress tensor can be computed from = J −1 P F T where the first Piola-Kirchhoff stress tensor is obtained from P = F . The nonlocal part of the energy density is taken from: 31 where c 1 is a nonlocal elasticity constant. With l = √ c 1 ∕ , the constant can be scaled to a length unit parameter. 10 In the following numerical examples the elasticity parameters are chosen with E = 500 MPa and = 0.3, unless stated otherwise. For the 2D analysis, the constitutive model is implemented under the plain strain assumption.

Unit square with smooth solution
For this test, the domain  is considered to be a 2D unit square with the dimensions 1 × 1 mm 2 (cf. Figure 1(A)), where also the element alignment for the initial mesh is depicted. A cross-pattern element patch has been chosen, since a mesh with a unidirectional element orientation may show inferior convergence results (cf. Reference 32). In order to construct an exact error measure, the large displacement smooth reference solution u of Reference 18(p. 12) is used. The reference solution u|  = 0 and its first-order gradient ∇u|  = 0 vanishes on the boundary, such that the strong form of (11) reads −Div P + Div Div G = f.

Element name Number of Gauss points AceGen ID
With (77), the volume load f is computed analytically for each Gauss point and then evaluated through r ext u of the discrete Equation (47). For each uniform mesh refinement step the error ||∇u − H h || L 2 () and ||∇ 2 u − ∇H h || L 2 () is measured by computing the numerical integrals, in which the exact solution and the interpolation of the finite element solution are evaluated at the Gauss points. The polynomial order of both pre-and postprocessing discretization is k + 1, where k is the polynomial order of the H h interpolation of the main step (cf. Table 4). The elements are compared to the P2 u -P1B H -P1 element of. 18 There, following the previously existing mixed approaches (cf. References 16 and 17), the displacement gradient is replaced by H only in the nonlocal energy nloc ∶= nloc (H) and loc ∶= loc (∇u) remains a functional of u. The discretization of the P2 u -P1B H -P1 element is in accordance to the naming convention used here (cf. Table 3). Moreover, the Lagrange multiplier is not additively decomposed and sought in L 2 . Therefore, the corresponding variational problem is not decoupled.

Influence of the stabilization term
For the results depicted in Figure 2 cubic with u h ∈  h ∩ P 3 , cf. Table 4). The convergence rate h 4 is reached even in the limit case of vanishing nonlocal contribution (l = 0). While the Taylor-Hood-type element P3 H -P2 behaves analogously with an approximate L 2 convergence rate of h 5 , the P1B H -P1 element experiences a shift from order h 3 to h 2 as l approaches zero. The tangent matrix of the element P2 u -P1B H -P1 becomes singular in this case. In Figure 2(B) the convergence behavior of the P2 H -P1 element depending on the numerical value of is depicted. While it can be observed that numerical values up to = 10 N∕mm 3 L do not significantly alter the convergence behavior, a higher value than = 10 N∕mm 3 L leads to slightly varying results. This may be due to the fact that in this case the rot-free constraint is dominated more by the stabilization term, which for high values of can be interpreted as penalty term. Moreover, for higher values of an increase of the condition number of the global tangent matrix and, for more extreme values, loss of convergence of the Newton-Raphson solution procedure is observed.

Computing efficiency
In order to compare the computational cost of the finite elements, Figure 3 shows a displacement convergence study at the point A = (0.5, 0.5) mm. In Figure 3(A) the displacement of the finite element solution relative to the exact solution u y (0.5, 0.5) = −0.1953 mm is plotted versus the computing time required for both assembly and solution of all iterative steps of the solution procedure in total at each refinement step. Since the scale of the problem is rather small, multiple simulations are made at each refinement step. The computing times of all simulations for each refinement step are averaged to take into account variations of the calculating capacity of the computer. Figure 3(B) depicts the average computing time of the elements corresponding to the converged state. The convergence criterion is |Δu| ≤ 0.001 mm and the corresponding refinement stage is marked by a bullet in Figure 3(A). While for all elements the computing time corresponding to the main problem (53) is considered, for element P3 H -P2 the total computing time including pre-and postprocessing is additionally taken into account (the corresponding plot in Figure 3(B) is marked with an asterisk). It becomes evident, that even for this small-sized problem, the proposed elements have a computational advantage over the three-field element P2 u -P1B H -P1 . Moreover, the contribution of the pre-and postprocessing step to the total computing time is relatively small even in the present case of P4 u,g elements (cf. Table 4).

Cook's membrane problem
The boundary of the Cook's membrane problem (cf. Figure 4 boundary conditions are prescribed. The nonlocal parameter is l = 1 mm and the stabilization parameter is = 10 L. The mesh (cf. Figure 4) is refined uniformly. We denote as the element average stress, where h is the Cauchy Stress obtained by postprocessing of the finite element solution.

Removal of singularities
In order to show the ability to avoid singularities, the gradient elasticity elements are compared to the quadratic local elasticity element P2 u (discretization of (11) with G = 0 and u h ∈ H 1 Γ D ∩ P 2 ). For this, the average stress ⟨ ⟩ e of the element adjacent to the singularity point B = (0, 44) mm is evaluated. In Figure 5, the element average stress is plotted versus the global number of degrees of freedom corresponding to each refinement step. It becomes evident, that for an increasing number of degrees of freedom, the element average stress of the gradient elasticity elements converges similarly toward a finite value, while the element average stress of the local P2 u element increases. This behavior is illustrated in Figure 6, in which contour plots of the average element stresses are shown for the sixth refinement step (cf. bullet marks in Figure 5).

3D Unit cube with smooth solution
For the investigation of the convergence behavior of the 3D elements the unit square domain is extended to 3D. Here,  has the dimensions 1 × 1 × 1 mm 3 (cf. Figure 7(A)). The smooth finite strain reference solution of Reference 18 (p. 14) is used and again, with the strong form (77) the right-hand side f is derived analytically. The convergence of the L 2 -error is depicted in Figure 7(A). The nonlocal parameter is l = 1 mm and the rot-rot stabilization term is switched off with = 0. It is observed that both investigated elements converge similarly. From Figure 7(B) it can be seen, that for both elements the L 2 -norm | |Rot H| | L 2 () of the constraint condition decreases approximately with the order h.

3D Cook's problem
In this subsection the 3D Cook's problem (cf. Figure 8

Influence of the nonlocal parameter on the displacement response
For a mesh refinement stage corresponding to 5120 elements (cf. Figure 8(B) through (D)), the displacement response at the corner point A = (48, 20, 60) mm is evaluated for various nonlocal parameters l relative to the length of the geometry L. The compared elements are the proposed P1FB H -RT0 -P0 and P1FB H -RT0 -P0 discretizations as well as the P2 u -P2 H -P1 element from Reference 18, which again is based on a non-decoupled approach for the Lagrange mutliplier and loc ∶= loc (∇u) remaining a functional of u (cf. explanation in Section 4.1). From Figure 9(A) it can be seen, that for l/L = 0 the solution of the proposed formulations coincides with the converged solution u z (A) = 21.41 mm of the local P u displacement formulation. Note, that in this case for the element P1FB H -RT0 -P0 (cf. Table 3) a stabilization parameter = 10 L yields sufficient results, while for the element P1FB H -RT0 -P0 a higher value = 500 L for the stabilization parameter is necessary to obtain results that are coinciding with the local converged solution. Moreover, in this case (l/L = 0) the Newton-Raphson solution procedure fails to converge for the P2 u -P2 H -P1 element. Nevertheless, for l/L ≥ 0.05 mm it can be seen from Figure 9(A), that for the three gradient elasticity elements the displacement response depending on l/L is similar. The dependency of the displacement response on the ratio l/L can be related to the modeling of size effects (cf. Reference 4). Similar to the 2D results of Section 4.2, the contour plot of Figure 8(C) also shows a reduction of the element average Cauchy stress for l = 0.0005 L at the singularity point B = (0, 0, 44) mm compared to the local element P2 u (cf. Figure 8(B)). The computing time, which corresponds to the simulation marked with a bullet in Figure 9(A) is comparable for both elements, see Figure 9(B). Again, as in the numerical evaluation of the proposed

CONCLUSION
A new C 0 -continuous mixed finite element formulation for finite strain gradient elasticity was introduced. The proposed approach was based on a split of the constraint term, which enforces compatibility between the mixed variables. As a result, a set of decoupled variational equations was obtained. In the linear setting, stability of the corresponding continuous formulation was proven. Robustness was shown even in the limit case of a vanishing nonlocal contribution and corresponding suitable finite element discretizations were introduced both in the 2D and the 3D case. Numerical tests of the nonlinear finite strain formulation showed appropriate convergence behavior and robustness of the proposed finite elements. The computing time was compared to the 2D three-field approach of Reference 18 (see also Reference 16) and a notable cost reduction was shown.
and similarly for d = 3, the relation is used.

A.2 Raviart-Thomas interpolation matrices
Similarly to the preceding explanations, the Raviart-Thomas interpolation matrices can be constructed from the following relations.
Here, the nodal degrees of freedom d I Φ correspond to the first moment of the element-surface normal direction of the solution variable. 27 Moreover, Ψ I depicts the lowest order Raviart-Thomas shape function in the physical coordinate space corresponding to the Ith mid-face element node. Opposed to the Lagrange interpolation operators, the Raviart-Thomas shape functions are additionally mapped via Piola transformation from the reference coordinate space to the physical coordinate space (for a description on the construction of Ψ I see Reference 32).

APPENDIX B. PROOFS
In this section proofs of the propositions and lemmas of the mathematical analysis of Section 3 are given. The proofs are listed in the order as the corresponding propositions etc. appear in the text. This proves the assertion for d = 2. Let now d = 3. Due to (Reference 34, proposition A.1) for each ∈  there is a ∈  with = rot and || || H 1 () ≤ c|| || L 2 () with a constant c < ∞. Since ||| • ||| ≲ || • || H 1 () as above, we obtain for the choice H = This completes the proof. ▪ Proposition 7 is the main ingredient of the following proof of proposition 1.
Together with (17) Therefore, u solves (33). The other direction follows from this equivalence and the existence of unique solutions to (33) and (17), (37), and (20). ▪ Proof of Proposition 3. To prove Proposition 3 we need an abstract lemma, which generalizes Brezzi's splitting lemma 35 from a system with one constraint to a mixed system with two constraints. The lemma is also used in the stability analysis of Proposition 6. Let , , and  be Hilbert spaces with norms || • ||  , || • ||  , and || • ||  and let a, b, and c be bilinear forms on  × ,  × ,  × , respectively. We consider the abstract problem: Find (H, , ) ∈  ×  ×  such that Thus, ã satisfies the inf-sup condition on Z(C) andc satisfies an inf-sup condition. An application of Brezzi's splitting lemma on ã andc yields the assertion. ▪ We proceed with the proof of Proposition 3. (Div , ) L 2 () || || H 1 () .

Proof of Proposition
The inf-sup condition for the Stokes equations, also known as Ladyzhenskaya lemma, 33 proves This proves the inf-sup condition for the bilinear form (Div •, •) L 2 () onQ × . The inf-sup condition for b is proved in Proposition 7. These two inf-sup conditions together with the continuity of all three involved bilinear forms and the coercivity of ã as in the proof of proposition 1 together with Lemma 1 yields the unique existence of solutions of (41).
If (H, , ) ∈  ×Q ×  is a solution to (41), then Div = 0 and therefore (H, ) ∈  ×  is a solution to (37). The other direction of the equivalence follows from this and the unique existence of solutions for both problems. ▪