Three‐field mixed finite element formulations for gradient elasticity at finite strains

Gradient elasticity formulations have the advantage of avoiding geometry‐induced singularities and corresponding mesh dependent finite element solution as apparent in classical elasticity formulations. Moreover, through the gradient enrichment the modeling of a scale‐dependent constitutive behavior becomes possible. In order to remain C0 continuity, three‐field mixed formulations can be used. Since so far in the literature these only appear in the small strain framework, in this contribution formulations within the general finite strain hyperelastic setting are investigated. In addition to that, an investigation of the inf sup condition is conducted and unveils a lack of existence of a stable solution with respect to the L2‐H1‐setting of the continuous formulation independent of the constitutive model. To investigate this further, various discretizations are analyzed and tested in numerical experiments. For several approximation spaces, which at first glance seem to be natural choices, further stability issues are uncovered. For some discretizations however, numerical experiments in the finite strain setting show convergence to the correct solution despite the stability issues of the continuous formulation. This gives motivation for further investigation of this circumstance in future research. Supplementary numerical results unveil the ability to avoid singularities, which would appear with classical elasticity formulations.

Another field of application is the modeling of very small structures. When the specimen of interest reaches a microscopic size, the influence of the microstructure on the constitutive behavior becomes significant. These so called size-effects cannot be modeled by classical local elasticity formulations. While molecular or even atomistic models can be used to capture these effects, corresponding simulations, however, quickly reach their computational limits, when the modeled nanostructure is more complex. Gradient elasticity formulations on the other hand capture the size effects through additional constitutive parameters imposed in the gradient enriched term of the elastic energy. Compared to atomistic models they are computationally much cheaper. Examples of the simulation of size-effects with a focus on the general formulation and numerical aspects can be found in Askes et al., [2,4,5] while investigations of gradient elasticity approaches with a focus on specific applications can be found in Aifantis et al. [6][7][8][9][10] Another application of gradient elasticity approaches is the instable elastic energy of domains with martensitic phase transformations. [2,3,11] Significant theoretical foundations of these various applications date back to the 1960s gradient enriched elastic continuum theories with notable works of Mindlin [12] and Toupin. [13,14] A noteworthy adaptation which proved to be a simplified special case of Mindlin's theory was established by Aifantis. [15] For a comprehensive historical overview the reader is referred to Askes and Aifantis. [2] The idea of gradient enrichment for elasticity problems gave also rise to numerous gradient enriched plasticity models (cf. Aifantis et al. [4,16,17] amongst others) and also to approaches for the regularization of damage formulations through gradient enrichment (cf. de Borst et al. [18,19] amongst others). Investigations with respect to evaluation of constitutive gradient elasticity constants are undertaken in Peerlings and Fleck. [20] While most of the gradient elasticity models are restricted to the regime of small strains, investigations of gradient elastic constitutive models for large deformations can be found in Beheshti [21] and Triantafyllidis et al. [22] Despite its before mentioned advantages over classical formulations, gradient elasticity formulations have not found a significant employment in general practical applications yet. One reason for this is a difficult numerical implementation due to their higher continuity requirements for corresponding finite element formulations. While the underlying equations of classical formulations are partial differential equations of second order, the corresponding equations of gradient elasticity formulations are of fourth order. Thus, discrete weak forms of classical elasticity formulations require 0 -continuous interpolating functions, while gradient elasticity formulations require 1 -continuous interpolations, whose numerical implementation, especially in the finite element framework is not straightforward.
An investigation of meshless discretization methods for the gradient elasticity formulation can be found in Askes and Aifantis, [5,23] whereas in Papanicolopulos et al. [24] a 1 -continuous finite element with a preprocessing smoothing based on a least squared minimization is presented. The drawback of the latter one is the restriction to structured meshes only. The works of Rudraraju et al. [3,11] consider isogeometric analysis (IGA) as discretization scheme. Although the realization of the 1 -continuity condition is straightforward in IGA, the well-known difficulties of meshing more complex structures and boundary conditions remain a research challenge.
Alternatively, by deriving mixed gradient elasticity formulations, standard 0 -continuous finite element approximation schemes can be obtained. Based on the idea of Ru et al., [25] in which the fourth-order gradient elasticity problem of Aifantis [15] is split into two second order problems, a 0 -continuous corresponding finite element formulation is possible. [26] While the latter is restricted to the special formulation of Aifantis, [15] in Shu et al. [27] and Zybell [28] a mixed approach is investigated, which enables 0 -continuous discretizations in the framework of the more general Mindlin-Toupin gradient elasticity theory.
The mixed approaches of Shu et al. [27] and Zybell [28] are promising with respect to their use in commercial software, because corresponding discrete forms are conform with well established standard finite element schemes, that is, the modeled structures can be arbitrarily complex since associated meshing procedures are available. However, due to the introduction of additional discretization variables, the computational cost is relatively high. Moreover, a mathematical analysis with respect to stability and the extension of specific formulations to nonlinear problems are still missing in the literature. This is the starting point of the present contribution, in which three-field mixed finite element formulations are proposed for the general nonlinear hyperelastic framework. These are then mathematically and numerically investigated with respect to stability.
In Section 2 a continuum mechanical overview of the gradient elasticity theory is given, followed by the introduction of the continuous three-field formulation for finite strains in Section 3 together with stability considerations in Section 4. The discretization procedure for the corresponding finite element formulations and the stability investigations of various approximations are given in Section 5. In Section 6 numerical examples are presented and a conclusion is given in Section 7.

Notation
For  = I R, I R or I R × , the space of tensor-valued functions whose components are square integrable over a body  is denoted by 2 (; ). The corresponding 2 -inner product and the 2 -norm for any tensor function Let Curl denote the Curl-operator, that is and let ∧ denote the cross product. For vector or tensor valued functions, the operators ∇, div and Curl are applied row-wise. The cross product ∧ ∶ I R 3×3 × I R 3 → I R 3×3 acts on tensor fields as for ∈ {1, 2, 3}. Let 1 (; ) denote the Sobolev space of 2 -functions with gradient in 2 and let 1 0 (; ) denote the space of 1 (; ) functions with vanishing trace, that is, wherein denotes the unit vector normal to .

Kinematic fundamentals
Let the bounded Lipschitz domain  ⊆ I R ∶ ∈ {2, 3} be the body of interest in the reference configuration and  ∈ I R ∶ ∈ {2, 3} be the body in the current configuration. For each material point defined by the position vector ∈  in the reference configuration there exists a corresponding position vector ∈  in the current configuration defined by the deformation map ∶  → . The deformation map can be described by ( , ) = + ( , ) in terms of the displacement function .
Denote ∇(•) = (•)∕ as the gradient of a vector-/tensor-valued function with respect to the position vector in the reference configuration. The deformation gradient is written as with 1 denoting the identity tensor of second order.

Gradient elasticity theory
According to the gradient elasticity theory, [12,14] the boundary of  can be decomposed in the following way: with Γ and Γ being the standard Dirichlet-and Neumann boundaries, respectively, and Γ and Γ being Dirichlet-and Neumann boundaries on which higher-order quantities are prescribed. The hyperelastic strain energy density, which in classical formulations is a function of the deformation gradient, is enriched by a dependency on the second-order deformation gradient, that is ⋅ ⋅ = ( ) (classical elasticity) ⇒ ⋅ ⋅ = ( , ∇ ) (gradient elasticity). (2.4) In order to retain frame invariance, the preferred choice for the non-gradient enriched elastic energy density are functions ∶= ( ), which can be expressed in terms of the right Cauchy tensor = . Corresponding investigations with respect to frame invariance of the gradient enriched elastic energy can be found in Triantafyllidis and Aifantis [22] and will be taken up in Section 6.
Denote D(•) = ∇(•) ⋅ as the derivative in normal direction defined on the surface . The gradient elasticity boundary value problem can be formulated as the minimization problem where Π ext is the external potential of volume load , surface traction and surface moment 1 , which is work conjugated to the displacement gradient D̃in normal direction. The minimizer , for which Π[̃]∕ | =0 = 0 holds with̃= + is sought in ∈ 2 (; I R ) with =̄on Γ and D =̄on Γ , where the essential boundary conditions consist of prescribed displacements̄and prescribed displacement gradients in normal direction̄. Because of the affiliation to 2 , corresponding discrete functions of have the requirement to be 1 -continuous, for which the design and implementation of conforming finite elements is not straightforward.

THREE-FIELD MIXED FORMULATION
In order to arrive at a formulation, which enables 0 -continuous discretizations, a mixed three-field approach is considered. Based on the three-field potential for small strains given in Shu et al., [27,28] we propose an associated formulation for finite strains as Here, in the gradient enriched part of , ∇ is replaced by the new independent variable ∈ 1 (; I R × ) and compatibility is enforced in an integral sense via Lagrange multiplier ∈ 2 (; I R × ). The introduction of enables a relaxation of the continuity requirement to 0 . We introduce where is known as the first Piola Kirchhoff stress tensor and denotes a higher order stress tensor. Consequently, the variational equation is obtained. For the mathematical analysis we consider the small strain setting and investigate the analogous formulation to (3.3) to (3.5) with  = Γ = Γ for the special case of quadratic local and nonlocal strain energies. Given ∈ 2 (; I R 2 ) find ∈ 2 0 (; R ) with Herein, denotes a material parameter associated with the higher-order stress response, C is the elasticity tensor of fourth order for linear elasticity, and denotes the linear strain tensor. Thereby, an energy function is considered which is additively split into a local (classical) quadratic part depending on and a nonlocal quadratic part depending on ∇ . The corresponding mixed formulation seeks ( , , where the bilinear forms and are defined by

MATHEMATICAL ANALYSIS OF THE CONTINUOUS FORMULATION
This section analyzes the problem (3.7). It proves in Lemma 1 the equivalence of this problem with (3.6) in the sense that the solutions of the two problems coincide. However, problem (3.7) is not stable with respect to the standard norms ‖(•)‖ 1 () , ‖(•)‖ 1 () and ‖(•)‖ 2 () for the three variables , , , respectively. This means that a small perturbation of the right-hand side, that is, the external potential, could lead to large deviations of the solutions in those norms. In particular, errors from, for example, numerical integration of the right-hand side could lead to large errors of the solution.
Since the bilinear form in (3.6) is continuous and coercive due to Poincaré's inequality, Lax-Milgram yields existence of a unique solution. Poincaré's inequality applied on (∇ , ∇ ) 2 () and a Korn's inequality applied on (C ( ), ( )) 2 () show that the bilinear form is coercive. Since the bilinear forms and are continuous, it would remain to proof an inf-sup condition to show that problem (3.7) is stable.
However, Proposition 1 shows that the inf-sup condition for the bilinear form (3.9) is not satisfied. Therefore the formulation (3.7) is not stable with respect to the standard norms and Brezzi's splitting lemma does not guarantee the unique solvability. .

Lemma 1 If is the solution of (3.6), there exists
, then is in 2 0 (; I R ) and a solution to (3.6).

FINITE ELEMENT FORMULATION
In the previous section it was shown that for the linear gradient elastic scenario at small strains stability with respect to the standard norms is not guaranteed. Since the investigation of continuity and coercivity of the nonlinear finite strain gradient elasticity problem (2.5) as well as applicability of Brezzi's splitting lemma to the corresponding mixed problem (5.2) is not straightforward, analogous mathematical analysis of the constitutively and geometrically nonlinear formulation is left for future research. Nevertheless, similar problems with respect to stability are expected in the nonlinear regime. In order to investigate this, numerical analysis will be used. Thus, several large strain three-field finite element formulations are introduced in the following.

Discretization
In the following the general procedure for constructing the discrete system of equations is presented. The derivation of the corresponding matrices is here based on automatic differentiation (AD), which is increasingly used in finite element research as the development time is accelerated. [29] For the finite element implementation and the numerical example computations, the software package AceGen/AceFEM has been used. Note that usage of AD is not required for the formulations discussed here and only applied as a technical tool. Starting point for the discretization is the total potential (3.1). From the element partition of the continuum body it follows Subsequently, the discrete element potential takes the form where the solution variables are discretized with the continuous Galerkin method. They are based in finite dimensional subspaces  ℎ ⊆ ,  ℎ ⊆  and  ℎ ⊆  and can be written in matrix notation as where , and are standard matrices containing polynomial shape functions and , and are vectors of the corresponding degrees of freedom associated with ℎ , ℎ , and ℎ , respectively. The individual element degrees of freedom can be collected in the complete element vector of degrees of freedom of the discrete system as The variation of (5.2) can then be written as A detailed description of the individual residual components can be found in Appendix 7. For the numerical solution of the nonlinear problem A = 0 using the Newton-Raphson procedure, the discrete variational equation (5.5) is linearized: Herein,̄represents the unknowns in the previous Newton iteration step, Δ depicts the linear increment of the nodal degrees of freedom and A is an appropriate assembly operator. Consequently, ∕ represents the element tangent matrix Again, a detailed description of the individual components of the tangent matrix is given in Appendix 7. Making use of the distributive property of the assembly operator yields the global system of equations All integrals are numerically evaluated with Gauss point integration over the isoparametric reference element. Note that with = ( ) = 0 and = 0 the global tangent matrix takes the form (5.10)

Finite element approximations
In this section exemplary finite element approximations are investigated. For the discrete system, we consider partitions  of  into simplices. We define , (), ( ) to be the edges, the inner edges and the boundary edges of  , respectively. We further define  ,  (),  ( ) to be the nodes, the inner nodes and the boundary nodes of  . We denote  Table 1 shows an overview of the finite element pairs considered in the following.

Unstable discretizations
This section discusses several discretizations which, at first glance, seem to be natural choices. However, it will be shown that they are in fact not stable with respect to the norms considered in Section 4. Therefore, Brezzi's splitting lemma [30] proves that a unique discrete solution does not exist in general. Lemma 2 proves that the P1-P1-P0 discretization in 2D is not stable.

Existence of a discrete solution and influence of the domain size
This section proves stability of a P2-P1B-P0 discretization (P1B denoting the space of P1 functions enriched by volume bubbles) for = 2 under a condition on the size of . The stability is proven with respect to a nonstandard mesh-dependent norm of . This implies that a unique discrete solution exists. Note that convergence of the discrete solution does not follow from Brezzi's splitting theorem due to the missing stability of the continuous problem (3.7).
Recall the definition of the Crouzeix-Raviart space CR 1 ( ; I R 2 ) from Lemma 2 and define

Remark 3
While infsup,Curl does only depend on the shape of , the Poincaré constant P depends on the diameter of . Therefore, the inequality 2 infsup,Curl − P > 0 is satisfied if the diameter of  is sufficiently small.

Remark 4 Lemma 5 and Proposition 2 imply the mesh-dependent inf-sup condition
This mesh-dependency can also be seen in the numerically computed singular values of the tangent matrix in Figure 1.
Further elements, which are covered in the numerical investigation of the next section are denoted in Table 2. The notation (; I R ) is introduced for polynomials corresponding to th order Lagrange quadrilaterals for = 2 and hexahedrals for = 3.

NUMERICAL EXAMPLES
In this section the proposed discretizations are numerically evaluated with regard to the finite element convergence as well as to stability. In order to keep the focus on the finite element discretization a comparatively simple material model is considered. The gradient enriched strain energy density consists of an additively decomposed local and nonlocal part Here, we choose a classical compressible Neo-Hooke energy for loc ( ) and thus, which is quadratic in ∇ . The frame indifference and isotropy of nloc was proven in Triantafyllidis and Aifantis. [22] The constitutive parameter ≥ 0 of the gradient enrichment can be linked to a microstructural length-parameter with = √ ∕ . [3,21] For the following numerical examples the nonlocal length-parameter is set to = 0.1 mm.

Unit square test
In order to analyze the convergence behavior for the case = 2, a unit square domain  with the dimensions 1 × 1 mm 2 (cf. Figure 1A) is considered. The constitutive model is implemented under the plain strain assumption. For the construction of an error estimator the displacement field is predefined as reference solution by which is obtained by variation of (2.5). By imposing as external loading to the discrete boundary value problem (5.9), the finite element solution is compared to the P2 interpolation of the reference solution in the relative L 2-norm, that is, = ‖I 2 ( ) − ℎ ‖ 2 () ∕‖I 2 ( )‖ 2 () with P2 interpolation I 2 ( ) of the exact solution. It can be seen from Figure 1B that for uniform mesh refinement convergence order of ℎ 2 is reached for the elements P2-P1B-P0 and P2-P1B-P1, whereas the order ℎ 4 is obtained by the element P2-P2-P1, with h being the element diameter. Since the P2 interpolation I 2 ( ) converges with order ℎ 3 to the exact solution, the convergence rates of the discrete solutions to the exact solution is ℎ 2 for the elements P2-P1B-P0 and P2-P1B-P1 and ℎ 3 for the element P2-P2-P1. Furthermore, the same convergence rates are obtained when measuring the error in the relative H1-seminorm = ‖I 2 (( ) − ℎ ‖ 2 () ‖I 2 (( )‖ 2 () . However, the element P2-P1B-P0 has evoked an interruption of the linear solver at the refinement stage corresponding to ℎ = 0.011 mm such that the last possible solution was obtained for the discretization marked as bullet in Figure 1B. This observation is in correspondence to Remark 3, stating that the stability of this element depends on the size of the computational domain . The fact, that this additional stability issue is uncovered already on this simple unit cube domain makes the P2-P1B-P0 element an unsuitable choice of discretization. To analyze this further, Figure 1C unveils an ℎ 2 -degeneration of minimal singular values of the submatrix = [ , ]. This corresponds to the considerations in Section 5.2 and may ultimately lead to a loss of invertibility of the global tangent matrix due to close to zero minimal singular values. For the performed refinement stages, yet, this has not occurred for elements Q2-Q2-Q1, P2-P1B-P1, and P2-P2-P1. Note that the results for elements P2-P1-P0 and Q2-Q1-P0 are not depicted because zero minimal singular values of the global tangent matrix occurred at all refinement stages precluding the solution of the linearized system of equations.

Unit cube test
To investigate the three dimensional elements corresponding to the considered approximations, the unit square problem of Section 6.1 is adapted to a three-dimensional unit cube with the dimensions 1 × 1 × 1 mm 3 . Consequently, the displacement field (6.4) is adjusted with , with ( ) = ( − 1)( − 1)( − 1) (6.6) and = −1000 mm −18 leading to a maximal reference displacement of (0.5, 0.5, 0.5) = −0.2441 mm. Convergence characteristics (cf. Figure 2B) of all investigated discretizations are in accordance to the two-dimensional case. Similar to Section 6.1, the element P2-P1B-P0 evokes an interruption of the linear solver at a mesh refinement stage corresponding to ℎ = 0.108 mm, which may again be in correspondence to Remark 3. While the minimal singular values of the tangent matrix of all elements degenerate as in the 2D-case, a loss of invertiblity of the global linearized system of equations corresponding to the P2-P1B-P1, P2-P2-P1, and Q2-Q2-Q1 discretization has again not been observed. It is notable, that the tangent matrix of the element Q2-Q1-P0, which is singular at all refinement stages, can be numerically stabilized by superposing it with a diagonal matrix 1 of the same dimension, with ∈ I R >0 and ≪ 1. [28] This may however induce an error, since the discrete system does not correspond to the initial formulation anymore. In Figure 2A, the relative displacement = ℎ (̄)∕ (̄) at the center point = (0.5, 0.5, 0.5) is plotted over the total computing time needed for the solution procedure. Summarizing, the elements Q2-Q2-Q1 and P2-P2-P1 have a favorable convergence behavior in this particular boundary value problem.

3D Cook's problem
In the 3D Cook's problem ( Figure 3A) the body is clamped on one side and loaded with surface traction on the other side. Thus, the assumption of essential boundary conditions applied to the whole boundary, made at the end of Section 3, is modified to | Γ = 0 and Γ = ∅ with a surface traction = (0, 0, 6.25 MPa) . In order to compare the convergence behavior of the elements, in Figure 3B the relative displacement in -direction, measured at the corner point = (48, 20, 60) mm, is plotted over the number of degrees of freedom. Since an analytical reference solution is not available for this problem, the reference displacement = 8.76 mm depicts the numerical solution at a high mesh refinement stage, in which the solutions of all considered elements are in concordance within an accuracy of three significant figures. In order to show the ability to avoid mesh dependent solutions due to singularities the nonlocal formulation using a P2-P2-P1 element is compared to a classical local formulation using a quadratic tetrahedral, that is, P2 ⋅ ⋅ = { ℎ = 2 ( ; I R 3 )}. The latter formulation is obtained if the terms In Figure 3D the absolute values of at the strain concentration point = (0, 0, 44) are plotted over the element size according to the mesh refinement stage. In order to pertain comparability is evaluated as postprocessing quantity of ℎ for both elements despite the fact that for element P2-P2-P1 could be directly evaluated from ℎ . It becomes evident that as ℎ decreases, the -values of the P2-element diverge to infinite absolute values, while those of the P2-P2-P1 element are bounded. This illustrates one of the advantages of the gradient elasticity approach over classical formulations also at finite strains.

CONCLUSION
The three-field mixed formulation for gradient elasticity (cf. Shu et al. [27,28] for the small strain framework) was generalized to the finite strain hyperelastic setting. Mathematical analysis in the linear framework of the continuous formulation has unveiled instability with respect to the standard norms and Brezzi's splitting lemma. Moreover, for the finite element space P1-P1-P0 in 2D and P2-P1-P0 in 3D the discrete inf-sup condition does not hold. Although instability of the linear formulation most likely implies instability of nonlinear formulations, the mathematical stability analysis is still missing and open problem for future research. However, three-field formulations similar to the ones investigated here have been proposed in the literature not only in the context of gradient elasticity at small strains. [27,28] Also for gradient damage at large strains similar mixed formulations have been proposed, [35,36] where any mathematical analysis (also of the continuous problem) is still missing. Therefore, a comprehensive numerical study on all elements has been conducted at finite strains in order to show potential problems resulting from the expected instabilities at large strains, which were indeed numerically verified. Interestingly, numerical experiments have however also shown that for some finite strain elements (eg, P2-P2-P1 and Q2-Q2-Q1) discrete instabilities have not occurred and finite element convergence could be observed. In addition, the numerically converged solutions of Section 6.1 and 6.2 are in concordance with the exact analytical solution. Also, the displacement convergence plot of the 3D Cook's problem of Section 6.3 unveils concurrent numerically converged solutions of the elements Q2-Q2-Q1, P2-P2-P1, P2-P1B-P1 and P2-P1B-P0. In summary, numerical convergence to different solutions could not be observed. This is mostly due to the fact, that although their quadratic order of decreasing minimal singular values with decreasing element size is unsatisfying (cf. Figure 1C), still suitable numerical results may be obtained as long as the minimal singular values do not approach zero before finite element convergence is reached. An overview of the observations regarding the investigated elements is given in Table 3. However, the fact that individual calculations converging to the solution were obtained even for the potentially unstable formulations should not be mistaken as indicator for a suitable FE formulation, since it is in general unclear if the converged solution is the solution to the problem. A mesh-sensitivity test on the 3D Cook's problem has shown the advantage of the gradient elasticity approach over a classical elasticity element even at finite strains. This motivates further research on robust and efficient finite element schemes, which avoid the observed instability issues.

APPENDIX A
Define the following matrix dimensions as ∈ {2, 3}, = dim , ℎ= dim , and = dim . (A.1) A suitable matrix notation for the derivatives is given by wherein ∈ I R 2 × and ∈ I R 3 ×ℎ are appropriate B-matrices including the first order derivatives of the shape functions belonging to and . Using these matrices, the analogous expressions for the variations and increments of the degrees of freedom are obtained. Let , , C loc , and C nloc denote analogous matrix notations for It becomes evident that the element tangent matrix is symmetrical.