A New C0‐Continuous FE‐Formulation for Finite Gradient Elasticity

Three field mixed C0‐continuous gradient elasticity finite element formulations contain both displacement and displacement gradients as solution variables. In this contribution a new formulation is proposed, in which the displacements are not part of the problem anymore, but only displacement gradients, leading to a reduced number of variables to be discretized. In numerical examples the proposed formulations are compared to finite strain adaptations of existing mixed three field approaches for small strains and unveil favorable computing time. For the 3D case a further reduction of variables through a penalty formulation is investigated.


Introduction
Classical elasticity formulations incorporate elastic energy functions, which depend on the first-order gradient of deformation F = ∇ϕ = 1 + ∇u, where ϕ : B → S maps points X ∈ B of the considered body in the reference configuration to points x ∈ S in the deformed configuration, ∇(•) = ∂(•)/∂X is the gradient operator and u is the displacement function. Although well established and sufficient for many applications, these formulations have limitations when it comes to specific problems. That is on the one hand the occurance of geometry-induced, non-physical stress singularities at e.g. sharp corners and on the other hand the lack of ability to capture the influence of the specimen size on the elastic behavior when modeling very small structures. Moreover, when using classical elasticity formulations, domains containing phase transformations may lead to non-convex elastic energy functions. Gradient elasticity formulations can represent a remedy in these cases by adding a functional dependency of the elastic energy W on the second order gradient of deformation ∇F .
W · · = W (F ) (classical elasticity) ⇒ W · · = W (F , ∇F ) (gradient elasticity) (1.1) Through the enrichment (1.1) the solution of the corresponding problem is smoothened and consequently stress singularities and the loss of ellipticity can be avoided. Furthermore, additional constitutive parameters take into account the size dependency of the elastic response when modeling small structures. The theoretical foundation for the small strain gradient-enriched continuum theory dates back to the works of [1] and [2] and the more recent adaptation of [3]. A historical overview is given in [4] and an investigation of finite strain models can be found in [5] and [6]. When finding suitable numerical solution schemes such as e.g. finite element formulations, the main challenge of gradient elasticity models is coping with the higher order continuity requirement that appears through the second-order gradient enrichment. The corresponding model problem reads where the displacement solution is sought in the Sobolev space u ∈ H 2 (B) of functions, whose gradients up to second order are L 2 (B) integrable. This leads to the requirement of C 1 -continuous interpolation functions, for which standard finite elements are not suitable. Isogeometric analysis e.g. as considered by [7] can satisfy the continuity requirement but has the kown difficulties when modeling complex structures and boundary conditions. According to gradient elasticity theory [1] the boundary ∂B of the body is decomposed into the standard Dirichlet and Neumann boundary ∂B = Γ D ∪ Γ N and higher order Dirichlet and Neumann boundary ∂B = Γ H ∪ Γ M . The external potential then reads with f and t being standard volume loads and surface tractions respectively, r being higher order tractions and n being the outwards unit normal vector. Alongside the essential boundary condition u =ū on Γ D on the higher order Dirichlet boundary displacement gradients are prescribed in normal direction with ∇u · n =h.

Finite element formulation 2.1 Three field mixed approach
In order to relax the continuity requirement for the discretization variables to C 0 and to make standard finite element approximations possible, a mixed three-field approach can be considered (c.f. [8] and [9] for the small strain framework). Following the large strain extension proposed in [10] the Lagrangian reads The idea is to replace ∇u in the gradient-enriched part of W by the new independent variable H ∈ H 1 (B) and enforce compatibility via the Lagrange multiplier λ. Through the introduction of H only first order derivatives appear in the corresponding weak formulation and thus, the continuity requirement is relaxed to C 0 . For the dimensions d ∈ {2, 3} (2D and 3D), the continuous function spaces of the solution variables are defined as follows [10]: In what follows for discrete functions u h and H h , the second order Lagrange approximation will be considered while for λ h first order Lagrange approximation is chosen (c.f. [8,10]). The corresponding finite elements will be denoted as P2 u -P2 H -P1 λ and Q2 u -Q2 H -Q1 λ for simplical and quadrilateral or hexahedral elements respectively.

Reformulation of the elastic potential
In order to reduce the number of degrees of freedom and to avoid possible stability issues, which may occur with the three field mixed approach (c.f. [10]), the displacement u can be decoupled from the problem in the spirit of [11,12]. The reformulation of (2.1) now reads making H and λ the only solution variables, while retaining the C 0 -continuity requirement. Moreover, H is sought in the standard solution space H ∈ G. The nature of vanishing rotations of gradient fields for simply connected domains with boundary conditions on connected subsets of ∂B is employed as constraint to make sure H is indeed a gradient field. It is notable that, since all stress and strain related quantities are measures of the displacement gradient, associated postprocessing contour plots can be directly obtained by evaluating the nodal solutions of (2.5). This is in contrast to classical elasticity approaches, where only the displacements (not the stresses/strains) are available as continous fields. Furthermore, in the 2D case (2.5) is equivalent to a Stokes-type problem [11] and the Lagrange multiplier is sought in the standard solution space λ ∈ Q. Suitable corresponding finite element approximations for this are well kown (see e.g. [13]) and in the following for H h and λ h the Taylor Hood type finite element discretizations P2 H -P1 λ and Q2 H -Q1 λ are considered. For solution spaces for the Lagrange multiplier and suitable discretizations in the 3D case the reader is referred to [11]. In order to determine the work quantity ∇g a preprocessing boundary value problem of the simple Laplace type is solved before the main problem, i.e.
For the discretization g h in (2.6) the straightforward H 1 -conforming Lagrange finite element interpolation can be employed. For plotting purposes, the displacements u can be determined from solving the postprocessing problem where H is the solution of (2.5). Problem (2.7) is again of the Laplace type for which the standard Lagrange discretization can be used. In what follows, for (2.6) and (2.7) second order Lagrange polynomials are chosen.

Penalty formulation
To further reduce the number of solution variables and consequently the number of degrees of freedom a straightforward penalty approach can be applied. The integrand of the constraint term associated with the Lagrange multiplier of (2.5) is then replaced by a penalty term 1 2 (∇ × H) 2 with 1. Consequently, H is the only discretization variable in this approach. The corresponding finite elements are named P2 H -pen for simplical meshes and Q2 H -pen for quadrilateral and hexahedral meshes. In a similar manner the integrand of the constraint term of (2.1) can be replaced by 1 2 (H − ∇u) 2 and corresponding finite elements are denoted by P2 u -P2 H -pen and Q2 u -Q2 H -pen. PAMM · Proc. Appl. Math. Mech. 19 (2019) Q2u-Q2H -Q1λ

Numerical examples
In this section the proposed finite element formulations are numerically evaluated with respect to convergence behavior and numerical efficiency. For this, the gradient-enriched strain energy W (F , ∇F ) = W loc (F ) + W nloc (∇F ) is chosen, in which W loc is a compressible Neo-Hooke energy function (c.f [14]) and W nloc (∇F ) = κ 2 ∇F · ∇F is a gradient enrichment term (c.f. [6]) with the additional constitutive parameter κ ≥ 0. In what follows the elasticity constants E = 500 MPa, ν = 0.3 and κ = 1.923 N are chosen.

Unit square convergence test
The first numerical example is a unit square domain B = (0, 1) × (0, 1) mm 2 for which the displacement solution field u = (u X , u Y ) T with u Y = a(XY (X − 1)(Y − 1)) 2 and u X = 0 is chosen. For this, the essential boundary conditions u = 0 and ∇u = 0 on Γ D = Γ H = ∂B hold. Note, that a = −50 mm −8 is a scaling factor so that the order of magnitude of the displacement field is within the nonlinear elastic range. The volume load f is computed analytically via the strong form obtained by variation of (1.3) under plain strain assumption. Error measures are then constructed from the finite element solution u h and H h obtained through evaluation of (2.1) through (2.7) and the analytical solution. As depicted in figure 1a,b convergence of order h 4 is observed for all elements both in the L 2 -norm of the displacement error and in the L 2 -norm of the displacement gradient error. While convergence rates are almost identical, the computing time (figure 1c) of the P2 H -P1 λ and Q2 H -Q1 λ elements is lower than the computing time of the P2 u -P2 H -P1 λ and Q2 u -Q2 H -Q1 λ elements at all refinement stages due to the reduced number of variables. Moreover, when adding the time of the pre and postprocessing problem (2.6) and (2.7) to the computing time of (2.5) (denoted by Q2 H -Q1 * λ ) the computing time only changes marginally.

Cook's membrane
In this numerical example, the proposed elements are tested on the Cook's membrane benchmark problem (figure 2a). As the domain is clamped at X = 0 mm and loaded at X = 48 mm the Dirichlet and Neumann boundary data is given with P2u-P2H -pen6

Unit cube convergence test
Similar to section 3.1 a unit cube domain B = (0, 1) 2 , u X = u Y = 0 and a = −1000 mm −18 . While the Q2 H -pen6, Q2 u -Q2 H -pen6 and P2 u -P2 H -pen6 show an almost identical convergence behavior for a penalty parameter 1/ = 10 6 , for the element P2 H -pen3 the choice 1/ = 10 3 yields suitable results. Again the elements derived from section 2.2 require a lower computing time than the elements derived from the three field mixed approach (figure 3c).

Conclusion
A C 0 continuous finite element formulation for gradient elasticity at finite strains with a reduced number of solution variables has been proposed and compared to the mixed three field approach as investigated in [10]. Numerical results in 2D show appropriate convergence rates and improved computing efficiency of the proposed formulations P2 H -P1 λ and Q2 H -Q1 λ using simple standard discretizations. In the 3D case the penalty approach, which further reduces the number of solution variables whilst enabling standard discretization, yields converging numerical results and again improved efficiency compared to the mixed three field approach. Nevertheless a critical view towards the correct choice of the magnitude of the penalty parameter and the possibility of badly conditioned systems gives motivation for further investigation.