# Harvesting Deformation Modes for Micromorphic Homogenization from
Experiments on Mechanical Metamaterials111The post-print version of this
article is published in _Int. J. Solids Struct._ ,
10.1016/j.ijsolstr.2024.112916. This manuscript version is made available
under the CC-BY 4.0 license.

S. Maraghechi<EMAIL_ADDRESS>O. Rokoš<EMAIL_ADDRESS>R.H.J. Peerlings
<EMAIL_ADDRESS>M.G.D. Geers<EMAIL_ADDRESS>J.P.M. Hoefnagels
<EMAIL_ADDRESS>Mechanics of Materials, Department of Mechanical
Engineering, Eindhoven University of Technology, P.O. Box 513, 5600 MB
Eindhoven, The Netherlands

###### Abstract

A micromorphic computational homogenization framework has recently been
developed to deal with materials showing long-range correlated interactions,
i.e. displaying patterning modes. Typical examples of such materials are
elastomeric mechanical metamaterials, in which patterning emerges from local
buckling of the underlying microstructure. Because pattern transformations
significantly influence the resulting effective behaviour, it is vital to
distinguish them from the overall deformation. To this end, the following
kinematic decomposition into three parts was introduced in the micromorphic
scheme:  (i) a smooth mean displacement field, corresponding to the slowly
varying deformation at the macro-scale,  (ii) a long-range correlated
fluctuation field, related to the buckling pattern at the meso-scale, and
(iii) the remaining uncorrelated local microfluctuation field at the micro-
scale.  The micromorphic framework has proven to be capable of predicting
relevant mechanical behaviour, including size effects and spatial as well as
temporal mixing of patterns in elastomeric metamaterials, making it a powerful
tool to design metamaterials for engineering applications. The long-range
correlated fluctuation fields need to be, however, provided _a priori_ as
input parameters. The main goal of this study is experimental identification
of the decomposed kinematics in cellular metamaterials based on the three-part
ansatz. To this end, a full-field micromorphic Integrated Digital Image
Correlation (IDIC) technique has been developed. The methodology is formulated
for finite-size cellular elastomeric metamaterial specimens deformed in  (i)
virtually generated images and  (ii) experimental images attained during in-
situ compression of specimens with millimetre sized microstructure using
optical microscopy.  The proposed IDIC method identifies the different
kinematic fields, both before and after the microstructural buckling, and
without any prior knowledge determines correctly the relevant patterning modes
required by the homogenization scheme. It is further argued that patterning
modes are independent of the unit cell size, the hole diameter to cell size
ratio, as well as local material properties, allowing for modelling and design
of (finite- and infinite-size) metamaterials and specimens with graded
microstructures in terms of geometry and/or material properties. It is shown
that the proposed methodology is also applicable to cellular metamaterials and
structures with different microstructural designs.

###### keywords:

Elastomeric Metamaterials , Homogenization , Integrated Digital Image
Correlation , Kinematic Decomposition , In-situ Testing

††journal: Int. J. Solids Struct.

##  1 Introduction

Mechanical metamaterials have recently gained substantial attention in the
literature due to their potential for various engineering applications such as
strain tunable photonic crystals [Krishnan and Johnson, 2009], tunable
acoustic properties [Zhao et al., 2015, Bilal et al., 2017, Wormser et al.,
2017, Lin et al., 2023], soft robotics [Yang et al., 2015, Mirzaali et al.,
2018b, Yu et al., 2018, Shen et al., 2023], programmable mechanical
metamaterials [Florijn et al., 2014, Bertoldi et al., 2017, Kolken and
Zadpoor, 2017, Li et al., 2023], biomedical prosthetics [Yavari et al., 2015,
Zadpoor, 2017], or functionally graded metamaterials [Francesconi et al.,
2019, Farzaneh et al., 2022]. Various aspects of such metamaterials with 2D
[Bertoldi et al., 2008] and 3D [Javid et al., 2016, Hedayati et al., 2019]
microstructures have been investigated. A number of studies focused on the
effect of the hole geometry and stacking [Bertoldi et al., 2010, Overvelde et
al., 2012, Overvelde and Bertoldi, 2014, Hu et al., 2013]. Several studies
reported in the literature made use of full-field displacement measurement
techniques to investigate cellular metamaterials, among which Slann et al.
[2015] used local Digital Image Correlation (DIC) to validate their numerical
study on auxetic metamaterials, Xu and Pasini [2016] used DIC in the design of
thermally expanding tunable 3D metamaterials, Easey et al. [2019] studied
different designs of auxetic dome-shaped structures using stereo DIC, and Tang
et al. [2015] compared DIC and finite element analysis results on highly
stretchable and reconfigurable metamaterials. Apart from displacement and
strain assessment, Mullin et al. [2007] made a qualitative comparison between
numerically and experimentally assessed stresses in cellular elastomers.
Besides the geometrical and morphological effects, also the type of material
is important. Mirzaali et al. [2018a] used 3D printing to manufacture and
study multi-material mechanical metamaterials, Wu et al. [2014] investigated
different pattern transformations as a result of swelling in cellular hydrogel
membranes, whereas Hu et al. [2018] investigated the effect of friction on
buckling pattern formation in cellular structures. A large number of works
also focused on experimental investigation of pantographic metamaterials,
e.g., [Auger et al., 2020, Hild et al., 2020, Barchiesi et al., 2021].

The common denominator to the listed applications are exotic properties of
metamaterials, typically emerging from a coordinated motion of their
microstructure, i.e., from the so-called patterning. The patterning usually
occurs upon reaching a critical level of compressive load, resulting in a
local buckling of the underlying microstructure and a significant change in
the overall mechanical properties. Fig. 1a shows an example of a cellular
elastomeric metamaterial specimen with a regular grid of circular holes, where
the local buckling results in a pattern of alternating elliptical holes (Fig.
1b). Here it can be observed that the patterning is locally constrained by
external boundaries, and depends on the size of the holes relative to the size
of the entire specimen (i.e., on the scale ratio), studied computationally in
detail by Ameen et al. [2018]. Since such irregularities in the pattern
directly influence the overall behaviour of the specimen—even on a larger
scale—, it is important from the engineering perspective to be able to
identify these distributions, independently of the global deformation, for a
proper assessment of and design for certain mechanical properties.

Apart from classical homogenization approaches mostly limited to linear
elastic behaviour, such as, e.g., [Toupin, 1962, Mindlin, 1964, 1965], an
effective tool for predicting mechanical behaviour of materials with complex
microstructures on an engineering scale and operating in a non-linear regime
is first-order computational homogenization [Kouznetsova et al., 2001, Miehe
and Koch, 2002]. To determine effective constitutive behaviour, computational
homogenization solves a microscopic boundary value problem associated with
each macroscopic integration point. Obtained results are used to compute
effective stresses, which are equilibrated at the macroscale. Because
separation of scales is assumed, such a scheme is able to properly model only
macroscopically local continuum. For better approximation, extensions to
second-order homogenization schemes were proposed [Kouznetsova et al., 2004,
dell’Isola et al., 2019, Giorgio et al., 2023], which use gradient of the
deformation gradient to account for non-local effects. Non-locality can also
be addressed using generalized micromorphic theories, e.g., [Forest, 2009,
Hütter, 2017, Biswas and Poh, 2017, Rokoš et al., 2019, 2020], in which
additional kinematic fields are introduced to account for the underlying
microfluctuations. For a recent overview and further approaches see, e.g.,
Findeisen et al. [2020].

(a) Undeformed microstructure

(b) Deformed microstructure

Figure 1: Elastomeric metamaterial specimen LABEL:sub@fig5:undef before and
LABEL:sub@fig5:def after the onset of microstructural buckling. The emerging
antisymmetric pattern is clearly visible, which is a result of a compression
in the horizontal direction.

Of particular interest in this contribution is the micromorphic computational
homogenization scheme based on recent theoretical and numerical studies [Rokoš
et al., 2019, 2020, van Bree et al., 2020], which relies on a general
decomposition of the displacement field according to the following ansatz

$\mathbf{u}(\mathbf{x})=\mathbf{v}_{0}(\mathbf{x})+\sum_{i=1}^{n}v_{i}(\mathbf{x})\boldsymbol{\varphi}_{i}(\mathbf{x})+\mathbf{w}(\mathbf{x}).$
(1)

Here, the total displacement $\mathbf{u}(\mathbf{x})$, which is a vector field
expressed as a function of the material coordinate in the undeformed
configuration $\mathbf{x}$, is divided into three parts:  (i) the mean smooth
displacement field, $\mathbf{v}_{0}(\mathbf{x})$, corresponding to the slow
variations at the macro-scale;  (ii) long-range correlated fluctuations at the
meso-scale, represented by the vector fields
$\boldsymbol{\varphi}_{i}(\mathbf{x})$, $i=1,\dots,n$, scaled by their spatial
magnitudes $v_{i}(\mathbf{x})$; and  (iii) a remaining microfluctuation field,
$\mathbf{w}(\mathbf{x})$, representing the uncorrelated fluctuations at the
micro-scale.  By means of the patterning modes, $\boldsymbol{\varphi}_{i}$,
the micromorhpic scheme is capable of capturing spatial non-local effects in
the microstructure, present typically because of the lack of separation of
scales and emerging due to effects such as boundary layers or temporal and
spatial mixing of modes. Non-local effects have been studied numerically as
well as experimentally within the context of metamaterial structures in, e.g.,
[Ameen et al., 2018, Maraghechi et al., 2020]. Using a variational formulation
towards a numerical approximation, [cf. van Bree et al., 2020], an effective
potential energy is minimized with respect to the kinematic decomposition of
Eq. (1), resulting in a micromorphic type of effective continuum, i.e., not of
the standard Cauchy–Born type. Although the standard balance of momentum
equation $\boldsymbol{\nabla}_{0}\cdot\mathbf{P}=\boldsymbol{0}$ governs the
evolution of the $\mathbf{v}_{0}(\mathbf{x})$ field, where $\mathbf{P}$ is the
homogenized first Piola–Kirchhoff stress tensor energetically conjugated with
the deformation gradient tensor
$\mathbf{F}=\mathbf{I}+(\boldsymbol{\nabla}_{0}\mathbf{v}_{0}(\mathbf{x}))^{\mathsf{T}}$,
an additional coupled micromorphic balance equation emerges. Herein,
$\boldsymbol{\nabla}_{0}$ denotes the gradient operator with respect to the
reference configuration. This equation governs the evolution of
$v_{i}(\mathbf{x})$, which is in the form of
$\boldsymbol{\nabla}_{0}\cdot\boldsymbol{\Lambda}-\Pi_{i}=0$, where
$\boldsymbol{\Lambda}$ is energetically conjugated with
$\boldsymbol{\nabla}_{0}v_{i}(\mathbf{x})$ while $\Pi_{i}$ is conjugated with
$v_{i}(\mathbf{x})$. Since $v_{i}(\mathbf{x})$ corresponds to the magnitude of
the patterning fluctuation field $\boldsymbol{\varphi}_{i}(\mathbf{x})$, its
kinematic interpretation is straightforward and essential boundary conditions
can be identified as $v_{i}(\mathbf{x})=0$ along any fixed boundaries. On the
contrary, $\boldsymbol{\Lambda}\cdot\boldsymbol{N}=0$ on the free boundaries
may be adopted (where $\boldsymbol{N}$ denotes the outer normal). We refer the
interested reader for further details and derivations to [Rokoš et al., 2020,
Section 3.3]. Note that unlike second-order computational homogenization
[e.g., Kouznetsova et al., 2004], which uses continuity of the gradient of the
deformation gradient tensor field $\mathbf{F}$ to account for non-locality and
size effects, micromorphic computational homogenization introduces an
additional independent kinematical field with a governing equation
complementing the standard balance of linear momentum. The homogenization
scheme can thus facilitate optimized designs of engineering systems, including
deeper understanding of their behaviour and development of new metamaterials.
Although the kinematic fields $\mathbf{v}_{0}$, $v_{i}$ and $\mathbf{w}$ are
an outcome of a numerical analysis, the microfluctuation patterning fields
$\boldsymbol{\varphi}_{i}$ are input parameters characteristic to the geometry
of the underlying microstructure that need to be provided _a priori_ , in
analogy to material constants of an underlying constitutive model. The
patterning modes can be estimated through an educated guess [cf. Rokoš et al.,
2019], or a numerical Bloch-type analysis [cf. Bertoldi et al., 2008]. It
would, however, be much better to measure the patterning modes directly from
experiments, as this would include the real material behaviour, all micro-
fabrication flaws, and errors in the correction of spatial distortions in the
geometrical images.

Therefore, the goal of this contribution is to develop a novel full-field
Integrated Digital Image Correlation (IDIC) framework for quantitative
experimental identification of the kinematic fields decomposed according to
Eq. (1), including the meso-scale long-range correlated fluctuation fields
$\boldsymbol{\varphi}_{i}$ for the micromorphic homogenization scheme. To this
end, high-resolution images captured at different stages of deformation during
an experimental test of a finite-sized specimen are analysed with this
micromorphic IDIC framework. IDIC is a powerful full-field displacement
measurement technique based on regularization of the kinematics according to
specific knowledge about the problem at hand, resulting in identification of
the optimal input model parameters minimizing the difference between
experimentally measured data and model predictions [Roux and Hild, 2006,
Neggers et al., 2015, Blaysat et al., 2015]. Knowledge of the nature of the
problem at hand can be integrated in the DIC scheme in the form of a
parametric analytical description, which is utilized here in the form of the
kinematic fields of interest. Specifically, the total displacement field is
regularized based on the ansatz of Eq. (1) by parametrizing all relevant
kinematic components, providing upon correlation the smooth mean field,
$\mathbf{v}_{0}$, the long-range correlated fluctuation modes carried by the
microstructure, $\boldsymbol{\varphi}_{i}$, as well as their spatial
distribution magnitudes, $v_{i}$.

The proposed methodology is tested on virtually generated and deformed images,
as well as optical microscopy images taken during an in-situ compression test
on a finite-size cellular elastomeric metamaterial with rectangular stacking
of circular holes, cf. Fig. 1. A spectral density analysis is shown to be an
effective tool for estimation, parametrization, and initialization of the
fluctuation mode for a specific metamaterial. It is demonstrated that the
micromorphic IDIC scheme leads to a proper decomposition of the kinematic
fields of interest, both before and after the emergence of the fluctuation
pattern. Considering geometrical arguments, the patterning mode is shown to be
independent of the considered periodic cell size, although effective
mechanical behaviour of such a material is still size dependent. Moreover, it
is furthermore shown that the patterning mode is only weakly dependent on the
hole diameter versus unit cell size ratio, which enables the design of
metamaterials with graded microstructures.

Although a square stacking of holes is used throughout to guide and validate
our developments, the proposed methodology is applicable for a general class
of cellular metamaterials including hexagonal stacking of holes and chiral or
re-entrant microstructures, as well as graded microstructures in terms of
geometry (stacking of holes) and material properties (e.g., stiffness). The
reason for focusing on the square stacking of holes herein is twofold. First,
our detailed analysis requires testing inside an in-situ micro-tensile stage
underneath an optical microscope, which means that the total specimen size
needs to fit within a small, rectangular space. Therefore, square stacking of
holes is preferred as it results in a thin, constant-thickness boundary layer
along the edge, in contrast to other microstructures that may result in a
boundary layer with a varying thickness along the edge. Second, a rectangular
specimen shape with square stacking of holes is easiest to manufacture with
high quality and apply a constant load along the boundary edge. Finally, the
remaining boundary layer effect of square-holed microstructure has been
studied in detail (for precisely the same rectangular specimen shape)
elsewhere [Ameen et al., 2018, Maraghechi et al., 2020], so this configuration
is known to work well. But, as indicated earlier, it will be demonstrated that
the herein proposed methodology is generally applicable also to various
(graded) geometries.

##  2 Micromorphic Integrated Digital Image Correlation

###  2.1 Integrated Digital Image Correlation

Integrated digital image correlation is a full-field displacement measurement
method based on the minimization of the residual field, which is the
difference of the reference image, $f$, and the deformed image, $g$, probed
at, respectively, reference coordinates and the deformed positions
corresponding to those coordinates based on the displacement field. The
problem is solved over the whole region of interest, where the residual field
is defined as

$r(\mathbf{x},\underline{\mathbb{a}})=f\left(\mathbf{x}\right)-g\left(\boldsymbol{\phi}_{M}(\mathbf{x},\underline{\mathbb{a}})\right),$
(2)

with
$\boldsymbol{\phi}_{M}(\mathbf{x},\underline{\mathbb{a}})=\mathbf{x}+\mathbf{u}\left(\mathbf{x},\underline{\mathbb{a}}\right)$
being the mapping function corresponding to the mechanical displacement field,
$\mathbf{u}$. The displacement field, in turn, is based on the regularization
according to a model with specific knowledge of the kinematics at hand,
thereby reducing the number of degrees of freedom (dof) contained in
$\underline{\mathbb{a}}$, and improving the robustness of the solution
[Ruybalid et al., 2016]. This model can be a numerical model in which the dofs
would be the material parameters or an analytical model describing the
deformation field. The optimal values of the degrees of freedom,
$\underline{\mathbb{a}}^{opt}$, are found as

$\underline{\mathbb{a}}^{opt}=\underset{\underline{\mathbb{a}}}{\arg\min}\left(\frac{1}{2}\int
r^{2}\left(\mathbf{x},\underline{\mathbb{a}}\right)d\mathbf{x}\right).$ (3)

This minimization problem is then solved by an iterative scheme such as
Gauss–Newton algorithm [Neggers et al., 2016].

###  2.2 Regularization of the Kinematics Based on the Micromorphic
Kinematical Ansatz

The kinematical ansatz introduced in Eq. (1) is used to experimentally obtain
the decomposed kinematic fields of cellular metamaterials. The local
microfluctuations are the uncorrelated part of the kinematics and by
definition unknown, therefore $\mathbf{w}(\mathbf{x})$ is not included in the
micromorphic IDIC parametrization for cellular metamaterials (in fact,
$\mathbf{w}(\mathbf{x})$ is directly related to the residual field
$r(\mathbf{x},\underline{\mathbb{a}})$ in Eq. (2) via the image gradient).
Accordingly, the displacement field for the IDIC problem is defined as

$\mathbf{u}(\mathbf{x},\underline{\mathbb{a}})=\mathbf{v}_{0}(\mathbf{x},\underline{\mathbb{a}}_{v_{0}})+\sum_{i=1}^{n}v_{i}(\mathbf{x},\underline{\mathbb{a}}_{v_{i}})\boldsymbol{\varphi}_{i}(\mathbf{x},\underline{\mathbb{a}}_{\varphi_{i}}),$
(4)

where $\underline{\mathbb{a}}=\left[\underline{\mathbb{a}}_{v_{0}}\
\underline{\mathbb{a}}_{v_{i}}\
\underline{\mathbb{a}}_{\varphi_{i}}\right]^{\mathsf{T}}$ is the column of
degrees of freedom, with $\underline{\mathbb{a}}_{v_{0}}$,
$\underline{\mathbb{a}}_{v_{i}}$ and $\underline{\mathbb{a}}_{\varphi_{i}}$
the sets of dofs for $\mathbf{v}_{0}$, $v_{i}$ and $\boldsymbol{\varphi}_{i}$.
Since each mode $\boldsymbol{\varphi}_{i}$ is periodic within an integer
multiple of the unit cell size [Rokoš et al., 2019], a truncated 2D Fourier
series is a natural choice for its parametrization. As will be shown below,
the selection of relevant sine waves can be easily done by means of a spectral
density analysis of one example for any family of cellular metamaterials and
loading conditions. Note that the number of modes that are activated by
different loading conditions in cellular metamaterials is often only one,
e.g., Fig. 1, and always limited to a few [Rokoš et al., 2020]. Each long-
range correlated fluctuation mode $\boldsymbol{\varphi}_{i}$ is defined such
that it is non-zero where the cellular microstructure exists and decays
linearly to zero inside a region of surrounding (non-cellular) bulk material
over a width of half the unit cell size. The boundary between bulk and
cellular material is a priori selected manually by visual inspection of the
reference image of the undeformed geometry; high accuracy in boundary
identification is not needed. Because $\mathbf{v}_{0}$ and $v_{i}$ are slow
macroscopic fields, smooth functions such as globally-supported polynomials
are the best choice for their parametrization. To this end, Chebyshev
polynomials of $5$th and $6$th order will be used hereafter to parametrize
$\mathbf{v}_{0}$ and $v_{1}$. Since the behaviour of these polynomials outside
the specimen’s domain is irrelevant (because no kinematics is considered
there), potentially large oscillations may be ignored in these regions.

Finally, in any in-situ test under optical visualization, in addition to the
mechanical displacement field, the images are affected by optical distortions
of the optical lenses. Such distortions may introduce large errors if
neglected. Therefore, the residual field of Eq. (2) is rewritten as

$r(\mathbf{x},\underline{\mathbb{a}})=f\left(\boldsymbol{\phi}_{S}(\mathbf{x})\right)-g(\boldsymbol{\phi}_{S}(\boldsymbol{\phi}_{M}(\mathbf{x},\underline{\mathbb{a}}))),$
(5)

where $\boldsymbol{\phi}_{S}(\mathbf{x})=\mathbf{x}+\mathbf{S}(\mathbf{x})$ is
the mapping function describing the spatial distortion [Maraghechi et al.,
2018]. The spatial distortion field, $\mathbf{S}(\mathbf{x})$, is determined a
priori in a calibration step by the method introduced in [Maraghechi et al.,
2019]. Hence, there are no additional degrees of freedom related to spatial
distortion that need to be identified in the micromorphic IDIC problem of Eq.
(3).

##  3 Results and Discussion

To validate the micromorphic IDIC framework, it is applied to both virtual
(Section 3.2) and real in-situ micro-compression experiments (Section 3.3) on
cellular elastomeric metamaterial. Its general applicability to various types
of metamaterial microstructures in terms of spectral analysis is first
demonstrated in Section 3.1 for square and hexagonal stacking of holes, as
well as for a chiral metamaterial, whereas considerations towards graded
microstructures are discussed in Section 3.4.

###  3.1 Parametrization and Initial Guess of the Patterning Modes
$\boldsymbol{\varphi}_{i}$

|  |  |
---|---|---|---
|  |  |
|  |  |

(a) square ref.

(b) square def.

(c) $S_{xx}$ [mm6] square

(d) $S_{yy}$ [mm6] square

(e) hexagonal ref.

(f) hexagonal def.

(g) $S_{xx}$ [mm6] hexagonal

(h) $S_{yy}$ [mm6] hexagonal

(i) chiral ref.

(j) chiral def.

(k) $S_{xx}$ [mm6] chiral

(l) $S_{yy}$ [mm6] chiral

Figure 2: Reference configuration, one particular deformed configuration, and
spectral densities $S_{xx}({\boldsymbol{\xi}})$ and
$S_{yy}({\boldsymbol{\xi}})$, normalized with respect to the unit cell size
$l$, expressed as a function of the wave vector
${\boldsymbol{\xi}}=[\xi_{1},\xi_{2}]$, corresponding to a square stacking of
holes LABEL:sub@fig:modesa–LABEL:sub@fig:modesd, hexagonal stacking of holes
LABEL:sub@fig:modese–LABEL:sub@fig:modesh, and a chiral microstructure
LABEL:sub@fig:modesi–LABEL:sub@fig:modesl. The most of the pattern energies
localize around a few particular wave vectors, allowing in all cases for
effective approximation through a truncated Fourier series.

In the underlying deformed or buckled microstructure, corresponding to an
infinitely large specimen or a periodic arrangement of unit cells, typically a
long-range correlated pattern emerges, implying non-local behaviour. Such an
order is a periodically repeating pattern which can well be represented with a
few lowest-frequency terms of a truncated Fourier series. To show this, we
perform a spectral density analysis [Grigoriu, 2002] on the reference
displacement field, available from numerical simulations obtained for infinite
periodic specimens and for which the affine part of the deformation has been
subtracted. Note that in real experiments the displacement fields will also be
available from regular DIC, as shown below in Section 3.3. Three typical
metamaterial microstructures of square and hexagonal stacking of holes and a
chiral metamaterial are considered. All three geometries are depicted in the
first column of Fig. 2, their corresponding deformed states for various types
of compression (uniaxial or biaxial) in the second column, whereas spectral
densities in the last two columns. The spectral densities are the energy
spectral densities of the displacement components in the horizontal and
vertical directions,
$S_{xx}=\hat{u}_{x}(\mathbf{\xi})\hat{u}_{x}^{*}(\mathbf{\xi})$ and
$S_{yy}=\hat{u}_{y}(\mathbf{\xi})\hat{u}_{y}^{*}(\mathbf{\xi})$, where $u_{x}$
and $u_{y}$ are the $x$\- and $y$-components of $\mathbf{u}(\mathbf{x})$, the
hat indicates the Fourier transformation, and $*$ complex conjugate.
Inspection of spectral density graphs indicates that most of the energy is
localized only around a few particular wave vectors, hence allowing to make a
specific choice on the parametrization and initialization of the long-range
correlated fluctuation field, see, e.g., Rokoš et al. [2019] and Eq. (6)
below, and this parametrization plus initialization is then applicable to a
broad range of metamaterials. Note that in the spectral density graphs there
are multiple secondary peaks present, corresponding to the microfluctuation
$\mathbf{w}$, which are (almost) invisible in the images by the naked eye due
to their small magnitudes.

###  3.2 Virtual Experiments

The micromorphic IDIC scheme will be first tested on virtually generated and
deformed speckle images. By this means the method is tested in a case where
there is no influence of image distortions and noise in the images, and there
are no imperfections in the specimen and the loading, etc. The geometry of the
specimen analysed virtually and later experimentally (although for a different
scale ratio) is depicted in Fig. 3. It has hole diameters $d=1.5$ mm and
centre to centre pitch $l=1.9$ mm (i.e., unit cell size $l$) with edges of 1.9
mm and 0.9 mm bulk material in loading and transverse directions. The length
of the specimen, excluding the bulk edges, is denoted $L$. The scale ratio of
the specimen, $L/l$, is equal to the number of holes in the loading direction.
Details on the specimen fabrication are discussed thoroughly in the PhD thesis
of Maraghechi [2019]. To obtain virtual data, a plane strain finite element
simulation, using quadratic isoparametric triangular elements with three-point
Gauss integration rule and large-strain Total Lagrangian formulation with
stability control and bifurcation, is performed for the compression of a
cellular elastomeric metamaterial specimen with scale ratio of $6$, i.e., $6$
and $10$ holes in loading and transverse directions, respectively. Dirichlet
boundary conditions are applied on the two vertical edges perpendicular to the
loading direction to apply overall compression to the specimen, while the two
remaining edges are free. Upon reaching the critical overall strain, the
imposed deformation triggers an instability in the underlying microstructure
leading to patterning. The material, i.e., PDMS, is modelled by a hyper-
elastic Ogden material model. The material considered, although rubber-like,
is of sufficiently high compressibility (i.e., the initial Poisson’s ratio
$<0.5$) to avoid any numerical problems with locking of elements. The
compressibility of the material has been verified numerically as well as
experimentally. More details on the numerical simulations are presented in PhD
thesis [Maraghechi, 2019]. The displacement fields obtained from this
simulation are used to deform a speckle pattern for two different time steps,
one before microstructural buckling and the other after the emergence of the
patterned fluctuation field, corresponding to $7.8\%$ and $12.8\%$ applied
nominal strain. A total formulation within the context of IDIC was employed,
in which two snapshots in a deformed state (pre-bifurcation and post-
bifurcation images) have been correlated with the reference undeformed image,
although an updated formulation correlating incremental states is equally
possible as well. To generate deformed images, pixel intensities of the
reference image are mapped by displacements, interpolated at nodal or pixel
positions, by inverting the iso-parametric mappings of the underlying FE
approximations. This step is performed to avoid any potential additional
biases due to interpolation errors. Fig. 4a depicts the deformed configuration
of the virtual speckle pattern after the buckling-induced pattern emerged,
which exhibits the pattern transformation more pronounced in the centre and
restricted at the edges, specifically the vertical edges. The employed speckle
pattern is obtained as a thresholded two-dimensional scalar and stationary
random Gaussian field, generated according to [Grigoriu, 2002, Section
5.3.1.2]. The spectral density is chosen as a radially-symmetric normal
probability density function with the mean zero and standard deviation $1/20$.
No additive image white noise has been added to such generated images. The
pattern is blurred using Gaussian filter with the standard deviation of $3$
pixels, resulting in the pattern shown in Fig. 4a.

Figure 3: Analysed cellular elastomeric metamaterial specimen, depicting the
total length of specimen excluding the bulk side edges ($L$: the large dashed
square), the size of the unit cells ($l=1.9$ mm: small dashed square), hole
diameter ($d=1.5$ mm), the size of the bulk side edges and the loading
direction.

As the virtual images are free of spatial distortions,
$\boldsymbol{\phi}_{S}(\mathbf{x})=\mathbf{x}$ is set for the virtual
experiments. Chebyshev polynomials of $5$th and $6$th order are used to
parametrize $\mathbf{v}_{0}$ and $v_{1}$, respectively.

(a) Deformed virtual image

(b) $S_{xx}$ [mm6]

(c) $S_{yy}$ [mm6]

Figure 4: LABEL:sub@fig5:def_virt Virtually generated and deformed speckle
pattern, corresponding to compression in horizontal direction after the onset
of pattern transformation in a cellular elastomeric metamaterial, revealing
the presence of one fluctuation mode inside the specimen. Energy spectral
density functions of the corresponding displacement field in
LABEL:sub@fig5:sxx_virt $x$ and LABEL:sub@fig5:syy_virt $y$ direction,
revealing the principal wave vectors related to the long-range correlated
fluctuation mode (four peaks approximately on the diagonals). The range of the
colour bars is reduced to highlight the frequencies related to the long-range
correlated fluctuation mode.

In order to attain a proper initial guess for
$\boldsymbol{\varphi}_{i}(\mathbf{x},\underline{\mathbb{a}}_{\varphi_{i}})$,
it is first realized that only one long-range correlated fluctuation mode is
triggered in the specimen, as can be easily seen in Fig. 4a, thus $n=1$ is set
in Eq. (1) with only $\boldsymbol{\varphi}_{1}$ considered. It may also be
clear that the triggered mode is periodic with a periodicity of two unit
cells, recall Figs. 2c–LABEL:sub@fig:modesd where the peaks are localized at
$\xi_{1}=\xi_{2}=\pm 1/2$. The horizontal and vertical lines in Figs. 4b and
4c (not present in Figs. 2c–LABEL:sub@fig:modesd) correspond to the
discontinuities in the mean deformation due to the edges of the specimen,
while the four peaks correspond to the correlated fluctuation mode, revealing
two sine functions roughly in diagonal directions of the Cartesian coordinate
system for the $x$ and $y$ components of $\boldsymbol{\varphi}_{1}$. The
purpose of the spectral density analysis is merely the rough identification of
the type of long-range correlated fluctuation modes, so only the corresponding
four peaks are considered here. Inspection of $\hat{u}_{x}(\mathbf{\xi})$ and
$\hat{u}_{y}(\mathbf{\xi})$, at the frequencies corresponding to the four
peaks of $S_{xx}$ and $S_{yy}$, suggests that the two sine waves are of the
same sign in the $x$ component, and of the opposite sign in the $y$ component
of the displacement field. This inspection also suggests that the two sine
waves are of the same amplitude, for both $x$ and $y$ components of the
displacement in agreement with the symmetry of the problem. Based on these
observations, $\boldsymbol{\varphi}_{1}$ is parametrized as the sum of two
sine functions for each component, i.e.,

$\displaystyle\boldsymbol{\varphi}_{1}(\mathbf{x},\underline{\mathbb{a}}_{\varphi_{1}})$
$\displaystyle=a_{1}\left[-\sin\frac{\pi}{l}(a_{2}x+a_{3}y+a_{6})-\sin\frac{\pi}{l}(a_{4}x+a_{5}y+a_{7})\right]\mathbf{e}_{x}$
(6)
$\displaystyle\>\quad+\left[\sin\frac{\pi}{l}(a_{2}x+a_{3}y+a_{6})-\sin\frac{\pi}{l}(a_{4}x+a_{5}y+a_{7})\right]\mathbf{e}_{y},$

where
$\underline{\mathbb{a}}_{\varphi_{1}}=[a_{1},a_{2},a_{3},a_{4},a_{5},a_{6},a_{7}]^{\mathsf{T}}$
is the column of degrees of freedom describing the mode. Note that the
direction of both sine waves and their wavelengths are free and independent,
which compensates for inaccuracies in the measurement of the unit cell size
due to processing uncertainties. This allows for an accurate identification of
the shape and magnitude of the fluctuation mode in the micromorphic IDIC
analysis, whereas the spectral density analysis is only used for a proper
parametrization of the fluctuation mode and for an approximate initial guess
of all these parameters. Note also that it is sufficient to perform the
spectral density analysis only once for a family of cellular metamaterials
triggering the same pattern. In cases where more than one mode is activated,
the spectral analysis should be done separately on each of the subregions
containing the individual modes, cf. Section 3.1.

The dofs describing $\boldsymbol{\varphi}_{1}$ are initialized as:
$\underline{\mathbb{a}}_{\varphi_{1}}=[1,1,1,1,1,0,0]^{\mathsf{T}}$. The
initial values for $a_{2}$ to $a_{5}$ are set to a wavelength of $2$ unit
cells and their signs and ratios are set such that the initial orientation of
the two sine waves is aligned with the two diagonals of the coordinate system.
Note that $a_{1}$ is only defining the ratio of the amplitude of the mode in
$x$ and $y$ direction, while $v_{1}(\mathbf{x})$ defines its absolute
amplitude according to Eq. (4). In order to assure a proper scaling of the
minimization problem, $\boldsymbol{\varphi}_{1}(\mathbf{x})$ is scaled such
that the initial maximum value of the mode is $0.04$. Normalization of the
microfluctuation field $\boldsymbol{\varphi}_{i}$ is an inherent choice of the
employed micromorphic computational homogenization formulation. We refer to
[Rokoš et al., 2020, Eq. (3)], for more details. The scaling
$\max_{\mathbf{x}}\|\boldsymbol{\varphi}_{i}(\mathbf{x})\|=0.04$ is chosen to
avoid poor conditioning of the resulting Hessian matrix during the solution
procedure using the Gauss–Newton algorithm. Moreover, since the sensitivity
functions of the mode $\boldsymbol{\varphi}_{1}$ used in the Gauss–Newton
optimization algorithm [Neggers et al., 2016] are a function of
$v_{1}(\mathbf{x})$, the initial guess for the latter needs to be non-zero to
avoid an ill-posed problem in the first iteration. To this end, coefficients
of linear combination of the Chebyshev polynomials corresponding to the
$v_{1}(\mathbf{x})$ field are initiated such that a small constant non-zero
spatial field results. The applied compressive global strain results in large
displacement on the edges of the specimen requiring an approximate initial
guess in order to ensure convergence. Thus, the smooth mean field
$\mathbf{v}_{0}(\mathbf{x})$ is initialized such that the first order term in
$x$ direction approximately accounts for the applied global strain. The IDIC
procedure overall identifies $77$ dofs, out of which $2\cdot 21=42$ correspond
to coefficients of linear combination $\underline{\mathbb{a}}_{v_{0}}$ of the
Chebyshev polynomials of 5th order describing the field
$\mathbf{v}_{0}(\mathbf{x})$, $28$ to coefficients of linear combination
$\underline{\mathbb{a}}_{v_{1}}$ of the Chebyshev polynomials of 6th order
describing the field $v_{1}(\mathbf{x})$, and $7$ corresponding to the degrees
of freedom $\underline{\mathbb{a}}_{\varphi_{1}}$ of the long-range correlated
fluctuation field $\boldsymbol{\varphi}_{1}(\mathbf{x})$ specified in Eq. (6).
This is a negligible amount of dofs compared to local DIC which needs at least
tens of thousands dofs to capture such locally-varying displacement fields.

Before microstructural buckling

(a) $\mathbf{v}_{0}.\mathbf{e}_{x}$ [$\mu$m]

(b) $\mathbf{v}_{0}.\mathbf{e}_{y}$ [$\mu$m]

(c) $v_{1}$ [$\mu$m]

(d) $\boldsymbol{\varphi}_{1}\leavevmode\nobreak\ {[-]}$

(e) $v_{1}\boldsymbol{\varphi}_{1}.\mathbf{e}_{x}$ [$\mu$m]

(f) $v_{1}\boldsymbol{\varphi}_{1}.\mathbf{e}_{y}$ [$\mu$m]

(g) residual [–]

(h) $||\mathbf{w}||$ [$\mu$m]

After microstructural buckling

(i) $\mathbf{v}_{0}.\mathbf{e}_{x}$ [$\mu$m]

(j) $\mathbf{v}_{0}.\mathbf{e}_{y}$ [$\mu$m]

(k) $v_{1}$ [$\mu$m]

(l) $\boldsymbol{\varphi}_{1}\leavevmode\nobreak\ {[-]}$

(m) $v_{1}\boldsymbol{\varphi}_{1}.\mathbf{e}_{x}$ [$\mu$m]

(n) $v_{1}\boldsymbol{\varphi}_{1}.\mathbf{e}_{y}$ [$\mu$m]

(o) residual [–]

(p) $||\mathbf{w}||$ [$\mu$m]

Figure 5: Results of the virtual experiments before (left) and after (right)
the onset of microstructural buckling: (LABEL:sub@fig5:virt1v0x,
LABEL:sub@fig5:virt1v0y, LABEL:sub@fig5:virt2v0x, LABEL:sub@fig5:virt2v0y) the
smooth mean displacement fields, $\mathbf{v}_{0}(\mathbf{x})$, in
(LABEL:sub@fig5:virt1v0x, LABEL:sub@fig5:virt2v0x) $x$ and
(LABEL:sub@fig5:virt1v0y, LABEL:sub@fig5:virt2v0y) $y$ directions;
(LABEL:sub@fig5:virt1v1, LABEL:sub@fig5:virt2v1) the pattern amplitudes,
$v_{1}(\mathbf{x})$, (LABEL:sub@fig5:virt1phi, LABEL:sub@fig5:virt2phi) the
deformed configuration of identified patterns/modes,
$\boldsymbol{\varphi}_{1}(\mathbf{x})$, scaled to assure the visibility of the
deformed shape, and (LABEL:sub@fig5:virt1flucx, LABEL:sub@fig5:virt1flucy,
LABEL:sub@fig5:virt2flucx, LABEL:sub@fig5:virt2flucy) the resulting fields,
(LABEL:sub@fig5:virt1flucx, LABEL:sub@fig5:virt2flucx) in $x$ and
(LABEL:sub@fig5:virt1flucy, LABEL:sub@fig5:virt2flucy) $y$ directions,
corresponding to the long-range correlated fluctuations. In
(LABEL:sub@fig5:virt1phi, LABEL:sub@fig5:virt2phi)
$\boldsymbol{\varphi}_{1}(\mathbf{x})$ are scaled to assure the visibility of
the deformed shape of the modes. The long-range correlated fluctuation fields,
$v_{1}(\mathbf{x})\boldsymbol{\varphi}_{1}(\mathbf{x})$, are plotted in the
deformed configuration corresponding to their contribution only. The residual
fields in normalized gray scales are in (LABEL:sub@fig5:virt1r,
LABEL:sub@fig5:virt2r), and (LABEL:sub@fig5:virt1w, LABEL:sub@fig5:virt2w)
show the Euclidean norm of the microfluctuation field,
$\mathbf{w}(\mathbf{x})$, plotted in deformed configuration based on the total
displacement.

Fig. 5 shows the results of the virtual experiments before (left) and after
(right) the emergence of the buckling pattern. The micromorphic IDIC results
in an accurate identification of the $\mathbf{v}_{0}(\mathbf{x})$,
$v_{1}(\mathbf{x})$ and $\boldsymbol{\varphi}_{1}(\mathbf{x})$, considering
the residual fields that are small everywhere both before and after buckling
occurs, see Figs. 5g and 5o. For sensible initial guesses the convergence is
consistent, but the procedure may be trapped in local minima during
minimization for poor initial guesses. For each case, the mean smooth
displacement field, $\mathbf{v}_{0}(\mathbf{x})$, in $x$ and $y$ directions,
and $v_{1}(\mathbf{x})$ are depicted in the first two rows, i.e., Figs.
5a–LABEL:sub@fig5:virt1v1 before and Figs. 5i–LABEL:sub@fig5:virt2v1 after
buckling occurs. The $\mathbf{v}_{0}$ fields show compression in $x$ direction
both before and after buckling (a linear profile). Before pattern
transformation, lateral expansion is visible in the $y$ direction, while
afterwards the auxetic effect, i.e., lateral contraction, is clearly visible
near the specimen’s centre. As expected, the $v_{1}$ field is close to zero
before buckling occurs, although its absolute value increases slightly in the
corners of the region with cellular microstructure (the area of the specimen
excluding the bulk edges), which is due to the fact that the fluctuation mode,
$\boldsymbol{\varphi}_{1}$, partially captures the deformation that occurs in
the corner holes where the edge effect is considerable. After the buckling
occurs, the $v_{1}$ field is close to zero on the edges and increases towards
the centre, as expected. Again, $v_{1}$ is not exactly zero at the edges and
the corners of the metamaterial, since the fluctuation mode partially
describes the deformation in the corners where the edge effect is
considerable. Note that the large values in $v_{1}(\mathbf{x})$ are due to the
scaling of the corresponding mode $\boldsymbol{\varphi}_{1}$ (maximum value of
0.04).

The identified long-range correlated mode,
$\boldsymbol{\varphi}_{1}(\mathbf{x})$, is depicted in Figs. 5d and 5l before
and after buckling, where their amplitudes are scaled to assure the visibility
of the deformed shape of the modes. The small $v_{1}$ prior to buckling
entails a low sensitivity of the mode parameters. As a result,
$\boldsymbol{\varphi}_{1}$ is significantly perturbed with respect to the
initial guess, see Fig. 5d where the periodicity does not match twice the unit
cell size any more. After buckling, however, the mode becomes representative,
and it resembles the initial guess more closely and matches the mode observed
in the simulations (Fig. 5l). The long-range correlated fluctuation fields,
$v_{1}\boldsymbol{\varphi}_{1}(\mathbf{x})$ in $x$ and $y$ directions, are
depicted in their deformed configurations prior to buckling in Figs. 5e and
5f, and beyond buckling in Figs. 5m and 5n. Prior to buckling,
$v_{1}\boldsymbol{\varphi}_{1}(\mathbf{x})$ is close to zero everywhere apart
from the four corners, where the identified mode, although significantly
perturbed with respect to the initial guess, partially captures the edge
effect described above. Post buckling, $\boldsymbol{\varphi}_{1}$ is
identified properly, resulting in small modifications to the wavelengths and
orientations of the sine waves compared to the initial guess based on spectral
density analysis.

The residual fields in normalized grey scales before and after buckling are
shown in Figs. 5g and 5o, in the deformed configuration based on the total
evaluated displacement $\mathbf{u}(\mathbf{x})$. The residual fields are
consistently small everywhere over the region of interest. However, the
microfluctuation fields, discarded in the micromorphic IDIC scheme which seeks
to extract the correlated fields only, are by definition reflected in the
residual fields. In order to evaluate the performance of the method, the
microfluctuation field $\mathbf{w}(\mathbf{x})$ is recovered as the difference
of the reference displacement field (here available from the simulation) and
the displacement field $\mathbf{u}(\mathbf{x})$ according to Eq. (4). The
Euclidean norm of the microfluctuation field, $||\mathbf{w}||$, is depicted in
Figs. 5h and 5p before and after the buckling, respectively. Note that the
micromorphic IDIC scheme is constrained to identify displacement fluctuations
with spatial frequencies bound to the highest spatial frequencies in the long-
range correlated fluctuation modes, $\boldsymbol{\varphi}_{i}$. Thus, the
microfluctuation field, $\mathbf{w}$, contains all parts of the field with
spatial frequencies higher than the characteristic frequency of
$\boldsymbol{\varphi}_{1}$, as well visible far from the edges of the specimen
in the post buckling case, see Fig. 5p. As expected, the microfluctuation
field, $\mathbf{w}$, is of a significantly smaller amplitude than the
correlated fluctuation field beyond buckling (compare Fig. 5p with Figs. 5m
and 5n). The opposite is true prior to buckling, due to the almost negligible
contribution of the long-range correlated fluctuation field (compare Fig. 5h
with Figs. 5e and 5f). Note the difference in scale bars for the figures of
the pre- and post-bifurcation state.

In order to assess the robustness of the method, the correlations of the
virtual experiments, after the onset of buckling, are repeated with different
initial guess values for the fluctuation mode. The values of
$\underline{\mathbb{a}}_{\varphi_{i}}$ are perturbed by a random value up to
$10\%$ of the nominal value stated above, resulting in random changes in the
wavelengths and orientations of the sine waves as well as the ratio of the $x$
and $y$ components of $\boldsymbol{\varphi}_{1}$. The perturbed correlations
are repeated 100 times, and the error in each case is evaluated as:
$\epsilon=100\times\frac{||\mathbf{u}||-||\mathbf{u}_{pert}||}{||\mathbf{u}||}$,
where $\mathbf{u}$ is the total displacement field evaluated with no
perturbation of initial guess, and $\mathbf{u}_{pert}$ is the perturbed total
displacement field. The error averaged over all the repetitions is $0.86\%$,
entailing only two cases with an erroneous assessment of the total
displacement field, i.e., $\epsilon>1.5\%$, confirming the $98\%$ robustness
of the method to initial guess. The average error $\epsilon$, ignoring the two
outlier cases, is $0.45\%$, confirming also the high accuracy of the
methodology.

###  3.3 Real Experiments

The results of an in-situ compression test on a cellular elastomeric
metamaterial specimen, manufactured using a customized mould as presented in
PhD thesis of Maraghechi [2019], are analysed with the micromorphic IDIC
method. The real experiment includes a range of additional complexities
compared to the virtual experiment, including the image distortions and image
capturing noise as well as the imperfections in the specimens and loading
conditions of the test. The specimen is $22.5$ mm thick and has 12 and 10
holes in loading and transverse directions, corresponding to a scale ratio 12.
The specimen is made of Sylgard 184, with a 10:1 mixing ratio. A two-layer
speckle pattern is deposited on the top surface (white powder sprayed to make
a fine-grained background and India ink sprayed using an air brush, see Fig.
6a), resulting in a speckle size between 30 and $80\leavevmode\nobreak\ \mu$m
(3 to 8 pixels). An in-situ compression test is performed with a Kammrath &
Weiss micro tensile/compression stage with a $50$ N load cell, while a Zeiss
V20 stereo microscope is used to acquire images at each load step at
$7.6\times$ magnification and $2751\times 2207$ pixels. The stereo microscope
and the objective was used for proper imaging of the specimen providing high-
enough resolution for the fine speckle pattern applied on the specimen with an
adequate level of detail. More details on the specimen processing and the in-
situ test can be found in PhD thesis of Maraghechi [2019]. Two images obtained
during the in-situ test, before and after the emergence of the buckling
pattern, corresponding to $8.2\%$ and $12.7\%$ applied nominal strain on the
metamaterial, are correlated to the undeformed reference configuration. The
deformed configuration after emergence of the buckling pattern is shown in
Fig. 6a.

A spatial distortion calibration step Maraghechi et al. [2018, 2019], is
performed to measure the spatial distortions of the imaging system at the
desired magnification. The resulting distortion field, which is as large as
$10\%$ of the actual displacements, is used to construct the distortion
mapping function $\mathbf{S}(\mathbf{x})$.

Chebyshev polynomials of 5th and 6th order are used to parametrize
$\mathbf{v}_{0}$ and $v_{1}$.

(a) Deformed Image

(b) $S_{xx}$ [mm6]

(c) $S_{yy}$ [mm6]

Figure 6: LABEL:sub@fig5:def_real Image of the deformed configuration of a
cellular elastomeric metamaterial with millimetre-sized circular holes in a
rectangular stacking, corresponding to a compressive load step after the onset
of pattern transformation, acquired during an in-situ compression test. Energy
spectral density functions of the displacement field, assessed by local DIC,
in LABEL:sub@fig5:sxx_real $x$ and LABEL:sub@fig5:syy_real $y$ direction,
revealing the principal frequencies related to the long-range correlated
fluctuation mode (four peaks approximately on the diagonals). The range of the
colour bars is reduced to highlight the frequencies related to the long-range
correlated fluctuation mode.

The same procedure as in the case of virtual experiments is performed for
estimating the reduced regularization of
$\boldsymbol{\varphi}_{1}(\mathbf{x})$ and the initialization of the
associated parameters. To this end, first, the post buckling kinematics is
determined using local DIC, in order to attain $S_{xx}$ and $S_{yy}$, shown in
Figs. 6b and 6c. Similar to the virtual experiments, vertical and horizontal
lines relate to discontinuities in the smooth mean deformations due to the
specimen’s edges, while the two sine waves in approximately diagonal
directions for each displacement component correspond to a single long-range
fluctuation mode. Based on these observations, the fluctuation mode is
assigned the same way as for the virtual experiments, i.e., Eq. (6) with
$\underline{\mathbb{a}}_{\varphi_{1}}=[1,1,1,1,1,0,0]^{\mathsf{T}}$ as initial
guess. In order to assure a proper scaling of the minimization problem,
$\boldsymbol{\varphi}_{1}(\mathbf{x})$ is scaled such that the initial maximum
value of the mode is $0.04$. In order to avoid an ill-posed problem in the
first iteration, $v_{1}(\mathbf{x})$ is again initiated with a small constant
value in space. The applied compressive global strain results in large
displacement on the edges of the specimen, which requires an approximate
initial guess in order to ensure convergence. Thus, the smooth mean field
$\mathbf{v}_{0}(\mathbf{x})$ is initialized such that the first order term in
$x$ direction approximately accounts for the applied global strain. Similarly
to the virtual experiment of Section 3.2, the IDIC procedure identifies in
total $77$ dofs $\underline{\mathbb{a}}$.

Before microstructural buckling

(a) $\mathbf{v}_{0}.\mathbf{e}_{x}$ [$\mu$m]

(b) $\mathbf{v}_{0}.\mathbf{e}_{y}$ [$\mu$m]

(c) $v_{1}$ [$\mu$m]

(d) $\boldsymbol{\varphi}_{1}\leavevmode\nobreak\ {[-]}$

(e) $v_{1}\boldsymbol{\varphi}_{1}.\mathbf{e}_{x}$ [$\mu$m]

(f) $v_{1}\boldsymbol{\varphi}_{1}.\mathbf{e}_{y}$ [$\mu$m]

(g) residual [–]

(h) $||\mathbf{w}||$ [$\mu$m]

After microstructural buckling

(i) $\mathbf{v}_{0}.\mathbf{e}_{x}$ [$\mu$m]

(j) $\mathbf{v}_{0}.\mathbf{e}_{y}$ [$\mu$m]

(k) $v_{1}$ [$\mu$m]

(l) $\boldsymbol{\varphi}_{1}\leavevmode\nobreak\ {[-]}$

(m) $v_{1}\boldsymbol{\varphi}_{1}.\mathbf{e}_{x}$ [$\mu$m]

(n) $v_{1}\boldsymbol{\varphi}_{1}.\mathbf{e}_{y}$ [$\mu$m]

(o) residual [–]

(p) $||\mathbf{w}||$ [$\mu$m]

Figure 7: Results of the real experiments, before (left) and after (right) the
onset of microstructural buckling:(LABEL:sub@fig5:real1v0x,
LABEL:sub@fig5:real1v0y, LABEL:sub@fig5:real2v0x, LABEL:sub@fig5:real2v0y) the
smooth mean displacement fields, $\mathbf{v}_{0}(\mathbf{x})$, in
(LABEL:sub@fig5:real1v0x, LABEL:sub@fig5:real2v0x) $x$ and
(LABEL:sub@fig5:real1v0y, LABEL:sub@fig5:real2v0y) $y$ directions;
(LABEL:sub@fig5:real1v1, LABEL:sub@fig5:real2v1) the amplitudes,
$v_{1}(\mathbf{x})$, (LABEL:sub@fig5:real1phi, LABEL:sub@fig5:real2phi) the
deformed configuration of identified modes,
$\boldsymbol{\varphi}_{1}(\mathbf{x})$, scaled to assure the visibility of the
deformed shape, and (LABEL:sub@fig5:real1flucx, LABEL:sub@fig5:real1flucy,
LABEL:sub@fig5:real2flucx, LABEL:sub@fig5:real2flucy) the resulting fields,
(LABEL:sub@fig5:real1flucx, LABEL:sub@fig5:real2flucx) in $x$ and
(LABEL:sub@fig5:real1flucy, LABEL:sub@fig5:real2flucy) $y$ directions,
corresponding to the long-range correlated fluctuations. The long-range
correlated fluctuation fields,
$v_{1}(\mathbf{x})\boldsymbol{\varphi}_{1}(\mathbf{x})$, are plotted in their
corresponding deformed configurations. (LABEL:sub@fig5:real1r,
LABEL:sub@fig5:real2r) are the residual fields in normalized gray scales; and
(LABEL:sub@fig5:real1w, LABEL:sub@fig5:real2w) the Euclidean norm of the
microfluctuation field, $\mathbf{w}(\mathbf{x})$, plotted in deformed
configuration based on the evaluated displacement, $\mathbf{u}(\mathbf{x})$.

For these real _in-situ_ compressive experiments, the micromorphic IDIC also
results in a proper and accurate identification of
$\mathbf{v}_{0}(\mathbf{x})$, $v_{1}(\mathbf{x})$ and
$\boldsymbol{\varphi}_{1}(\mathbf{x})$, whereby the residual fields are small
everywhere both before and after buckling, see Figs. 7g and 7o. Convergence is
consistent. Fig. 7 shows the results of applying micromorphic IDIC on real
experiments before (left) and after (right) the emergence of the buckling
pattern. For each case, the mean smooth displacement field,
$\mathbf{v}_{0}(\mathbf{x})$, in $x$ and $y$ directions, and
$v_{1}(\mathbf{x})$ are depicted in the first two rows, i.e., Figs.
7a–LABEL:sub@fig5:real1v1 before and Figs. 7i–LABEL:sub@fig5:real2v1 after the
onset of buckling. The $\mathbf{v}_{0}$ field shows compression in $x$
direction both before and after buckling. Similar results as in the virtual
experiments are observed in $\mathbf{v}_{0}$ and $v_{1}$. Note that the large
values in $v_{1}(\mathbf{x})$ are due to the scaling of the corresponding mode
$\boldsymbol{\varphi}_{1}$ (maximum value of 0.04). The auxetic effect is more
pronounced in the real experiments (Figs. 7j,LABEL:sub@fig5:real2r) compared
to the virtual experiments (Figs. 5j,LABEL:sub@fig5:virt2r), due to the larger
specimen size, i.e., the larger scale ratio. In a specimen with a larger scale
ratio the buckling pattern in the centre is less affected by the two
stiffening lateral edges and, thus, it is more pronounced and results in
larger auxetic effect.

The long-range correlated mode, $\boldsymbol{\varphi}_{1}(\mathbf{x})$, as
well as the long-range correlated fluctuation fields,
$v_{1}\boldsymbol{\varphi}_{1}(\mathbf{x})$, in $x$ and $y$ directions, are
depicted in their deformed configurations in Figs.
7d–LABEL:sub@fig5:real1flucy and Figs. 7l–LABEL:sub@fig5:real2flucy, before
and after buckling, respectively. In Figs. 7d and 7l,
$\boldsymbol{\varphi}_{1}(\mathbf{x})$ is scaled to assure the visibility of
the deformed shape of the modes. As expected, small values for $v_{1}$ prior
to buckling entail a low sensitivity to the fluctuation mode, i.e., a slightly
different mode is identified at this stage. Indeed, the long-range correlated
fluctuation mode is as good as absent before the emergence of the buckling
pattern, where it is not yet representative. The edge effect prior to
buckling, which is most pronounced at the corner unit cells, is partially
captured by the identified mode. This explains the non-zero values of $v_{1}$
and thus $v_{1}(\mathbf{x})\boldsymbol{\varphi}_{1}(\mathbf{x})$ in these
areas.

The residual fields in normalized grey scales are shown in Figs. 7g and 7o, in
the deformed configuration based on the evaluated displacement
$\mathbf{u}(\mathbf{x})$ from Eq. (4), before and after buckling. In order to
extract the microfluctuation field $\mathbf{w}(\mathbf{x})$, the displacement
fields obtained by local DIC on the same images are used as the reference.
Note that the local DIC procedure is not a part of the methodology introduced
here, and is used to assess the microfluctuation field (with errors from local
DIC) in order to evaluate the performance of the method. Although local DIC
gives the displacement field with minimum kinematical constraints, yet with
significant statistical and systematic errors, it requires tens of thousands
of dofs to fit the displacement field compared to the $77$ dofs for the
micromorphic IDIC framework. Figs. 7h and 7p show the Euclidean norm of the
microfluctuation field, $||\mathbf{w}||$, before and after buckling. The
microfluctuation field $\mathbf{w}$, is of a higher spatial frequency than the
characteristic frequency of $\boldsymbol{\varphi}_{1}$, which is more
pronounced far from edges of the specimen in the post-buckling regime (Fig.
7p). Yet, it can be observed prior to buckling as well (Fig. 7h). Beyond
buckling, the microfluctuation field $\mathbf{w}$ has a significantly smaller
amplitude than the correlated fluctuation field. This is not the case prior to
buckling, considering the almost zero contribution of the long-range
correlated fluctuation field. The residual field is larger in certain areas,
which corresponds directly to higher values in the microfluctuation field,
both before and after buckling.

All the identified fields, i.e., $\mathbf{v}_{0}(\mathbf{x})$,
$v_{1}(\mathbf{x})$ and $\boldsymbol{\varphi}_{1}(\mathbf{x})$, satisfy the
expected mechanical behaviour, and the correlation residual, reflected also in
the evaluated microfluctuation field $\mathbf{w}(\mathbf{x})$, confirms the
good performance of the micromorphic IDIC method.

###  3.4 Graded Microstructures and Size Independence of the Modes

As long as the constitutive law of the base elastomeric material is size
independent (e.g., the hyperelastic Ogden material used in Section 3.2,
usually used for constitutive modelling of metamaterials), the modes
$\boldsymbol{\varphi}_{i}$ are independent of the unit cell size $l$. To show
that the modes are independent also of $d/l$ ratio (i.e., hole diameter to
cell size ratio), an infinitely large specimen with a square stacking of holes
is considered, cf. Figs. 2a–LABEL:sub@fig:modesb. In such a case the buckling
analysis carried out on a $2\times 2$ periodic cell provides the reference
(exact) mode—exact within the accuracy of the adopted numerical model—,
against which the approximation of $\boldsymbol{\varphi}_{1}$ given by Eq. (6)
can be compared. Results are obtained for a larger and smaller $d/l$ ratio of,
respectively, 0.95 and 0.25, compared to the reference value of 0.79
identified by the micromorphic IDIC procedure, see Fig. 8. The deformed
configuration of the reference mode is shown in black colour, whereas the
deformed shape of the experimentally identified mode is shown in light grey.
Note that only one quarter of the $2\times 2$ periodic cell is shown due to
various symmetries [see e.g. Bertoldi et al., 2008]. Fig. 8a suggests that the
overall shape of the mode is captured sufficiently accurately for $d/l=0.95$
(recall also Fig. 5), although local discrepancies can be observed in higher-
frequency components. Such components are not visible in Figs.
4b–LABEL:sub@fig5:syy_virt and 6b–LABEL:sub@fig5:syy_real, since their
contributions to spectral densities are very small. More interestingly,
comparable errors are observed also in Fig. 8b for $d/l=0.79$, meaning that
identified mode coefficients $\underline{\mathbb{a}}_{\varphi_{1}}$ do not
change substantially as a function of small variations of $d/l$. For a very
different ratio, $d/l=0.25$, however, larger discrepancies between the two
modes are observed in Fig. 8c, suggesting that the experimental mode should be
identified again. Although this insight, i.e., that the modes are independent
for moderate variations in the hole diameter to cell size ratio, has been
demonstrated here for a square stacking of holes, it is a general feature of
cellular metamaterials.

(a) $d/\ell=0.95$

(b) $d/\ell=0.79$ (reference)

(c) $d/\ell=0.25$

Figure 8: Patterning modes for an infinite microstructure with a square
stacking of holes for a unit cell size $l=1.9$ mm and a varying hole diameter
$d\in\\{1.81,1.50,0.48\\}$ mm. Numerically computed modes are shown in black,
whereas experimentally obtained modes according to Eq. (6) in Section 3.3 are
shown in light grey. Large hole diameter to cell size ratio $d/l=0.95$ is
shown in LABEL:sub@fig7:gradeda, $d/l=0.79$ corresponding to the reference
geometry analysed in detail in Section 3 is shown in LABEL:sub@fig7:gradedb,
while small hole diameter $d/l=0.25$ is shown in LABEL:sub@fig7:gradedc.

This brings us to the realization that the introduced micromorphic IDIC
technique, combined with the previously developed micromorphic computational
homogenization scheme, provides a general tool for the design of engineering
systems at two scales.  (i) Size effects and microstructural changes are
communicated through the patterning modes to the macroscale, where such
information can be utilized for the design of engineering systems.  (ii)
Microstructures with graded properties smoothly varying in space (e.g.,
$d(\mathbf{x})$) or material property (e.g., Young’s modulus, $E(\mathbf{x})$,
or Poission’s ratio, $\nu(\mathbf{x})$), can be optimally adjusted and
understood at the micro- or meso-scale level. For more details on the
modelling and application of functionally graded materials the reader is
referred to the extensive literature on this topic, see, e.g., [Anthoine,
2010, El-Helou and Harne, 2019, Francesconi et al., 2019, Farzaneh et al.,
2022, Morris et al., 2023].  It is important to realize that considerable
design freedom is available even for geometrically simple microstructures,
such as the herein studied square stacking of holes. Adopting different types
of microstructural morphologies (e.g., hexagonal stacking of holes, chiral or
re-entrant microstructures) expands the design space further.

##  4 Summary and Conclusions

Cellular elastomeric metamaterials generally exhibit a pattern transformations
based on microstructural buckling at a critical level of compressive load.
These pattern transformations are correlated over long ranges in the specimen
and result in considerable changes in the overall mechanical response after
their emergence. A micromorphic Integrated Digital Image Correlation
(micromorphic IDIC) method, was here developed to experimentally identify and
quantify the decomposed kinematics for cellular metamaterials. The
decomposition follows a recently introduced micromorphic computational
homogenization scheme based on a kinematical ansatz consisting of three parts:
(i) mean smooth displacements, corresponding to the slow-scale material level;
(ii) long-range correlated fluctuation fields, representing the buckling
patterns; and  (iii) a microfluctuation field, which captures the remaining
uncorrelated part of the fluctuations.  This kinematical ansatz was integrated
within the digital image correlation scheme, enabling a direct identification
of the smooth displacement, and the long-range correlated fluctuation pattern
along with its spatial distribution field. On the basis of spectral analysis,
a 2D Fourier series approximation of the long-range fluctuation mode with a
proper parametrization was performed. The identified patterning modes serve as
input parameters for the micromorphic computational homogenization scheme,
which is capable of providing mechanical predictions of engineering systems
including size effects as well as temporal and spatial mixing of patterning
modes.

The methodology is validated on both virtual and real experiments for a
specific case of elastomeric metamaterials with rectangular stacking of
millimetre-sized circular holes. A standard spectral density analysis is used
to attain a regularization with a reduced number of degrees of freedom and to
initialize the associated parameters of the long-range fluctuation mode for
the considered metamaterial. The method successfully decomposes the
displacement field of cellular elastomers under compressive load, both before
and after local microstructural buckling leading to long-range correlated
pattern transformations. The performance of the method is assessed by
identifying the microfluctuation fields through a comparison with local DIC
data. In all cases, the microfluctuation field is of higher spatial frequency
than the long-range fluctuations and of smaller amplitude after the emergence
of the buckling pattern, thereby validating the kinematical ansatz as well as
the micromorphic IDIC identification routine. An initial guess robustness
study, performed by random perturbations up to $10\%$ in fluctuation mode
parameters, resulted in $98\%$ robustness of the correlations and negligible
errors ($0.45\%$ on average) in the results, confirming the robustness and the
accuracy of the proposed method.

The micromorphic IDIC methodology may be considered novel and unique in the
following aspects:

  1. (i)

The methodology aims at and succeeds in identifying long-range correlated
fluctuation fields, which is not existing yet in commercial DIC packages.

  2. (ii)

The decomposition of the displacement field is performed in a single
minimization step, identifying the long-range correlated fluctuation mode, its
amplitude in space, and the mean smooth displacement field;

  3. (iii)

The methodology can be readily extended to any cellular metamaterial by making
proper initialization choices of the fluctuation modes;

  4. (iv)

It is easy to implement based on minimal modifications to conventional global
DIC formulations.

The proposed micromorphic IDIC technique, in combination with the previously
developed micromorphic computational homogenization scheme, thus provides an
integrated experimental-computational tool for optimizing the design of
metamaterial-based engineering systems, including materials with various
metamaterial microstructures and/or with graded properties in space, such as
the hole size, or material properties, such as the Young’s modulus or
Poisson’s ratio.

## Acknowledgements

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].

## References

  * Ameen et al. [2018] Ameen, M.M., Rokoš, O., Peerlings, R.H.J., Geers, M.G.D., 2018.  Size effects in nonlinear periodic materials exhibiting reversible pattern transformations.  Mechanics of Materials 124, 55–70.  doi:10.1016/j.mechmat.2018.05.011.
  * Anthoine [2010] Anthoine, A., 2010.  Second-order homogenisation of functionally graded materials.  International Journal of Solids and Structures 47, 1477 – 1489.  doi:https://doi.org/10.1016/j.ijsolstr.2010.02.004.
  * Auger et al. [2020] Auger, P., Lavigne, T., Smaniotto, B., Spagnuolo, M., Dell’Isola, F., Hild, F., 2020.  Poynting effects in pantographic metamaterial captured via multiscale DVC.  Journal of Strain Analysis for Engineering Design 56, 462–477.
  * Barchiesi et al. [2021] Barchiesi, E., Dell’Isola, F., Hild, F., 2021.  On the Validation of Homogenized Modeling for Bi-Pantographic Metamaterials via Digital Image Correlation.  International Journal of Solids and Structures , 49–62doi:10.1016/j.ijsolstr.2020.09.036.
  * 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.  doi:10.1016/j.jmps.2008.03.006.
  * Bertoldi et al. [2010] Bertoldi, K., Reis, P.M., Willshaw, S., Mullin, T., 2010.  Negative Poisson’s ratio behavior induced by an elastic instability.  Advanced Materials 22, 361–366.  doi:10.1002/adma.200901956.
  * Bertoldi et al. [2017] Bertoldi, K., Vitelli, V., Christensen, J., van Hecke, M., 2017.  Flexible mechanical metamaterials.  Nature Reviews Materials 2, 17066.  doi:10.1038/natrevmats.2017.66.
  * Bilal et al. [2017] Bilal, O.R., Foehr, A., Daraio, C., 2017.  Reprogrammable phononic metasurfaces.  Advanced materials 29, 1700628.  doi:10.1002/adma.201700628.
  * Biswas and Poh [2017] Biswas, R., Poh, L., 2017.  A micromorphic computational homogenization framework for heterogeneous materials.  Journal of the Mechanics and Physics of Solids 102, 187–208.  doi:10.1016/j.jmps.2017.02.012.
  * Blaysat et al. [2015] Blaysat, B., Hoefnagels, J.P.M., Lubineau, G., Alfano, M., Geers, M.G.D., 2015.  Interface debonding characterization by image correlation integrated with double cantilever beam kinematics.  International Journal of Solids and Structures 55, 79 – 91.  doi:https://doi.org/10.1016/j.ijsolstr.2014.06.012.
  * 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.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0045782520305181, doi:10.1016/j.cma.2020.113333. arXiv: 2008.12850.
  * dell’Isola et al. [2019] dell’Isola, F., Seppecher, P., Alibert, J.J., Lekszycki, T., Grygoruk, R., Pawlikowski, M., Steigmann, D., Giorgio, I., Andreaus, U., Turco, E., Gołaszewski, M., Rizzi, N., Boutin, C., Eremeyev, V.A., Misra, A., Placidi, L., Barchiesi, E., Greco, L., Cuomo, M., Cazzani, A., Corte, A.D., Battista, A., Scerrato, D., Eremeeva, I.Z., Rahali, Y., Ganghoffer, J.F., Müller, W., Ganzosch, G., Spagnuolo, M., Pfaff, A., Barcz, K., Hoschke, K., Neggers, J., Hild, F., 2019.  Pantographic metamaterials: an example of mathematically driven design and of its technological challenges.  Continuum Mech. Thermodyn. 31, 851–884.  URL: https://doi.org/10.1007/s00161-018-0689-8, doi:10.1007/s00161-018-0689-8.
  * Easey et al. [2019] Easey, N., Chuprynyuk, D., Musa, W.W., Banks, A., Dobah, Y., Shterenlikht, A., Scarpa, F., 2019.  Dome-shape auxetic cellular metamaterials: manufacturing, modelling and testing.  Frontiers in Materials 6, 86.  doi:10.3389/fmats.2019.00086.
  * El-Helou and Harne [2019] El-Helou, C., Harne, R.L., 2019.  Exploiting functionally graded elastomeric materials to program collapse and mechanical properties.  Advanced Engineering Materials 21, 1900807.  doi:10.1002/adem.201900807, arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/adem.201900807.
  * Farzaneh et al. [2022] Farzaneh, A., Pawar, N., Portela, C.M., Hopkins, J.B., 2022.  Sequential metamaterials with alternating Poisson’s ratios.  Nature Communications 13.  doi:10.1038/s41467-022-28696-9.
  * Findeisen et al. [2020] Findeisen, C., Forest, S., Hohe, J., Gumbsch, P., 2020.  Discrete and continuum modelling of size effects in architectured unstable metamaterials.  Continuum Mechanics and Thermodynamics doi:10.1007/s00161-020-00870-8.
  * Florijn et al. [2014] Florijn, B., Coulais, C., Van Hecke, M., 2014.  Programmable mechanical metamaterials.  Physical Review Letters 113, 1–5.  doi:10.1103/PhysRevLett.113.175503.
  * Forest [2009] Forest, S., 2009.  Micromorphic Approach for Gradient Elasticity, Viscoplasticity, and Damage.  Journal of Engineering Mechanics 135, 117–131.  doi:10.1061/(ASCE)0733-9399(2009)135:3(117).
  * Francesconi et al. [2019] Francesconi, L., Baldi, A., Liang, X., Aymerich, F., Taylor, M., 2019.  Variable Poisson’s ratio materials for globally stable static and dynamic compression resistance.  Extreme Mechanics Letters 26, 1–7.  doi:10.1016/j.eml.2018.11.001.
  * Giorgio et al. [2023] Giorgio, I., dell’Isola, F., Steigmann, D.J., 2023.  Second-grade elasticity of three-dimensional pantographic lattices: theory and numerical experiments.  Continuum Mech. Thermodyn. URL: https://doi.org/10.1007/s00161-023-01240-w, doi:10.1007/s00161-023-01240-w.
  * Grigoriu [2002] Grigoriu, M., 2002.  Stochastic Calculus: Applications in Science and Engineering.  Birkhäuser.
  * Hedayati et al. [2019] Hedayati, R., Jedari Salami, S., Li, Y., Sadighi, M., Zadpoor, A.A., 2019.  Semianalytical geometry-property relationships for some generalized classes of pentamodelike additively manufactured mechanical metamaterials.  Physical Review Applied 11.  doi:10.1103/PhysRevApplied.11.034057.
  * Hild et al. [2020] Hild, F., Ching, W.Y., Misra, A., Dell’Isola, F., 2020.  Multiscale DIC applied to Pantographic Structures.  Experimental Mechanics 61, 431–443.  doi:10.1007/s11340-020-00636-y.
  * Hu et al. [2013] Hu, J., He, Y., Lei, J., Liu, Z., 2013.  Novel mechanical behavior of periodic structure with the pattern transformation.  Theoretical and Applied Mechanics Letters 3, 054007.  doi:10.1063/2.1305407.
  * Hu et al. [2018] Hu, J., Zhou, Y., Liu, Z., 2018.  The Friction Effect on Buckling Behavior of Cellular Structures Under Axial Load.  International Journal of Applied Mechanics 10, 1850013.  doi:10.1142/s1758825118500138.
  * Hütter [2017] Hütter, G., 2017.  Homogenization of a Cauchy continuum towards a micromorphic continuum.  Journal of the Mechanics and Physics of Solids 99, 394–408.  doi:10.1016/j.jmps.2016.09.010.
  * Javid et al. [2016] Javid, F., Liu, J., Shim, J., Weaver, J.C., Shanian, A., Bertoldi, K., 2016.  Mechanics of instability-induced pattern transformations in elastomeric porous cylinders.  Journal of the Mechanics and Physics of Solids 96, 1–17.  doi:10.1016/j.jmps.2016.06.015.
  * Kolken and Zadpoor [2017] Kolken, H.M., Zadpoor, A., 2017.  Auxetic mechanical metamaterials.  RSC Advances 7, 5111–5129.  doi:10.1039/c6ra27333e.
  * Kouznetsova et al. [2001] Kouznetsova, V., Brekelmans, W.A.M., Baaijens, F.P.T., 2001.  An approach to micro-macro modeling of heterogeneous materials.  Computational Mechanics 27, 37–48.  doi:10.1007/s004660000212.
  * Kouznetsova et al. [2004] Kouznetsova, V., Geers, M., Brekelmans, W., 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.  doi:10.1016/j.cma.2003.12.073.
  * Krishnan and Johnson [2009] Krishnan, D., Johnson, H.T., 2009.  Optical properties of two-dimensional polymer photonic crystals after deformation-induced pattern transformations.  Journal of the Mechanics and Physics of Solids 57, 1500–1513.  doi:10.1016/j.jmps.2009.05.012.
  * Li et al. [2023] Li, L., Yao, H., Mi, S., 2023.  Magnetically Driven Modular Mechanical Metamaterials with High Programmability, Reconfigurability, and Multiple Applications 15, 3486–3496.  doi:10.1021/acsami.2c19679.
  * Lin et al. [2023] Lin, Q., Zhou, J., Wang, K., Xu, D., Wen, G., Wang, Q., 2023.  Three-dimensional quasi-zero-stiffness metamaterial for low-frequency and wide complete band gap 307, 116656–116656.  doi:10.1016/j.compstruct.2022.116656.
  * Maraghechi [2019] Maraghechi, S., 2019.  Experimental full-field assessment of kinematical and imaging fluctuation patterns.  Ph.D. thesis. Department of Mechanical Engineering.  Proefschrift.
  * Maraghechi et al. [2018] Maraghechi, S., Hoefnagels, J.P.M., Peerlings, R.H.J., Geers, M.G.D., 2018.  Correction of scan line shift artifacts in scanning electron microscopy: An extended digital image correlation framework.  Ultramicroscopy 187, 144–163.  doi:10.1016/j.ultramic.2018.01.002.
  * Maraghechi et al. [2019] Maraghechi, S., Hoefnagels, J.P.M., Peerlings, R.H.J., Rokoš, O., Geers, M.G.D., 2019.  Correction of scanning electron microscope imaging artifacts in a novel digital image correlation framework.  Experimental Mechanics 59, 489–516.  doi:10.1007/s11340-018-00469-w.
  * Maraghechi et al. [2020] Maraghechi, S., Hoefnagels, J.P.M., Peerlings, R.H.J., Rokoš, O., Geers, M.G.D., 2020.  Experimental full-field analysis of size effects in miniaturized cellular elastomeric metamaterials.  Materials & Design 193, 108684.  doi:10.1016/j.matdes.2020.108684.
  * Miehe and Koch [2002] Miehe, C., Koch, A., 2002.  Computational micro-to-macro transitions of discretized microstructures undergoing small strains.  Archive of Applied Mechanics (Ingenieur Archiv) 72, 300–317.  doi:10.1007/s00419-002-0212-2.
  * Mindlin [1964] Mindlin, R.D., 1964.  Micro-structure in linear elasticity.  Arch. Rational Mech. Anal. 16, 51–78.  URL: https://doi.org/10.1007/BF00248490, doi:10.1007/BF00248490.
  * 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.  URL: https://www.sciencedirect.com/science/article/pii/0020768365900065, doi:10.1016/0020-7683(65)90006-5.
  * Mirzaali et al. [2018a] Mirzaali, M.J., Caracciolo, A., Pahlavani, H., Janbaz, S., Vergani, L., Zadpoor, A.A., 2018a.  Multi-material 3D printed mechanical metamaterials: Rational design of elastic properties through spatial distribution of hard and soft phases.  Applied Physics Letters 113.  doi:10.1063/1.5064864.
  * Mirzaali et al. [2018b] Mirzaali, M.J., Janbaz, S., Strano, M., Vergani, L., Zadpoor, A.A., 2018b.  Shape-matching soft mechanical metamaterials.  Scientific Reports 8, 1–7.  doi:10.1038/s41598-018-19381-3.
  * Morris et al. [2023] Morris, J., Wang, W., Plaisted, T., Hansen, C., Amirkhizi, A., 2023.  Optimization of Graded Metamaterials for Control of Energy Transmission Using a Genetic Algorithm ARXIV_ID: 2306.10628.
  * Mullin et al. [2007] Mullin, T., Deschanel, S., Bertoldi, K., Boyce, M.C., 2007.  Pattern transformation triggered by deformation.  Physical review letters 99, 084301.  doi:10.1103/PhysRevLett.99.084301.
  * Neggers et al. [2016] Neggers, J., Blaysat, B., Hoefnagels, J.P.M., Geers, M.G.D., 2016.  On image gradients in digital image correlation.  International Journal for Numerical Methods in Engineering 105, 243–260.  doi:10.1002/nme.4971.
  * Neggers et al. [2015] Neggers, J., Hoefnagels, J.P.M., Geers, M.G.D., Hild, F., Roux, S., 2015.  Time-resolved integrated digital image correlation.  International Journal for Numerical Methods in Engineering 103, 157–182.  doi:10.1002/nme.4882.
  * Overvelde and Bertoldi [2014] Overvelde, J.T.B., Bertoldi, K., 2014.  Relating pore shape to the non-linear response of periodic elastomeric structures.  Journal of the Mechanics and Physics of Solids 64, 351–366.  doi:10.1016/j.jmps.2013.11.014.
  * Overvelde et al. [2012] Overvelde, J.T.B., Shan, S., Bertoldi, K., 2012.  Compaction through buckling in 2D periodic, soft and porous structures: Effect of pore shape.  Advanced Materials 24, 2337–2342.  doi:10.1002/adma.201104395.
  * 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.  doi:https://doi.org/10.1016/j.jmps.2018.08.019. the N.A. Fleck 60th Anniversary Volume.
  * Rokoš et al. [2020] Rokoš, O., Ameen, M.M., Peerlings, R.H.J., Geers, M.G.D., 2020.  Extended micromorphic computational homogenization for mechanical metamaterials exhibiting multiple geometric pattern transformations.  Extreme Mechanics Letters 37, 100708.  doi:https://doi.org/10.1016/j.eml.2020.100708.
  * Roux and Hild [2006] Roux, S., Hild, F., 2006.  Stress intensity factor measurements from digital image correlation: post-processing and integrated approaches.  International Journal of Fracture 140, 141–157.  doi:10.1007/s10704-006-6631-2.
  * Ruybalid et al. [2016] Ruybalid, A.P., Hoefnagels, J.P.M., van der Sluis, O., Geers, M.G.D., 2016.  Comparison of the identification performance of conventional FEM updating and integrated DIC.  International Journal for Numerical Methods in Engineering 106, 298–320.  doi:10.1002/nme.5127.
  * Shen et al. [2023] Shen, J., Garrad, M., Zhang, Q., Leao, O., Pirrera, A., Groh, R., 2023.  Active reconfiguration of multistable metamaterials for linear locomotion 107.  doi:10.1103/physrevb.107.214103.
  * Slann et al. [2015] Slann, A., White, W., Scarpa, F., Boba, K., Farrow, I., 2015.  Cellular plates with auxetic rectangular perforations.  Physica Status Solidi (b) 252, 1533–1539.
  * Tang et al. [2015] Tang, Y., Lin, G., Han, L., Qiu, S., Yang, S., Yin, J., 2015.  Design of hierarchically cut hinges for highly stretchable and reconfigurable metamaterials with enhanced strength.  Advanced Materials 27, 7181–7190.
  * Toupin [1962] Toupin, R.A., 1962.  Elastic materials with couple-stresses.  Arch. Rational Mech. Anal. 11, 385–414.  URL: https://doi.org/10.1007/BF00253945, doi:10.1007/BF00253945.
  * Wormser et al. [2017] Wormser, M., Wein, F., Stingl, M., Körner, C., 2017.  Design and additive manufacturing of 3D phononic band gap structures based on gradient based optimization.  Materials 10, 1125.
  * Wu et al. [2014] Wu, G., Xia, Y., Yang, S., 2014.  Buckling, symmetry breaking, and cavitation in periodically micro-structured hydrogel membranes.  Soft Matter 10, 1392–1399.  doi:10.1039/c3sm51640g.
  * Xu and Pasini [2016] Xu, H., Pasini, D., 2016.  Structurally efficient three-dimensional metamaterials with controllable thermal expansion.  Scientific reports 6, 34924.
  * Yang et al. [2015] Yang, D., Mosadegh, B., Ainla, A., Lee, B., Khashai, F., Suo, Z., Bertoldi, K., Whitesides, G.M., 2015.  Buckling of Elastomeric Beams Enables Actuation of Soft Machines.  Advanced Materials 27, 6323–6327.  doi:10.1002/adma.201503188.
  * Yavari et al. [2015] Yavari, S.A., Ahmadi, S.M., Wauthle, R., Pouran, B., Schrooten, J., Weinans, H., Zadpoor, A.A., 2015.  Relationship between unit cell type and porosity and the fatigue behavior of selective laser melted meta-biomaterials.  Journal of the mechanical behavior of biomedical materials 43, 91–100.
  * Yu et al. [2018] Yu, X., Zhou, J., Liang, H., Jiang, Z., Wu, L., 2018.  Mechanical metamaterials associated with stiffness, rigidity and compressibility: A brief review.  Progress in Materials Science 94, 114–173.  doi:10.1016/j.pmatsci.2017.12.003.
  * Zadpoor [2017] Zadpoor, A., 2017.  Design for additive bio-manufacturing: From patient-specific medical devices to rationally designed meta-biomaterials.  International journal of molecular sciences 18, 1607.
  * Zhao et al. [2015] Zhao, D.G., Li, Y., Zhu, X.F., 2015.  Broadband lamb wave trapping in cellular metamaterial plates with multiple local resonances.  Scientific reports 5, 9376.