Reduced integration schemes in micromorphic computational homogenization of elastomeric mechanical metamaterials

Exotic behaviour of mechanical metamaterials often relies on an internal transformation of the underlying microstructure triggered by its local instabilities, rearrangements, and rotations. Depending on the presence and magnitude of such a transformation, eﬀective properties of a metamaterial may change signiﬁcantly. To capture this phenomenon accu-rately and eﬃciently, homogenization schemes are required that reﬂect microstructural as well as macro-structural instabilities, large deformations, and non-local eﬀects. To this end, a micromorphic computational homogenization scheme has recently been developed, which employs the particular microstructural transformation as a non-local mechanism, magnitude of which is governed by an additional coupled partial diﬀerential equation. Upon discretizing the resulting problem it turns out that the macroscopic stiﬀness matrix requires integration of macro-element basis functions as well as their derivatives, thus calling for a higher-order integration rules. Because evaluation of constitutive law in multiscale schemes involves an expensive solution of a non-linear boundary value problem, computational eﬃciency can be improved by reducing the number of integration points. Therefore, the goal of this paper is to investigate reduced-order schemes in computational homogenization, with emphasis on the stability of the resulting elements. In particular, arguments for lowering the order of integration from the expensive mass-matrix to a cheaper stiﬀness-matrix equivalent are ﬁrst outlined. An eﬃcient one-point integration quadrilateral element is then introduced and proper hourglass stabilization discussed. Performance of the resulting set of elements is ﬁnally tested on a benchmark bending example, showing that we achieve accuracy comparable to the full quadrature rules, whereas computational costs decrease proportionally to the reduction in the number of quadrature points used.


Introduction
Mechanical metamaterials have recently received a great amount of attention in the engineering literature, aiming at applications ranging from acoustics to soft robotics [27,18,19,26]. In the particular case of elastomeric mechanical metamaterials, effective properties often depend on the internal state of the microstructure, dictated by the so-called pattern transformation occurring upon microstructural buckling, see e.g. [6], example of which is shown in Fig. 1a. Such patterning results in abrupt changes in effective material and structural behaviour, from which non-local effects emerge.
For engineering design and optimization of mechanical metamaterials, efficient numerical modelling tools are required. One of the options available is computational homogenization [14], which replaces effective constitutive behaviour of the macroscopic material with a mechanical system, specified by a Representative Volume Element (RVE). Advantages of such an approach are that all microstructural features are taken into account, including microstructural morphology, non-linear material behaviour, or local microstructural buckling. However, the macroscopic continuum is treated as local, i.e. it ignores any communication among neighbouring points (a consequence of the assumption of the formulation on separation of scales). To alleviate this limitation, various enriched theories have been proposed in the literature, including higher-order and micromorphic computational homogenization schemes, see e.g. [15,16,7]. In both first-and higher-order homogenization schemes, the evaluation of the effective constitutive law represents the most computationally intensive operation: it involves a separate non-linear FE analysis on the RVE level, from which average quantities such as homogenized stress and constitutive tangent stiffness are identified. Multiple approaches can be used to reduce the computational effort, including reducedorder modelling at the level of each RVE or considering equivalent surrogate models, see e.g. [28,20,25]. In this contribution, we focus on a yet another approach to reduce the computational effort in the context of the micromorphic computational homogenization [22] by decreasing the number of macroscopic integration points. We expect that computing times will scale linearly with the number of macroscopic integration points used. Furthermore, we systematically test and review performance of several element types. In standard FE technology, in addition to the reduction in computational costs, under-integration is frequently used to remove artificial stiffening that can be especially severe in the limit of incompressibility and bending-dominated simulations (i.e. volumetric and shear locking). Within the context of the micromorphic computational homogenization, however, the efficiency is the primal motivation, as employed microstructures are typically heterogeneous and porous, and hence no locking occurs.
Because the stiffness matrix featuring in the discretized version of the micromorphic formulation contains integrals of basis shape functions as well as their derivatives, full integration entails expensive integration rules corresponding to mass matrix equivalents. On the basis of linear analysis it can be conjectured that integration rules accurate for basis functions' derivatives are sufficient to maintain the proper rank of element stiffness matrices. Furthermore, we demonstrate that a one-integration point quadrilateral, which suffers from the well-known hourglassing, see e.g. [5] and Fig. 1b, can be enhanced with standard stabilization procedures and used in micromorphic computational homogenization.
Throughout the manuscript, scalars are denoted a, vectors a, second-order tensors σ, fourth-order tensors C, column matrices a, and matrices A. Position vector (in two dimensions) is denoted X = X 1 e 1 + X 2 e 2 , gradient operator ∇ a = ∂a j ∂X i e i e j , single contraction A · B = A ik B kj e i e j , and double contraction A : B = A ij B ji . Summation over repeated indices is assumed for tensor operations, unless indicated otherwise.

Kinematic decomposition and governing equations
The micromorphic computational homogenization relies on a kinematic decomposition of the displacement field u( X, X m ) in the close vicinity of a point X in the following form Two scales are present, a macroscopic scale, spanned by X ∈ Ω, and a local microscopic scale, spanned by X m ∈ Ω m , where Ω denotes the macroscopic domain and Ω m a microscopic RVE associated with each macroscopic point X. The unknown fields to be computed are the macroscopic mean solution v 0 , micromorphic fields v i , i = 1, . . . , n, that act as magnitudes of patterning modes ϕ i , and a microfluctuation field w. Individual patterning modes ϕ i , i = 1, . . . , n, are assumed to be known a priori, identified either experimentally or through a buckling or Floquet-Bloch analysis. An example of a patterning mode in an infinite microstructure with a square stacking of holes is shown in Fig. 1a.
The mean field v 0 obeys the following balance equation of linear momentum where ∇ = ∂ ∂X i denotes the macroscopic gradient operator, N is the outer unit normal to the boundary of the macroscopic domain ∂Ω, and Γ N ⊆ ∂Ω is the part of the macroscopic domain where zero tractions are prescribed. The magnitude v i of each patterning field is governed by the following scalar partial differential equation In Eqs. (2)-(5), the homogenized stress quantities read as where P (F ) = ∂W (F ) ∂F T is the microscopic first Piola-Kirchhoff stress tensor with W being an elastic energy density function, ∇ m = ∂ ∂X m,i denotes microscopic gradient operator, and |Ω m | stands for the RVE area.
The microfluctuation fields w( X, X m ), X m ∈ Ω m , associated with each macroscopic point X, satisfy the following set of equations in Ω m , In the governing equation (9), µ, ν i , and η i are Lagrange multipliers associated with orthogonality constraints of Eqs. (12)- (14), λ is the Lagrange multiplier (equivalent to RVE boundary tractions) enforcing the periodicity constraint of Eq. (10), and w( X, X m ) = w( X, ∂Ω + m ) − w( X, ∂Ω − m ) denotes the jump of the field w( X, X m ) across the RVE boundary ∂Ω m = ∂Ω + m ∪ ∂Ω − m split into two parts, image Ω + m and mirror Ω − m (distinguished in Fig. 5b in green and blue colour).
Full details on the derivation of macro-and micro-scopic governing equations and energy considerations are available in [22,23], whereas comparison of the micro-morphic scheme with first-and second-order computational homogenization appears in [24].

Solution procedure and discretization
The macroscopic governing equations are discretized using the standard FE procedures, resulting in the following iterative system of macroscopic equations where N • are the standard matrices of the macroscopic shape interpolation functions, B • stand for their derivatives, and w ig are integration weights with the corresponding Jacobians J ig . The column matrices Θ and Λ i store components of the homogenized stress quantities defined in Eqs. (6) and (8).
The macroscopic single-element stiffness matrix can be expressed as where the first term on the right-hand side reflects stiffness quantities obtained by volume averaging of microscopic constitutive tangent whereas the second term corresponds to the Schur complement of the microfluctuation field w solved for at the microscale and thus coupled to all macroscopic quantities. Note that D ww = 1 |Ωm| H ww , where H ww is defined in Eq. (21) below, and D iw , i = 0, . . . , n, are coupling terms between the macroscopic and microfluctuation fields not discussed here in detail for brevity. The stiffness density term associated with the mean field v 0 has the standard form, that is where C stands for the matrix representation of the constitutive tangent C. The micromorphic stiffness densities have the following structure Averaged constitutive tangents, E ij , F ij , and G ij , are more involved weighted averages of the microscopic constitutive tensor C over the RVE domain, again skipped for brevity. Note that the matrices of basis interpolation functions N i as well as their derivatives B i , i = 0, . . . , n, depend only on the macroscopic spatial variable X, whereas all quantities related to material properties (that is Θ, Λ i , Π i , C, E ij , F ij , and G ij ) depend in a non-linear manner on the micro-fluctuations w( X, X m ). Finally, at each macroscopic quadrature point, i g , a microscopic problem for w is itera- where K ww is the microscopic stiffness matrix, f w stores microscopic internal forces, λ collects all discretized Lagrange multipliers, and A is a constraint matrix reflecting all equality constraints required for w, listed in Eqs. (11)- (14). Further details on discretization and numerical implementation of the micromorphic computational homogenization scheme can be found in [9,8].

Full integration
It may appear from the structure of D ij in Eq. (20) that for the full integration of the micromorphic element stiffness matrix the mass-matrix equivalent integration rule should be employed, recalled for convenience in Tab. 1 for a few typical element types. Note that the listed quadrature rules are accurate for polynomial functions. The integrated quantities in FEM are, however, rational functions in general. Such integration rules thus might represent a sub-optimal choice, yet they are often employed in numerical simulations for non-linear problems and isoparametric elements, see e.g. [1,Section 5.5.5].
Because evaluation of the macroscopic constitutive law is a rather expensive operation involving a separate solution of a non-linear FEM system (21), it is desirable to reduce the number of integration points. This may be achieved in two ways: (i) using a lower-order interpolation scheme for the micromorphic fields v i than for the mean field v 0 (and thus requiring a lower order of integration); (ii) under-integrating the element stiffness expression while keeping the same interpolation scheme for v 0 and v i .
a In standard two-dimensional FEM, two spurious singular modes requiring hourglass control are observed. b Induces one spurious singular mode in standard FEM; the mode is non-communicable, i.e. it typically does not occur in an element patch and hence no stabilization is usually required, cf. To achieve high approximation quality of all involved fields, especially in the regime of small separation of scales, we opt for the latter alternative, although from the standard asymptotic homogenization theory one might expect that the magnitude of the fluctuation v i would scale with the gradient of v 0 in the limit of infinite scale separation, and hence lower approximation space might suffice.

Reduced integration and spurious singular modes
To conjecture that the order of integration can be lowered from the mass-matrix type of quadrature to the stiffness matrix equivalent without introducing any additional spurious singular modes, let us start with a set of supporting arguments carried out for a single element stiffness matrix considered in the reference configuration (i.e. in the limit of small deformations).
For the sake of argument we assume the penalty method being used instead of the equivalent Lagrange multiplier formulation, merely to avoid indefiniteness of D ww in Eq. (17). Assuming further that the employed mechanical system is stable in the reference configuration, i.e. the corresponding stiffness matrix is positive definite (upon suppressing proper physical singular modes, i.e. Rigid Body Modes, abbreviated as RBM in what follows), one can argue by means of [12,Eq. (7.7.5)], that the Schur complement in Eq. (17) will not introduce additional spurious singular modes, hence it can be neglected in further considerations. From central symmetry of the employed RVE, Ω m , and the associated patterning modes, ϕ i , it follows (or can be numerically verified) that quantities E ij , F ij , and G ij , i, j = 0, . . . , n, i = j, in Eqs. (17)- (20) are zero for square as well as for hexagonal stacking of holes. The matrix H is thus block-diagonal and its rank, and therefore also the rank of K e M , reads where n rbm is the number of RBM associated with v 0 , n is the number of micromorphic fields, recall Eq. (1), and order (A) = m for a symmetric m×m matrix A. To eliminate any spurious singular modes in the system, D ii is required to have a full rank, i.e. order (D ii ) = rank (D ii ).
This condition will be satisfied for both integration rules listed in Tab. 1 (i.e. for stiffnessas well as mass-equivalent, but not for their under-integrated versions in brackets). This is because the constant vector 1 is in all cases considered an element of the kernel of U ii , whereas it is not in the kernel of V ii (recall Eq. (20) for the definition of U ii and V ii ). For the stiffness-eqivalent integration rule we thus conclude that the matrix D ii maintains its full rank.
To further support the above considerations and to extend their validity to the non-linear regime, we perform a numerical test on a metamaterial with a square stacking of holes, shown in Fig. 1a and 5b. In such a case, only one patterning mode occurs (i.e. n = 1), and the analytical expression of the mode ϕ 1 can be found e.g. in [22, Eq. (7)]. A single element in the reference configuration ( v 0 = 0, v 1 = 0) and in a deformed configuration ( ∇ v 0 = −0.05 e 2 e 2 , v 1 ( X) = (X 1 + X 2 + 3)/2, X i ∈ [−0.5, 0.5], i = 1, 2) is considered. Errors in the internal forces (of Eq. (16)) and stiffness matrix (of Eq. (17)) relative to those obtained with the highest order of integration (chosen as n g = 13 for triangles and n g = 49 for quadrilaterals) is shown in Fig. 2 as a function of the number of integration points. These results show that the error quickly drops to zero, proper integration rules coincide with those of Tab. 1, and that the stiffness-equivalent integration rule provides error as low as 0.5 % for triangles and 5 % for quadrilaterals. Integration rules for which spurious singular modes occur are emphasized by red circles, and correspond to T6G1, and Q4G1, Q8G1, Q8G4, where the notation TiGj (or QiGj) stands for a triangular (or quadrilateral) element with i nodes and j integration points. Spurious singular modes for Q4G1 and Q8G4 quadrilaterals are shown in Fig. 4.
The eight-node quadrilateral with four integration points, Q8G4, shows a non-communicable mode present only in the v 0 field. This mode is observed also in the standard FE formulations, but is usually not of any concern as it typically appears for a single element only and not in an element patch, cf. Fig. 4a, [1, Section 5.5.7.], and [10].

Q4G1 element
A question now arises whether the Q4G1 element could be stabilized to obtain an efficient element with a single Gauss integration point only for the micromorphic computational homogenization. As it may be seen from Figs. 4b-4d, all observed spurious singular modes correspond to the standard hourglassing mode, see Fig. 1b and e.g. [11]. In the context of the micromorphic homogenization, the hourglassing mode occurs twice in the mean field v 0 , and once for each micromorphic field v i , thus inducing 2 + n rank deficiency of the resulting stiffness matrix in addition to proper deficiency of the three RBM.
Following standard stabilization procedures, see e.g. [11,3,17,5,2,13], stabilized micromorphic Q4G1 element stiffness matrix can be written as where K e M is obtained from Eq. (17) with n g = 1. The stabilization matrix K stab has the  (17); elements from Tab. 1 are used. Configuration of the element corresponds to homogeneous overall deformation ∇ v 0 = −0.05 e 2 e 2 with an affine micromorphic field v 1 ( X) = (X 1 + X 2 + 3)/2, X i ∈ [−0.5, 0.5], i = 1, 2 (microstructure with a square stacking of holes used, i.e., n = 1). Figure 4: Spurious singular mode observed for Q8G4 element in (a), and three spurious singular modes for Q4G1 element in (b)-(d) (typical hourglassing modes emerging also in Fig. 1b). Figures (a)-(c) correspond to deformed configurations for which the micromorphic field vanishes, i.e. v 1 ( X) = 0, whereas (d) shows a surface plot of v 1 ( X) for which the mean field vanishes, i.e. v 0 ( X) = 0. Notation QiGj designates a quadrilateral element with i nodes and j integration points. In all cases, bottom left node is fixed in both directions, whereas the right bottom node is fixed in the vertical direction to eliminate rigid body modes. A microstructure with a square stacking of holes used, i.e., n = 1.
following form where x, y store nodal coordinates of a Q4 element in the reference configuration, b x , b y are discrete versions of the gradient operator, and α 1 , α 2 , and β i , i = 1, . . . , n, are stabilization constants for which the explicit values can be found e.g. in [2,Tab. 1]. For simplicity, all mutual coupling terms between v 0 and v i fields have been neglected, although a more in-depth derivation would yield off-diagonal stabilization terms as well. Such considerations are, however, outside the scope of this study and are left for future considerations.
In what follows, the formulation by Reese [21] is adopted with small amendments. This formulation was originally developed for the stabilization of quadrilateral elements in largedeformation thermo-mechanically coupled problems. The mechanical part of the stabilization matrix reads where the first term on the right-hand side corresponds to hourglass stabilization, while the Schur complement accounts for the contribution of the incompatible (enhanced) assumed strain. The constitutive stabilization tangent (denoted as A 0 in [21]) is taken as C (in Eq. (18)), and updated only once at the beginning of the simulation. The structure of the micromorphic stabilization matrices reads where L hg = [η, ξ, η, ξ] T with ξ, η ∈ [−1, 1] local coordinates on the parent element, j 0 is a matrix storing entries of the Jacobian of the isoparametric mapping evaluated at the centre of the element, E ii is taken as E ii and again updated only once at the beginning of the simulation, whereas the integral is taken over an element's volume Ω e . For a detailed definition of all quantities involved in Eqs. (26)- (27) and for further details we refer to [21].

Performance evaluation
In this section, performance of individual elements is tested in terms of convergence properties and computational efficiency. A simple plane strain bending example shown in Fig. 5 is considered, having a microstructure consisting of a square stacking of circular holes in an otherwise homogeneous hyperelastic matrix. Constitutive behaviour of the base material is specified by the following elastic energy density function where F = I + ( ∇ u) T is the deformation gradient, J = det F , I is the second-order identity tensor, and I 1 = tr C and I 2 = 1 2 [(tr C) 2 − tr (C 2 )] are the first two invariants of the right Cauchy-Green deformation tensor C = F T · F . Constitutive parameters are chosen according to [6] as a 1 = 0.55 MPa, a 2 = 0.3 MPa, and K = 55 MPa. The adopted RVE is shown in Fig. 5b, having cell size = 1 and diameter d = 0.85 , whereas the only patterning mode ϕ 1 (i.e. n = 1) is shown in Fig. 1a (for analytical expression we again refer to [22, Eq. (7)]). Distributed external loading, applied at the right vertical edge of the cantilever specimen, was prescribed in the form where t is parametrization time. Note that a similar problem was also considered in [4,Section 8.8.2]. Four element types with various integration rules are used for the discretization of the macroscopic domain Ω, namely linear triangles T3G1 and T3G3, quadratic triangles T6G3 and T6G6, bi-linear quadrilaterals Q4G4 and the stabilized quadrilateral Q4G1, and serendipity quadratic quadrilaterals Q8G9 and Q8G4. For all element types, five element sizes h ∈ {1, 2, 4, 8, 16} are used, discretized with right-angled triangles and rectangular quadrilaterals of good aspect ratios. For the discretization of the RVE problem, quadratic triangular elements of a typical size /5 with three Gauss integration points are used.
The numerical solution of the macroscopic problem obtained for Q8G9 elements with h = 4 at t = 1 appears in Fig. 6. From the micromorphic field v 1 shown in Fig. 6b we can clearly   see that the patterning mode is triggered only in the compressive region at the bottom-left part of the specimen domain, and decays rapidly in the close vicinity of prescribed displacements (the left vertical edge). Although the obtained results correspond to expectations, no comparison with a direct numerical simulation fully accounting for the underlying microstructure is available. Instead, the solution corresponding to the T6G6 element with h = 1 is considered as the reference solution against which the accuracy of other discretizations is measured. Fig. 7 shows convergence plots of the relative error for individual element types, , which is defined as where g denotes a quantity of interest ( v 0 , v 1 , or E), and g ref is its corresponding reference value represented by the T6G6, h = 1 , solution. In particular, Fig. 7a shows error in the mean field v 0 as a function of the element size h evaluated at the end of the loading history, i.e. at t = 1. Similarly, Fig. 7b shows the error evaluated for the micromorphic field, whereas Fig. 7c captures the error in energy evolutions evaluated as a function of time.
Observed trends for all three figures are very similar: the linear and bi-linear elements (T3, Q4) converge slower than the quadratic elements (T6, Q8), ordered according to their level of interpolation. Note that decreased accuracy of Q8G9 element for h = 1 relative to h = 2 might be explained by the fact that exact solution is not known, and is only approximated. The under-integrated element Q4G1 reaches a similar rate of convergence as T3G1, T3G3, and Q4G4, although it is slightly less accurate, especially in the energy norm. Its corresponding convergence in the reaction moment is shown in Fig. 7d, where a good agreement for the two finest approximations h = 1 and h = 2 is observed. The hourglass energy required for stabilization of all elements, H, relative to the elastic energy of the entire system, E, is summarized in Tab. 2. The entries presented suggest that due to the bending character of the deformation and loading, the stabilization energy reaches considerable values for coarse discretizations, being acceptable only for as fine discretizations as h = 2 .
Computing times T relative to the fastest discretization (Q4G1) as well as number of macroscopic elements n e , Gauss integration points n g , and nodes n node , are summarized for h = 2 in Tab. 3. Here we conclude that although Q4G1 achieves comparable accuracy as T3G1 and Q4G4, it is significantly faster. Another interesting combination is the element Q8G4, which is roughly four times slower in comparison with Q4G1, but achieves a significantly higher accuracy (recall Fig. 7).
Let us finally note that when too small stabilization factors are chosen or when they vanish altogether for v 0 (i.e. α 1 = α 2 = 0), spurious instabilities pollute the mean solution. Such a situation is depicted in Fig. 8a, where a deformed configuration is shown along with the horizontal component of the mean field v 0 (in colour), featuring clear hourglass oscillations. When such oscillations become too severe, the Newton solver even fails to converge. Upon choosing too weak or vanishing stabilization in the micromorphic field v 1 while stabilizing the mean solution v 0 (i.e. β 1 = 0), a buckling mode, needed to initialize the Newton solver towards a proper equilibrium path upon physical instability, is also distorted by the spurious hourglassing mode clearly seen in Fig. 8b. Spurious oscillations may occur even despite the proper stabilization, i.e. for sufficiently large α 1 , α 2 , and β i ; it is then recommended to switch (locally or globally) to the full four-point integration scheme.

Summary and conclusions
In this paper, proper choices of elements and integration rules for the micromorphic computational homogenization scheme have been discussed. It has been shown that although the micromorphic part of the stiffness matrix necessitates a mass-matrix equivalent integration rule, this requirement can be relaxed with only a negligible drop in accuracy and no artificial spurious modes. An efficient one-integration point quadrilateral element was further discussed, which offers the advantage of being computationally efficient, while maintaining a reasonable accuracy and convergence. Because standard hourglassing modes are observed for this type of element based on the eigenvalue analysis of the stiffness matrix, a suitable stabilization technique was introduced and discussed, based on methods available in the literature. A dedicated stabilization technique might be developed to further improve the performance of the under-integrated element, which lies, however, outside of the scope of the current study. Using a benchmark numerical example it was further shown that the achieved performance of the one-integration point element is satisfactory in comparison to the full integration rule, requiring, nevertheless, approximately only one fourth of the computational effort in comparison to the fully integrated four-node quadrilateral. Another interesting option for the micromorphic computational homogenization scheme is the eightnoded quadrilateral with four integration points, which achieves good performance and high accuracy.