# Enriched Computational Homogenization Schemes Applied to Pattern-
Transforming Elastomeric Mechanical Metamaterials

S.O. Sperling<EMAIL_ADDRESS>T. Guo<EMAIL_ADDRESS>R.H.J. Peerlings
<EMAIL_ADDRESS>V.G. Kouznetsova<EMAIL_ADDRESS>M.G.D. Geers
<EMAIL_ADDRESS>O. Rokoš<EMAIL_ADDRESS>Mechanics of Materials, Department
of Mechanical Engineering, Eindhoven University of Technology, P.O. Box 513,
5600 MB Eindhoven, The Netherlands  Institute for Complex Molecular Systems,
Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, The
Netherlands  Centre for Analysis, Scientific Computing and Applications,
Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, The
Netherlands

###### Abstract

Elastomeric mechanical metamaterials exhibit unconventional mechanical
behaviour owing to their complex microstructures. A clear transition in the
effective properties emerges under compressive loading, which is triggered by
local instabilities and pattern transformations of the underlying cellular
microstructure. Such transformations trigger a non-local mechanical response
resulting in strong size effects. For predictive modelling of engineering
applications, the effective homogenized material properties are generally of
interest. For mechanical metamaterials, these can be obtained in an expensive
manner by ensemble averaging of the direct numerical simulations for a series
of translated microstructures, applicable especially in the regime of small
separation of scales. To circumvent this expensive step, computational
homogenization methods are of benefit, employing volume averaging instead.
Classical first-order computational homogenization, which relies on the
standard separation of scales principle, is unable to capture any size and
boundary effects. Second-order computational homogenization has the ability to
capture strain gradient effects at the macro-scale, thus accounting for the
presence of non-localities. Another alternative is micromorphic computational
homogenization scheme, which is tailored to pattern-transforming metamaterials
by incorporating prior kinematic knowledge. In this contribution, a systematic
study is performed, assessing the predictive ability of computational
homogenization schemes in the realm of elastomeric metamaterials. Three
representative examples with distinct mechanical loading are employed for this
purpose: uniform compression and bending of an infinite specimen, and
compression of a finite specimen. Qualitative and quantitative analyses are
performed for each of the load cases where the ensemble average solution is
set as a reference. It is shown that although the second-order computational
homogenization provides good qualitative predictions, it fails in accurately
capturing bifurcation strains at which pattern-transformations occur, and
slightly over-predicts the homogenized response. Overall, the micromorphic
method provides the most accurate prediction of the effective material
behaviour, although the presence of soft boundary layers induce large errors
for small scale ratios.

###### keywords:

Mechanical metamaterials , cellular solids , computational homogenization ,
strain gradient continuum , micromorphic continuum

††journal: Springer

##  1 Introduction

Mechanical metamaterials are characterized by their exotic material properties
owing to their engineered structure. Supported by extensive research in recent
years, various types of metamaterials hold potential for dedicated
applications in distinct fields of engineering, including biomedical,
automotive, and aerospace industries, cf. [Prawoto, 2012, Krushynska et al.,
2023]. An interesting subgroup of such materials, applicable especially in
soft robotics and receiving significant scientific attention, is called
elastomeric metamaterials. Their extraordinary behaviour relies on the
transformation of recurrent patterns within the microstructure during
compression. The pattern reconfiguration yields a significant change in
effective material properties such as Young’s modulus or Poisson’s ratio, see,
e.g., Hu et al. [2013], Florijn et al. [2014], Wu et al. [2015], or Coulais
[2016] and Dykstra et al. [2023]. One of the earlier works by Mullin et al.
[2007] already showed that the pattern transformation and corresponding change
in effective properties is initiated by buckling of the ligaments which
separate individual cells.

Pattern transformation of hexagonally shaped cellular structures has been
investigated by Gao et al. [2018], where two distinct re-configurations were
identified. It was shown that the wall thickness to cell size ratio dictates
the post-buckling pattern and associated effective characteristics. Wu et al.
[2015] evaluated the effect of microstructural buckling on the effective
stress–strain path for a triangular structure built from the kagome lattices.
Here, the ligaments constructed from mirrored triangles were pre-bent,
allowing a smooth structural re-configuration and avoiding abrupt changes in
the mechanical response. The in-plane and out-of-plane buckling behaviour of
square cellular plates was studied by Nikman and Akbarzadeh [2018], where it
was concluded that the slenderness ratio, defined as the ratio between the
length and the thickness of the specimen, governs the buckling mode and
therefore determines pattern switching.

A rigorous numerical study on pattern-transforming mechanical metamaterials is
provided by Ameen et al. [2018b]. The hyper-elastic specimens considered in
this study contain a structure of periodically distributed circular holes,
which upon compression transform into an arrangement of alternating ellipses.
It was found that the post-transformation material behaviour is strongly
influenced by the specimen size when the ratio between the specimen height and
cell size is relatively small. These size effects are induced by constraints
on the pattern transformation close to the sample’s fixed edges, which trigger
stiff boundary layers. In metamaterial design, as well as predictive modelling
involving mechanical metamaterials, one is interested in the effective
material properties. Moreover, engineering simulations of components or large
structures cannot be done with full account of the microstructure. A
homogenized approach is the most efficient solution method for engineering
purposes. To retrieve the actual effective material and kinematic behaviour
corresponding to a single specimen, Ameen et al. [2018b] employed ensemble
averaging. To acquire the ensemble averaged solution, the material behaviour
for a family of translated microstructures is required, meaning that a number
of boundary value problems needs to be solved. Since such a brute force
approach is computationally expensive, more efficient homogenization
techniques are called for. A rich variety of options has been presented in the
literature. For a well-characterized microstructure, however, computational
homogenization schemes are typically assumed to be among the most accurate
approaches [Geers et al., 2010]. In spite of its high accuracy, standard
first-order computational homogenization cannot capture size effects, since
the scale separation assumption does not allow for non-local behaviour [Ameen
et al., 2018a]. Enriched computational homogenization schemes are thus
required to achieve a higher accuracy, especially for specimens with small
differences between the length of the microstructural features and considered
macrostructural specimen size.

In order to reflect for non-locality and size effects, second-order
computational homogenization [Kouznetsova et al., 2002, 2004] incorporates a
generalized continuum formulation at the macro-scale, which makes it a serious
candidate for elastomeric mechanical metamaterials. Anthoine [2010] applied
second-order homogenization to functionally graded materials to quantify the
effect of microstructural grading. The intrinsic length scale which influences
the higher-order stiffness can be directly related to the graded
microstructure of the material. Forest and Trinh [2010] used the same method
to investigate the size effect of Representative Volume Elements (RVEs) for
elastic and elastoplastic composites subjected to quadratic Dirichlet boundary
conditions. The research conducted by Noels and Nguyen [2014] studied the
behaviour of elastic and elastoplastic materials containing a cellular
hexagonal microstructure through second-order computational homogenization.
Via comparison to full numerical simulations, it was concluded that this
approach accurately predicts the behaviour of structures consisting of
repeated hexagonal patterns, but the method fails for high degrees of
macroscopic imperfections. Generalization of second-order homogenization
method to include body forces avoiding spurious effects were carried out in
[Luscher et al., 2010, Blanco et al., 2016, Yvonnet et al., 2020], and in the
context of non-linear cellular metamaterials in [Wu et al., 2023]. An overview
comparison of second-order models has been presented by Lopes and Pires
[2022a, b].

Another suitable option for homogenization of elastomeric metamaterials is
micromorphic homogenization scheme introduced by Rokoš et al. [2019], which
incorporates kinematic coupling between individual cells through an additional
field variable. The essence of the adopted approach lies in the decomposition
of the displacement field into slow and fast varying components. In addition,
a predefined fluctuation field corresponding to the pattern transformation is
introduced, representing prior knowledge about the most important
microstructural deformation modes. The method has been compared to ensemble
averaged results for uniformly compressed structures by Ameen et al. [2018b],
where the accuracy of the method has been assessed for a macroscopically one-
dimensional problem. The methodology has been generalized to multiple
micromorphic fields in [Rokoš et al., 2020a], and multiscale buckling
addressed in [van Bree et al., 2020].

The aim of this paper is to compare the performance of the two above-mentioned
enriched computational homogenization schemes (i.e., the second-order and
micromorphic) in the context of elastomeric mechanical metamaterials
exhibiting pattern transformations against the reference ensemble-averaged
solution on three examples. Both methods include an enhanced continuum
description at the macro-scale, allowing to account for relevant size and
boundary effects. For completeness, results corresponding to the scale-
independent first-order computational homogenization scheme are included as
well.

The remainder of this paper is organized as follows. First, Section 2
specifies the elastomeric metamaterial considered, through its geometry,
constitutive laws, and basic mechanical behaviour. Section 3 introduces the
reference solution represented by the ensemble average of Direct Numerical
Simulations (DNS) corresponding to a family of translated microstructures. The
classical first-order computational homogenization method is then briefly
reviewed, followed by a description of the second-order and micromorphic
homogenization schemes. In Section 4, the considered load cases are presented
and the results obtained with the different homogenization methods are
discussed. Finally, overall summary and conclusions are provided in Section 5.

Throughout this paper, the following notational conventions are adopted:
scalars $a$, vectors $\vec{a}=a_{i}\vec{e}_{i}$, second-order tensors
$\boldsymbol{A}=A_{ij}\vec{e}_{i}\vec{e}_{j}$, third-order tensors
${}^{3}\boldsymbol{A}=A_{ijk}\vec{e}_{i}\vec{e}_{j}\vec{e}_{k}$, scalar
product $\vec{a}\cdot\vec{b}=a_{i}b_{i}$, single contraction
$\boldsymbol{A}\cdot\vec{b}=A_{ij}b_{j}\vec{e}_{i}$, double contraction
$\boldsymbol{A}:\boldsymbol{B}=A_{ij}B_{ji}$, transpose
$\boldsymbol{A}^{\mathrm{T}},A_{ij}^{\mathrm{T}}=A_{ji}$, gradient operator
$\vec{\nabla}\vec{a}=\frac{\partial a_{j}}{\partial
X_{i}}\vec{e}_{i}\vec{e}_{j}$, and divergence operator
$\vec{\nabla}\cdot\vec{a}=\frac{\partial a_{i}}{\partial X_{i}}$. Unless
indicated otherwise, summation over repeated indices is assumed for tensor
operations.

##  2 Metamaterial Specification

The geometry considered in this study is inspired by Specimen 1, as described
in [Bertoldi et al., 2008], which was adopted also in the study of Ameen et
al. [2018b]. A schematic representation of the cellular structure is provided
in Fig. 1a, illustrating the underlying square unit cell, its size $l$, and
hole diameter $d$. In this study, the adopted geometrical values are set as
$l=9.97$ mm and $d=8.67$ mm. The specimen’s width and height are denoted by
$W$ and $H$ and depend on the amount of cells in the horizontal or vertical
direction. Upon application of a compressive load, typically a bi-linear
stress–strain diagram is obtained, as shown in Fig. 1b. The diagram splits
into two parts by the vertical dashed line, indicating the bifurcation strain
of the microstructure. At this strain, individual circular voids transform
abruptly into ellipses with alternating horizontal and vertical major axes,
resulting in a significant drop of the effective stiffness. The pattern
transformation and resulting mechanical response have been thoroughly
investigated, e.g., by Mullin et al. [2007], Bertoldi et al. [2008], Hu et al.
[2013], or Coulais [2016]. Fig. 1c illustrates a $4\times 4$ part of an
infinite specimen, subjected to $7.5\%$ overall vertical compressive strain.
At that strain, the specimen is in the post-bifurcation regime, i.e., the
pattern is in the transformed state and the material’s mechanical response is
now related to the new microstructural setting (cf. Fig. 1b). When replacing
the periodicity conditions imposed on the $4\times 4$ subsection by prescribed
displacements, the fluctuating edges resulting from the pattern transformation
will be constrained. This explains the mechanism behind the formation of stiff
boundary layers, dictating the mechanical response for specimens containing a
low ratio between the specimen height and cell length ($H/l$). This ratio will
be referred to as the scale ratio in what follows.

(a) cellular structure

(b) stress–strain diagram

(c) deformation pattern

Figure 1: (a) Sketch of the two-dimensional elastomeric metamaterial,
including the geometrical parameters, and (b) corresponding stress–strain
response in compression. (c) Deformation pattern of a $4\times 4$ periodic
specimen subjected to $7.5\%$ of overall compressive strain; the colour
indicates the norm of the Green–Lagrange strain $\boldsymbol{E}$, defined as
$\epsilon=\sqrt{\boldsymbol{E}:\boldsymbol{E}}$.

The elastomeric base material, employed throughout all examples discussed in
Section 4, is a hyperelastic material. The strain energy density is described
by a two-term $I_{1}$-based compressible Mooney–Rivlin model, which in a two-
dimensional setting reads as

$W(\boldsymbol{F}(\vec{X}))=c_{1}(I_{1}-2)+c_{2}(I_{1}-2)^{2}-2c_{1}\mathrm{log}J+\frac{K}{2}(J-1)^{2},$
(1)

where $I_{1}$ indicates the first invariant of the right Cauchy–Green tensor
$\boldsymbol{C}=\boldsymbol{F}^{\mathrm{T}}\cdot\boldsymbol{F}$, and
$J=\det{({\boldsymbol{F}})}$ represents the volume change ratio. Constants
$c_{1}$, $c_{2}$, and $K$ are constitutive parameters, set to $c_{1}=0.55$
MPa, $c_{2}=0.3$ MPa, and $K=55$ MPa, according to the experimental data of
Bertoldi et al. [2008].

##  3 Homogenization Methodologies

###  3.1 Reference Solution Based on Ensemble Averaging

In the case of heterogeneous materials, the DNS results typically reveal
highly oscillatory kinematic fields, which tend to shift periodically with the
microstructural position relative to the macroscopic sample. To obtain a
representative reference solution to assess the accuracy of computational
homogenization methods, an ensemble averaged DNS solution of a family of
translated microstructures is established. Averaging the deformation fields
over all translated microstructures filters out the fast oscillations and
results in a macroscopically smooth solution. For constructing this reference
solution, each translated microstructure is assumed to have equal probability
of occurrence, as employed by, e.g., Drugan and Willis [1996], Smyshlyaev and
Cherednichenko [2000], or Ameen et al. [2018a], relying on the assumption that
the exact positioning of the microstructure is not known a-priori.

In Fig. 2a, a schematic representation of the macroscopic specimen in the
reference configuration is given, with domain $\Omega$, outer boundary
$\Gamma$, and position vector $\vec{X}$. The vector $\vec{\zeta}$ is related
to the translated positions of the microstructure with respect to its
reference position. For the ensemble averaged (reference) solution, one needs
to average over the recurrent periodic pattern, which consists of $2\times 2$
cells in the considered case. Consequently, the domain including the required
microstructural translations spans $\vec{\zeta}\in Q=[-l,l]\times[-l,l]$. The
reference solution is then obtained as

$\vec{\overline{u}}(\vec{X})=\frac{1}{|Q|}\int_{Q}\vec{u}(\vec{X},\vec{\zeta})\,\mathrm{d}\vec{\zeta},$
(2)

where $\vec{u}(\vec{X},\vec{\zeta})$ denotes the displacement field
corresponding to the microstructure translated by $\vec{\zeta}$. For
convenience, the reference homogenized solution $\vec{\overline{u}}(\vec{X})$
is computed numerically by discretizing the integral of Eq. (2) into a set of
grid points spanning $Q$, and calculating the DNS solution
$\vec{u}(\vec{X},\vec{\zeta})$ for each realization.

(a) microstructural translation

(b) first-order computational homogenization

Figure 2: Sketch of a macroscopic problem domain $\Omega$ with boundary
$\Gamma$; Dirichlet and Neumann boundary conditions are applied on boundary
segments $\Gamma_{\mathrm{D}}$ and $\Gamma_{\mathrm{N}}$. (a) Illustration of
a microstructural translation $\vec{\zeta}$ (in blue) with respect to the
reference configuration (black). (b) A sketch of the computational
homogenization scheme; the microscale problem in the reference (bottom right)
and deformed (bottom left) configuration considered at a macroscopic point of
interest $\vec{X}$.

###  3.2 First-Order Computational Homogenization

Computational homogenization is a multiscale technique, which incorporates
detailed microstructural information by calculating the microscale and macro-
scale boundary value problem in a nested manner [Geers et al., 2010]. This
method provides a powerful numerical tool for macroscopically homogeneous and
microscopically heterogeneous materials, for which the macroscopic
constitutive law is not available in a closed form, but implicitly through the
solution of an underlying microstructural problem. The macroscopic first
Piola–Kirchhoff stress tensor can be expressed as

$\boldsymbol{P}=\boldsymbol{f}\big{(}\vec{X},\boldsymbol{P}_{\mathrm{m}}\big{(}\boldsymbol{F}(\vec{X})\big{)}\big{)},\quad\vec{X}\in\Omega,$
(3)

where $\boldsymbol{F}$ denotes the macroscopic deformation gradient and
$\boldsymbol{P}_{\mathrm{m}}$ the microscopic first Piola–Kirchhoff stress
tensor. From Eq. (3) it is clear that the macroscopic stress is related to the
microscopic stress resulting from the macroscopic deformation. In Fig. 2b a
schematic representation of the computational homogenization concept applied
to elastomeric mechanical metamaterials is given. Note that two nested
boundary value problems are involved: (i) a macroscopic problem (Fig. 2b top)
governed by the classical balance equation, i.e.,

$\vec{\nabla}\cdot\boldsymbol{P}^{\mathrm{T}}(\vec{X})=\vec{0},\quad\vec{X}\in\Omega,$
(4)

and (ii) a microscopic problem (Fig. 2b bottom) with a similar balance
equation

$\vec{\nabla}_{\mathrm{m}}\cdot\boldsymbol{P}_{\mathrm{m}}^{\mathrm{T}}(\vec{X}_{\mathrm{m}})=\vec{0},\quad\vec{X}_{\mathrm{m}}\in\Omega_{\mathrm{m}}.$
(5)

In Eq. (5), $\vec{X}_{\mathrm{m}}$ denotes the microscopic spatial position
vector and $\Omega_{\mathrm{m}}$ the microscopic domain. Recall that,
throughout this paper, the macroscopic $\vec{\nabla}$ and microscopic
$\vec{\nabla}_{\mathrm{m}}$ differential operators are defined with respect to
the undeformed configurations.

To obtain the macroscopic first Piola–Kirchhoff stress tensor $\boldsymbol{P}$
at a point of interest $\vec{X}$, the corresponding macroscopic deformation
gradient tensor $\boldsymbol{F}({\vec{X}})$ is imposed on the microscopic
problem. The relation between the microscopic deformation field and
macroscopic deformation gradient is given by

$\Delta\vec{x}_{\mathrm{m}}=\boldsymbol{F}\cdot\Delta\vec{X}_{\mathrm{m}}+\vec{w}(\Delta\vec{X}_{\mathrm{m}}),$
(6)

which corresponds to a Taylor expansion of the macroscopic deformation field
truncated after its leading term. In Eq. (6), $\Delta\vec{x}_{\mathrm{m}}$ and
$\Delta\vec{X}_{\mathrm{m}}$ are the finite line elements in the deformed and
reference configuration of the microscopic problem, whereas $\vec{w}$
represents the local zero-mean micro-fluctuation correction field. The
boundary conditions imposed on the microscopic problem, also referred to as
Representative Volume Element (RVE), are a combination of periodicity
conditions,

$\vec{w}(\Gamma_{\mathrm{m}_{\mathrm{T}}})=\vec{w}(\Gamma_{\mathrm{m}_{\mathrm{B}}})\quad\text{and}\quad\vec{w}(\Gamma_{\mathrm{m}_{\mathrm{L}}})=\vec{w}(\Gamma_{\mathrm{m}_{\mathrm{R}}}),$
(7)

and prescribed displacements resulting from the macroscopic deformation
gradient,

$\vec{x}_{\mathrm{m},i}=\boldsymbol{F}\cdot\vec{X}_{\mathrm{m},i}\quad\text{for}\quad
i=1,2,4,$ (8)

see Fig. 3. In Eq. (7), $\Gamma_{\mathrm{m}_{\bullet}}$ denotes individual RVE
boundary segments, and $\bullet$ represents right (R), left (L), top (T), and
bottom (B) RVE boundaries. Spatial positions $\vec{X}_{\mathrm{m},i}$,
$i=1,2,4$, denote in Eq. (8) three RVE corner points, as indicated in Fig. 3
where a schematic representation of a rectangular RVE domain is shown. Note
that equivalently to fixing the three control points in Eq. (8), $\vec{w}$ may
be considered periodic on top of $\boldsymbol{F}\cdot\vec{X}_{\mathrm{m}}$
while being orthogonal with respect to a constant inside the RVE domain
$\Omega_{\mathrm{m}}$ to remove rigid body motions [see, e.g., Miehe and
Bayreuther, 2007].

After solving the microscopic RVE problem according to Eq. (5), the
macroscopic stress at a material point $\vec{X}$ is obtained by volume
averaging the microscopic stress field over $\Omega_{\mathrm{m}}$, unlike
ensemble averaging considered for the underlying solution of Section 3.1,
i.e.,

$\boldsymbol{P}(\vec{X})=\frac{1}{|\Omega_{\mathrm{m}}|}\int_{\Omega_{\mathrm{m}}}\boldsymbol{P}_{\mathrm{m}}\,\mathrm{d}\Omega_{\mathrm{m}},\quad\vec{X}\in\Omega,$
(9)

which is consistent with the Hill–Mandel condition. This process can be
repeated for each macroscopic point of interest. Generally, the micro- and
macro-boundary value problems are conveniently solved using the Finite Element
Method, where each macroscopic integration point on $\Omega$ is related to a
RVE. The macroscopic problem of Eq. (4) can be solved using a standard Newton
iterative algorithm, resulting in an integrated scheme in which the micro- and
macro-boundary value problems are solved in a nested manner.

###  3.3 Second-Order Computational Homogenization

The second-order homogenization method may be regarded as an extension of its
first-order counterpart, relying on a Taylor series expansion of the
deformation field truncated after the second term. The local deformation field
within a RVE then takes the form (cf. Eq. (6))

$\Delta\vec{x}_{\mathrm{m}}=\boldsymbol{F}\cdot\Delta\vec{X}_{\mathrm{m}}+\frac{1}{2}\Delta\vec{X}_{\mathrm{m}}\cdot{}^{3}\boldsymbol{G}\cdot\Delta\vec{X}_{\mathrm{m}}+\vec{w}(\Delta\vec{X}_{\mathrm{m}}),$
(10)

where the second term on the right-hand side accounts for non-local strain
gradient effects through the gradient of the deformation gradient,
${}^{3}{\boldsymbol{G}}=\vec{\nabla}{\boldsymbol{F}}$. Such a formulation
requires an extended continuum framework, which is here given by Mindlin’s
strain gradient continuum at the macro level, cf., e.g., Mindlin [1965] or
Mindlin and Eshel [1968]. Following the lines of reasoning presented by
Kouznetsova et al. [2004], a classical continuum of Eq. (5) can be adopted at
the microscale level.

Figure 3: Schematic representation of the adopted square RVE with domain
$\Omega_{\mathrm{m}}$, consisting of $2\times 2$ unit cells. The RVE has four
corner points, $1,\dots,4$, four boundary segments
$\Gamma_{\mathrm{m}_{\bullet}}$ with associated outer normals
$\vec{N}_{\mathrm{m}_{\bullet}}$, where $\bullet$ represents right (R), left
(L), top (T), and bottom (B) edges, and is parametrized by two normalized
coordinates $\xi\in[-1,1]$ and $\eta\in[-1,1]$.

Introducing the higher-order stress, ${}^{3}\boldsymbol{Q}$, which is work
conjugate to ${}^{3}\boldsymbol{G}$, the non-classical balance equation reads

$\vec{\nabla}\cdot\big{(}\boldsymbol{P}^{\mathrm{T}}-(\vec{\nabla}\cdot{}^{3}\boldsymbol{Q})^{\mathrm{T}}\big{)}=\vec{0},\quad\vec{X}\in\Omega.$
(11)

In Eq. (11), body forces as well as Neumann boundary conditions have been
neglected for simplicity. It may be clear from Eq. (11) that $C^{1}$
continuity is required over the entire domain $\Omega$ to obtain the gradient
of the deformation gradient at every point [Lesičar et al., 2014]. From a
finite element perspective such a requirement poses an additional
complication, since the conventional shape functions typically possess $C^{0}$
continuity only. In what follows, this is addressed by resorting to a mixed
finite element formulation, which introduces a so-called relaxed deformation
gradient field $\widehat{{\boldsymbol{F}}}$. The gradient of the displacement
and the relaxed deformation gradient field are coupled via Lagrange
multipliers, requiring that
$\widehat{{\boldsymbol{F}}}={\boldsymbol{I}}+(\vec{\nabla}\vec{u})^{\mathrm{T}}$
holds in the weak sense; for more details see Xia and Hutchinson [1996].

In addition to the macroscopic first Piola–Kirchhoff stress tensor, specified
in Eqs. (3) and (9), the higher-order macroscopic stress tensor is computed
from each RVE solution as

${}^{3}\boldsymbol{Q}=\frac{1}{2|\Omega_{\mathrm{m}}|}\int_{\Omega_{\mathrm{m}}}\big{(}\boldsymbol{P}^{\mathrm{T}}_{\mathrm{m}}\vec{X}_{\mathrm{m}}+\vec{X}_{\mathrm{m}}\boldsymbol{P}_{\mathrm{m}}\big{)}\,\mathrm{d}\Omega_{\mathrm{m}},$
(12)

which does not involve higher-order microscopic stresses since the RVE obeys a
classical continuum description. This choice effectively circumvents the need
for a higher-order characterization of the microstructural problem.

However, to enforce the gradient effects occurring at the macroscopic scale on
the micro problem, specific RVE boundary conditions are required, as
elaborated in [Kouznetsova et al., 2004]. The starting point for deriving the
microscopic RVE boundary conditions is Eq. (10), which provides the
microscopic deformation gradient in the form

$\boldsymbol{F}_{\mathrm{m}}=(\vec{\nabla}_{\mathrm{m}}\Delta\vec{x}_{\mathrm{m}})^{\mathrm{T}}=\boldsymbol{F}+\Delta\vec{X}_{\mathrm{m}}\cdot{}^{3}\boldsymbol{G}+(\vec{\nabla}_{\mathrm{m}}\vec{w})^{\mathrm{T}}.$
(13)

To conform to the classical averaging theorem given by

$\boldsymbol{F}=\frac{1}{|\Omega_{\mathrm{m}}|}\int_{\Omega_{\mathrm{m}}}\boldsymbol{F}_{\mathrm{m}}\,\mathrm{d}\Omega_{\mathrm{m}},$
(14)

two requirements have to be satisfied. First,

$\frac{1}{|\Omega_{\mathrm{m}}|}\int_{\Omega_{\mathrm{m}}}\Delta\vec{X}_{\mathrm{m}}\,\mathrm{d}\Omega_{\mathrm{m}}=\vec{0},$
(15)

which is guaranteed by simply shifting the coordinate system to the RVE
centre. The second requirement reads

$\frac{1}{|\Omega_{\mathrm{m}}|}\int_{\Omega_{\mathrm{m}}}(\vec{\nabla}\vec{w})^{\mathrm{T}}\,\mathrm{d}\Omega_{\mathrm{m}}=\frac{1}{|\Omega_{\mathrm{m}}|}\int_{\Gamma_{\mathrm{m}}}\vec{w}\vec{N}_{\mathrm{m}}\,\mathrm{d}\Gamma_{\mathrm{m}}=\boldsymbol{0},$
(16)

where the divergence theorem has been used to transform the volume integral
into a contour integral along $\Gamma_{\mathrm{m}}$, with
$\vec{N}_{\mathrm{m}}$ denoting its outward unit normal, cf. Fig. 3. Assuming
${}^{3}{\boldsymbol{G}}={}^{3}{\boldsymbol{0}}$ yields
$\vec{N}_{\mathrm{m}_{\mathrm{T}}}(\Gamma_{\mathrm{m}_{\mathrm{T}}})=-\vec{N}_{\mathrm{m}_{\mathrm{B}}}(\Gamma_{\mathrm{m}_{\mathrm{B}}})$
and
$\vec{N}_{\mathrm{m}_{\mathrm{R}}}(\Gamma_{\mathrm{m}_{\mathrm{R}}})=-\vec{N}_{\mathrm{m}_{\mathrm{L}}}(\Gamma_{\mathrm{m}_{\mathrm{L}}})$,
which in turn reduces the requirement of Eq. (16) to the standard periodicity
condition of Eq. (7). When ${}^{3}\boldsymbol{G}\neq{}^{3}\boldsymbol{0}$,
however, a non-periodic deformation results, requiring periodicity in a
generalized sense, defined as

$\displaystyle\vec{x}_{\mathrm{m}_{\mathrm{R}}}$
$\displaystyle=\vec{x}_{\mathrm{m}_{\mathrm{L}}}+2l\boldsymbol{F}\cdot\vec{N}_{\mathrm{m}_{\mathrm{R}}}+2l^{2}\eta\vec{N}_{\mathrm{m}_{\mathrm{R}}}\cdot{}^{3}\boldsymbol{G}\cdot\vec{N}_{\mathrm{m}_{\mathrm{T}}},$
(17) $\displaystyle\vec{x}_{\mathrm{m}_{\mathrm{T}}}$
$\displaystyle=\vec{x}_{\mathrm{m}_{\mathrm{B}}}+2l\boldsymbol{F}\cdot\vec{N}_{\mathrm{m}_{\mathrm{T}}}+2l^{2}\xi\vec{N}_{\mathrm{m}_{\mathrm{T}}}\cdot{}^{3}\boldsymbol{G}\cdot\vec{N}_{\mathrm{m}_{\mathrm{R}}},$
(18)

where $\xi$ and $\eta$ are normalized RVE coordinates, cf. Fig. 3. The four
corner points are displaced according to Eq. (10) in which
$\vec{w}(\vec{X}_{\mathrm{m},i})=\vec{0}$, $i=1,\dots,4$, similarly to Eq.
(8). Since this requirement may lead to stress concentrations close to RVE
corner points [Lesičar et al., 2014, Lopes and Pires, 2022a], prescribed
corner point displacements may be relaxed and replaced by a condition
requiring orthogonality of $\vec{w}$ with respect to a constant [Luscher et
al., 2010]. In the examples that follow, the former option according to
Kouznetsova et al. [2002] has been adopted. Finally, to relate the microscopic
variables to the macroscopic gradient of the deformation gradient, the
following conditions must hold

$\int_{\Gamma_{\mathrm{m_{L}}}}\vec{w}\,\mathrm{d}\Gamma_{\mathrm{m}}=\vec{0}\quad\text{and}\quad\int_{\Gamma_{\mathrm{m_{B}}}}\vec{w}\,\mathrm{d}\Gamma_{\mathrm{m}}=\vec{0},$
(19)

which are consistent with the generalized periodicity of $\vec{w}$ on the
right and top boundaries, i.e., Eqs. (17) and (18).

###  3.4 Micromorphic Computational Homogenization

A micromorphic continuum theory is based on the introduction of additional
kinematic variables capturing the microstructural displacement field, which
include non-local effects on the macro level. The Cosserat continuum [Cosserat
and Cosserat, 1909] is a micromorphic continuum where the additional kinematic
variables take the microstructural rotations into account. In [Rokoš et al.,
2019], the micromorphic theory is used to construct a homogenization framework
tailored to cellular elastomeric metamaterials. Here, the characteristic
deformation mode, corresponding to the pattern transformation, is extracted
from the displacement field, and an additional variable controlling its
magnitude is introduced.

The following decomposition ansatz of the kinematic field then holds

$\vec{u}(\vec{X},\vec{\zeta})=\vec{v_{0}}(\vec{X})+\sum^{n}_{i=1}v_{i}(\vec{X})\vec{\varphi_{i}}(\vec{X},\vec{\zeta})+\vec{w}(\vec{X},\vec{\zeta}),$
(20)

in which the vector function $\vec{v}_{0}(\vec{X})$ corresponds to the mean
effective displacement field, $v_{i}(\vec{X})$ are scalar fields representing
the magnitudes of the predetermined long-range correlated patterning modes
$\vec{\varphi}_{i}(\vec{X},\vec{\zeta})$, and $\vec{w}(\vec{X},\vec{\zeta})$
is the remaining microfluctuation field. The notion of ensemble averaging is
adopted, recall Eq. (2) and the discussion therein, which is reflected in Eq.
(20) through the dependence of the micro-fluctuating quantities on the
microstructural translation $\vec{\zeta}\in Q$. For the problem considered
here, only a single spatially correlated mode is considered (i.e., $n=1$,
shown in Fig. 1c), which is considered as an input quantity (prior kinematic
knowledge), and can be obtained through a Bloch-type analysis [Bertoldi et
al., 2008], estimated analytically from full-scale numerical simulations
[Rokoš et al., 2019], or identified experimentally [Maraghechi et al., 2024].

To establish a computational homogenization scheme, a Taylor expansion of the
coarse fields $\vec{v}_{0}$ and $v_{1}$ is performed in analogy to the first-
and second-order computational methods (cf. Eqs. (6) and (10)), yielding

$\Delta\vec{x}_{\mathrm{m}}={\boldsymbol{F}}\cdot\Delta\vec{X}_{\mathrm{m}}+\big{[}v_{1}+\Delta\vec{X}_{\mathrm{m}}\cdot\vec{\nabla}v_{1}\big{]}\vec{\varphi}_{1}(\Delta\vec{X}_{\mathrm{m}})+\vec{w}(\Delta\vec{X}_{\mathrm{m}}),$
(21)

where the translation variable $\vec{\zeta}\in Q$ has been replaced by the
spatial variable $\Delta\vec{X}_{\mathrm{m}}\in\Omega_{\mathrm{m}}$,
effectively replacing the expensive ensemble averaging by substantially
cheaper volume averaging. Note that since the patterning mode
$\vec{\varphi}_{1}$ is defined with respect to the reference microscopic
configuration $\vec{X}_{\mathrm{m}}$, macroscopic rotation needs to be
factored out especially in cases it is non-negligible. This can be achieved
through polar decomposition of the macroscopic deformation gradient
$\boldsymbol{F}=\boldsymbol{R}\cdot\boldsymbol{U}$, where $\boldsymbol{R}$ is
the macroscopic rotation tensor and $\boldsymbol{U}$ the macroscopic stretch
tensor. This effectively means that the deformation gradient $\boldsymbol{F}$
is replaced with $\boldsymbol{U}$ in Eq. (21), see also [van Bree et al.,
2020, Section 2.1]. Unlike the second-order method, the non-local information
of the micromorphic model is contained in the magnitude of the patterning mode
$v_{1}$ and its spatial variation $\vec{\nabla}v_{1}$. For more details see
also [Rokoš et al., 2019, Section 5.1]. The kinematic approximation of Eq.
(21) is considered over a local neighbourhood of a macroscopic point
$\vec{X}$, spanned by the adopted RVE, over which the micro-fluctuation field
$\vec{w}$ can be computed. Inserting the decomposition of Eq. (21) into the
internal virtual work, two macroscopic balance equations result,

$\displaystyle\vec{\nabla}\cdot\boldsymbol{\Theta}^{\mathrm{T}}$
$\displaystyle=\vec{0},\quad\vec{X}\in\Omega,$ (22)
$\displaystyle\vec{\nabla}\cdot\vec{\Lambda}-\Pi$
$\displaystyle=0,\quad\vec{X}\in\Omega,$ (23)

along with the following homogenized stress quantities

$\displaystyle\boldsymbol{\Theta}(\vec{X})$
$\displaystyle=\frac{1}{|\Omega_{\mathrm{m}}|}\int_{\Omega_{\mathrm{m}}}\boldsymbol{P}_{\mathrm{m}}\,\mathrm{d}\Omega_{\mathrm{m}},$
(24) $\displaystyle\Pi(\vec{X})$
$\displaystyle=\frac{1}{|\Omega_{\mathrm{m}}|}\int_{\Omega_{\mathrm{m}}}\boldsymbol{P}_{\mathrm{m}}:\vec{\nabla}_{\mathrm{m}}\vec{\varphi}_{1}\,\mathrm{d}\Omega_{\mathrm{m}},$
(25) $\displaystyle\vec{\Lambda}(\vec{X})$
$\displaystyle=\frac{1}{|\Omega_{\mathrm{m}}|}\int_{\Omega_{\mathrm{m}}}\boldsymbol{P}_{\mathrm{m}}^{\mathrm{T}}\cdot\vec{\varphi}_{1}+\vec{X}_{\mathrm{m}}[\boldsymbol{P}_{\mathrm{m}}:\vec{\nabla}_{\mathrm{m}}\vec{\varphi}_{1}]\,\mathrm{d}\Omega_{\mathrm{m}}.$
(26)

At the microscale the classical balance law of Eq. (5) applies, which requires
in addition to periodicity of $\vec{w}$ in Eq. (7) also orthogonality of
$\vec{w}$ with respect to a constant field, fluctuation mode
$\vec{\varphi}_{1}$, and linear component
$\vec{\varphi}_{1}\vec{X}_{\mathrm{m}}$, to guarantee uniqueness of the
solution. Unlike the second-order homogenization scheme, the four RVE corner
points are free, and rigid body motions are restricted by a combination of
periodicity conditions and the orthogonality with respect to a constant field.
For more details on implementation and theory see van Bree et al. [2020] and
Rokoš et al. [2020a].

###  3.5 Concluding Remarks

Each of the presented homogenization schemes relies on a classical continuum
at the microscale level. The first-order method employs a classical continuum
description also at the macro-scale level, while the second-order method
extends the balance equation by including higher-order stresses arising from
strain gradients. The micromorphic method, on the contrary, introduces an
additional scalar balance law governing the macroscopic magnitude of the
fluctuation pattern. Whereas the first-order and second-order methods require
no prior information on the considered kinematics, the micromorphic scheme
assumes the typical long-range correlated patterning field to be known in
advance for the adopted kinematic ansatz.

##  4 Numerical Examples and Comparison

This section provides a performance analysis of the homogenization
methodologies discussed in Section 3, applied to pattern transforming
elastomeric metamaterials. The accuracy of the enriched homogenization
solutions is assessed relative to the reference ensemble-averaged solution,
whereas the first-order homogenization solution is considered as the current
standard (most widely used). Three distinct loading cases are evaluated,
comparing mechanical responses obtained by each of the methods. To this end,
first an infinitely wide specimen loaded in compression is considered. Second,
a specimen subjected to bending is evaluated. A finite rectangular specimen,
subjected to a compressive load with two free boundaries experiencing an
auxetic effect, is analysed as the final example.

###  4.1 Uniform Compression

The first example considers an infinite layer of elastomeric material loaded
in compression along the vertical direction. The specimen’s height is
$H=n_{v}l$, where $n_{v}$ is the number of unit cells in the vertical
direction. The simulation domain spans the region $2l\times H$, as shown in
Fig. 4, whereby only two unit cells are modelled in the horizontal direction,
which is sufficient to capture the complete pattern transformation. The
infinite structure is mimicked through periodicity along the two vertical
edges, AD and BC. The bottom edge AB is fixed, whereas the top edge CD is
displaced downwards by $u^{*}$. As a result of these kinematic constraints,
the macroscopic problem can be modelled as a one-dimensional continuum along
$\vec{e}_{2}$, with the following essential boundary conditions. For the
first-order computational homogenization:

$u_{2}(0)=0,\,u_{2}(H)=u^{*},$ (27)

for the second-order method:

$u_{2}(0)=0,\,u_{2}(H)=u^{*},\quad F_{22}(0)=F_{22}(H)=1,$ (28)

and for the micromorphic approach:

$v_{0,2}(0)=0,\,v_{0,2}(H)=u^{*},\quad v_{1}(0)=v_{1}(H)=0.$ (29)

Note that the choice of $F_{22}=1$ in Eq. (28) is motivated by physical
insights of no deformation at the top and bottom part of the domain boundary.
Although not presented for conciseness, numerical investigations showed high
influence of prescribed values of $F_{22}$ on achieved accuracy in comparison
with DNS. For all methods, the macroscopic displacement fields are discretized
by piece-wise quadratic one-dimensional elements with a two-point Gauss
integration rule. For the second-order method, the Lagrange multiplier field
and relaxed deformation gradient field are discretized using piece-wise linear
shape functions. Scale ratios up to $H/l=10$ are calculated using $10$
macroscopic one-dimensional elements, whereas simulations at higher scale
ratios have as many elements as unit cells, i.e., $H/l$. The refined
discretization for scale ratios $H/l<10$ is adopted to ensure sufficiently
smooth and converged kinematic fields.

Figure 4: Uniform compression case. A sketch of the considered simulation
domain and boundary conditions for a specimen with scale ratio $H/l=4$.
Periodic boundary conditions are enforced along the two vertical edges AD and
BC, representing an infinitely wide structure.

(a) vertical displacement

(b) vertical deformation

Figure 5: Uniform compression case. Comparison of (a) vertical displacement
component $u_{2}(X_{2})$, and (b) deformation gradient $F_{22}(X_{2})$ over
the specimen’s height. Overall applied compressive strain corresponds to $1\%$
and $7.5\%$, for a specimen with scale ratio $H/l=6$.

|
---|---
|

(a) DNS

(b) first-order

(c) second-order

(d) micromorphic

Figure 6: Uniform compression case. Deformed configuration of an infinitely
wide specimen with scale ratio $H/l=6$, subjected to $7.5\%$ nominal
compressive strain. (a) DNS solution corresponding to zero microstructural
translation $\vec{\zeta}=\vec{0}$. $F_{22}$ component of the macroscopic
deformation gradient plotted against the vertical coordinate in the reference
configuration, including deformed RVEs of the first three integration points,
corresponding to the (b) first-order, (c) second-order, and (d) micromorphic
computational homogenization schemes.

(a) $H/l=6$

(b) $H/l=12$

Figure 7: Uniform compression case. Stress–strain diagram for scale ratios (a)
$H/l=6$ and (b) $H/l=12$, obtained via DNS and the first-order, second-order,
and micromorphic homogenization schemes.

(a) $u^{*}/H=0.01$

(b) $u^{*}/H=0.075$

Figure 8: Uniform compression case. Nominal compressive stress in the (a) pre-
bifurcation and (b) post bifurcation regime as a function of the scale ratio
$H/l\in[4,36]$.

(a) $u^{*}/H=0.01$

(b) $u^{*}/H=0.075$

Figure 9: Uniform compression example. Tangent stiffness in the (a) pre-
bifurcation and (b) post-bifurcation regime as a function of the scale ratio
$H/l\in[4,36]$.

The macroscopic kinematic quantities in the pre- and post-bifurcation regimes
for a fixed scale ratio $H/l=6$ are compared in Fig. 5. Here, the reference
solution indicates that during pre-bifurcation, no strain gradients are
present, which is accurately captured by the micromorphic and first-order
method. The second-order scheme exhibits, however, artificial gradient effects
close to the two horizontal boundaries. Post-bifurcation, the micromorphic
method shows a better deformation trend compared to the second-order scheme,
which again exhibits too strong gradient effects, which are mainly caused by
the kinematic constraints on $F_{22}$ in Eq. (28). Note that
$F_{22}(0)=F_{22}(H)\approx 0.95$ in Fig. 5b for the reference solution,
implying that the kinematic constraints in the second-order scheme are too
restrictive. However, relaxing these boundary conditions introduces large
errors in the post-bifurcation regime, resulting in an overall response that
is too compliant.

Fig. 6 presents a deformed DNS configuration along with three RVEs positioned
close to the specimen’s bottom edge for each homogenization method. The
gradient effect on the deformation of the individual RVEs is clearly visible
close to the two boundaries. The first-order method (Fig. 6b) does not include
strain gradients and provides the same compressed RVE pattern at each
integration point, even close to the boundary. The RVE deformations obtained
by the second-order method (Fig. 6c), show a gradual compression from the
bottom up. However, the pattern is not representative for what is observed in
the DNS, especially for the bottom RVE near the boundary. The micromorphic
RVEs, shown in Fig. 6d, reproduce the DNS deformation adequately.

The obtained nominal stress–strain diagrams are provided in Fig. 7 for scale
ratios $6$ and $12$. The reference solution shows that the bifurcation strain
(indicated by the black dashed vertical line) decreases with increasing scale
ratio, and that the bifurcation strain is captured accurately by the
micromorphic homogenization method (the error in bifurcation stress is
approximately $2.8\%$ and $0.82\%$ for scale ratios $6$ and $12$). The first-
order method, which is insensitive to the macroscopic size, provides a
bifurcation strain $u^{*}/H\approx 0.0225$. Interestingly, the second-order
method exhibits an opposite effect, in which the bifurcation strain increases
with increasing scale ratio. The final buckling strain corresponds to
$u^{*}/H\approx 0.0225$, which is expected since the first-order solution
should be recovered for $H/l\rightarrow\infty$.

The effect of the macroscopic size on the stress response is shown in Fig. 8,
where the homogenized stresses at fixed nominal strains $u^{*}/H=0.01$ and
$u^{*}/H=0.075$ are plotted against the scale ratio $H/l\in[4,36]$. Before the
pattern transformation occurs, no non-local effects (and thus almost no size
effects) can be observed in the reference, first-order, and micromorphic
homogenized solutions (Fig. 8a). This is not the case, however, for the
second-order homogenization scheme, which shows an excessively stiff response
in the pre-bifurcation regime due to the kinematic constraint $F_{22}=1$. In
the post-bifurcation regime, the second-order and micromorphic homogenization
schemes capture the overall trend, but both underestimate the effective
stress. For the smallest scale ratio considered, i.e., $H/l=4$, the first-
order homogenization entails approximately $30\%$, second-order $11.5\%$, and
micromorphic $7.5\%$ relative error.

The effective stiffnesses, considered again at $u^{*}/H=0.01$ and
$u^{*}/H=0.075$ of the overall nominal strain and expressed as a function of
the scale ratio $H/l,$ are shown in Fig. 9. Pre-bifurcation, similar trends as
for the nominal stresses are observed (cf. Fig. 8a), whereas post-bifurcation
(Fig. 9b), observed differences are approximately equal for both non-local
schemes. For scale ratios $H/l>10$, nevertheless, the second-order method
follows the reference trend more accurately compared to the micromorphic
method.

###  4.2 Uniform Bending

Similarly to the previous example, a semi-infinite domain is considered, which
is subjected to uniform bending. The simulation domain again spans a region
$2l\times H$, where an integer number of unit cells along the vertical
direction is considered, cf. Fig. 10. The two vertical edges, AD and BC,
deform periodically on top of the uniform bending contribution. The coordinate
system is shifted to the midpoint of the simulation domain such that the
neutral axis corresponds to $X_{2}=0$. The neutral axis is furthermore fixed
in the vertical direction by constraining horizontal and vertical
displacements at the centre left and right nodes, whereas the top and bottom
edges are considered as free surfaces. Mathematically, the applied boundary
conditions are expressed as follows. For the first-order method:

$\vec{u}(\Gamma_{\mathrm{AD}})=\vec{u}(\Gamma_{\mathrm{BC}})+\theta
X_{2}\vec{e}_{1},$ (30)

for the second-order method:

$\vec{u}(\Gamma_{\mathrm{AD}})=\vec{u}(\Gamma_{\mathrm{BC}})+\theta
X_{2}\vec{e}_{1},\quad
F_{11}(\Gamma_{\mathrm{AD}})=F_{11}(\Gamma_{\mathrm{BC}}),\,F_{22}(\Gamma_{\mathrm{AD}})=F_{22}(\Gamma_{\mathrm{BC}}),$
(31)

and for the micromorphic method:

$\vec{v}_{0}(\Gamma_{\mathrm{AD}})=\vec{v}_{0}(\Gamma_{\mathrm{BC}})+\theta
X_{2}\vec{e}_{1},\quad
v_{1}(\Gamma_{\mathrm{AD}})=v_{1}(\Gamma_{\mathrm{BC}}).$ (32)

The bending angle $\theta$ is prescribed such that the nominal strain at the
top and bottom edges (expressed as $\theta H/(4l)$) does not exceed $15\%$.
The macroscopic discretization consists of $8$-node two-dimensional square
quadratic quadrilateral elements with edge size $l$, integrated by a $2\times
2$ Gauss integration rule. The relaxed deformation gradient field
$\widehat{{\boldsymbol{F}}}$, present in the second-order homogenization
method, is discretized by $4$-node bi-linear quadrilaterals. More details and
performance analysis of various two-dimensional element types for micromorphic
scheme can be found in Rokoš et al. [2020b], and for second-order
computational homogenization, e.g., in Kouznetsova et al. [2004].

Figure 10: Uniform bending case. A sketch of the simulation domain and
boundary conditions for a specimen with scale ratio $H/l=4$. Superposed
periodic boundary conditions are enforced along the vertical edges AD and BC
to mimic an infinitely wide structure.

|
---|---
|

(a) DNS

(b) first-order

(c) second-order

(d) micromorphic

Figure 11: Uniform bending case. An infinitely wide specimen with scale ratio
$6$ bent up to $\theta H/(4l)=0.1125$. (a) DNS solution with zero
microstructural shift $\vec{\zeta}=\vec{0}$. The macroscopic deformed
structures including deformed RVEs corresponding to three macroscopic
integration points for the (b) first-order, (c) second-order, and (d)
micromorphic homogenization method. The dotted blue line indicates the
deformed contour of the ensemble-averaged solution.

(a) stress–strain diagram

(b) homogenized stress $P_{11}$

Figure 12: Uniform bending case corresponding to scale ratio $12$. (a)
Stress–strain diagram; a close-up of the bifurcation point in the inset. (b)
$P_{11}$ stress component along the vertical reference coordinate $X_{2}$ at
$\theta H/(4l)=0.1125$.

(a) $\theta H/(4l)=0.015$

(b) $\theta H/(4l)=0.1125$

Figure 13: Uniform bending case. Nominal stress in the (a) pre-bifurcation and
(b) post-bifurcation regime as a function of the scale ratio $H/l\in[4,36]$.

The deformed configurations, obtained for all methods, are shown in Fig. 11
for an applied nominal strain $\theta H/(4l)=0.1125$ (i.e., in the post-
bifurcation regime). Since only the bottom half of the specimen is subjected
to compression, the pattern transformation is limited to the region situated
below the neutral axis, cf. Fig. 11a. The macroscopic deformation along with
deformed RVEs corresponding to three macroscopic integration points are shown
in Figs. 11b–11d. Here a clear distinction in the RVE deformations for the
different methods emerges. In particular, the bending effect of the applied
higher-order gradient $\partial F_{11}/\partial X_{2}=G_{211}$ is visible for
the second-order homogenization method, whereas the effect of the pattern
gradient $\partial v_{1}/\partial X_{2}$ is present in the micromorphic
scheme.

The nominal bending stress, defined as $6M/H^{2}$, where $M$ denotes the
moment acting on the vertical specimen edge, is plotted as a function of the
nominal strain $\theta H/(4l)$ for a scale ratio $12$ in Fig. 12a. Even though
both enriched homogenization schemes provide similar stress–strain results,
they systematically overestimate the stress in the post-bifurcation regime.
The micromorphic method exhibits a clear kink when moving from pre- to post-
bifurcation regime, whereas the first- and second-order methods show a more
gradual bifurcation transition. The corresponding homogenized stress component
of the first Piola–Kirchhoff stress tensor, $P_{11}(X_{2})$, is shown in Fig.
12b as a function of the vertical coordinate $X_{2}$ for scale ratio $12$ and
a nominal strain $\theta H/(4l)=0.1125$. Due to the pattern transformation,
the bottom half of the specimen displays a distinct plateau, whereas the top
part shows a linear profile. All homogenization methods considered provide
relatively accurate predictions.

The scale separation plots corresponding to strains $\theta H/(4l)=0.015$ and
$\theta H/(4l)=0.1125$, given in Fig. 13, reflect the scale ratio effect in
the pre- and post-bifurcation regime. In both cases, the micromorphic
prediction is closer to the reference solution compared to the second-order
method. The size effect is, however, not as strong as for the previously
discussed compression case, through which the first-order method (neglecting
size effects) provides the most accurate solution in the pre- and post-
bifurcation regimes. Note that although the ensemble averaged solution is not
varying much with the scale ratio, individual microstructural translations do
result in a large variance (cf. Fig. 13).

###  4.3 Compression of a Finite Specimen

The final case considers a finite rectangular domain spanning $n\times n$ unit
cells subjected to compression, cf. Fig. 14. The applied compression acts
along the vertical direction, and is induced by fixing the bottom horizontal
edge AB while displacing the top horizontal edge CD downwards by $u^{*}$.
Unlike the uniform compression case of Section 4.1, the two vertical edges AD
and BC are now free surfaces, allowing for compliant boundary layers and
auxetic effects. Mathematically, the applied boundary conditions are expressed
for the first-order homogenization method as:

$\vec{u}(\Gamma_{\mathrm{AB}})=\vec{0},\,\vec{u}(\Gamma_{\mathrm{CD}})=-u^{*}\,\vec{e}_{2},$
(33)

for the second-order method:

$\vec{u}(\Gamma_{\mathrm{AB}})=\vec{0},\,\vec{u}(\Gamma_{\mathrm{CD}})=-u^{*}\,\vec{e}_{2},\quad
F_{22}(\Gamma_{\mathrm{AB}})=F_{22}(\Gamma_{\mathrm{CD}})=1,$ (34)

and for the micromorphic method:

$\vec{v}_{0}(\Gamma_{\mathrm{AB}})=\vec{0},\,\vec{v}_{0}(\Gamma_{\mathrm{CD}})=-u^{*}\,\vec{e}_{2},\quad
v_{1}(\Gamma_{\mathrm{AB}})=v_{1}(\Gamma_{\mathrm{CD}})=0.$ (35)

In order to avoid global buckling, which has been addressed elsewhere [cf. van
Bree et al., 2020, Wu et al., 2023], and to focus on boundary layers and non-
local behaviour, only one half of the specimen is modelled, introducing a
vertical symmetry axis at the middle of the specimen along which the
horizontal displacements are restricted while the vertical displacements are
left free. The same discretization and integration rules as those in Section
4.2 are used.

Figure 14: Finite specimen compression case. Illustration of the simulation
$4l\times 4l$ domain subjected to a vertical compressive load. The two lateral
edges, AD and BC, are free surfaces.

|
---|---
|

(a) DNS

(b) first-order

(c) second-order

(d) micromorphic

Figure 15: Finite specimen compression case. Right half of a deformed
$10l\times 10l$ specimen obtained by (a) DNS subjected to $u^{*}/H=0.075$
compressive strain, corresponding to a zero microstructural translation
$\vec{\zeta}=\vec{0}$. The right half of the macroscopic deformed structure,
including deformed RVEs positioned at three macroscopic integration points,
corresponding to the (b) first-order, (c) second-order, and (d) micromorphic
homogenization method. The dotted blue line indicates the deformed contour of
the ensemble-averaged solution.

(a) $H/l=4$

(b) $H/l=10$

Figure 16: Finite specimen compression case. Stress–strain diagrams
corresponding to scale ratios (a) $4$, and (b) $10$, obtained by the first-
order, second-order, and micromorphic homogenization schemes, as compared to
the reference ensemble averaged solution.

(a) $u^{*}/H=0.01$

(b) $u^{*}/H=0.05$

Figure 17: Finite specimen compression case. Nominal compressive stress in the
(a) pre-bifurcation and (b) post-bifurcation regime as a function of the scale
ratio $H/l\in[4,36]$.

(a) $u^{*}/H=0.01$

(b) $u^{*}/H=0.05$

Figure 18: Finite specimen compression case. Tangent stiffness in the (a) pre-
bifurcation and (b) post-bifurcation regime as a function of the scale ratio
$H/l\in[4,36]$.

Deformed configurations of a half specimen (due to symmetry), corresponding to
scale ratio $H/l=10$ and subjected to $7.5\%$ of overall compressive strain
(i.e., post-bifurcation), are shown in Fig. 15. The auxetic effect resulting
from the pattern transformation is clearly visible and captured by all
methods. The two enriched homogenization schemes (Figs. 15c and 15d) reproduce
the auxetic behaviour more accurately than the first-order method (Fig. 15b),
judging from the outer contours of the ensemble-averaged solution denoted by
the blue dotted line. A significant drawback of the first-order method is the
presence of a kink in the displacements at the centres of the two lateral
edges. The overall shape of the compliant boundary and the deformed RVEs
corresponding to the top and bottom Gauss integration points are most
accurately captured by the micromorphic method. The investigated integration
points close to the middle height of the specimen are similar for both
enriched methods, whereas the RVEs corresponding to the first-order method
deform severely close to the kinks.

The stress–strain diagrams, corresponding to $4l\times 4l$ and $10l\times 10l$
specimens, are provided in Fig. 16. Similarly to the uniform compression case
of Section 4.1, the bifurcation strain is not captured accurately by the
second-order homogenization method, whereas it is captured adequately by the
micromorphic scheme. In the case of the $4l\times 4l$ specimen, the post-
bifurcation stiffness (i.e., the slope of the stress–strain diagram), although
comparable for both enriched homogenization schemes, significantly
overestimates the expected value, which even reveals softening (a negative
slope). The observed softening results from the presence of the two vertical
compliant boundaries, which for such a small scale ratio have a larger
influence on the overall response than the stiff boundary layers induced by
the two horizontal edges. Although this effect rapidly vanishes with
increasing scale ratio, neither of the two enriched homogenization schemes is
able to capture this phenomenon properly. In the case of the $10l\times 10l$
specimen, the post-bifurcation stiffness is reproduced accurately by both the
second-order as well as the micromorphic method. The first-order method gives,
however, the most accurate prediction of the plateau behaviour despite its
large error in the predicted bifurcation strain, cf. also Fig. 18b.

Graphs of the applied nominal stress are shown in Fig. 17 as a function of
scale ratio for two different nominal strains $u^{*}/H=0.01$ (pre-bifurcation)
and $u^{*}/H=0.05$ (post-bifurcation). From the reported results it may be
concluded that the micromorphic homogenization method provides the most
accurate or balanced results, although an overestimation of the stress is
still observed in Fig. 17b. In the same figure, the second-order method
provides overall more accurate result, which is partially due to the choice of
the nominal strain at which the nominal stress is evaluated. Pre-bifurcation,
a significant overestimation similar to the uniform compression case is
observed for second-order scheme.

(a) DNS

(b) first-order

(c) second-order

(d) micromorphic

(e) DNS

(f) first-order

(g) second-order

(h) micromorphic

Figure 19: Finite specimen compression case. Resulting field plots for
$u^{*}/H=0.075$ and $H/l=10$, for the right half of the domain. The top row
corresponds to the horizontal displacement component $\vec{u}\cdot\vec{e}_{1}$
[mm] for (a) DNS, (b) first-order, (c) second-order, and (d) micromorphic
simulation. The bottom row corresponds to the $F_{22}-1$ [-] component of the
macroscopic deformation for (e) DNS, (f) first-order, (g) second-order, and
(g) micromorphic simulation. All fields are plotted over the reference
configuration. DNS ensemble averages are computed for $21$ shifts
$\vec{\zeta}$ in each direction, showing high fluctuations due to low
sampling.

Fig. 18 shows the tangent stiffness at nominal strains $u^{*}/H=0.01$ and
$u^{*}/H=0.05$ for all adopted scale ratios. In particular, Fig. 18b confirms
that both enriched homogenization methods fail to capture the softening
behaviour observed for the lowest scale ratios considered (cf. Fig. 16a), and
that tangent stiffnesses are significantly overestimated. The micromorphic
method is, nevertheless, closest to the reference solution within the
evaluated range. Fig. 18a is similar to Fig. 17a, hence the same discussion
and conclusions apply.

Corresponding field plots of the effective horizontal displacements are shown
in Figs. 19a–19d, whereas field plots of $F_{22}-1$ component in Figs.
19e–19h. In all cases, the state corresponds to applied nominal strain
$u^{*}/H=0.075$ in the post-bifurcation regime. The field plots clearly show
that the DNS results are under-sampled and exhibit strong fluctuations, in
spite of using $21$ shifts in both horizontal and vertical directions.
Nevertheless, a clear auxetic effect can be observed in Fig. 19a, whereas the
average vertical strain component is almost uninformative, although boundary
layers of size $3l$ and lower values on average close to the mid-height of the
specimen can be observed. For first-order computational homogenization, a
clear localization due to separation of scales (i.e., locality) assumption can
be observed, in contrast to the two enriched homogenization schemes. Here, the
second-order approach shows a slightly more extensive and less localized
region of negative displacements as compared to the micromorphic scheme. For
the $F_{22}-1$ component, second-order homogenization shows a significantly
more localized boundary layer, similarly to the conclusions of Section 4.1 in
Fig. 5b, which is caused by the prescribed value of the deformation gradient
$F_{22}=1$ at the two horizontal edges.

###  4.4 Concluding Remarks

Table 1: General performance overview of tested computational homogenization methods applied to elastomeric mechanical metamaterials. property/method | first-order | second-order | micromorphic
---|---|---|---
pre-bifurcation (no size effect) | ✓ | ✗ | ✓
post-bifurcation (size effect) | ✗ | ✓ | ✓
shift of bifurcation point | ✗ | ✗ | ✓
no prior knowledge | ✓ | ✓ | ✗

The overall performance of all considered homogenization methods applied to
elastomeric mechanical metamaterials, as obtained from the numerical examples,
is summarized in Tab. 1. In general, it may be concluded that the micromorphic
computational homogenization scheme is able to capture most of the present
size effects, relying on prior information on the patterning mode
$\vec{\varphi}_{1}$. The second-order homogenization scheme is less accurate
(in particular, it exhibits too strong size effects, especially in the pre-
bifurcation regime), but on the other hand is fully general, requiring no
prior knowledge.

##  5 Summary and Conclusions

In this paper, a thorough comparison of two enhanced computational
homogenization schemes applied to pattern-transforming elastomeric mechanical
metamaterials has been performed. The first method considered is second-order
computational homogenization scheme, which accounts for strain gradients at
the macro-scale, and thus includes non-local effects induced by pattern
transformations. The second method, a micromorphic computational
homogenization scheme, introduces the magnitude of the emerging pattern as an
additional macroscopic field, which captures kinematic coupling between
individual cells. Whereas the second-order computational homogenization
requires solution of a higher-order continuum problem at the macro-scale, the
micromorphic scheme adds one coupled scalar partial differential equation next
to the classical macroscopic balance equation. Both methodologies have been
compared against a reference, ensemble-averaged solution, which is obtained
through Direct Numerical Simulations (DNS) by including a series of translated
microstructures. The first-order homogenization method is included as well to
reveal the limitations of this current standard in homogenization solutions.
Three distinct loading cases were evaluated to demonstrate the accuracy of the
enriched methods: uniform compression, bending, and compression of finite
specimens. The most important conclusions can be summarized as follows:

  1. 1.

Both enhanced homogenization schemes provide a reasonable estimate of the
effective behaviour exhibited by the elastomeric mechanical metamaterials
adopted in this study.

  2. 2.

The micromorphic homogenization scheme provides more accurate results in terms
of the kinematics as well as the stress quantities as compared to second-order
computational homogenization scheme.

  3. 3.

The deformed Representative Volume Elements (RVEs) obtained by both methods
correspond well with the DNS results due to the presence of the macroscopic
strain gradient or a micromorphic field.

  4. 4.

Due to very small size effects in the case of bending, both enhanced
homogenization methods overestimate the effective stress in the post-
bifurcation regime. In the worst case, i.e., for scale ratio $4$, the second-
order and micromorphic method are approximately $38\%$ and $15\%$ above the
corresponding reference whereas the first-order homogenization method leads to
about $5\%$ of error.

  5. 5.

The auxetic effect, present in the compression of the finite specimen example,
is captured accurately by both enhanced homogenization schemes. The shape of
the soft boundary, however, is more accurate for the micromorphic method.

  6. 6.

The compression of finite specimens proves to be a challenging task for all
computational homogenization schemes, especially in terms of the nominal
stress–strain response. Whereas the micromorphic method is still able to
capture the bifurcation point with reasonable accuracy, the second-order
computational homogenization introduces an error of the order of $50\%$ for
$4l\times 4l$ specimen. The post-bifurcation stiffness for the same specimen
is significantly overestimated, but converges rapidly towards the correct
value with increasing scale ratio.

From the presented results, it may be concluded that the micromorphic
computational homogenization method proves to be the most suitable numerical
tool for homogenization of elastomeric mechanical metamaterials experiencing
pattern-transformations, although second-order computational homogenization
holds promise as well. Its main drawback is the overly strong non-local effect
in the absence of patterning, often leading to overestimation of effective
properties.

## Acknowledgments

The research leading to these results has received funding from the European
Research Council under the European Union’s Seventh Framework Programme
(FP7/2007-2013)/ERC grant agreement № [339392]. The authors would furthermore
like to thank Maqsood M. Ameen for numerous fruitful discussions in the early
stages of the project, related to micromorphic computational homogenization.

## References

  * Ameen et al. [2018a] Ameen, M.M., Peerlings, R.H.J., Geers, M.G.D., 2018a.  A quantitative assessment of the scale separation limits of classical and higher-order asymptotic homogenization.  European Journal of Mechanics /A Solids 71, 89–100.
  * Ameen et al. [2018b] Ameen, M.M., Rokoš, O., Peerlings, R.H.J., Geers, M.G.D., 2018b.  Size effects in nonlinear periodic materials exhibiting reversible pattern transformations.  Mechanics of Materials 124, 55–70.
  * Anthoine [2010] Anthoine, A., 2010.  Second-order homogenisation of functionally graded materials.  International Journal of Solids and Structures 47, 1477–1489.
  * Bertoldi et al. [2008] Bertoldi, K., Boyce, M.C., Deschanel, S., Prange, S.M., Mullin, T., 2008.  Mechanics of deformation-triggered pattern transformations and superelastic behavior in periodic elastomeric structures.  Journal of the Mechanics and Physics of Solids 56, 2642–2668.
  * Blanco et al. [2016] Blanco, P.J., Sánchez, P.J., de Souza Neto, E., Feijóo, R., 2016.  The method of multiscale virtual power for the derivation of a second order mechanical model.  Mechanics of Materials 99, 53–67.
  * Cosserat and Cosserat [1909] Cosserat, E., Cosserat, F., 1909.  Théorie des corps déformables.  A. Hermann et fils.
  * Coulais [2016] Coulais, C., 2016.  Periodic cellular materials with nonlinear elastic homogenized stress-strain response at small strains.  International Journal of Solids and Structures 98, 226–238.
  * Drugan and Willis [1996] Drugan, W.J., Willis, J.R., 1996.  A micromechanics-based nonlocal constitutive equation and estimates of representative volume element size for elastic composites.  Journal of the Mechanics and Physics of Solids 44, 497–524.
  * Dykstra et al. [2023] Dykstra, D.M., Lenting, C., Masurier, A., Coulais, C., 2023.  Buckling metamaterials for extreme vibration damping.  Advanced Materials n/a, 2301747.
  * Florijn et al. [2014] Florijn, B., Coulais, C., van Hecke, M., 2014.  Programmable mechanical metamaterials.  Physical Review Letters 113, 175503.
  * Forest and Trinh [2010] Forest, S., Trinh, D.K., 2010.  Generalized continua and non-homogeneous boundary conditions in homogenization methods.  Zeitschrift für Angewandte Mathematik und Mechanik 91, 90–109.
  * Gao et al. [2018] Gao, C., Slesarenko, V., Boyce, M.C., Rudykh, S., Li, Y., 2018.  Instability-induced pattern transformation in soft metamaterial with hexagonal networks for tunable wave propagation.  Scientific Reports , 8:11834.
  * Geers et al. [2010] Geers, M.G.D., Kouznetsova, V.G., Brekelmans, W.A.M., 2010.  Multi-scale computational homogenization: trends and challenges.  Journal of Computational and Applied Mathematics 234, 2175–2182.
  * Geuzaine and Remacle [2009] Geuzaine, C., Remacle, J.F., 2009.  Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities.  International journal for numerical methods in engineering 79, 1309–1331.
  * Hu et al. [2013] Hu, J., He, Y., Lei, J., Lui, Z., 2013.  Novel mechanical behavior of periodic structure with the pattern transformation.  Theoretical and Applied Mechanics Letters 3.
  * Kouznetsova et al. [2002] Kouznetsova, V.G., Geers, M.G.D., Brekelmans, W.A.M., 2002.  Multi‐scale constitutive modelling of heterogeneous materials with a gradient-enhanced computational homogenization scheme.  International Journal for Numerical Methods in Engineering 54, 1235–1260.
  * Kouznetsova et al. [2004] Kouznetsova, V.G., Geers, M.G.D., Brekelmans, W.A.M., 2004.  Multi-scale second-order computational homogenization of multi-phase materials: a nested finite element solution strategy.  Computer methods in applied Mechanics and Engineering 193, 5525–5550.
  * Krushynska et al. [2023] Krushynska, A.O., Torrent, D., Aragón, A.M., Ardito, T., Bilal, O.R., Bonello, B., Bosia, F., Chen, Y., Christensen, J., Colombi, A., Cummer, S.A., Djafari-Rouhani, B., Fraternali, F., Galich, P.I., Garcia, P.D., Groby, J.P., Guenneau, S., Haberman, M.R., Hussein, M.I., Janbaz, S., Jiménez, N., Khelif, A., Laude, V., Mirzaali, M.J., Packo, P., Palermo, A., Pennec, Y., Picó, R., López, M.R., Rudykh, S., Serra-Garcia, M., Torres, C.M.S., A., S.T., Tournat, V., Wright, O.B., 2023.  Emerging topics in nanophononics and elastic, acoustic, and mechanical metamaterials: an overview.  Nanophotonics 12, 659–686.
  * Lesičar et al. [2014] Lesičar, T., Tonković, Z., Sorić, J., 2014.  A second-order two-scale homogenization procedure using $C^{1}$ macrolevel discretization.  Computational Mechanics 54, 425–441.
  * Lopes and Pires [2022a] Lopes, I.A.R., Pires, F.M.A., 2022a.  An assessment of multi-scale models based on second-order computational homogenisation.  Computers & Structures 259, 106679.
  * Lopes and Pires [2022b] Lopes, I.A.R., Pires, F.M.A., 2022b.  Unlocking the potential of second-order computational homogenisation: An overview of distinct formulations and a guide for their implementation.  Archives of Computational Methods in Engineering 29, 1339–1393.
  * Luscher et al. [2010] Luscher, D.J., McDowell, D.L., Bronkhorst, C.A., 2010.  A second gradient theoretical framework for hierarchical multiscale modeling of materials.  International Journal of Plasticity 26, 1248–1275.  Special Issue In Honor of Lallit Anand.
  * Maraghechi et al. [2024] Maraghechi, S., Rokoš, O., Hoefnagels, J.P.M., Peerlings, R.H.J., Geers, M.G.D., 2024.  Harvesting deformation modes for micromorphic homogenization from experiments on mechanical metamaterials.  Submitted , 1–24.
  * Miehe and Bayreuther [2007] Miehe, C., Bayreuther, C.G., 2007.  On multiscale FE analyses of heterogeneous structures: from homogenization to multigrid solvers.  International Journal for Numerical Methods in Engineering 71, 1135–1180.
  * Mindlin [1965] Mindlin, R.D., 1965.  Second gradient of strain and surface-tension in linear elasticity.  International Journal of Solids and Structures 1, 417–438.
  * Mindlin and Eshel [1968] Mindlin, R.D., Eshel, N.N., 1968.  On first strain-gradient theories in linear elasticity.  International Journal of Solids and Structures 4, 109–124.
  * Mullin et al. [2007] Mullin, T., Deschanel, S., Bertoldi, K., Boyce, M.C., 2007.  Pattern transformation triggered by deformation.  Physical Review Letters 99.
  * Nikman and Akbarzadeh [2018] Nikman, H., Akbarzadeh, A.H., 2018.  In-plane and out-of-plane buckling of architected cellular plates: Numerical and experimental study.  Composite Structures 206, 739–749.
  * Noels and Nguyen [2014] Noels, L., Nguyen, V.D., 2014.  Computational homogenization of cellular materials.  International Journal of Solids and Structures 51, 2183–2203.
  * Prawoto [2012] Prawoto, Y., 2012.  Seeing auxetic materials from the mechanics point of view: A structural review on the negative Poisson’s ratio.  Computational Materials Science 58, 140–153.
  * Rokoš et al. [2019] Rokoš, O., Ameen, M.M., Peerlings, R.H.J., Geers, M.G.D., 2019.  Micromorphic computational homogenization for mechanical metamaterials with patterning fluctuation fields.  Journal of the Mechanics and Physics of Solids 123, 119–137.
  * Rokoš et al. [2020a] Rokoš, O., Ameen, M.M., Peerlings, R.H.J., Geers, M.G.D., 2020a.  Extended micromorphic computational homogenization for mechanical metamaterials exhibiting multiple geometric pattern transformations.  Extreme Mechanics Letters 37, 100708.
  * Rokoš et al. [2020b] Rokoš, O., Zeman, J., Doškář, M., Krysl, P., 2020b.  Reduced integration schemes in micromorphic computational homogenization of elastomeric mechanical metamaterials.  Under review , 1–19.
  * Smyshlyaev and Cherednichenko [2000] Smyshlyaev, V.P., Cherednichenko, K.D., 2000.  On rigorous derivation of strain gradient effects in the overall behaviour of periodic heterogeneous media.  Journal of the Mechanics and Physics of Solids 48, 1325–1357.
  * van Bree et al. [2020] van Bree, S.E.H.M., Rokoš, O., Peerlings, R.H.J., Doškář, M., Geers, M.G.D., 2020.  A Newton solver for micromorphic computational homogenization enabling multiscale buckling analysis of pattern-transforming metamaterials.  Computer Methods in Applied Mechanics and Engineering 372, 113333.
  * Wu et al. [2015] Wu, G., Cho, Y., Choi, I., Ge, D., Li, J., Nam Han, H., Lubensky, T., Yang, S., 2015.  Directing the deformation paths of soft metamaterials with prescribed asymmetric units.  Advanced Materials 27, 2747–2752.
  * Wu et al. [2023] Wu, L., Mustafa, M., Segurado, J., Noels, L., 2023.  Second-order computational homogenisation enhanced with non-uniform body forces for non-linear cellular materials and metamaterials.  Computer Methods in Applied Mechanics and Engineering 407, 115931.
  * Xia and Hutchinson [1996] Xia, Z.C., Hutchinson, J.W., 1996.  Crack tip fields in strain gradient plasticity.  Journal of the Mechanics and Physics of Solids 44, 1621–1648.
  * Yvonnet et al. [2020] Yvonnet, J., Auffray, N., Monchiet, V., 2020.  Computational second-order homogenization of materials with effective anisotropic strain-gradient behavior.  International Journal of Solids and Structures 191-192, 434–448.