Lagrange and H(curl,B)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H({\text {curl}},{{\mathcal {B}}})$$\end{document} based finite element formulations for the relaxed micromorphic model

Modeling the unusual mechanical properties of metamaterials is a challenging topic for the mechanics community and enriched continuum theories are promising computational tools for such materials. The so-called relaxed micromorphic model has shown many advantages in this field. In this contribution, we present significant aspects related to the relaxed micromorphic model realization with the finite element method (FEM). The variational problem is derived and different FEM-formulations for the two-dimensional case are presented. These are a nodal standard formulation H1(B)×H1(B)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H^1({{\mathcal {B}}}) \times H^1({{\mathcal {B}}})$$\end{document} and a nodal-edge formulation H1(B)×H(curl,B)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H^1({{\mathcal {B}}}) \times H({\text {curl}}, {{\mathcal {B}}})$$\end{document}, where the latter employs the Nédélec space. In this framework, the implementation of higher-order Nédélec elements is not trivial and requires some technicalities which are demonstrated. We discuss the computational convergence behavior of Lagrange-type and tangential-conforming finite element discretizations. Moreover, we analyze the characteristic length effect on the different components of the model and reveal how the size-effect property is captured via this characteristic length parameter.


. Introduction
Metamaterials are receiving tremendous attention in both academia and industry due to their unconventional mechanical properties.These are not solely governed by the bulk mechanical properties but also by the geometry of the unit cells which can be designed to attain the desired functionality, see [14,19,52,54,55].Moreover, the recent advances of the additive manufacturing (AM, or 3D printing) techniques are empowering the fabrication process of three-dimensional architected metamaterials, c.f. [16,20,32,42].To simplify the design process of novel metamaterials, suitable computational tools are needed to capture their unprecedented effective mechanical properties.The classical Cauchy-Boltzmann theory and the first-order homogenization procedures often fail to describe the mechanical macroscopic behavior of mechanical metamaterials since they exhibit the size-effect phenomenon, i.e. small specimens are stiffer than big specimens, and therefore other generalized theories are needed such as the classical Mindlin-Eringen micromorphic theory [11,12,21,22,29,51], the Cosserat theory [7,35,36], gradient elasticity [2,13,30] or others.

arXiv:2112.00382v1 [math.NA] 1 Dec 2021
The relaxed micromorphic model, which we adopt in this work, has been introduced recently in [15,37,38].It keeps the full kinematics of the micromorphic theory but employs the matrix Curl operator of a non-symmetric second-order micro-distortion field for the curvature measurement.The relaxed micromorphic model reduces the complexity of the classical micromorphic theory by decreasing the number of material parameters and has shown many advantages such as the separation of the material parameters into scaledependent and scale-independent ones, see for example [9].Furthermore, it has already been used to obtain the main mechanical characteristics (stiffness, anisotropy, dispersion) of the targeted metamaterials for many well-posed dynamical and statical problems, e.g.[1,3,4,23,24,25,26,27,28,39].Recently, the scale-independent short range elastic parameters in the relaxed micromorphic model were determined for artificial periodic microstructures in [40] which are used to capture the band-gaps as a dynamical property of mechanical metamaterial in [9].Analytical solutions of the relaxed micromorphic compared to the solutions of other generalized continua for some essential boundary value problems, i.e. pure shear, bending, torsion and uniaxial tension, are discussed in [44,45,46,47], emphasizing the validity of the relaxed micromorphic model for small sizes (bounded stiffness), where most of the other generalized continua exhibit unphysical stiffness properties.
As a result of employing the matrix Curl operator of the micro-distortion field for the curvature measurement, the relaxed micromorphic model seeks the solution of the microdistortion in H(curl, B), while the displacement solution is still in H 1 (B).The appropriate finite elements of such case must be conforming in H(curl, B) (tangentially conforming).The first formulation of edge elements was presented in [43].In fact, the name "edge" elements was used because the degree of freedoms (dofs) are associated only with edges for the first-order approximation.H(curl, B)-conforming finite elements of first kind were introduced in [33] and second kind in [34], which are comparable with H(div, B)-conforming elements of first kind in [43] and second kind in [6].An extension to elements with curved edges, based on covariant projections, was developed by [8].A general implementation of Nédélec elements of first kind is presented in [41] and a detailed review about H(div, B)and H(curl, B)-conforming finite elements is available in [17] and [48].Furthermore, hierarchical H(curl)-conforming finite elements are used to solve Maxwell boundary and eigenvalue problems in [49].A H 1 (B) × H(curl, B) finite element formulation for a simplified anti-plane shear case of the relaxed micromorphic model utilizing a scalar displacement field and a vectorial micro-distortion field is available in [50].
In our work, we demonstrate the main technologies related to the finite element realization of the theoretically-sound relaxed micromorphic model.The proper finite element approximation of the micro-distortion field is the Nédélec space which utilizes tangentialconforming vectorial shape functions.We provide a comprehensive description of the construction of H 1 (B) × H(curl, B) elements with Nédélec formulation of first kind on triangular and quadrilateral meshes.Six finite elements are built, which are different in the approximation space of the micro-distortion: two triangular elements with firstand second-order Nédélec formulation, two quadrilateral elements with first-and secondorder Nédélec formulation, and two nodal triangular elements with standard first-and second-order Lagrangian formulation.The paper is organized as follows.In Section 2, we introduce the relaxed micromorphic model and derive the variational problem with the resulting strong forms and the associated boundary conditions which are modulated in a physical point of view by the so-called consistent coupling condition.We cover in Section 3 the main components of the implementation of standard nodal and nodal-edge elements.Two numerical example are introduced in Section 3. The first numerical example is designed to check the convergence behavior of the different finite elements when the solution is discontinuous in the micro-distortion field.We investigate the influence of the characteristic length in a second example which covers the size-effect property.We conclude the paper in Section 5.

. The relaxed micromorphic model
The relaxed micromorphic model is a continuum model which describes the kinematics of a material point using a displacement vector u : B ⊆ R 3 → R 3 and a non-symmetric micro-distortion field P : B ⊆ R 3 → R 3×3 .Both are defined for the static case by the minimization of potential with (u, P ) ∈ H 1 (B) × H(curl, B).The vector f and the tensor M describe, respectively, the given body force and the body moment, while t is the traction vector acting on the boundary ∂B t ⊂ ∂B.The elastic energy density W reads (2) Here, C micro , C e are fourth-order positive definite standard elasticity tensors, C c is a fourthorder positive semi-definite rotational coupling tensor, L is a fourth-order tensor, L c is a non-negative parameter describing the characteristic length scale and µ is a typical shear modulus.The characteristic length parameter plays a significant role in the relaxed micromorphic model.This parameter is related to the size of the microstructure and determines its influence on the macroscopic mechanical behavior.A relation of the relaxed micromorphic model to the classical Cauchy model was shown in [40] for limiting values of L c , which we can also observe in our numerical examples.
The variation of the potential with respect to the displacement field, i.e. δ u Π = 0, with leads after integration by parts and applying the divergence theorem to the weak form where σ denotes the non-symmetric force stress tensor.The associated strong form with the related boundary conditions reads satisfying ∂B u ∩∂B t = ∅ and ∂B u ∪∂B t = ∂B and n is the outward normal on the boundary.In a similar way, the variation of the potential with respect to the micro-distortion field, i.e. δ P Π = 0, with yields after integration by parts and applying Stokes' theorem where σ micro and m are the micro-and moment stresses, respectively, and m i and δP i denote the row vectors of the associated second-order tensors.Using the identity of the scalar triple vector product allows for the reformulation The associated strong form reads with the related boundary conditions where ∂B P ∩∂B m = ∅ and ∂B P ∪∂B m = ∂B.A dependency between the displacement field and the micro-distortion field on the boundary was proposed by [40] and subsequently considered in [10,45,46,50].This so-called consistent coupling condition is defined by where τ is a tangential vector on the Dirichlet boundary.This condition relates the projection of the displacement gradient on the tangential plane of the boundary to the respective parts of the micro-distortion.
The first strong form in Equation ( 5) represents a generalized balance of linear momentum (force balance) while the second strong form in Equation ( 10) outlines a generalized balance of angular momentum (moment balance).The generalized moment balance invokes the Cosserat theory with the Curl Curl operator rising from the matrix Curl operator of the second-order moment stress m.In comparison to the classical micromorphic model, see [11,37], the relaxed micromorphic model uses the same kinematical measures but employs a curvature measure form the Cosserat theory, see [36].The micro-distortion field has the following general form for the three-dimensional case where P i are the row vectors of P .We let the Curl operator act on the row vectors of the micro-distortion field P , i.e., For the two-dimensional case, the micro-distortion field and its Curl operator are reduced to 3 .Approximation spaces We introduce the formulation of a standard nodal element utilizing Lagrange-type shape functions for both displacement and micro-distortion field, see for example [53].Let us assume that there are n u nodes in each element for the discretization of the displacement field u and n P nodes for micro-distortion field P in two dimensions.Geometry and displacement field are approximated employing the related Lagrangian shape functions N u I defined in the parameter space with the natural coordinates ξ = {ξ, η} by where X I are the coordinates of the displacement node I and d u I are its displacement degrees of freedom.The deformation gradient is obtained in the physical space by where J = ∂X ∂ξ is the Jacobian, ∇ and ∇ ξ denote the gradient operators with respect to X and ξ, respectively.The micro-distortion field P for the 2D case is approximated using the relevant scalar shape functions where d P 1 I and d P 2 I are the micro-distortion row vectors degrees of freedom of node I.In order to calculate the Curl of P , the gradient of the row vectors in the physical space can be calculated by and the rotation of the vector P i h reads curl 3.2.Nodal-edge elements (u, The here presented formulation uses different spaces to describe the micro-distortion field. The geometry and the displacement field are approximated in the standard Lagrange space as in Equation ( 16).For the micro-distortion field, its solution is in H(curl, B) and the suitable finite element space is known as Nédélec space, see [33,34].In this work, we choose the Nédélec space of first kind.For more details the reader is referred to [5,17,31,48].Nédélec formulations use vectorial shape functions which satisfy the tangential continuity at element interfaces.The lowest-order two-dimensional Nédélec elements are depicted in Figure 1.(1,-1) where IP k−1 is the linear space of polynomials of degree k − 1 or less and Ĩ P k is the linear space of homogeneous polynomials of degree k.Equivalently, the space can be characterized by The dimension of this linear space is k(k + 2).Quadrilateral Nédélec elements of order k are based on the linear space with dim N D 2 k = 2k(k + 1).The vectorial shape functions v k in the parametric space are obtained by constructing a linear system of equations based on a set of inner and outer dofs.For the 2D case, the outer dofs of an edge e i are defined by the integral where r j is a polynomial IP k−1 along the edge e i and t i is the normalized tangential vector of the edge e i .The inner dofs are introduced for triangular elements by while they are given for quadrilateral elements by The scalar-valued and vectorial functions r j and q i are linearly independent polynomials which are chosen as Lagrange polynomials in our work.For lowest-order element (k = 1), only outer dofs occur.For higher-order elements (k ≥ 2), the number of outer dofs are increased and additional inner dofs are introduced.E.g. for the N D 2 2 with a dimension 8, we have 6 outer dofs and 2 inner ones.The derivations of the H(curl, B)-conforming vectorial shape functions is shown in Appendix A. Mapping the vectorial shape functions v k I from the parametric space to ψk I in the physical space must conserve the tangential continuity property.This is guaranteed by using the covariant Piola transformation, see for example [48], which reads For our implementation of the H(curl, B)-conforming elements, we modify the mapping to enforce the required orientation of the degrees of freedom at the inter-element boundaries and to attach a direct physical interpretation to the Neumann-type boundary conditions.Hence, two additional parameters, α and β, appear for the vectorial shape functions associated with edge dofs where α I = ±1 is the orientation consistency function which ensures that on an edge, belonging to two neighboring finite elements, a positive tangential flux direction is defined.Therefore, a positive tangential direction is defined based on a positive x-coordinate.A tangential component pointing in negative x-direction is multiplied by a value α I = −1 to obtain the overall positive tangential flux on each edge.If the tangential direction has no projection on x-axis, then the same procedure is employed on y-direction.Figure 2 illustrates an example of calculating the orientation parameter values of two neighboring elements.
The normalization parameter β I enforces that the sum of the vectorial shape functions ψ k I at a common edge scalar multiplied with the associated tangential vector has to be  equal one in the physical space.Furthermore, the sum of the shape functions belonging to one edge scalar multiplied with the tangential vector of the other edges must vanish.These conditions are reflected by Here, J ψ k is the sum of shape vectors related to outer dofs of an edge E J evaluated on the edge E I and τ I is the normalized tangential vector of an edge E I where E denotes the edges in the physical space.Based on Equations (28) 1 and (29) 1 , we compute straightforward the parameters β I .In detail we get for the first-and second-order elements respectively, where L I denotes the length of the edge E I in the physical space.For the 2D case, the rotation of the vectorial shape functions only has one active component out of the plane which reads The micro-distortion field P is approximated by the vectorial dofs d P I representing its tangential components at the location I = 1, ..., n P .The micro-distortion field and its Curl are interpolated as The non-vanishing components of the Curl operator of the micro-distortion field for the 2D case are obtained by

Implemented finite elements
In this work, we present four nodal-edge elements based on the formulation in Section 3.2 and two standard nodal elements based on Section 3.1.All implemented finite elements employ scalar quadratic shape functions of Lagrange-type for the displacement field approximation with the notation T2 for triangles and Q2 for quadrilaterals.The micro-distortion field is approximated using different formulations introduced in Sections 3.1 and 3.2.For the standard nodal elements, Lagrange-type ansatz functions are used resulting in the element types T2T1 (linear ansatz for P ) and T2T2 (quadratic ansatz for P ).Different nodal-edge elements are built utilizing first-and second-order Nédélec formulations with tangential-conforming shape functions denoted as NT1 and NT2 for triangular elements and QT1 and QT2 for quadrilateral elements.The micro-distortion dofs in the standard nodal elements T2T1 and T2T2 are tensorial with 2 × 2 entries while the nodal-edge elements use vectorial dofs for the micro-distortion field which represent the tangential components.The full micro-distortion tensor is restored based on Equation (32).The used finite elements are depicted in the parameter space in Figure 3.

Numerical examples
For our numerical examples, we neglect the body forces and moments, i.e. f = 0 and M = 0, and assume isotropic material behavior which can be described by the set of material parameters λ micro , µ micro , λ e , µ e , µ c , µ and L c , where λ * and µ * denote the Lamé coefficients.Furthermore, we assume L as a fourth-order identity tensor II.Throughout all examples, we consider the Cosserat modulus µ c = 0, cf.[35,38], leading to the symmetry of the force stress tensor.The simulations presented in this paper are performed within AceGen and AceFEM programs, which are developed and maintained by Jože Korelc (University of Ljubljana).The interested reader is referred to [18].

Discontinuous solution: an interface between two different materials
In the first boundary value problem (bvp), we consider a rectangular domain B with length l = 2 and height h = 1 which consists of a side-by-side arrangement of two different materials, see Figure 4.The bottom edge is fixed in both directions, u = (0, 0) T , and we apply a displacement to the upper edge, u = (0.01, 0.01) T , while the left and right edges are subjected to displacement given by u = (0.01y 2 , 0.01y 2 ) T .For the micro-distortion field, the consistent coupling boundary condition P • τ = ∇u • τ is enforced on the entire boundary ∂B.The material parameters are given in Table 1.Due to the different material parameters of the domains, a discontinuous solution in the micro-distortion is expected which allows us to compare and evaluate the behavior of the implemented finite elements.In Figure 5, the displacement and micro-distortion fields obtained for a discretization

Characteristic length analysis: pure shear problem
We introduce a second boundary value problem, see Figure 8, consisting of a circular domain B with a radius r o = 25 and a circular hole at its center with a radius r i = 2.We fix the displacement field u = 0 on the inner boundary ∂B i and we rotate the outer boundary ∂B o counter clockwise with ū = (− ∆ ro y, ∆ ro x) T where ∆ = 0.01.For the microdistortion field, we apply the consistent coupling boundary condition (P • τ = ∇u • τ ) on all boundaries ∂B = ∂B i ∪ ∂B o .Two different cases are discussed in the following.For case A, a single material is assumed whereas for case B two materials are considered.The second material is located as a ring with an outer radius r m = 10 and an inner radius r i = 2.The material parameters are shown in Table 2.For the analysis of the influence of the characteristic length L c , the characteristic length will be varied.The problem results in a rotationally-symmetric solution where only the shear components (u r,θ , u θ,r , P rθ , P θr = 0) are non-vanishing.The convergence behavior of the different elements is investigated for case B and L c = 5 using three different mesh densities (410, 3044 and 30620 triangular elements and 448, 3040 and 30256 quadrilateral elements).Since the micro-distortion field is in H(curl, B), the tangential shear component P rθ has to be continuous while the radial shear component P θr exhibits a jump, see Figure 9, where the Q2NQ1 element is used.Similar to Section 4.1, the H 1 (B) × H 1 (B) elements are unable to capture this discontinuity in P θr , which is shown in Figure 10.Actually, H 1 (B) × H 1 (B) elements approximate the discontinuous H(curl, B) solution only when a very fine discretization is used.The discontinuous solution of the micro-distortion field is demonstrated in Figure 11 using the H 1 (B) × H(curl, B) elements.The higher-order Nédélec formulation in T2NT2 and Q2NQ2 elements exhibit very satisfactory results already with the coarse mesh.The case L c → ∞ resembles an infinite zoom into the material, where an equivalence to linear elasticity with C micro can be derived, cf.[40].In the latter case, it can be shown that it holds that P = ∇u.In our study, we approximate the limiting cases by L c = 10 −3 and L c = 10 3 , respectively.Figure 12 shows the elastic energy W along the radius and Figure 13 illustrates the non-vanishing  components of P together with the respective displacement gradient components using 30624 Q2NQ2 elements.In Figure 14, we plot the total potential of the relaxed micromorphic model varying the characteristic length parameter L c .The figures clearly show the above described behavior.The bounded behavior of the relaxed micromorphic model for small sizes, i.e.L c → 0, is an important advantage which most other generalized continua miss.For the linear elasticity model, a standard T2 nodal element is implemented.17 is confined between the linear elasticity stress with elasticity tensor C micro from above and the one with C macro from below for large and small values of the characteristic length, respectively.Summarizing the previous findings shortly, increasing the characteristic length diminishes force stress and raises the microand moment stresses, while both force and moment stresses vanish for large and small values of the characteristic length, respectively, the micro-stress is always present.

. Conclusions
The relaxed micromorphic model is a generalized continuum model which can suitably reproduce the macroscopic effective properties of mechanical metamaterials.First, we derived the variational problem with the relevant weak and strong forms and the associated boundary conditions.We put together the main components of standard nodal and nodaledge finite element formulations of the relaxed micromorphic model.The standard nodal elements H 1 (B) × H 1 (B) are incapable to achieve satisfactory results for a discontinuous solution unlike H 1 (B) × H(curl, B) elements which capture the jumps of the normal components of the micro-distortion field and therefore converges efficiently.We reveal the role of the characteristic length which governs the scale-dependency property of the relaxed micromorphic model.For L c → 0, the model is equivalent to the standard Cauchy linear elasticity model with C macro defined as the Reuss lower-limit of elasticity tensors C e and C micro while the model is corresponding to Cauchy linear elasticity model with C micro with P = ∇u for L c → ∞.Furthermore, we have shown the dependency of the different stress measurements on the characteristic length.The force stress is at maximum for L c → 0 and it vanishes for L c → ∞ but the moment stress behaves in the opposite way.The micro-stress varies between Cauchy linear elasticity stresses with C micro and C macro for L c → ∞ and L c → 0, respectively.where a i , i = 1, .., 12 are coefficients yet to be defined based on the dofs.Starting from Equations ( 24) and ( 26), the explicit functions r j and q i are set as edge 1: The edge and inner dofs are calculated according to Equations ( 24) and ( 26) considering that the tangential vector and the coordinates coloration are same as in NQ1 element

Figure 1 :
Figure 1: Lowest-order (k = 1) Nédélec elements: triangle N D 2 1 (right) and quadrilateral N D 2 1 (left).Definition of the individual edges e i .The red arrows indicate the orientation of the tangential flux.

Figure 2 :
Figure 2: Example of assembling of two neighboring elements with satisfying the orientation consistency via the orientation parameter α I .

Figure 3 :
Figure 3: The implemented finite elements in the parameter space.Black dots represent the displacement nodes while red squares stand for micro-distortion field nodes associated with tensorial dofs.Red arrows and crosses indicate the edge and inner vectorial dofs, respectively, of the micro-distortion field used in Nédélec formulation.

Figure 4 :
Figure 4: 2D rectangular bvp with two different materials.Inspection line, see Figures 6and 7, can be seen in red color.
with 1600 N2NT2 finite elements are shown.Note that the elements solution is plotted in this work without the usual averaging or smoothing on the element edges in order to investigate the possible discontinuities.The solution of the displacement field is continuous through the interface between the two domains.The tangential components of the microdistortion P •e 2 = (P 12 , P 22 ) T are continuous on the interface while the normal components P • e 1 = (P 11 , P 21 ) T exhibit discontinuities.

Figure 5 :
Figure 5: The displacement and the micro-distortion fields of the first boundary problem.

Figure 6 :Figure 7 :
Figure 6: Illustration of P 21 along the inspection line y = 0.5 using the nodal elements with different mesh densities.

Figure 8 :
Figure 8: The geometry of the second boundary value problem.

Figure 9 :
Figure 9: The non-vanishing micro-distortion components of the second boundary problem using 30256 Q2NQ2 elements for L c = 5.

Figure 10 :
Figure 10: The non-vanishing components of P along the radius for the H 1 (B) × H 1 (B) elements using three mesh densities and L c = 5.

Figure 11 :
Figure 11: The non-vanishing components of P along the radius for H 1 (B) × H(curl, B) elements using three mesh densities and L c = 5.

Figure 12 :
Figure 12: Elastic energy W along the radius.

PFigure 13 :Figure 14 :
Figure 13: The non-vanishing components of P and ∇u along the radius.∇u is not influenced by the value of the characteristic length L c .

v 1 4 Figure 20 :
Figure 20: Tangential-conforming vectorial shape functions of NQ1 element.Blue circles indicate the position where the dofs are defined.

Table 1 :
Material parameters of the first boundary problem, see Figure4.

Table 2 :
Material parameters of the second boundary problem, see Figure8.