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, effective properties of a metamaterial may change significantly. To capture this phenomenon accurately and efficiently, homogenization schemes are required that reflect microstructural as well as macro-structural instabilities, large deformations, and non-local effects. 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 differential equation. Upon discretizing the resulting problem it turns out that the macroscopic stiffness matrix requires integration of macro-element basis functions as well as their derivatives, thus calling for higher-order integration rules. Because evaluation of a constitutive law in multiscale schemes involves an expensive solution of a non-linear boundary value problem, computational efficiency of the micromorphic scheme 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 expensive mass-matrix to a cheaper stiffness-matrix equivalent are outlined first. An efficient one-point integration quadrilateral element is then introduced and a proper hourglass stabilization is discussed. Performance of the resulting set of elements is finally tested on a benchmark bending example, showing that we achieve accuracy comparable to the full quadrature rules, whereas computational cost decreases proportionally to the reduction in the number of quadrature points used.

Mechanical metamaterials have recently received a great amount of attention in the engineering literature, aiming at applications ranging from acoustics to soft robotics [1,2,3,4]. 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. [5], an example of which is shown in Fig. 1a. Such patterning results in abrupt changes in the effective material and structural behaviour, from which non-local effects emerge.

a Typical patterning mode in mechanical metamaterial consisting of a square stacking of holes, cf. [5]. b Hourglassing observed in under-integrated elements in standard finite element technologies

Full size image

For engineering design and optimization of mechanical metamaterials, efficient numerical modelling tools are required. One of the available options is computational homogenization [6], which replaces the 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. [7,8,9]. 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 Finite Element (FE) analysis on the RVE level, from which average quantities such as the homogenized stress and constitutive tangent stiffness are identified. Multiple approaches can be used to reduce the computational effort, including reduced-order modelling at the level of each RVE or considering equivalent surrogate models, see e.g. [10,11,12].

In this contribution, we focus on a yet another approach to reduce the computational effort in the context of the micromorphic computational homogenization [13] by decreasing the number of macroscopic integration points. We expect that computing times will scale approximately linearly with the number of macroscopic integration points used, since most of the computing effort in a multiscale computational homogenization analysis is spent on evaluating the macro-scale constitutive law from the solutions of individual RVE boundary value problems. Depending on the bandwidth and size of the discretized macroscopic stiffness matrix, a certain overhead due to a repetitive solution of the resulting system of linear equations occurs as well. Furthermore, we systematically test and review performance of several element types. In standard FE technology, in addition to the reduction in computational costs, the 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 occurring 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 a 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. [14] and Fig. 1b, can be enhanced with standard stabilization procedures and used in micromorphic computational homogenization.

Throughout the manuscript, scalars are denoted a, vectors \(\vec {a}\), second-order tensors \({\varvec{\sigma }}\), fourth-order tensors \({\mathbb {C}}\), column matrices \({\underline{a}}\), and matrices \({\varvec{\mathsf {A}}}\). A position vector in the reference configuration (in two dimensions) is denoted \( \vec {X} = X_1\vec {e}_1 + X_2\vec {e}_2 \), a gradient operator is defined as \( \vec {\nabla } \vec {a} = \textstyle \frac{\partial a_j}{\partial X_i} \vec {e}_i \vec {e}_j \), single contraction \({\varvec{A}}\cdot {\varvec{B}} = A_{ik}B_{kj}\vec {e}_i\vec {e}_j \), double contraction \({\varvec{A}}:{\varvec{B}} = A_{ij}B_{ji}\), and right transposition of a fourth-order tensor \({\mathbb {C}} = C_{ijkl}\) as \({\mathbb {C}}^\mathsf {RT} = C_{ijlk}\). Summation over repeated indices is assumed for tensor operations, unless indicated otherwise.

Kinematic decomposition and governing equations

The micromorphic computational homogenization, developed originally in [13] and extended later on in [15], relies on a multiscale kinematic decomposition of the displacement field \(\vec {u}(\vec {X},\vec {X}_\mathrm {m})\) in the close vicinity of a macroscopic point \(\vec {X}\) in the following form

Two scales are present, a macroscopic scale, spanned by \(\vec {X} \in \Omega \), and a local microscopic scale, spanned by \(\vec {X}_\mathrm {m} \in \Omega _\mathrm {m}\), where \(\Omega \) denotes the macroscopic domain and \(\Omega _\mathrm {m}\) a microscopic RVE associated with each macroscopic point \(\vec {X}\). The unknown fields to be computed are the macroscopic mean solution \(\vec {v}_0\), micromorphic fields \(v_i\), \(i = 1,\ldots ,n\), that act as magnitudes of patterning modes \(\vec {\varphi }_i\), and a microfluctuation field \(\vec {w}\). Individual patterning modes \(\vec {\varphi }_i\), \(i = 1, \ldots , n\), are assumed to be known a priori, identified either experimentally [16] or through a buckling or Floquet–Bloch analysis [5]. An example of a patterning mode in an infinite microstructure with a square stacking of holes is shown in Fig. 1a.

Assuming a hyperelastic behaviour of the elastomeric base material, specified by an elastic energy density function W, an averaged total potential energy \({\mathcal {E}}\) can be defined as

Dropping for simplicity of the explanation any additional constraints imposed on the \(\vec {w}\) term (specified later in Eqs. (13)–(16)), the first variation of \({\mathcal {E}}\) reads

where \({\varvec{P}}({\varvec{F}}) = \frac{\partial W({\varvec{F}})}{\partial {{\varvec{F}}}^{\mathsf {T}}}\) is the microscopic first Piola–Kirchhoff stress tensor, \(\vec {\nabla }_\mathrm {m} = \frac{\partial }{\partial X_{\mathrm {m},i}}\) denotes microscopic gradient operator with respect to the reference configuration, \(\vec {u}\) is a function of all unknown field quantities (i.e. \(\vec {v}_0\), \(v_i\), and \(\vec {w}\)), and \(|\Omega _\mathrm {m}|\) stands for the RVE area. Making use of the kinematic decomposition of Eq. (1), a set of macroscopic and microscopic governing equations can be derived following the standard rules of the calculus of variations, see [13, 15] for further details.

The mean field \(\vec {v}_0\) is found to be governed by the following balance equation of linear momentum

where \(\vec {\nabla } = \frac{\partial }{\partial X_i}\) denotes the macroscopic gradient operator, \(\vec {N}\) is the outer unit normal to the boundary of the macroscopic domain \(\partial \Omega \), and \(\Gamma _\mathrm {N} \subseteq \partial \Omega \) 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. (4)–(7), the homogenized stress quantities read as

The microfluctuation fields \(\vec {w}(\vec {X},\vec {X}_\mathrm {m})\), \(\vec {X}_\mathrm {m} \in \Omega _\mathrm {m}\), associated with each macroscopic point \(\vec {X}\), satisfy the following set of equations

In the governing equation (11), \(\mu \), \(\nu _i\), and \(\vec {\eta }_i\) denote the Lagrange multipliers associated with the orthogonality constraints imposed on \(\vec {w}\), Eqs. (14)–(16), \(\vec {\lambda }\) is the Lagrange multiplier (equivalent to RVE boundary tractions) enforcing the periodicity constraint of Eq. (12), and \(\llbracket \vec {w}(\vec {X},\vec {X}_\mathrm {m}) \rrbracket = \vec {w}(\vec {X},\partial \Omega _\mathrm {m}^+) - \vec {w}(\vec {X},\partial \Omega _\mathrm {m}^-)\) denotes the jump of the field \(\vec {w}(\vec {X},\vec {X}_\mathrm {m})\) across the RVE boundary \(\partial \Omega _\mathrm {m} = \partial \Omega _\mathrm {m}^+ \cup \partial \Omega _\mathrm {m}^-\) split into two parts, image \(\Omega _\mathrm {m}^+\) and mirror \(\Omega _\mathrm {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 [13, 15], whereas comparison of the micro-morphic scheme with first- and second-order computational homogenization is provided in [17].

Solution procedure and discretization

Using the standard FE procedures as described in [18], the macroscopic governing equations can be discretized, resulting in the iterative system of macroscopic equations

where \(\Delta {\underline{v}} = [\Delta {\underline{v}}_0,\dots ,\Delta {\underline{v}}_n]^{\mathsf {T}}\) denotes the increment of macroscopic Degrees Of Freedom (DOFs), \({\underline{f}}_\mathrm {ext}\) represents an external loading applied to \({\underline{v}}_0\), is the vector of internal macroscopic forces, is the global macroscopic stiffness matrix, and  denotes the assembly operator.

Internal forces for each macroscopic element \(e = 1,\dots ,n_e\) are computed as

where \({\varvec{\mathsf {N}}}_\bullet \) are the standard matrices of the macroscopic shape interpolation functions, \({\varvec{\mathsf {B}}}_\bullet \) stand for their derivatives, and \(w_{i_\mathrm {g}}\) are integration weights with the corresponding Jacobians \(J_{i_\mathrm {g}}\). The column matrices \({\underline{\Theta }}\) and \({\underline{\Lambda }}_i\) store components of the homogenized stress quantities defined in Eqs. (8) and (10).

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 \({\mathbb {C}}({\varvec{F}}) = \frac{\partial {\varvec{P}}({\varvec{F}})}{\partial {{\varvec{F}}}^{\mathsf {T}}}\) weighted by \(\vec {\varphi }_i\), \(\vec {\nabla }_\mathrm {m}\vec {\varphi }_i\), or \(\vec {X}_\mathrm {m}\vec {\nabla }_\mathrm {m}\vec {\varphi }_i\), which can be derived from the second variation of the averaged total potential energy \({\mathcal {E}}\), i.e.

while making use of the kinematic decomposition of Eq. (1). The second term of Eq. (21) corresponds to the Schur complement of the microfluctuation field \(\vec {w}\) solved for at the microscale, and thus coupled to all macroscopic quantities. Note that \({\varvec{\mathsf {D}}}_{ww} = \frac{1}{|\Omega _\mathrm {m}|}{\varvec{\mathsf {H}}}_{ww}\), where \({\varvec{\mathsf {H}}}_{ww}\) is defined in Eq. (31), and \({\varvec{\mathsf {D}}}_{iw}\), \(i = 0, \dots , 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 \(\vec {v}_0\) has the standard form, that is

where \({\varvec{\mathsf {C}}}\) stands for the matrix representation of the constitutive tangent \({\mathbb {C}}\). The micromorphic stiffness densities have the following structure

where no summation on repeated indices i and j is implied. The averaged constitutive tangents \(\overline{{\varvec{\mathsf {E}}}}_{ij}\), \(\overline{{\underline{F}}}_{ij}\), and \({\overline{G}}_{ij}\) are obtained as more involved weighted averages of the microscopic constitutive tensor \({\mathbb {C}}\) over the RVE domain, recall Eq. (22), expressed explicitly in tensor notation as

see [18] for more details. Note that the matrices of basis interpolation functions \({\varvec{\mathsf {N}}}_i\) as well as their derivatives \({\varvec{\mathsf {B}}}_i\), \(i = 0,\dots ,n\), in Eqs. (23)–(25) depend only on the macroscopic spatial variable \(\vec {X}\), whereas all quantities related to material properties (that is \({\underline{\Theta }}\), \({\underline{\Lambda }}_{i}\), \(\Pi _i\), \(\overline{{\varvec{\mathsf {C}}}}\), \(\overline{{\varvec{\mathsf {E}}}}_{ij}\), \(\overline{{\underline{F}}}_{ij}\), and \({\overline{G}}_{ij}\)) depend in a non-linear manner on the micro-fluctuations \(\vec {w}(\vec {X},\vec {X}_\mathrm {m})\).

Finally, at each macroscopic quadrature point, \(i_\mathrm {g}\), a microscopic problem for \({\underline{w}}\) is iteratively solved from

where \({\varvec{\mathsf {K}}}_{ww}\) is the microscopic stiffness matrix, \({\underline{f}}_w\) stores microscopic internal forces, \({\underline{\lambda }}\) collects all discretized Lagrange multipliers, and \({\varvec{\mathsf {A}}}\) is a constraint matrix reflecting all equality constraints required for \(\vec {w}\), listed in Eqs. (13)–(16).

The overall solution procedure is summarized as follows. For each macroscopic Gauss integration point \(i_\mathrm {g}\), an associated non-linear microscopic RVE problem of Eq. (31) is solved for \(\vec {w}\). When a converged \(\vec {w}\) is achieved, the microscopic first Piola–Kirchhoff stress tensor \({\varvec{P}}\) and constitutive tangent \({\mathbb {C}}\) can be established in the RVE domain, and the homogenized stresses (8)–(10) along with corresponding stiffness (23)–(30) follow by their averaging. The macroscopic internal forces (18) and stiffness matrices (21) are subsequently integrated, and the resulting macroscopic system (17) assembled and solved. Since all microfluctuation fields \({\underline{w}}\) are condensed out (i.e. solved at the Gauss integration point level), all macroscopic quantities \({\underline{v}}_0\) and \({\underline{v}}_i\) directly follow from Eq. (17). Further details on discretization and numerical implementation of the micromorphic computational homogenization scheme can be found in [18, 19].

Full integration

It may appear from the structure of \({\varvec{\mathsf {D}}}_{ij}\) in Eq. (25) that for the full integration of the macroscopic element stiffness matrix associated with the micromorphic field a mass-matrix equivalent integration rule should be employed, i.e. a rule exactly integrating the \(\int _\Omega {\varvec{\mathsf {N}}}^{\mathsf {T}} {\varvec{\mathsf {N}}} \, \mathrm {d}\vec {X}\) term for an undistorted nor mapped element, recalled for convenience in Table 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. [20, 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 (31), it is desirable to reduce the number of macroscopic integration points. This may be achieved in two ways:

1. (i)

   using a lower-order interpolation scheme for the micromorphic fields \(v_i\) than for the mean field \(\vec {v}_0\) (and thus requiring a lower order of integration);

2. (ii)

   under-integrating the element stiffness expression while keeping the same interpolation scheme for \(\vec {v}_0\) and \(v_i\).

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 \(\vec {v}_0\) in the limit of infinite scale separation, and hence lower approximation space might suffice. Note that because a standard Cauchy continuum is used at the microscale, a standard FE discretization with a full integration scheme is used to discretize the microscopic balance equations (11)–(16).

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 \({\varvec{\mathsf {D}}}_{ww}\) in Eq. (21). 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 [23, Eq. (7.7.5)] that the Schur complement in Eq. (21) will not introduce additional spurious singular modes, hence it can be neglected in further considerations. From central symmetry of the employed RVE, \(\Omega _\mathrm {m}\), and the associated patterning modes, \(\vec {\varphi }_i\), it follows (or can be numerically verified) that quantities \(\overline{{\varvec{\mathsf {E}}}}_{ij}\), \(\overline{{\underline{F}}}_{ij}\), and \({\overline{G}}_{ij}\), \(i,\,j = 0,\dots ,n\), \(i\ne j\), in Eqs. (21)–(25) are zero for square as well as for hexagonal stacking of holes [15]. The matrix \({\varvec{\mathsf {H}}}\) is thus block-diagonal and its rank, and therefore also the rank of \({\varvec{\mathsf {K}}}_\mathrm {M}^e\), reads

where \(n_\mathrm {rbm}\) is the number of RBM associated with \(\vec {v}_0\), n is the number of micromorphic fields, recall Eq. (1), and \({\mathrm {order\,{(}}}{\varvec{\mathsf {A}}}) = m\) for a symmetric \(m \times m\) matrix \({\varvec{\mathsf {A}}}\). To eliminate any spurious singular modes in the system, \({\varvec{\mathsf {D}}}_{ii}\) is required to have a full rank, i.e. \({\mathrm {order\,{(}}}{\varvec{\mathsf {D}}}_{ii}) = {\mathrm {rank\,{(}}}{\varvec{\mathsf {D}}}_{ii})\). This condition will be satisfied for both integration rules listed in Table 1 (i.e. for stiffness- as well as mass-equivalent, but not for their under-integrated versions in brackets). This is because the constant vector \({\underline{1}}\) is in all cases considered an element of the kernel of \({\varvec{\mathsf {U}}}_{ii}\), whereas it is not in the kernel of \({\varvec{\mathsf {V}}}_{ii}\) (recall Eq. (25) for the definition of \({\varvec{\mathsf {U}}}_{ii}\) and \({\varvec{\mathsf {V}}}_{ii}\)). For the stiffness-equivalent integration rule we thus conclude that the matrix \({\varvec{\mathsf {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 Figs. 1a and 5b. In such a case, only one patterning mode occurs (i.e. \(n = 1\)), and the analytical expression of the mode \(\vec {\varphi }_1\) can be found e.g. in [13, Eq. (7)]. A single element in the reference configuration (\(\vec {v}_0 = \vec {0}\), \(v_1 = 0\)) and in a deformed configuration (\(\vec {\nabla }\vec {v}_0 = - 0.05\vec {e}_2\vec {e}_2\), \(v_1(\vec {X}) = (X_1+X_2+3)/2\), \(X_i \in [-0.5,0.5]\), \(i = 1,2\)) is considered. Errors in the internal forces (of Eq. (18)) and stiffness matrix (of Eq. (21)) relative to those obtained with the highest order of integration (chosen as \(n_\mathrm {g} = 13\) for triangles and \(n_\mathrm {g} = 49\) for quadrilaterals) is shown in Fig. 2 for the reference and in Fig. 3 for the deformed configuration as a function of the number of integration points. The results of Fig. 2 show that the error quickly drops to zero, proper integration rules coincide with those of Table 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 a quadrilateral) element with i nodes and j integration points. Spurious singular modes for Q4G1 and Q8G4 quadrilaterals are shown in Fig. 4. The results of Fig. 3 show a more complicated dependence, although observed errors drop quickly with increasing number of integration points for all tested element types.

The eight-node quadrilateral with four integration points, Q8G4, shows a non-communicable mode present only in the \(\vec {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, [20, Section 5.5.7.], and [22].

\(L_2\) errors of a single element stiffness matrix for triangular (a) and quadrilateral (b) elements as a function of the number of Gauss integration points. The accuracy is measured relative to the stiffness matrix obtained for the maximum number of integration points used (\(n_\mathrm {g} = 13\) for triangles and \(n_\mathrm {g} = 49\) for quadrilaterals). Solid lines correspond to \({\varvec{\mathsf {K}}}_{00}^e\) and dashed lines to \({\varvec{\mathsf {K}}}_{11}^e\) sub-matrices of \({\varvec{\mathsf {K}}}^e_\mathrm {M}\) in Eq. (21), whereas \({\varvec{\mathsf {K}}}_{01}^e = ({\varvec{\mathsf {K}}}_{10}^e)^{\mathsf {T}} = {\varvec{\mathsf {0}}}\); elements from Table 1 are used. Configuration of elements correspond to the reference undeformed state (microstructure with a square stacking of holes used, i.e., \(n=1\)). Integration rules for which spurious singular modes occur are denoted with red circles

Full size image

\(L_2\) errors of a single element internal forces (a) and stiffnesses (b) as a function of the number of Gauss integration points. The accuracy is measured relative to the force or stiffness matrix obtained for the maximum number of integration points used (\(n_\mathrm {g} = 13\) for triangles and \(n_\mathrm {g} = 49\) for quadrilaterals). Solid lines correspond to \({\underline{f}}_0^e\) or \({\varvec{\mathsf {K}}}_{00}^e\), dashed lines to \({\underline{f}}_1^e\) or \({\varvec{\mathsf {K}}}_{11}^e\), and dash-dot lines to \({\varvec{\mathsf {K}}}_{01}^e = ({\varvec{\mathsf {K}}}_{10}^e)^{\mathsf {T}}\) sub-matrices of \({\underline{f}}_\mathrm {M}^e\) or \({\varvec{\mathsf {K}}}^e_\mathrm {M}\) in Eqs. (18) and (21); elements from Table 1 are used. Configuration of the element corresponds to homogeneous overall deformation \(\vec {\nabla }\vec {v}_0 = - 0.05\vec {e}_2\vec {e}_2\) with an affine micromorphic field \(v_1(\vec {X}) = (X_1+X_2+3)/2\), \(X_i \in [-0.5,0.5]\), \(i = 1,2\) (microstructure with a square stacking of holes used, i.e., \(n=1\))

Full size image

A spurious singular mode observed for Q8G4 element in a, and three spurious singular modes for Q4G1 element in b–d (typical hourglassing modes are shown also in Fig. 1b). a–c correspond to deformed configurations for which the micromorphic field vanishes, i.e. \(v_1(\vec {X}) = 0\), whereas d shows a surface plot of \(v_1(\vec {X})\) for which the mean field vanishes, i.e. \(\vec {v}_0(\vec {X}) = \vec {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\)

Full size image

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 Fig. 4b–d, all observed spurious singular modes correspond to the standard hourglassing mode, see Fig. 1b and e.g. [24]. In the context of the micromorphic homogenization, the hourglassing mode occurs twice in the mean field \(\vec {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. [14, 24,25,26,27,28], stabilized micromorphic Q4G1 element stiffness matrix can be written as

where \({\varvec{\mathsf {K}}}^e_\mathrm {M}\) is obtained from Eq. (21) with \(n_\mathrm {g} = 1\). The stabilization matrix \({\varvec{\mathsf {K}}}^\mathrm {stab}\) has the following form

where

\({\underline{x}}\), \({\underline{y}}\) store nodal coordinates of a Q4 element in the reference configuration, \({\underline{b}}_x\), \({\underline{b}}_y\) are discrete versions of the gradient operator, and \(\alpha _1\), \(\alpha _2\), and \(\beta _i\), \(i = 1,\dots ,n\), are stabilization constants. An example of explicit stabilization values for standard continuum can be found e.g. in [27, Table 1]. For simplicity, all mutual coupling terms between \(\vec {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 [29] is adopted with small amendments. This formulation was originally developed for the stabilization of quadrilateral elements in large-deformation thermo-mechanically coupled problems, in which a scalar temperature field is considered in addition to a vector displacement field. A direct link between temperature and the micromorphic fields \(v_i\) can be thus established. The mechanical part of the stabilization matrix from (34) reads

where the first term on the right-hand side \({\varvec{\mathsf {K}}}_{uu,\mathrm {hg}}\) corresponds to hourglass stabilization, while the Schur complement accounts for the contribution of the incompatible (enhanced) assumed strain, not elaborated here in more detail. The constitutive stabilization tangent (denoted as \({\varvec{A}}_0^\star \) in [29]) is taken as \(\overline{{\varvec{\mathsf {C}}}}\) (defined in Eq. (23)), and updated only once at the beginning of the simulation. Although the stabilization matrix \({\varvec{\mathsf {K}}}^\mathrm {stab}_{00}\) can be expressed by a pair of stabilization constants \(\alpha _1\) and \(\alpha _2\) [27], we do not provide here their explicit forms, and rather define them implicitly through the second expression on the right hand side of Eq. (37) due to the presence of the Schur complement term \({\varvec{\mathsf {K}}}_{uv}{\varvec{\mathsf {K}}}_{vv}^{-1}{\varvec{\mathsf {K}}}_{vu}\). The structure of the micromorphic stabilization matrices reads

where the stabilization coefficient \(\beta _i\) can be expressed in an explicit form as

and where \({\underline{L}}_\mathrm {hg} = [\eta ,\xi ,\eta ,\xi ]^{\mathsf {T}}\), with \(\xi ,\,\eta \in [-1,1]\) being local coordinates on the parent element, \({\varvec{\mathsf {j}}}_0\) is a matrix storing entries of the Jacobian of the isoparametric mapping evaluated at the centre of the element, \({\varvec{\mathsf {E}}}_{ii}^\star \) is taken as \(\overline{{\varvec{\mathsf {E}}}}_{ii}\) (introduced in Eq. (25), see also Eq. (28)) and again updated only once at the beginning of the simulation, whereas the integral is taken over an element’s volume \(\Omega _e\) using the full four-point Gauss integration rule. For a thorough definition of all quantities involved in Eqs. (37)–(39) and for further details we refer to [29].

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 \({\varvec{F}} = {\varvec{I}} + (\vec {\nabla }\vec {u})^{\mathsf {T}}\) is the deformation gradient, \(J = \det {{\varvec{F}}}\), \({\varvec{I}}\) is the second-order identity tensor, and \(I_1 = {\mathrm {tr\,{{\varvec{C}}}}}\) is the first invariant of the right Cauchy–Green deformation tensor \({\varvec{C}} = {\varvec{F}}^{\mathsf {T}} \cdot {\varvec{F}}\). Constitutive parameters are chosen according to [5] as \( a_1 = 0.55~\) MPa, \( a_2 = 0.3\) MPa, and \(\kappa = 55\) MPa. The adopted RVE is shown in Fig. 5b, having cell size \(\ell = 1\) and diameter \(d = 0.85\,\ell \), whereas the only patterning mode \(\vec {\varphi }_1\) (i.e. \(n = 1\)) is shown in Fig. 1a (for analytical expression we again refer to [13, 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 [30, Section 8.8.2].

Four element types with various integration rules are used for the discretization of the macroscopic domain \(\Omega \), 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 \in \{1, 2, 4, 8, 16\} \, \ell \) 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 \(\ell /5\) with three Gauss integration points are used.

a Considered benchmark example of a cantilever beam subjected to bending. b Representative volume element associated with each macroscopic Gauss integration point (the corresponding patterning mode is shown in Fig. 1a)

Full size image

a Deformed mesh of \(\vec {v}_0(\vec {X})\); colour plot shows \(\vec {v}_0\cdot \vec {e}_2\). b Corresponding field plot of \(v_1(\vec {X})\), close-up on the left quarter of the specimen. In both cases, \(t = 1\), and the third finest discretization (i.e. \(h = 4\ell \)) using Q8G9 elements for the bending example shown in Fig. 5a is used

Full size image

The numerical solution of the macroscopic problem obtained for Q8G9 elements with \(h = 4\ell \) 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\ell \) is considered as the reference solution against which the accuracy of other discretizations is measured.

Figure 7 shows convergence plots of the relative error for individual element types, \(\epsilon \), which is defined as

where g denotes a quantity of interest (\(\vec {v}_0\), \(v_1\), or \({\mathcal {E}}\)), and \(g_\mathrm {ref}\) is its corresponding reference value represented by the T6G6, \(h = 1\ell \), solution. In particular, Fig. 7a shows error in the mean field \(\vec {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.

Convergence of a\(L_2\) displacement error in \(\vec {v}_0(\vec {X})\), \(t = 1\), b\(L_2\) micromorphic error in \(v_1(\vec {X})\), \(t = 1\), and c error in energy evolutions \({\mathcal {E}}(t)\), \(t \in [0,1]\). d Reaction bending moment M(t) for Q4G1 and various element sizes h in comparison with the reference solution. Data points for T6G6, \(h=1\ell \), are missing, because this solution is considered as the reference solution and the corresponding error thus vanishes

Full size image

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\ell \) and \(h = 2\ell \) is observed. The hourglass energy required for stabilization of all elements, \({\mathcal {H}}\), relative to the elastic energy of the entire system, \({\mathcal {E}}\), is summarized in Table 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\ell \).

Computing times T relative to the fastest discretization (Q4G1) as well as the number of macroscopic elements \(n_e\), Gauss integration points \(n_\mathrm {g}\), and nodes \(n_\mathrm {node}\), are summarized for \(h = 2\ell \) in Table 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 note finally that unlike the definitions of Eqs. (37)–(39) which perform optimally with respect to incompatible enhanced strains, the stabilization factors \(\alpha _1\), \(\alpha _2\), and \(\beta _i\), \(i = 1, \dots , n\), can also be chosen manually. When their values are chosen too small or when they vanish altogether for \(\vec {v}_0\) (i.e. \(\alpha _1 = \alpha _2 = 0\)), spurious instabilities pollute the mean solution. For comparison, such a situation is depicted in Fig. 8a, where a deformed configuration is shown along with the horizontal component of the mean field \(\vec {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 \(\vec {v}_0\) (i.e. \(\beta _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 \(\alpha _1\), \(\alpha _2\), and \(\beta _i\); it is then recommended to switch (locally or globally) to the full four-point integration scheme.

Example of wrong choices of stabilization parameters for Q4G1, \(h = 4\ell \). a Artificial instability polluting the mean field \(\vec {v}_0\); deformed configuration at \(t = 0.5\) obtained for \(\alpha _1 = \alpha _2 = 0\), where colour plot shows \(\vec {v}_0\cdot \vec {e}_1\). b Artificial instability in the micromorphic field polluting a buckling mode obtained for \(\beta _1 = 0\). Recall the structure of \({\varvec{\mathsf {K}}}^\mathrm {stab}\) in Eq. (34)

Full size image

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 eight-noded quadrilateral with four integration points, which achieves good performance and high accuracy.

Not applicable.

• 29 May 2020

  Following publication of the original article [1], the authors reported the errors in the equations.

RVE:

Representative volume element

FE:

Finite element

FEM:

Finite element method

RBM:

Rigid body modes

TiGj (or QiGj):

A triangular (or a quadrilateral) element with i nodes and j integration points

DOF:

Degree of freedom

The research leading to these results has received funding from the Czech Science Foundation (GAČR) under Grant Agreement No. (19-26143X). PK would like to thank the Ministry of Education, Youth and Sports of the Czech Republic for financial support of his stay at the Czech Technical University in Prague in the framework of International Mobility of Researchers, project No. CZ.02.2.69/0.0/0.0/16_027/0008465. Petr Krysl—Currently Visiting Professor at Czech Technical University in Prague.

Authors and Affiliations

1. Department of Mechanics, Faculty of Civil Engineering, Czech Technical University in Prague, Thákurova 7, 166 29, Prague 6, Czech Republic

   Ondřej Rokoš, Jan Zeman & Martin Doškář

2. Structural Engineering Department, University of California, La Jolla, San Diego, USA

   Petr Krysl

Contributions

OR carried out most of the study, performed numerical simulations, and drafted the manuscript. MD helped with implementation and numerical issues. All authors developed the methodology, conceived of the study, and participated in its design, coordination, and critical review of the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Ondřej Rokoš.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

The original version of this article was revised: The errors in the equations have been corrected.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions