# Single-test evaluation of directional elastic properties of anisotropic
structured materials

Jagannadh Boddapati111These authors contributed equally. Moritz
Flaschel222These authors contributed equally. Siddhant Kumar  Laura De
Lorenzis  Chiara Daraio<EMAIL_ADDRESS>Division of Engineering and
Applied Science, California Institute of Technology, Pasadena, CA 91125, USA
Department of Mechanical and Process Engineering, ETH Zürich, 8092 Zürich,
Switzerland  Department of Materials Science and Engineering, Delft University
of Technology, 2628 CD Delft, The Netherlands

###### Abstract

When the elastic properties of structured materials become direction-
dependent, the number of their descriptors increases. For example, in two-
dimensions, the anisotropic behavior of materials is described by up to 6
independent elastic stiffness parameters, as opposed to only 2 needed for
isotropic materials. Such high number of parameters expands the design space
of structured materials and leads to unusual phenomena, such as materials that
can shear under uniaxial compression. However, an increased number of
properties descriptors and the coupling between shear and normal deformations
render the experimental evaluation of material properties more challenging. In
this paper, we propose a methodology based on the virtual fields method to
identify six separate stiffness tensor parameters of two-dimensional
anisotropic structured materials using just one tension test, thus eliminating
the need for multiple experiments, as it is typical in traditional methods.
The approach requires no stress data and uses full-field displacement data and
global force data. We show the accuracy of our method using synthetic data
generated from finite element simulations as well as experimental data from
additively manufactured specimens.

###### keywords:

Anisotropy , Shear-normal coupling , Virtual fields method , Metamaterial
design , Data-driven identification , Inverse problems

††journal: Composites Part B: Engineering

##  1 Introduction

The advent of additive manufacturing has allowed the design and engineering of
a new class of materials known as metamaterials, or structured/architected
materials. Mechanical metamaterials are a special branch of metamaterials that
derive special functionalities from their peculiar deformation, dynamic motion
and/or elastic energy distribution (Lee et al., 2012; Christensen et al.,
2015; Zadpoor, 2016; Bertoldi et al., 2017; Surjadi et al., 2019).
Metamaterials derive their effective properties from both the micro- and meso-
structure and their constitutive material properties. They often exhibit
mechanical properties that deviate from those of their constituent materials,
showing unusual behaviors, such as negative Poisson’s ratios (Greaves et al.,
2011), vanishing shear moduli (Kadic et al., 2012), and negative refractive
indices (Kaina et al., 2015).

By carefully selecting the geometry of the micro- and meso-structures with
varying symmetries (Milton and Cherkaev, 1995; Kadic et al., 2012; Wu et al.,
2019; Kulagin et al., 2020; Mao et al., 2020; Bastek et al., 2022),
metamaterial designers can explore novel anisotropy classes in the material
responses. In turn, the presence of rich anisotropy expands the materials’
functionality space, by exploiting coupled-deformation mechanisms that are
non-existent in symmetric structures. Examples include metamaterials that
twist under compression (Frenzel et al., 2017; Chen et al., 2018; Wu et al.,
2019; Yuan et al., 2021), shear under thermal loading (Ni et al., 2019) and
shape-morph (Guseinov et al., 2020; Risso et al., 2021; Agnelli et al., 2022).
In the dynamic regime, anisotropy allows observing phenomena like conical
refraction (Ahn et al., 2017) and control of broadband elastic waves (Zheng et
al., 2019; Yang et al., 2019; Zheng et al., 2020).

In a two-dimensional continuum, the elastic behavior of an anisotropic
material is described using six independent elastic parameters (Ting and Chi-
Tsai, 1996). In experiments, characterizing these many independent elastic
parameters is quite complex. Indeed, the presence of shear-normal coupling
makes it hard to measure even one of the six parameters from a single
experiment. Prior work suggested different approaches to experimentally
measure the elastic parameters for different anisotropy classes (Schittny et
al., 2013; Gras et al., 2015; Lee et al., 2016; Kim et al., 2020; Agnelli et
al., 2021). However, most of these approaches focus on measuring the stiffness
tensor components when the off-diagonal, shear-normal coupling, components are
absent. In addition, several of these approaches require multiple experimental
steps. For example, techniques based on the detection of different acoustic
wave speeds along different material directions involve multiple tests and
assume a certain material symmetry in predicting elastic parameters (Every and
Sachse, 1990; François et al., 1998). To date, there are no experimental
methods that can measure the stiffness parameters of fully anisotropic
structures from a single experiment.

Traditional material parameter identification methods rely on single-load
experimental setups with homogeneous (constant) strain distributions within
the tested specimen, which allow the derivation of closed-form stress-strain
relations. However, the amount of data that can be acquired through a one
dimensional tension test, for example, is limited (e.g., one stress-strain
data pair for each measurement). When characterizing complex materials,
multiple experimental setups with different loading conditions are needed.
Full-field identification methods allow extracting additional information from
single-load experiments. Measuring the full displacement field, e.g., through
Digital Image Correlation (DIC), of arbitrarily shaped specimens under loading
maximizes the amount of data generated from a single experimental test. Such
data can then be used to characterize the material by applying inverse
identification methods such as, among others, Finite Element Model Updating,
the Equilibrium Gap Method or the Virtual Fields Method (VFM), see (Avril et
al., 2008; Roux and Hild, 2020; Pierron, 2023) for a review.

These methods have in common that they are used to calibrate the parameters of
an a priori chosen material model, i.e., the mathematical functions and
operations that describe the material response need to be fixed by means of
the intuition or modeling experience of the user. However, the selection of
inappropriate a priori assumptions about the model and its underlying
mathematical structure can introduce errors. Recent research used full-field
data to train machine-learning-models, whose versatile ansatz spaces promise
to mitigate modeling errors. Flaschel et al. (2021), for example, proposed the
method EUCLID (Efficient Unsupervised Constitutive Law Identification and
Discovery) that uses sparse regression (Tibshirani, 1996) informed by full-
field displacement data and net reaction force data, to automatically select
interpretable material models from a potentially large predefined set of
candidate material models. EUCLID has been applied to hyperelasticity
(Flaschel et al., 2021), elastoplasticity (Flaschel et al., 2022),
viscoelasticity (Marino et al., 2023), and generalized standard materials
(Flaschel et al., 2023), see Flaschel (2023) for an overview. Further, EUCLID
was formulated in a Bayesian setting by Joshi et al. (2022) to simultaneously
perform model selection and quantification of uncertainty in the material
parameters. In contrast to selecting interpretable material models through
sparse regression, full-field data may also be used to train black-box
material model surrogates like neural networks, as shown by Man and Furukawa
(2011); Huang et al. (2020); Liu et al. (2020) for small strain elasticity and
by Thakolkaran et al. (2022) for hyperelasticity. In the present work, it is
assumed that the material response does not leave the realm of elasticity at
infinitesimal strains. Thus, the material model can be assumed to be known a
priori, and its parameters are calibrated with the VFM.

The VFM, originally proposed by Grédiac (1989) (see also Grédiac et al.
(2008); Pierron and Grédiac (2012)), employs the balance of linear momentum in
its weak form, to identify unknown material parameters. The VFM method assumes
that the kinematic fields in the specimen, as well as the reaction forces at
the boundaries, are known from experiments. As such, material parameters
remain the only unknowns in the balance equations and can be calculated using
standard linear or nonlinear solvers. In essence, the VFM describes the
inverse problem to the classical Finite Element Method (FEM). The method has
been applied in various cases, such as small-strain elasticity, elasto-
plasticity (Grédiac and Pierron, 2006), and hyperelasticity (Promma et al.,
2009), among others.

The accuracy of the VFM in identifying unknown material parameters and its
sensitivity to noise are highly dependent on the choice of the functions for
which the weak linear momentum balance is tested, also known as the virtual
displacement fields. A distinction can be made between global virtual fields
that are defined over the whole specimen domain, such as polynomials, and
local virtual fields with compact support, such as in the Bubnov-Galerkin
discretization with piecewise polynomial shape functions. As the choice of the
virtual fields is arbitrary and user-dependent, several attempts have been
made to automate and optimize it (Avril et al., 2004; Pierron et al., 2010;
Marek et al., 2017).

In this article, full-field measurement based identification, and in
particular the VFM, is explored in the context of anisotropic structured
materials and compared to traditional identification methods. We focus in
particular on the identification of shear-normal coupling parameters,
notoriously complex to extract from conventional experiments. The rest of the
paper is organized as follows. In Section 2, we discuss the theory of
anisotropic linear elasticity and introduce our model setup used for parameter
identification. In Section 3, we present our virtual fields method. In Section
4, we describe our experimental and numerical data acquisition methods. In
Section 5, we discuss our results, including experimental validation, and we
draw our conclusions in Section 6.

##  2 Material model and geometry

In this section, we review the fundamental equations of linear elasticity at
infinitesimal strains and introduce our model setup used to identify the
governing material parameters of anisotropic metamaterials.

###  2.1 Anisotropic linear elasticity

Under the small strain assumption, the constitutive law for a general
anisotropic solid, which relates the Cauchy stress tensor $\bm{\sigma}$ and
the infinitesimal strain tensor $\bm{\varepsilon}$, is given by the
generalized Hooke’s law (Rychlewski, 1984; Ting and Chi-Tsai, 1996),

$\bm{\sigma}=\bm{C\varepsilon}\hskip 10.0pt\text{or}\hskip
10.0pt(\sigma_{ij}=C_{ijkl}\varepsilon_{kl}),$ (1)

where $C$ is a fourth-order tensor, known as the elasticity tensor or the
stiffness tensor, and Einstein’s notation for summation over repeated indices
is followed. For a two-dimensional anisotropic solid, under plane stress
conditions, Eq. 1 can be written using Voigt notation as

$\left[\begin{array}[]{l}\sigma_{11}\\\ \sigma_{22}\\\
\sigma_{12}\end{array}\right]=\left[\begin{array}[]{lll}C_{1111}&C_{1122}&C_{1112}\\\
C_{1122}&C_{2222}&C_{2212}\\\ C_{1112}&C_{2212}&C_{1212}\\\
\end{array}\right]\left[\begin{array}[]{c}\varepsilon_{11}\\\
\varepsilon_{22}\\\ 2\varepsilon_{12}\end{array}\right],$ (2)

where $C_{1111},\ C_{1122},\ C_{2222},\ C_{1112},\ C_{2212},\ C_{1212}$ are
the elasticity tensor parameters in a given reference frame,
$\varepsilon_{11},\varepsilon_{22}$ are the axial strains, $\varepsilon_{12}$
is the shear strain, $\sigma_{11},\sigma_{22}$ are the axial stresses, and
$\sigma_{12}$ is the shear stress. For readability, we combine the pair of
indices as follows:
$()_{11}\rightarrow()_{1},()_{22}\rightarrow()_{2},()_{12}\rightarrow()_{6}$
and write Eq. 2 as

$\left[\begin{array}[]{l}\sigma_{{1}}\\\ \sigma_{{2}}\\\
\sigma_{{6}}\end{array}\right]=\left[\begin{array}[]{lll}C_{11}&C_{12}&C_{16}\\\
C_{12}&C_{22}&C_{26}\\\ C_{16}&C_{26}&C_{66}\\\
\end{array}\right]\left[\begin{array}[]{c}\varepsilon_{{1}}\\\
\varepsilon_{{2}}\\\ 2\varepsilon_{{6}}\end{array}\right].$ (3)

Our objective is to identify these six material parameters
$C_{11},C_{12},C_{22},C_{16},C_{26},C_{66}$ from experimental measurements
while fulfilling certain constraints. From thermodynamic constraints, the
elasticity tensor has to be positive definite, which implies

$\displaystyle C_{11}>0,\quad C_{22}>0,\quad C_{66}>0,$ (4a) $\displaystyle
C_{11}C_{22}-C_{12}^{2}>0,\quad C_{11}C_{66}-C_{16}^{2}>0,\quad
C_{22}C_{66}-C_{26}^{2}>0.$ (4b)

The stiffness parameter $C_{12}$ represents the extension-to-extension
deformation coupling. The stiffness parameters $C_{16},C_{26}$ represent the
extension-to-shear coupling, also known as shear-normal coupling, which
induces shear stress from axial strains, and axial stresses from shear
strains. Shear-normal coupling has been explored in the context of structured
materials by Karathanasopoulos et al. (2020); Dos Reis and Karathanasopoulos
(2022). As a result of these anisotropy-induced couplings, the experimental
identification of the material parameters becomes non-trivial because a
constant state of strain is hard to achieve, even in a standard uniaxial
tension test.

Note that the parameters $C_{16}$ and $C_{26}$ will be zero if the material
has symmetry planes along the $x_{1}$ and $x_{2}$ axes. Thus, the existence of
shear-normal coupling and the maximum number of independent stiffness tensor
parameters depend on the symmetries associated with the material microscopic
topology (Ting and Chi-Tsai, 1996; Podestá et al., 2019). In plane elasticity,
stiffness tensors are categorized into four symmetry classes. They are denoted
as $O(2)$ for Isotropic, $D_{4}$ for Tetragonal, $D_{2}$ for Orthotropic and
$Z_{2}$ for Digonal (fully anisotropic) with 2, 3, 4 and 6 independent
parameters respectively. This categorization is based on the invariants of the
stiffness tensor (Forte and Vianello, 2014; Auffray and Ropars, 2016).
However, in our methods of parameter identification, we do not consider any
prior information on the material symmetries or the number of independent
material parameters.

###  2.2 Model setup

Without loss of generality, we study two-dimensional structured solids,
obtained from finite periodic tessellation of square unit cells (Fig. 1)333Our
methods are easily extendable to non-square unit cells.. We focus on
identifying the effective anisotropic material parameters of these composite
assemblies, as linear elastic continua.

To design unit cells, we follow an approach inspired by Cahn’s method of
generating Gaussian random fields by superposing plane waves of fixed
wavelength but random in phase and direction (Cahn, 1965; Soyarslan et al.,
2018; Kumar et al., 2020). We first define a function $f(x_{1},x_{2})$, as a
linear superposition of cosine periodic functions:

$f(x_{1},x_{2})=\sum_{m,n}A_{mn}\cos\left(2\pi(mx_{1}+nx_{2})\right),\quad\forall(x_{1},x_{2})\in[-0.5,0.5],\quad\forall
m,n\in[-3,-2,-1,0,1,2,3],$ (5)

where $m,n$ are spatial frequencies, and $A_{mn}$ are the corresponding cosine
function weights. The function is then thresholded at a value $\xi$, to
generate a binary image which represents a unit cell, as shown in Fig. 1,
panels a, b. Each unit cell is pixelated and discretized with a 100 $\times$
100 square mesh. In this pixelated representation, the gray phase represents a
stiffer material and the black phase represents a softer material (see Section
4.3.1).

The periodicity is ensured from the choice of the cosine functions directly.
We randomly sample the weights $A_{mn}$ and the threshold value $\xi$ to
generate a small database of unit cells (about 100), from which we pick four
unit cells to study in this paper. The four unit cells are chosen such that
they are diverse in anisotropic properties and suitable for additive
manufacturing (see Section 4.1). We consider a unit cell as suitable for
manufacturing if the stiff phase is connected in the finite periodic
tessellation with a minimum feature size of 5 pixels, matching the resolution
of our chosen additive manufacturing approach.

A schematic of our setup is shown in Fig. 1c. A two-dimensional square
anisotropic structured solid with $10\times 10$ unit cell tessellation, with
side length L, is subjected to a displacement-controlled tension test. The
boundary conditions are such that the bottom end is fixed, while a
displacement of $\bm{u}=[0,u_{p}]^{T}$ is prescribed at the top end. The
reaction force components measured at the fixed end are denoted as
$F_{1},F_{2}$.

Figure 1: a) Design of an anisotropic unit cell geometry by thresholding a
periodic function $f(x_{1},x_{2})$. b) A two-phase unit cell geometry
consisting of a stiffer (gray) and a softer phase (black). c) A two-
dimensional anisotropic metamaterial created by tessellating the unit cell
geometry (shown in the inset) ten times along both $x_{1}$\- and $x_{2}$\-
axes.

##  3 Virtual fields method for anisotropic metamaterials

Many parameter identification methods rely on conducting multiple experiments,
which are time consuming, complex and require specialized equipment. To
circumvent these drawbacks, we explore a material characterization method
based on the VFM that solely relies on full-field displacements and net
reaction force measurements from a single experimental test. In this section,
after discussing the assumptions underlying the adoption of the VFM for
metamaterials, we outline all the components of the proposed method.

###  3.1 Basic assumptions

The VFM (Grédiac, 1989; Grédiac et al., 2008; Pierron and Grédiac, 2012)
exploits the weak formulation of linear momentum balance, i.e., the principle
of virtual work, as a constraint on the material parameter space. Since the
full displacement field over the specimen and the net reaction forces at the
specimen boundaries are known, testing the weak formulation for a suitable set
of test functions (also known as virtual fields) results in a system of
equations that can be solved for the unknown material parameters. By choosing
the test functions as not constant in space, the linear momentum balance is
tested in different regions of the considered specimen domain. As such, the
VFM takes advantage of the local strain data, as opposed to global methods for
parameter identification.

In the following, the VFM is used to characterize the mechanical behavior of
metamaterials. However, it should be noted that – due to the non-homogeneous
nature of the metamaterials – the application of identification methods based
on full-field measurements is not trivial. Full-field measurement techniques
such as DIC measure the kinematic fields locally, i.e., at several points on
the considered specimen surface. The studied metamaterials are not expected to
behave at these local points as their homogenized counterparts, especially
when the number of repeating unit cells is low in comparison to the size of
the specimen. To give an example, in Section 5.1.1 the deformation of a
heterogeneous metamaterial specimen will be compared to that of an equally-
dimensioned homogeneous body, whose stiffness is set to the homogenized
stiffness of the metamaterial. Under the same loading conditions, the two
specimens exhibit different local displacements, which is likely caused by
local size effects and the different boundary conditions that are assumed
during the loading of the macroscopic structure and the homogenization of the
microscopic unit cell. It is observed that deviations between the kinematic
fields are predominant at the boundary and in particular at the corners of the
domain. This agrees with theoretical studies on heterogeneous metamaterials,
which suggest the usage of non-local – e.g., higher-order strain-gradient
based – theories as proposed by (Mindlin and Eshel, 1968), to model size
effects and wedge forces appearing at corners of non-homogeneous bodies
(Fischer et al., 2011; Andreaus et al., 2016; Yang et al., 2021). Within this
work, such theories are avoided for the sake of simplicity and to keep a
reasonably low number of material parameters. Hence, the assumption is made
that the global material behavior of the metamaterials can be characterized
based on local kinematic measurements within a local constitutive theory. As
we will see later, this assumption will introduce errors in the identification
procedure, which are, however, below a practically relevant level. During the
development of the VFM, we found that the locally measured kinematic data must
be treated with care, especially at the boundary and the corners of the
specimen. We will later introduce specifically designed virtual fields that
reduce the influence of data acquired at the specimen boundary and corners
(see Section 3.5 for details).444We note at this point that reducing the
influence of data acquired at the specimen boundary and corners may be
beneficial not only when studying heterogeneous materials. Even for
homogeneous specimens, the acquisition of kinematic data at the specimen
boundary via DIC is known to be difficult.

###  3.2 Required data

To identify the unknown parameters, the VFM needs diverse local strain data,
i.e., strain fields that are not homogeneous. Therefore, data that serve as
input for the VFM are usually generated by testing complex specimen geometries
under complex loading conditions. For our purposes we will show that, due to
the anisotropy of the material, a clamped square plate under uniaxial tension
produces a sufficiently heterogeneous strain field. We hence consider a
displacement-controlled uniaxial tension experiment of a square-shaped
specimen that consists of $n_{c}\times n_{c}$ repeating square unit cells of
the considered metamaterial (Fig. 1). At the fixed boundary of the specimen, a
load cell measures the net reaction force. Further, the full-field deformation
of the specimen is tracked through DIC, which measures the local displacements
of the solid material. After preprocessing the data, the VFM takes as input
the displacement measurements at the $(n_{c}+1)\times(n_{c}+1)$ unit cell
corners and the net reaction forces. A quadrilateral finite element mesh is
generated such that each of the $n_{c}\times n_{c}$ elements corresponds to
one unit cell and the element nodes correspond to the unit cell corners with
experimentally known displacement values. The continuous displacement field
$\bm{u}(\bm{x})$ is hence approximated by

$\displaystyle\bm{u}(\bm{x})=\sum_{a=1}^{n_{n}}N^{a}(\bm{x})\bm{u}^{a},$ (6)

where $n_{n}=(n_{c}+1)^{2}$ denotes the number of nodes in the finite element
mesh and $\bm{u}^{a}$ are the known nodal displacements, while $N^{a}(\bm{x})$
are the standard ansatz functions of bilinear quadrilateral finite elements.
The infinitesimal strain field is then obtained as the symmetric gradient of
the displacement field, i.e.,
$\bm{\varepsilon}(\bm{x})=\frac{1}{2}\left(\nabla\bm{u}(\bm{x})+\left(\nabla\bm{u}(\bm{x})\right)^{T}\right)$.

###  3.3 Weak formulation of linear momentum balance

We denote the specimen domain and its boundary as $\Omega$ and
$\partial\Omega$, respectively, and the surface traction force acting on
$\partial\Omega$ as $\bm{t}$. Assuming no inertia and body forces, the weak
form of linear momentum balance reads

$\displaystyle\int_{\Omega}\bm{\sigma}(\bm{x})\colon\nabla\bm{v}(\bm{x})\
\mathrm{d}A-\int_{\partial\Omega}\bm{t}\cdot\bm{v}(\bm{x})\ \mathrm{d}s=0,$
(7)

which has to hold true for all admissible, i.e., sufficiently regular, test
functions $\bm{v}(\bm{x})$. Note that we are not introducing the classical
distinction between Dirichlet and Neumann portions of the boundary;
accordingly, we are not requiring admissible test functions to vanish
anywhere.

###  3.4 Discretization

The weak form of linear momentum balance has to hold true for any chosen set
of admissible test functions. Here, we adopt the standard (Bubnov-Galerkin)
approach and express the test functions as a linear combination of the same
shape functions $N^{a}(\bm{x})$ used to interpolate the displacement data

$\displaystyle\bm{v}(\bm{x})=\sum_{a=1}^{n_{n}}N^{a}(\bm{x})\bm{v}^{a}.$ (8)

Inserting the test function ansatz into the weak form of linear momentum
balance results in

$\displaystyle\sum_{a=1}^{n_{n}}\bm{v}^{a}\cdot\left[\underbrace{\int_{\Omega}\bm{\sigma}\nabla
N^{a}(\bm{x})\
\mathrm{d}A}_{\bm{F}_{\text{int}}^{a}}-\underbrace{\int_{\partial\Omega}\bm{t}N^{a}(\bm{x})\
\mathrm{d}S}_{\bm{F}_{\text{ext}}^{a}}\right]=0,$ (9)

where the first and second integral are the nodal internal forces
$\bm{F}_{\text{int}}^{a}$ and nodal external forces $\bm{F}_{\text{ext}}^{a}$,
respectively. By employing the constitutive relation Eq. 3, the nodal internal
forces may be written as

$\displaystyle\bm{F}_{\text{int}}^{a}$
$\displaystyle=\int_{\Omega}\bm{\sigma}\nabla N^{a}\ \mathrm{d}A,$
$\displaystyle=\int_{\Omega}\begin{bmatrix}\sigma_{{1}}N_{,x}^{a}+\sigma_{{6}}N_{,y}^{a}\\\
\sigma_{{6}}N_{,x}^{a}+\sigma_{{2}}N_{,y}^{a}\\\ \end{bmatrix}\ \mathrm{d}A,$
$\displaystyle=\int_{\Omega}\begin{bmatrix}C_{11}\varepsilon_{{1}}N_{,x}^{a}+C_{12}\varepsilon_{{2}}N_{,x}^{a}+2C_{16}\varepsilon_{{6}}N_{,x}^{a}+C_{16}\varepsilon_{{1}}N_{,y}^{a}+C_{26}\varepsilon_{{2}}N_{,y}^{a}+2C_{66}\varepsilon_{{6}}N_{,y}^{a}\\\
C_{16}\varepsilon_{{1}}N_{,x}^{a}+C_{26}\varepsilon_{{2}}N_{,x}^{a}+2C_{66}\varepsilon_{{6}}N_{,x}^{a}+C_{12}\varepsilon_{{1}}N_{,y}^{a}+C_{22}\varepsilon_{{2}}N_{,y}^{a}+2C_{26}\varepsilon_{{6}}N_{,y}^{a}\\\
\end{bmatrix}\ \mathrm{d}A,$
$\displaystyle=\int_{\Omega}\begin{bmatrix}\varepsilon_{{1}}N_{,x}^{a}&\varepsilon_{{2}}N_{,x}^{a}&0&2\varepsilon_{{6}}N_{,x}^{a}+\varepsilon_{{1}}N_{,y}^{a}&\varepsilon_{{2}}N_{,y}^{a}&2\varepsilon_{{6}}N_{,y}^{a}\\\
0&\varepsilon_{{1}}N_{,y}^{a}&\varepsilon_{{2}}N_{,y}^{a}&\varepsilon_{{1}}N_{,x}^{a}&\varepsilon_{{2}}N_{,x}^{a}+2\varepsilon_{{6}}N_{,y}^{a}&2\varepsilon_{{6}}N_{,x}^{a}\\\
\end{bmatrix}\ \mathrm{d}A\ \bm{C}_{\text{vec}},$ (10)

where the elasticity tensor parameters $\bm{C}_{\text{vec}}=[C_{11}\ C_{12}\
C_{22}\ C_{16}\ C_{26}\ C_{66}]^{T}$ are assumed to be constant in space.

###  3.5 Choice of test functions

Choosing a test function in the form of (8) and evaluating (9) results in two
linear equations with the material parameters as unknowns. As the weak linear
momentum balance has to hold true for any test function this provides an
infinite supply of linear equations. Hence, the problem at hand is
overdetermined and different choices of test functions will yield different
solutions for the unknown material parameters.

As discussed in Section 3.1, the deformation of a heterogeneous specimen and
that of its homogenized counterpart under the same loading conditions are
locally different, a phenomenon that is best observed at the boundary and at
the corners of the specimen where local effects are especially pronounced. In
the following, this special characteristic of the problem at hand motivates a
special choice of the test functions that avoids evaluations of the linear
momentum balance in the boundary regions of the specimen.

First, we define test functions that are constant at the nodes corresponding
to one finite element, i.e., one unit cell, and zero at all other nodes. To
this end, we define $\mathcal{C}=\\{1,\dots,n_{c}^{2}\\}$ as the set of all
unit cells and $\mathcal{D}^{c}$ as the set of all nodes corresponding to the
unit cell $c\in\mathcal{C}$, and define a set of test functions as

$\displaystyle\mathcal{V}=\left\\{\bm{v}(\bm{x})=\frac{1}{n_{nc}}\sum_{a\in\mathcal{D}^{c}}N^{a}(\bm{x})\bm{e}_{i}\
|\ c\in\mathcal{C},\ i\in\\{1,2\\}\right\\},$ (11)

where $\bm{e}_{i}$ are the unit vectors in the corresponding $x$\- and
$y$-direction. Note that the test functions are normalized by dividing by the
number of nodes corresponding to the unit cell $n_{nc}$ (equal to 4 in our
case).

Using the test functions in $\mathcal{V}$ to test weak linear momentum balance
would cause two problems. First, at elements adjacent to the loaded and to the
restrained portions of the boundary, the external force contributions
$\bm{F}_{\text{ext}}^{a}$ in (9) are unknown, leading to equations that could
not be solved for the unknown material parameters. And second, we want to
avoid using data at the specimen boundary due to the reasons discussed
earlier. Therefore, we modify (11) such that

$\displaystyle\mathcal{V}^{\text{int}}=\left\\{\bm{v}(\bm{x})=\frac{1}{n_{nc}}\sum_{a\in\mathcal{D}^{c}}N^{a}(\bm{x})\bm{e}_{i}\
|\ c\in\mathcal{C}^{\text{int}},\ i\in\\{1,2\\}\right\\},$ (12)

where $\mathcal{C}^{\text{int}}\subset\mathcal{C}$ denotes a reduced set of
unit cells that does not include unit cells close to the boundary. We found
that ignoring two rows of unit cells at the top and bottom boundary as well as
two columns of unit cells at the left and right boundary are a good
compromise, and we kept this choice constant throughout all tests. As the
fields in $\mathcal{V}^{\text{int}}$ depend on $\bm{e}_{i}$, each field is
zero in either $x$\- or $y$-direction. The non-zero component of an exemplary
virtual field in $\mathcal{V}^{\text{int}}$ is shown in Fig. 2 (left).

Evaluating Eq. 9 for this set of functions leads to

$\displaystyle\frac{1}{n_{nc}}\sum_{a\in\mathcal{D}^{c}}\bm{F}_{\text{int}}^{a}=\bm{0},\quad\forall
c\in\mathcal{C}^{\text{int}}.$ (13)

Hence, this choice of virtual fields can be interpreted physically as
enforcing that the sum of internal forces over one unit cell should vanish.

Figure 2: Non-zero component of a virtual field in $\mathcal{V}^{\text{int}}$
(left) and non-zero component of a virtual field in
$\mathcal{V}^{\text{center}}$ (right).

Equations (13) are not sufficient to identify the unknown material parameters,
as the trivial solution $\bm{C}_{\text{vec}}=\bm{0}$ fulfills (13). To obtain
a well-posed problem, the measured reaction forces need to be incorporated. At
the same time, we want to avoid using displacement data at the specimen
boundary. Therefore, we consider the free-body diagram of the lower half of
the domain as depicted in Fig. 2 (right). Denoting the half-body domain as
$\Omega^{*}=\\{\bm{x}\ |\ 0\leq x_{1}\leq L,\ \ 0\leq x_{2}\leq\frac{L}{2}\\}$
and its boundary as $\partial\Omega^{*}$, the weak form of linear momentum
balance for this domain reads

$\displaystyle\int_{\Omega^{*}}\bm{\sigma}(\bm{x})\colon\nabla\bm{v}(\bm{x})\
\mathrm{d}A-\int_{\partial\Omega^{*}}\bm{t}\cdot\bm{v}(\bm{x})\
\mathrm{d}s=0.$ (14)

Inserting the test function ansatz leads to

$\displaystyle\sum_{a=1}^{n_{n}}\bm{v}^{a}\cdot\left[\underbrace{\int_{\Omega^{*}}\bm{\sigma}\nabla
N^{a}(\bm{x})\
\mathrm{d}A}_{\bm{F}_{\text{int}}^{*a}}-\underbrace{\int_{\partial\Omega^{*}}\bm{t}N^{a}(\bm{x})\
\mathrm{d}S}_{\bm{F}_{\text{ext}}^{*a}}\right]=0.$ (15)

We define $\mathcal{D}^{\text{center}}=\\{a\ |\ y^{a}=\frac{L}{2}\\}$ as the
set of nodes in the center of the specimen. If the tessellated geometry
consists of an odd number of unit cells in each spatial direction, i.e., there
are no nodes at $y^{a}=\frac{L}{2}$, we consider instead
$\mathcal{D}^{\text{center}}=\\{a\ |\ y^{a}=\frac{L}{2}+\frac{L}{2n_{c}}\\}$.
We choose a set of virtual fields $\mathcal{V}^{\text{center}}$ that are
constant along $\mathcal{D}^{\text{center}}$ and zero at all other nodes

$\displaystyle\mathcal{V}^{\text{center}}=\left\\{\bm{v}(\bm{x})=\frac{1}{n_{nc}}\sum_{a\in\mathcal{D}^{\text{center}}}N^{a}(\bm{x})\bm{e}_{i}\
|\ i\in\\{1,2\\}\right\\}.$ (16)

Evaluating (15) for these particularly chosen test functions results in

$\displaystyle\sum_{a\in\mathcal{D}^{\text{center}}}\bm{F}_{\text{int}}^{*a}=\int_{\partial\Omega^{\text{center}}}\bm{t}\
\mathrm{d}S=\bm{R},$ (17)

where $\partial\Omega^{\text{center}}$ is the top boundary of $\Omega^{*}$.
Note that due to the specific choice of the test functions, the surface
integral simplifies in such a way that it equals the global reaction force
$\bm{R}$, meaning that the sum of the internal forces at
$\partial\Omega^{\text{center}}$ must equal the net reaction force.

###  3.6 Deterministic parameter identification

After choosing the virtual fields and considering (3.4), the linear equations
in (13) can be assembled in a system of equations

$\bm{A}^{\text{int}}\bm{C}_{\text{vec}}=\bm{0},$ (18)

and the linear equations in (17) can be rewritten as

$\bm{A}^{\text{center}}\bm{C}_{\text{vec}}=\bm{R},$ (19)

where $\bm{A}^{\text{int}}$ and $\bm{A}^{\text{center}}$ are in general non-
symmetric matrices. The system formed by the linear equations (18) and (19) is
overdetermined, i.e., it consists of more equations than unknown parameters.
Assuming that the equations in the overdetermined system are not linearly
dependent (which is a valid assumption as every equation is perturbed by noise
when considering experimental data), there is no unique solution that
satisfies all equations. Instead, we obtain an approximate solution of the
overdetermined system by minimizing the sum of squared residuals

$\bm{C}_{\text{vec}}^{\text{opt}}=\operatornamewithlimits{arg\
min}_{\bm{C}_{\text{vec}}}\left(\|\bm{A}^{\text{int}}\bm{C}_{\text{vec}}\|^{2}+\lambda_{r}\|\bm{A}^{\text{center}}\bm{C}_{\text{vec}}-\bm{R}\|^{2}\right),$
(20)

where $\|\cdot\|$ is the Euclidean norm and $\lambda_{r}>0$ is a weighting
parameter that scales the different contributions to the minimization problem.
As there are less equations in the system (19) than in (18), the weighting
parameter should be chosen sufficiently larger than one ($\lambda_{r}>>1$).
Following previous works Flaschel et al. (2021, 2022, 2023), we choose
$\lambda_{r}=100$ and keep it constant throughout this work. Based on our
experience, the choice of $\lambda_{r}$ is not crucial for the success of the
method (see also Joshi et al. (2022); Thakolkaran et al. (2022); Marino et al.
(2023)). The necessary condition for a minimum is

$\bar{\bm{A}}\bm{C}_{\text{vec}}^{\text{opt}}=\bar{\bm{R}},\qquad\text{with}\qquad\bar{\bm{A}}=\left(\bm{A}^{\text{int}}\right)^{T}\bm{A}^{\text{int}}+\lambda_{r}\left(\bm{A}^{\text{center}}\right)^{T}\bm{A}^{\text{center}},~{}\bar{\bm{R}}=\lambda_{r}\left(\bm{A}^{\text{center}}\right)^{T}\bm{R},$
(21)

which leads to a determined system of equations that can be solved for
$\bm{C}_{\text{vec}}^{\text{opt}}$. The minimization problem in Eq. 20 can
alternatively be written as

$\bm{C}_{\text{vec}}^{\text{opt}}=\operatornamewithlimits{arg\
min}_{\bm{C}_{\text{vec}}}\|\bm{A}\bm{C}_{\text{vec}}-\bm{B}\|^{2},$ (22)

where we have defined

$\bm{A}=\begin{bmatrix}\bm{A}^{\text{int}}\\\
\sqrt{\lambda_{r}}\bm{A}^{\text{center}}\\\
\end{bmatrix},\qquad\bm{B}=\begin{bmatrix}\bm{0}\\\
\sqrt{\lambda_{r}}\bm{R}\\\ \end{bmatrix}.$ (23)

The necessary condition for a minimum then reads

$\bm{A}^{T}\bm{A}\bm{C}_{\text{vec}}^{\text{opt}}=\bm{A}^{T}\bm{B}.$ (24)

###  3.7 Bayesian inference

Besides the previously introduced deterministic approach, we further study the
problem from a stochastic perspective. To this end, we construct a Bayesian
linear regression model, for which we assume no intercept and a diffuse prior,
as implemented in the Matlab ® built-in function bayeslm.

We denote the number of rows in $\bm{A}$ as $n_{eq}$ and we define
$\bm{A}_{i}$ with $i\in\\{1,\dots,n_{eq}\\}$ as the $i^{\text{th}}$ row of
$\bm{A}$. For each equation in the overdetermined system of equations
$\bm{A}\bm{C}_{\text{vec}}=\bm{B}$, we assume the likelihood of obtaining
$B_{i}$ as a Gaussian likelihood with mean
$\bm{A}_{i}\cdot\bm{C}_{\text{vec}}$ and standard deviation $\sigma>0$, i.e.,

$p(B_{i}|\bm{A}_{i},\bm{C}_{\text{vec}},\sigma^{2})=\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left[-\frac{\left(B_{i}-\bm{A}_{i}\cdot\bm{C}_{\text{vec}}\right)^{2}}{2\sigma^{2}}\right],$
(25)

where $\bm{C}_{\text{vec}}$ and $\sigma^{2}$ are treated as random variables.
Assuming further that the likelihoods are conditionally independent, we define
the joint likelihood as

$p(\bm{B}|\bm{A},\bm{C}_{\text{vec}},\sigma^{2})=\prod_{i=1}^{n_{eq}}p_{i}(B_{i}|\bm{A}_{i},\bm{C}_{\text{vec}},\sigma^{2}).$
(26)

Assuming here a diffuse prior for the joint prior distribution of
$\bm{C}_{\text{vec}}$ and $\sigma^{2}$, i.e.,

$p(\bm{C}_{\text{vec}},\sigma^{2})\propto\frac{1}{\sigma^{2}},$ (27)

the marginal posterior distributions of $\bm{C}_{\text{vec}}$ and $\sigma^{2}$
are analytically tractable and implemented in the Matlab ® function bayeslm.

##  4 Data acquisition

In this section, we first discuss the unit cell geometries considered for
identification of the material parameters. Then, we describe our numerical and
experimental data acquisition methods, including details on fabrication,
experimental setup, testing and DIC.

###  4.1 Design and choice of unit cell geometries

We pick four unit cells with distinct/diverse effective stiffness tensor
parameters (all with six non-zero stiffness parameters). Table 1 shows the
unit cells along with their symmetry class and homogenized stiffness tensor.
Geometry #1 has $C_{22}$ as the largest stiffness parameter with $C_{16}$
almost comparable to $C_{12}$ and $C_{26}>C_{16}$. While geometry #2 has
$C_{11}$ as the largest stiffness parameter with $C_{16}>C_{26}$, geometry #3
has negative values for all of the off-diagonal parameters. Geometry #4 has
four independent stiffness parameters with $C_{66}$ as one of the largest
values among other stiffness parameters, along with $C_{11}=C_{22}$ and
$C_{16}=C_{26}$. The fill fraction of the stiff phase for all the unit cells
lies between 60 and 70 %.

Unit Cell Geometry |  Name |  Homogenized Stiffness Tensor $(\bm{C}^{\text{H}})$ [MPa] |  Elastic Symmetry Class
---|---|---|---
|  Geometry #1 |  ${\left[\begin{array}[]{rrr}131.62&61.98&63.58\\\ 61.98&198.38&83.87\\\ 63.58&83.87&95.30\end{array}\right]}$ |  $Z_{2}$
|  Geometry #2 |  ${\left[\begin{array}[]{rrr}127.14&59.50&73.42\\\ 59.50&105.83&55.16\\\ 73.42&55.16&110.15\end{array}\right]}$ |  $Z_{2}$
|  Geometry #3 |  ${\left[\begin{array}[]{rrr}44.70&-9.42&-12.52\\\ -9.42&107.19&-20.71\\\ -12.52&-20.71&105.35\end{array}\right]}$ |  $Z_{2}$
|  Geometry #4 |  ${\left[\begin{array}[]{rrr}65.74&40.36&18.95\\\ 40.36&65.74&18.95\\\ 18.95&18.95&86.47\end{array}\right]}$ |  $D_{2}$
Table 1: Unit cell geometries considered in this study along with their
mechanical and symmetry properties.

###  4.2 Numerical data generation

Simulations: We use synthetic data generated using the FEM to verify our
methods and aid our analysis before performing the experiments. Each pixel is
discretized using a four-node plane-stress bilinear quadrilateral element. For
tessellation, we vary the number of unit cells $n_{c}$ between 5 and 25.

Homogenization: We compute the effective mechanical properties of the unit
cells using the theory of homogenization implemented using the FEM (as in
Andreassen and Andreasen (2014)).

###  4.3 Experimental data generation

####  4.3.1 Fabrication

As specimens with a large number of unit cells are difficult to fabricate, we
here pick 10 $\times$ 10 tessellations to perform experimental validations. We
use a commercial multi-material polyjet technology based 3D printer, Stratasys
Objet500 Connex, to fabricate all the specimens. The dimensions of the
specimen are 75 × 75 × 5 mm excluding the portion that goes into the grips. We
use Stratasys’ proprietary material DM8530 for the stiff phase and TangoBlack
for the soft phase. The material properties (DM8530: Young’s modulus E = 1000
$\pm$ 90 MPa and Poisson’s ratio $\nu$= 0.35, TangoBlack: Young’s modulus E =
0.7 MPa and Poisson’s ratio $\nu$= 0.49) are experimentally measured following
the ASTM D638-14 standard test method and the same values are used in the
numerical computations.

####  4.3.2 Experimental setup and testing

We subject the additively manufactured specimens to displacement-controlled
tension tests using a universal testing machine, Instron E3000, mounted with a
multi-axis force-torque sensor (ATI Mini85) as shown in Fig. 3. The force-
torque sensor is acquired from ATI Industrial Automation. We apply a vertical
displacement of 1.5 mm at the top boundary at a rate of 0.5 mm/min resulting
in a global axial strain of $\tilde{\varepsilon}_{22}=0.02$ and a global
strain rate of $1.1\times 10^{-4}$s-1. Custom designed grips are fabricated
out of aluminum and are serrated to hold the specimens firmly and prevent any
lateral slipping. We use the same strain rate while measuring the constitutive
material properties of the individual phases.

We use DIC, an image-based optical technique, to measure the full-field
displacements (Sutton et al., 2009). We capture images at a frequency of 1 Hz
using a Nikon D750 camera equipped with a Nikon AF-S NIKKOR 24-120mm f/4G ED
VR zoom lens. We use manual mode at an exposure rate of 1/640 sec, an ISO
setting of 1250 and an aperture setting of F8. The camera has a 6016 by 4016
square pixel resolution and the region of interest we studied is about 3060 by
3060 pixels. We place a ring light between the additively manufactured
specimen and the camera to illuminate the surface uniformly and we place the
camera lens at a distance of about 35-40 cm from the specimen plane.

Figure 3: Experimental setup for displacement-controlled uniaxial testing of
an anisotropic metamaterial.

####  4.3.3 Digital image correlation

Given a reference image $f$ and and a deformed image $g$, the correlation
algorithm aims at minimizing the sum of squared differences over the
considered domain $\Omega$

$\mathcal{T}=\int_{\Omega}\left(g(\bm{x}+\bm{u}(\bm{x}))-f(\bm{x})\right)^{2}d\bm{x},$
(28)

where $\bm{x}$ is the position in the reference image and $\bm{u}(\bm{x})$ is
the displacement field which is interpolated as

$\bm{u}(\bm{x})=\sum u_{n}\bm{\phi}_{n}(\bm{x}),$ (29)

where $\bm{\phi}_{n}$ are a set of shape functions and $u_{n}$ the associated
degrees of freedom. There are two approaches to determine the unknowns
$u_{n}$, local DIC and global DIC (Hild and Roux, 2012). In the local
approach, the region of interest ($\Omega$) is divided into several sub-images
known as subsets and the mean displacement of each subset is computed
independently while minimizing the objective Eq. 28. In the global approach,
shape functions defined through a finite element mesh over the whole region of
interest are used (Besnard et al., 2006). The global approach assumes
continuity of displacements over the entire region of interest which is well
suited when the structure is heterogeneous. Moreover, the global approach
provides the displacement information at the boundaries, which is hard to
obtain using the local approach. The displacement data at the boundaries are
an important input for the VFM. Hence, we follow the global approach to
perform the correlation in this study.

We perform DIC using piece-wise linear shape functions defined on a triangular
mesh to compute the displacements (as in Agnelli et al. (2021)). We choose an
edge length of 18 pixels ($\sim$ 0.44 mm) to construct the triangular mesh. We
observe a noise floor of the order of 0.04 mm in the displacement data which
is obtained from correlation performed on static images. The data provided by
the DIC correspond to the nodes that might not always align with the unit cell
corners. To obtain the displacements of the unit cell corners, we further
average the displacement data from the nodes that fall within 1 mm radius of a
unit cell corner.

##  5 Results and discussion

In this section, we discuss the data generation from both numerical
simulations and experiments. Afterwards, we apply the proposed deterministic
parameter identification method to the data and discuss the results. Finally,
at the end of the section, we apply the Bayesian method to the data.

###  5.1 Generation of full-field displacement data

####  5.1.1 Synthetic data

In the following, we investigate the synthetically generated displacement data
for a heterogeneous structure in comparison to the computed displacement field
of a homogeneous body, whose stiffness is equal to the homogenized stiffness
of the heterogeneous structure. To simulate the displacement of a homogeneous
body, we assume a $10\times 10$ bilinear quadrilateral finite element mesh.
The displacement of the heterogeneous body is computed on a much finer mesh
with $1000\times 1000$ elements. To allow for a comparison with the
displacement field of the homogeneous body, the computed displacements at the
unit cell corners of the heterogeneous body (i.e. the data of interest for the
VFM) are extracted and interpolated with a bilinear polynomial for each unit
cell. It can be seen in Fig. 4 that there is a good qualitative agreement
between the two displacement fields for geometry #1 (see Fig. S1, Fig. S2 and
Fig. S3 for the other geometries). However, there are quantitative differences
due to local effects in the heterogeneous structure, which appear to be
dominant at the boundary and corners of the specimen.

Figure 4: Comparison between the displacement fields obtained from finite
element simulations of a homogeneous specimen (left) and a heterogeneous
structure made of geometry #1 (center). For the homogeneous specimen a finite
element simulation using $10\times 10$ bilinear quadrilateral elements was
executed. The heterogeneous specimen was simulated using $1000\times 1000$
bilinear quadrilateral elements. Afterwards, the displacement data at the unit
cell corners were extracted and interpolated with a bilinear polynomial for
each unit cell, to allow for a comparison with the homogeneous specimen. The
difference between the fields is shown on the right.

####  5.1.2 Comparison between experimental and synthetic data

In Fig. 5, we compare the full-field displacement and strain fields between
the numerical and experimental data on the heterogeneous structure for
geometry #1. (See Figs. S4, S5 and S6 for the other geometries). We observe
very good agreement between the numerical and experimental data, especially
for the variables $u_{2},\varepsilon_{22}$. However, the experimentally
measured $u_{1}$ appears to be slightly higher than the numerical data, by
about 0.1 mm, for all the geometries. Also the two $\varepsilon_{11}$ fields
are in good qualitative agreement, but experimental strains are larger. As
expected, most of the strain is localized in the softer phase, although the
applied global strain ($\tilde{\varepsilon}_{22}$) is 0.02.

Figure 5: Comparison between numerical (top) and experimentally measured
(bottom) full-field displacement and strain field data for the 10 unit cell
tessellation of geometry #1 subjected to displacement-controlled uniaxial
tension test.

Further, a comparison of the displacement fields after postprocessing the
synthetic and experimental data, i.e., after extracting and interpolating the
displacements at the unit cell corners for all the geometries for geometry #2
are shown in Fig. S7, Fig. S8, Fig. S9, and Fig. S10. All displacements are in
good agreement. An exception is observed for geometry #3 (see Fig. S5 and Fig.
S9), for which the experimentally measured horizontal displacement $u_{1}$
does not compare well to the corresponding finite element results. The unit
cell architecture of geometry #3 leads to highly nonlinear mechanical behavior
(see Fig. 10), which is not captured well in the simulations.

###  5.2 Parameter identification based on synthetic data

Since the homogenization theory assumes length scale separation and periodic
boundary conditions in identifying the effective material parameters, it is
important to understand the continuum behavior of the heterogeneous structures
as the number of unit cells change. For this, we apply the VFM described in
Section 3.6 on the synthetic data to identify material parameters as the
number of unit cells are varied. Further, we also use synthetic data to
identify parameters using multiple tests (as in the conventional approach). A
discussion on this conventional approach is provided in A. The relative error
is defined as

$\text{LSE}_{\|\cdot\|_{2}}=\frac{\left\|\bm{C}^{\mathrm{H}}_{\text{vec}}-\bm{C}^{\mathrm{M}}_{\text{vec}}\right\|_{2}}{\left\|\bm{C}^{\mathrm{H}}_{\text{vec}}\right\|_{2}}\quad\text{
with }\bm{C}_{\text{vec}}\in\mathbb{R}^{6},$ (30)

where $\bm{C}^{\mathrm{H}}_{\text{vec}}$ is the vectorized homogenized
stiffness tensor obtained from computational homogenization and
$\bm{C}^{\mathrm{M}}_{\text{vec}}$ is the vectorized stiffness tensor
identified using the VFM and the conventional methods.

We compare the relative error in parameter identification when performing
multiple tests (as in the conventional approach) and when using the VFM (Fig.
6). Since we exclude two rows and columns of boundary unit cells in the
proposed VFM, the number of unit cells available to form the system of
equations is guaranteed only when there are at least 7 unit cells and the
results are shown starting with this number. For geometry #1, as the number of
unit cells increases, the error calculated for the conventional method based
on multiple tests decreases monotonically from 13.4% at 5 unit cell
tessellation to 2.3% at 25 unit cell tessellation. Similarly, the error for
the VFM decreases monotonically from 13.1% at 7 unit cell tessellation to 2.4%
at 25 unit cell tessellation. This shows that the parameters identified using
our VFM are as good as those obtained by performing multiple tests, as long as
there are at least ten repeated unit cells in the domain of interest. We found
this general conclusion to hold for most of the considered geometries. The
only exception is geometry #3, for which the error remains at 7.0% (for
multiple tests) and 11.8% (for the VFM) after 10 unit cell tessellation. In
the case of geometry #3, a major portion of the error lies in just two of the
parameters $C_{12}$ and $C_{16}$. These two parameters are quite small
relative to the rest of the parameters and hence, they are hard to accurately
estimate in comparison to the others.

Figure 6: Variation of least square error between homogenized stiffness tensor
and stiffness tensor identified using the VFM and the conventional methods as
the number of unit cells are varied.

As the number of unit cells increases, the ratio of the number of unit cells
along the boundary to the number of unit cells in the interior decreases. As a
result, the boundary effects described in Section 5.1.1 diminish and the
behavior of the structure approaches the continuum equivalent. In Table 2, we
summarize the parameters identified for geometry #1 from both the methods
against homogenization for 25 unit cell tessellation (see Table 3, Table 4,
Table 5 for the other geometries).

Method |  $C_{11}$ (MPa) |  $C_{12}$ (MPa) |  $C_{22}$ (MPa) |  $C_{16}$ (MPa) |  $C_{26}$ (MPa) |  $C_{66}$ (MPa)
---|---|---|---|---|---|---
Homogenization | 131.62 | 61.98 | 198.38 | 63.58 | 83.87 | 95.30
VFM | 125.97 | 62.57 | 196.15 | 61.35 | 82.39 | 93.48
Multiple tests | 129.86 | 67.93 | 199.37 | 64.01 | 85.47 | 95.54
Table 2: Comparison of stiffness tensor parameters identified for geometry #1
with 25 unit cell tessellation based on synthetic data using the VFM and the
conventional methods against the computational homogenization.

###  5.3 Parameter identification based on experimental data

Fig. 7 summarizes the material parameters identified by the VFM using the
simulated and experimental data for 10 unit cell tessellations in comparison
to the homogenized stiffness. The parameters identified using synthetic data
compare well with the homogenized properties for all the geometries. Further,
a good qualitative agreement is observed for the parameters identified using
experimental data. For some of the parameters, such as $C_{12},C_{22}$,
$C_{16}$, $C_{26}$, the experimentally determined parameters match the
expectations quantitatively. In contrast, there is a larger discrepancy in the
values of $C_{11},C_{66}$, for almost all the geometries. These discrepancies
are related to the fact that the experimentally measured displacement $u_{1}$
appears higher than in the simulations, i.e., about 0.1 mm, leading to an
under-prediction of the stiffness in the lateral directions. An interesting
observation is made for geometry #4. Based on the numerical data, we know that
$C_{11}=C_{22}$ and $C_{16}=C_{26}$. However, we observe experimentally that
$C_{11}<C_{22}$ and $C_{16}<C_{26}$. In geometry #3, the discrepancy may be
caused by the architecture itself. The structure has thin and sharp features
in the soft phase. As it is well known, the behavior of materials in the
vicinity of such sharp discontinuities is quite different from a linear
elastic continuum (Rosakis and Ravi-Chandar (1986)). In fact, the load-
displacement data for geometry #3 show nonlinear behavior (see Fig. 10). In
such micro-structures, our linear elastic model assumption fails.

Figure 7: Comparison of material parameters identified using the VFM from
numerical and experimental data of 10 unit cell tessellations.

We finalize the study by applying the Bayesian method described in Section 3.7
to the experimental data. The resulting marginal posterior probability
distributions of the material parameters are shown in Fig. 8. It is observed
that the computed mean values of the marginal posteriors are similar to the
deterministic results shown in Fig. 7. Beyond that, the standard deviation of
the marginal posteriors indicate (un)certainty in the parameter predictions.
Matching our expectations, the parameter $C_{22}$ is identified with the
highest certainty, while for example the identification of the parameter
$C_{11}$ shows a high uncertainty. Additionally, it is noteworthy that the
marginal posteriors of the parameters identified when the Bayesian method is
applied to the numerical data show low standard deviations as the data is not
affected by the experimental noise (see Fig. 11).

Figure 8: Marginal posterior probability distributions of the material
parameters obtained through Bayesian linear regression on the experimental
data. The red lines indicate the mean of the marginal posterior distributions.
The blue boxes indicate the standard deviation from the mean, i.e., the $68\%$
probability interval. The black intervals indicate three times the standard
deviation from the mean, i.e., the $99\%$ probability interval.

We note that, for geometry #3, the marginal posterior probability
distributions of the parameters exceed the thermodynamically admissible range,
e.g., the marginal posterior of $C_{11}$ is partially negative. This must be
considered when interpreting the results. In this work, no measure was taken
to enforce thermodynamic admissibility in the Bayesian method, which thus
remains a future objective.

##  6 Conclusions

In this paper, we present an approach to identify the 6 independent elastic
material parameters of plane anisotropic elasticity from a single experiment,
using the virtual fields method. This approach allows identifying shear-normal
coupling parameters experimentally, a task that had remained challenging so
far. We first demonstrate the effectiveness of our method using numerically
generated data from a single tension test. We then experimentally validate the
method on additively manufactured specimens, by measuring full-field
displacement data and traction forces. We show that our method is effective
for materials that include at least 10 repeated unit cells in their structure,
to satisfy homogenenization conditions. We calculate the uncertainity in the
identification estimation of the material parameters using Bayesian linear
regression. In the future, to further refine the experimental parameter
identification, it is necessary to optimize the shape of the specimens to
ensure strong contributions of strains from different stiffness tensor
components. The proposed approach has potential for measurements of elasticity
parameters of complex, anisotropic, three-dimensional structured materials and
composites with shear-shear couplings, and for the study of their nonlinear
behavior. A further potential application of the method could be for parameter
identification of constitutive tensors corresponding to different types of
coupled behavior, such as generalized piezoelectric, flexoelectric and
piezomagnetic tensors.

## CRediT authorship contribution statement

Jagannadh Boddapati: Conceptualization, Investigation, Software, Formal
analysis, Writing - original Draft. Moritz Flaschel: Conceptualization,
Investigation, Software, Formal analysis, Writing - original Draft. Siddhant
Kumar:. Conceptualization, Writing - review and editing. Laura De Lorenzis:
Conceptualization, Writing - review and editing. Chiara Daraio:.
Conceptualization, Writing - review and editing, Supervision, Funding
acquisition.

## Acknowledgements

We thank Pierre Margerit (École Polytechnique) for the discussion on digital
image correlation, and Jihoon Ahn (Caltech) and Perry Samimy (Caltech) for
their help on designing the experiments. C.D. and J.B. acknowledge support
from the US National Science Foundation (NSF), grant number 1835735. M.F. and
L.D.L. acknowledge support from the Swiss National Science Foundation (SNF),
project number 200021_204316.

##  Appendix A Parameter identification based on multiple tests

In this section, we explore a method of parameter identification that involves
multiple tests (as in the conventional approach) in the context of anisotropic
metamaterials. We subject the metamaterial to three different tests namely
Test A, Test B, and Test C as shown in Fig. 9. Test A and Test C are tension
tests along $x_{2}$ and $x_{1}$ axis respectively, and Test B is a simple
shear test. We assume that the average strains
$\tilde{\bm{\varepsilon}}_{ij}^{A,B,C}$ are known experimentally from full-
field measurements. In addition, the reaction forces at the fixed end are
known experimentally from load sensor measurements. We will show that the
material parameters can be identified from the average strains and the net
reaction forces from these three tests.

From Gauss’ divergence theorem, the average stresses $\tilde{\bm{\sigma}}$ are
related to the tractions $\bm{t}$ at the fixed end as

$t_{i}=\tilde{\sigma}_{ij}n_{j},$ (31)

where $\bm{n}$ is the unit outward normal. For Test A, the unit outward normal
$\bm{n}$ at the fixed end is $[0,-1]^{T}$. Using Eqs. 31 and 3, and assuming
homogenized effective continuum behavior for the structured solid, we get

$\displaystyle F_{1}^{A}/\mathcal{A}$
$\displaystyle=\tilde{\sigma}_{6}^{A}=C_{16}\tilde{\varepsilon}_{11}^{A}+C_{26}\tilde{\varepsilon}_{2}^{A}+C_{66}\left(2\tilde{\varepsilon}_{6}^{A}\right),$
(32a) $\displaystyle F_{2}^{A}/\mathcal{A}$
$\displaystyle=\tilde{\sigma}_{2}^{A}=C_{12}\tilde{\varepsilon}_{1}^{A}+C_{22}\tilde{\varepsilon}_{2}^{A}+C_{26}\left(2\tilde{\varepsilon}_{6}^{A}\right),$
(32b)

where $\tilde{\sigma}_{12}^{A},\tilde{\sigma}_{22}^{A}$ are the average stress
components, $F_{1}^{A},F_{2}^{A}$ are the reaction force components at the
fixed end from Test A and $\mathcal{A}$ is the cross sectional area of the
fixed end.

Similarly, from Test B and Test C, we get

$\displaystyle F_{1}^{B}/\mathcal{A}$
$\displaystyle=\tilde{\sigma}_{6}^{B}=C_{16}\tilde{\varepsilon}_{1}^{B}+C_{26}\tilde{\varepsilon}_{2}^{B}+C_{66}\left(2\tilde{\varepsilon}_{6}^{B}\right),$
(33a) $\displaystyle F_{2}^{B}/\mathcal{A}$
$\displaystyle=\tilde{\sigma}_{2}^{B}=C_{12}\tilde{\varepsilon}_{1}^{B}+C_{22}\tilde{\varepsilon}_{2}^{B}+C_{26}\left(2\tilde{\varepsilon}_{6}^{B}\right),$
(33b) $\displaystyle F_{1}^{C}/\mathcal{A}$
$\displaystyle=\tilde{\sigma}_{1}^{C}=C_{11}\tilde{\varepsilon}_{1}^{C}+C_{12}\tilde{\varepsilon}_{2}^{C}+C_{16}\left(2\tilde{\varepsilon}_{6}^{C}\right),$
(33c) $\displaystyle F_{2}^{C}/\mathcal{A}$
$\displaystyle=\tilde{\sigma}_{6}^{C}=C_{16}\tilde{\varepsilon}_{1}^{C}+C_{26}\tilde{\varepsilon}_{2}^{C}+C_{66}\left(2\tilde{\varepsilon}_{6}^{C}\right).$
(33d)

Rearranging Eqs. 32a, 32b, 33a, 33b, 33c and 33d into a matrix form, we obtain
a system of linear equations,

$\left[\begin{array}[]{cccccc}0&0&2\tilde{\varepsilon}_{6}^{A}&\tilde{\varepsilon}_{1}^{A}&\tilde{\varepsilon}_{2}^{A}&0\\\
0&\tilde{\varepsilon}_{2}^{A}&0&0&2\tilde{\varepsilon}_{6}^{A}&\tilde{\varepsilon}_{1}^{A}\\\
0&0&2\tilde{\varepsilon}_{6}^{B}&\tilde{\varepsilon}_{1}^{B}&\tilde{\varepsilon}_{2}^{B}&0\\\
0&\tilde{\varepsilon}_{2}^{B}&0&0&2\tilde{\varepsilon}_{6}^{B}&\tilde{\varepsilon}_{1}^{B}\\\
\tilde{\varepsilon}_{1}^{C}&0&0&2\tilde{\varepsilon}_{6}^{C}&0&\tilde{\varepsilon}_{2}^{C}\\\
0&0&2\tilde{\varepsilon}_{6}^{C}&\tilde{\varepsilon}_{1}^{C}&\tilde{\varepsilon}_{2}&0\end{array}\right]\left[\begin{array}[]{l}C_{11}\\\
C_{22}\\\ C_{66}\\\ C_{16}\\\ C_{26}\\\
C_{12}\end{array}\right]=\frac{1}{\mathcal{A}}\left[\begin{array}[]{c}F_{1}^{A}\\\
F_{2}^{A}\\\ F_{1}^{B}\\\ F_{2}^{B}\\\ F_{1}^{C}\\\
F_{2}^{C}\end{array}\right].$ (34)

For readability, Eq. 34 is written as

$\mathcal{A}\tilde{\bm{\varepsilon}}\bm{C}_{\text{vec}}=\bm{F}_{\text{vec}},$
(35)

where $\tilde{\bm{\varepsilon}}$ is a non-symmetric square matrix of size 6
containing average strain components from all of the tests and
$\bm{F}_{\text{vec}}$ is a vector containing net reaction force components
from all of the tests. Then the material parameters
$\bm{C}^{\text{opt}}_{\text{vec}}$ can be obtained as a solution to the least
squares minimization problem,

$\bm{C}^{\text{opt}}_{\text{vec}}=\operatornamewithlimits{arg\
min}_{\bm{C}_{\text{vec}}}{\|\mathcal{A}\tilde{\bm{\varepsilon}}\bm{C}_{\text{{vec}}}-\bm{F}_{\text{vec}}\|^{2}}.$
(36)

It should be noted that we use this method for material parameter
identification, only using the numerical data. We did not experimentally
validate this method, since shear testing is non-trivial and requires
dedicated setups, such as a hexapod machine (Dalemat, 2019).

Figure 9: Parameter identification of an anisotropic metamaterial by
performing three different tests.

##  Appendix B Additional data

Method |  $C_{11}$ (MPa) |  $C_{12}$ (MPa) |  $C_{22}$ (MPa) |  $C_{16}$ (MPa) |  $C_{26}$ (MPa) |  $C_{66}$ (MPa)
---|---|---|---|---|---|---
Homogenization | 127.14 | 59.50 | 105.83 | 73.42 | 55.16 | 110.15
VFM | 129.92 | 58.84 | 102.72 | 77.43 | 55.25 | 113.73
Multiple tests | 123.82 | 57.44 | 102.35 | 72.74 | 54.68 | 110.12
Table 3: Comparison of stiffness tensor parameters identified for geometry #2 with 25 unit cell tessellation based on synthetic data using the VFM and the conventional methods against the computational homogenization. Method |  $C_{11}$ (MPa) |  $C_{12}$ (MPa) |  $C_{22}$ (MPa) |  $C_{16}$ (MPa) |  $C_{26}$ (MPa) |  $C_{66}$ (MPa)
---|---|---|---|---|---|---
Homogenization | 44.70 | -9.42 | 107.19 | -12.52 | -20.71 | 105.35
VFM | 33.95 | -5.84 | 106.48 | -4.44 | -20.78 | 92.65
Multiple tests | 45.78 | -19.27 | 112.01 | -11.30 | -20.30 | 105.65
Table 4: Comparison of stiffness tensor parameters identified for geometry #3 with 25 unit cell tessellation based on synthetic data using the VFM and the conventional methods against the computational homogenization. Method |  $C_{11}$ (MPa) |  $C_{12}$ (MPa) |  $C_{22}$ (MPa) |  $C_{16}$ (MPa) |  $C_{26}$ (MPa) |  $C_{66}$ (MPa)
---|---|---|---|---|---|---
Homogenization | 65.74 | 40.36 | 65.74 | 18.95 | 18.95 | 86.47
VFM | 65.45 | 40.94 | 65.43 | 17.93 | 17.76 | 82.19
Multiple tests | 66.66 | 41.58 | 66.66 | 18.75 | 18.75 | 86.86
Table 5: Comparison of stiffness tensor parameters identified for geometry #4
with 25 unit cell tessellation based on synthetic data using the VFM and the
conventional methods against the computational homogenization.

Figure 10: Axial and shear load-displacement data for all the experimentally
tested specimens. Figure 11: Marginal posterior probability distributions of
the material parameters obtained through Bayesian linear regression on the
numerical data.

## References

  * Agnelli et al. (2021) Agnelli, F., Margerit, P., Celli, P., Daraio, C., Constantinescu, A., 2021.  Systematic two-scale image analysis of extreme deformations in soft architectured sheets.  International Journal of Mechanical Sciences 194, 106205.  URL: https://www.sciencedirect.com/science/article/pii/S0020740320343101, doi:10.1016/j.ijmecsci.2020.106205.
  * Agnelli et al. (2022) Agnelli, F., Nika, G., Constantinescu, A., 2022.  Design of thin micro-architectured panels with extension–bending coupling effects using topology optimization.  Computer Methods in Applied Mechanics and Engineering 391, 114496.  URL: https://www.sciencedirect.com/science/article/pii/S0045782521007040, doi:10.1016/j.cma.2021.114496.
  * Ahn et al. (2017) Ahn, Y.K., Lee, H.J., Kim, Y.Y., 2017.  Conical Refraction of Elastic Waves by Anisotropic Metamaterials and Application for Parallel Translation of Elastic Waves.  Scientific Reports 7, 10072\.  URL: http://www.nature.com/articles/s41598-017-10691-6, doi:10.1038/s41598-017-10691-6.
  * Andreassen and Andreasen (2014) Andreassen, E., Andreasen, C.S., 2014\.  How to determine composite material properties using numerical homogenization.  Computational Materials Science 83, 488–495.  URL: https://www.sciencedirect.com/science/article/pii/S0927025613005302, doi:10.1016/j.commatsci.2013.09.006.
  * Andreaus et al. (2016) Andreaus, U., dell’Isola, F., Giorgio, I., Placidi, L., Lekszycki, T., Rizzi, N.L., 2016\.  Numerical simulations of classical problems in two-dimensional (non) linear second gradient elasticity.  International Journal of Engineering Science 108, 34–50.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0020722516308187, doi:10.1016/j.ijengsci.2016.08.003.
  * Auffray and Ropars (2016) Auffray, N., Ropars, P., 2016\.  Invariant-based reconstruction of bidimensional elasticity tensors.  International Journal of Solids and Structures 87, 183–193.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0020768316000676, doi:10.1016/j.ijsolstr.2016.02.013.
  * Avril et al. (2008) Avril, S., Bonnet, M., Bretelle, A.S., Grédiac, M., Hild, F., Ienny, P., Latourte, F., Lemosse, D., Pagano, S., Pagnacco, E., Pierron, F., 2008.  Overview of Identification Methods of Mechanical Parameters Based on Full-field Measurements.  Experimental Mechanics 48, 381–402.  URL: http://link.springer.com/10.1007/s11340-008-9148-y, doi:10.1007/s11340-008-9148-y.
  * Avril et al. (2004) Avril, S., Grédiac, M., Pierron, F., 2004.  Sensitivity of the virtual fields method to noisy data.  Computational Mechanics 34, 439–452.  URL: https://doi.org/10.1007/s00466-004-0589-6, doi:10.1007/s00466-004-0589-6.
  * Bastek et al. (2022) Bastek, J.H., Kumar, S., Telgen, B., Glaesener, R.N., Kochmann, D.M., 2022.  Inverting the structure–property map of truss metamaterials by deep learning.  Proceedings of the National Academy of Sciences 119, e2111505119.  URL: https://www.pnas.org/doi/abs/10.1073/pnas.2111505119, doi:10.1073/pnas.2111505119. publisher: Proceedings of the National Academy of Sciences.
  * Bertoldi et al. (2017) Bertoldi, K., Vitelli, V., Christensen, J., van Hecke, M., 2017\.  Flexible mechanical metamaterials.  Nature Reviews Materials 2, 1–11.  URL: https://www.nature.com/articles/natrevmats201766, doi:10.1038/natrevmats.2017.66. number: 11 Publisher: Nature Publishing Group.
  * Besnard et al. (2006) Besnard, G., Hild, F., Roux, S., 2006.  “Finite-Element” Displacement Fields Analysis from Digital Images: Application to Portevin–Le Châtelier Bands.  Experimental Mechanics 46, 789–803.  URL: https://doi.org/10.1007/s11340-006-9824-8, doi:10.1007/s11340-006-9824-8.
  * Cahn (1965) Cahn, J.W., 1965.  Phase Separation by Spinodal Decomposition in Isotropic Systems.  The Journal of Chemical Physics 42, 93–99.  URL: https://aip.scitation.org/doi/10.1063/1.1695731, doi:10.1063/1.1695731. publisher: American Institute of Physics.
  * Chen et al. (2018) Chen, W., Ruan, D., Huang, X., 2018.  Optimization for twist chirality of structural materials induced by axial strain.  Materials Today Communications 15, 175–184.  URL: https://www.sciencedirect.com/science/article/pii/S235249281830031X, doi:10.1016/j.mtcomm.2018.03.010.
  * Christensen et al. (2015) Christensen, J., Kadic, M., Wegener, M., Kraft, O., Wegener, M., 2015.  Vibrant times for mechanical metamaterials.  MRS Communications 5, 453–462.  URL: http://link.springer.com/10.1557/mrc.2015.51, doi:10.1557/mrc.2015.51.
  * Dalemat (2019) Dalemat, M., 2019.  Une expérimentation réussie pour l’identification de la réponse mécanique sans loi de comportement : Approche data-driven appliquée aux membranes élastomères.  Theses. École centrale de Nantes.  URL: https://tel.archives-ouvertes.fr/tel-02506891.
  * Dos Reis and Karathanasopoulos (2022) Dos Reis, F., Karathanasopoulos, N., 2022\.  Inverse metamaterial design combining genetic algorithms with asymptotic homogenization schemes.  International Journal of Solids and Structures 250, 111702.  URL: https://www.sciencedirect.com/science/article/pii/S0020768322002104, doi:10.1016/j.ijsolstr.2022.111702.
  * Every and Sachse (1990) Every, A.G., Sachse, W., 1990\.  Determination of the elastic constants of anisotropic solids from acoustic-wave group-velocity measurements.  Physical Review B 42, 8196–8205.  URL: https://link.aps.org/doi/10.1103/PhysRevB.42.8196, doi:10.1103/PhysRevB.42.8196. publisher: American Physical Society.
  * Fischer et al. (2011) Fischer, P., Klassen, M., Mergheim, J., Steinmann, P., Müller, R., 2011.  Isogeometric analysis of 2D gradient elasticity.  Computational Mechanics 47, 325–334.  URL: http://link.springer.com/10.1007/s00466-010-0543-8, doi:10.1007/s00466-010-0543-8.
  * Flaschel (2023) Flaschel, M., 2023.  Automated Discovery of Material Models in Continuum Solid Mechanics.  Ph.D. thesis. ETH Zurich.  URL: http://hdl.handle.net/20.500.11850/602750, doi:10.3929/ETHZ-B-000602750.
  * Flaschel et al. (2021) Flaschel, M., Kumar, S., De Lorenzis, L., 2021.  Unsupervised discovery of interpretable hyperelastic constitutive laws.  Computer Methods in Applied Mechanics and Engineering 381, 113852.  doi:10.1016/j.cma.2021.113852.
  * Flaschel et al. (2022) Flaschel, M., Kumar, S., De Lorenzis, L., 2022.  Discovering plasticity models without stress data.  npj Computational Materials 8, 91\.  URL: https://www.nature.com/articles/s41524-022-00752-4, doi:10.1038/s41524-022-00752-4.
  * Flaschel et al. (2023) Flaschel, M., Kumar, S., De Lorenzis, L., 2023.  Automated discovery of generalized standard material models with EUCLID.  Computer Methods in Applied Mechanics and Engineering 405, 115867.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0045782522008234, doi:10.1016/j.cma.2022.115867.
  * Forte and Vianello (2014) Forte, S., Vianello, M., 2014\.  A unified approach to invariants of plane elasticity tensors.  Meccanica 49, 2001–2012.  URL: http://link.springer.com/10.1007/s11012-014-9916-y, doi:10.1007/s11012-014-9916-y.
  * François et al. (1998) François, M., Geymonat, G., Berthaud, Y., 1998.  Determination of the symmetries of an experimentally determined stiffness tensor: Application to acoustic measurements.  International Journal of Solids and Structures 35, 4091–4106.  URL: https://www.sciencedirect.com/science/article/pii/S002076839700303X, doi:10.1016/S0020-7683(97)00303-X.
  * Frenzel et al. (2017) Frenzel, T., Kadic, M., Wegener, M., 2017.  Three-dimensional mechanical metamaterials with a twist.  Science 358, 1072–1074.  URL: https://www.sciencemag.org/lookup/doi/10.1126/science.aao4640, doi:10.1126/science.aao4640.
  * Gras et al. (2015) Gras, R., Leclerc, H., Hild, F., Roux, S., Schneider, J., 2015.  Identification of a set of macroscopic elastic parameters in a 3D woven composite: Uncertainty analysis and regularization.  International Journal of Solids and Structures 55, 2–16.  URL: https://www.sciencedirect.com/science/article/pii/S0020768313004964, doi:10.1016/j.ijsolstr.2013.12.023.
  * Greaves et al. (2011) Greaves, G.N., Greer, A.L., Lakes, R.S., Rouxel, T., 2011\.  Poisson’s ratio and modern materials.  Nature Materials 10, 823–837.  URL: https://www.nature.com/articles/nmat3134, doi:10.1038/nmat3134. number: 11 Publisher: Nature Publishing Group.
  * Grédiac (1989) Grédiac, M., 1989.  Principle of virtual work and identification.  Comptes Rendus de L Academie des Sciences Serie Ii , 1–5.
  * Grédiac and Pierron (2006) Grédiac, M., Pierron, F., 2006\.  Applying the Virtual Fields Method to the identification of elasto-plastic constitutive parameters.  International Journal of Plasticity 22, 602–627.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0749641905000896, doi:10.1016/j.ijplas.2005.04.007.
  * Grédiac et al. (2008) Grédiac, M., Pierron, F., Avril, S., Toussaint, E., 2008\.  The Virtual Fields Method for Extracting Constitutive Parameters From Full-Field Measurements: a Review.  Strain 42, 233–253.  URL: http://doi.wiley.com/10.1111/j.1475-1305.2006.tb01504.x, doi:10.1111/j.1475-1305.2006.tb01504.x.
  * Guseinov et al. (2020) Guseinov, R., McMahan, C., Pérez, J., Daraio, C., Bickel, B., 2020.  Programming temporal morphing of self-actuated shells.  Nature Communications 11, 237\.  URL: https://www.nature.com/articles/s41467-019-14015-2, doi:10.1038/s41467-019-14015-2. number: 1 Publisher: Nature Publishing Group.
  * Hild and Roux (2012) Hild, F., Roux, S., 2012.  Comparison of Local and Global Approaches to Digital Image Correlation.  Experimental Mechanics 52, 1503–1519.  URL: http://link.springer.com/10.1007/s11340-012-9603-7, doi:10.1007/s11340-012-9603-7.
  * Huang et al. (2020) Huang, D.Z., Xu, K., Farhat, C., Darve, E., 2020\.  Learning constitutive relations from indirect observations using deep neural networks.  Journal of Computational Physics 416, 109491.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0021999120302655, doi:10.1016/j.jcp.2020.109491.
  * Joshi et al. (2022) Joshi, A., Thakolkaran, P., Zheng, Y., Escande, M., Flaschel, M., De Lorenzis, L., Kumar, S., 2022.  Bayesian-EUCLID: Discovering hyperelastic material laws with uncertainties.  Computer Methods in Applied Mechanics and Engineering 398, 115225.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0045782522003681, doi:10.1016/j.cma.2022.115225.
  * Kadic et al. (2012) Kadic, M., Bückmann, T., Stenger, N., Thiel, M., Wegener, M., 2012.  On the practicability of pentamode mechanical metamaterials.  Applied Physics Letters 100, 191901\.  URL: https://aip.scitation.org/doi/10.1063/1.4709436, doi:10.1063/1.4709436. publisher: American Institute of Physics.
  * Kaina et al. (2015) Kaina, N., Lemoult, F., Fink, M., Lerosey, G., 2015\.  Negative refractive index and acoustic superlens from multiple scattering in single negative metamaterials.  Nature 525, 77–81.  URL: https://www.nature.com/articles/nature14678, doi:10.1038/nature14678. number: 7567 Publisher: Nature Publishing Group.
  * Karathanasopoulos et al. (2020) Karathanasopoulos, N., Dos Reis, F., Diamantopoulou, M., Ganghoffer, J.F., 2020\.  Mechanics of beams made from chiral metamaterials: Tuning deflections through normal-shear strain couplings.  Materials & Design 189, 108520\.  URL: https://www.sciencedirect.com/science/article/pii/S0264127520300538, doi:10.1016/j.matdes.2020.108520.
  * Kim et al. (2020) Kim, C., Kim, J.H., Lee, M.G., 2020.  A virtual fields method for identifying anisotropic elastic constants of fiber reinforced composites using a single tension test: Theory and validation.  Composites Part B: Engineering 200, 108338.  URL: https://www.sciencedirect.com/science/article/pii/S1359836820333874, doi:10.1016/j.compositesb.2020.108338.
  * Kulagin et al. (2020) Kulagin, R., Beygelzimer, Y., Estrin, Y., Schumilin, A., Gumbsch, P., 2020.  Architectured Lattice Materials with Tunable Anisotropy: Design and Analysis of the Material Property Space with the Aid of Machine Learning.  Advanced Engineering Materials 22, 2001069.  URL: https://onlinelibrary.wiley.com/doi/10.1002/adem.202001069, doi:10.1002/adem.202001069.
  * Kumar et al. (2020) Kumar, S., Tan, S., Zheng, L., Kochmann, D.M., 2020.  Inverse-designed spinodoid metamaterials.  npj Computational Materials 6.  URL: https://doi.org/10.1038/s41524-020-0341-6, doi:10.1038/s41524-020-0341-6.
  * Lee et al. (2016) Lee, H.J., Lee, H.S., Ma, P.S., Kim, Y.Y., 2016.  Effective material parameter retrieval of anisotropic elastic metamaterials with inherent nonlocality.  Journal of Applied Physics 120, 104902\.  URL: https://aip.scitation.org/doi/10.1063/1.4962274, doi:10.1063/1.4962274. publisher: American Institute of Physics.
  * Lee et al. (2012) Lee, J.H., Singer, J.P., Thomas, E.L., 2012.  Micro-/Nanostructured Mechanical Metamaterials.  Advanced Materials 24, 4782–4810.  URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/adma.201201644, doi:10.1002/adma.201201644. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/adma.201201644.
  * Liu et al. (2020) Liu, X., Tao, F., Du, H., Yu, W., Xu, K., 2020.  Learning Nonlinear Constitutive Laws Using Neural Network Models Based on Indirectly Measurable Data.  Journal of Applied Mechanics 87, 081003\.  URL: https://asmedigitalcollection.asme.org/appliedmechanics/article/doi/10.1115/1.4047036/1083320/Learning-Nonlinear-Constitutive-Laws-Using-Neural, doi:10.1115/1.4047036.
  * Man and Furukawa (2011) Man, H., Furukawa, T., 2011\.  Neural network constitutive modelling for non-linear characterization of anisotropic materials.  International Journal for Numerical Methods in Engineering 85, 939–957.  URL: https://onlinelibrary.wiley.com/doi/10.1002/nme.2999, doi:10.1002/nme.2999.
  * Mao et al. (2020) Mao, H., Rumpler, R., Gaborit, M., Göransson, P., Kennedy, J., O’Connor, D., Trimble, D., Rice, H., 2020\.  Twist, tilt and stretch: From isometric Kelvin cells to anisotropic cellular materials.  Materials & Design 193, 108855\.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0264127520303890, doi:10.1016/j.matdes.2020.108855.
  * Marek et al. (2017) Marek, A., Davis, F.M., Pierron, F., 2017.  Sensitivity-based virtual fields for the non-linear virtual fields method.  Computational Mechanics 60, 409–431.  URL: http://link.springer.com/10.1007/s00466-017-1411-6, doi:10.1007/s00466-017-1411-6.
  * Marino et al. (2023) Marino, E., Flaschel, M., Kumar, S., De Lorenzis, L., 2023\.  Automated identification of linear viscoelastic constitutive laws with EUCLID.  Mechanics of Materials 181, 104643\.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0167663623000893, doi:10.1016/j.mechmat.2023.104643.
  * Milton and Cherkaev (1995) Milton, G.W., Cherkaev, A.V., 1995\.  Which Elasticity Tensors are Realizable?  Journal of Engineering Materials and Technology 117, 483–493.  URL: https://doi.org/10.1115/1.2804743, doi:10.1115/1.2804743.
  * 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.  URL: https://www.sciencedirect.com/science/article/pii/002076836890036X, doi:10.1016/0020-7683(68)90036-X.
  * Ni et al. (2019) Ni, X., Guo, X., Li, J., Huang, Y., Zhang, Y., Rogers, J.A., 2019.  2D Mechanical Metamaterials with Widely Tunable Unusual Modes of Thermal Expansion.  Advanced Materials 31, 1905405\.  URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/adma.201905405, doi:10.1002/adma.201905405. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/adma.201905405.
  * Pierron (2023) Pierron, F., 2023.  Material testing 2.0: A brief review.  Strain n/a, e12434.  URL: https://onlinelibrary.wiley.com/doi/abs/10.1111/str.12434, doi:10.1111/str.12434.
  * Pierron et al. (2010) Pierron, F., Avril, S., Tran, V.T., 2010.  Extension of the virtual fields method to elasto-plastic material identification with cyclic loads and kinematic hardening.  International Journal of Solids and Structures 47, 2993–3010.  URL: https://www.sciencedirect.com/science/article/pii/S002076831000243X, doi:10.1016/j.ijsolstr.2010.06.022.
  * Pierron and Grédiac (2012) Pierron, F., Grédiac, M., 2012\.  The Virtual Fields Method.  Springer New York, New York, NY.  URL: http://link.springer.com/10.1007/978-1-4614-1824-5, doi:10.1007/978-1-4614-1824-5.
  * Podestá et al. (2019) Podestá, J., Méndez, C., Toro, S., Huespe, A., 2019\.  Symmetry considerations for topology design in the elastic inverse homogenization problem.  Journal of the Mechanics and Physics of Solids 128, 54–78.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0022509618305210, doi:10.1016/j.jmps.2019.03.018.
  * Promma et al. (2009) Promma, N., Raka, B., Grédiac, M., Toussaint, E., Le Cam, J.B., Balandraud, X., Hild, F., 2009.  Application of the virtual fields method to mechanical characterization of elastomeric materials.  International Journal of Solids and Structures 46, 698–715.  URL: https://www.sciencedirect.com/science/article/pii/S0020768308003946, doi:10.1016/j.ijsolstr.2008.09.025.
  * Risso et al. (2021) Risso, G., Sakovsky, M., Ermanni, P., 2021.  Instability-driven shape forming of fiber reinforced polymer frames.  Composite Structures 268, 113946\.  URL: https://www.sciencedirect.com/science/article/pii/S0263822321004062, doi:10.1016/j.compstruct.2021.113946.
  * Rosakis and Ravi-Chandar (1986) Rosakis, A.J., Ravi-Chandar, K., 1986\.  On crack-tip stress state: An experimental evaluation of three-dimensional effects.  International Journal of Solids and Structures 22, 121–134.  URL: https://www.sciencedirect.com/science/article/pii/0020768386900028, doi:10.1016/0020-7683(86)90002-8.
  * Roux and Hild (2020) Roux, S., Hild, F., 2020.  Optimal procedure for the identification of constitutive parameters from experimentally measured displacement fields.  International Journal of Solids and Structures 184, 14–23.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0020768318304542, doi:10.1016/j.ijsolstr.2018.11.008.
  * Rychlewski (1984) Rychlewski, J., 1984.  On Hooke’s law.  Journal of Applied Mathematics and Mechanics 48, 303–314.  URL: https://www.sciencedirect.com/science/article/pii/0021892884901370, doi:10.1016/0021-8928(84)90137-0.
  * Schittny et al. (2013) Schittny, R., Bückmann, T., Kadic, M., Wegener, M., 2013\.  Elastic measurements on macroscopic three-dimensional pentamode metamaterials.  Applied Physics Letters 103, 231905\.  URL: https://aip.scitation.org/doi/10.1063/1.4838663, doi:10.1063/1.4838663. publisher: American Institute of Physics.
  * Soyarslan et al. (2018) Soyarslan, C., Bargmann, S., Pradas, M., Weissmüller, J., 2018\.  3d stochastic bicontinuous microstructures: Generation, topology and elasticity.  Acta Materialia 149, 326–340.  URL: https://doi.org/10.1016/j.actamat.2018.01.005, doi:10.1016/j.actamat.2018.01.005.
  * Surjadi et al. (2019) Surjadi, J.U., Gao, L., Du, H., Li, X., Xiong, X., Fang, N.X., Lu, Y., 2019\.  Mechanical Metamaterials and Their Engineering Applications.  Advanced Engineering Materials 21, 1800864.  URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/adem.201800864, doi:10.1002/adem.201800864. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/adem.201800864.
  * Sutton et al. (2009) Sutton, M.A., Orteu, J.J., Schreier, H., 2009.  Image correlation for shape, motion and deformation measurements: basic concepts, theory and applications.  Springer Science & Business Media.
  * Thakolkaran et al. (2022) Thakolkaran, P., Joshi, A., Zheng, Y., Flaschel, M., De Lorenzis, L., Kumar, S., 2022\.  NN-EUCLID: Deep-learning hyperelasticity without stress data.  Journal of the Mechanics and Physics of Solids 169, 105076.  URL: https://linkinghub.elsevier.com/retrieve/pii/S0022509622002538, doi:10.1016/j.jmps.2022.105076.
  * Tibshirani (1996) Tibshirani, R., 1996.  Regression Shrinkage and Selection via the Lasso.  Journal of the Royal Statistical Society: Series B (Methodological) 58, 267–288.  URL: http://doi.wiley.com/10.1111/j.2517-6161.1996.tb02080.x, doi:10.1111/j.2517-6161.1996.tb02080.x.
  * Ting and Chi-Tsai (1996) Ting, Chi-Tsai, T., 1996.  Anisotropic elasticity: theory and applications.  45, Oxford University Press on Demand.
  * Wu et al. (2019) Wu, W., Hu, W., Qian, G., Liao, H., Xu, X., Berto, F., 2019.  Mechanical design and multifunctional applications of chiral mechanical metamaterials: A review.  Materials & Design 180, 107950\.  URL: https://www.sciencedirect.com/science/article/pii/S0264127519303880, doi:10.1016/j.matdes.2019.107950.
  * Yang et al. (2021) Yang, H., Timofeev, D., Abali, B.E., Li, B., Müller, W.H., 2021.  Verification of strain gradient elasticity computation by analytical solutions.  ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik URL: https://onlinelibrary.wiley.com/doi/10.1002/zamm.202100023, doi:10.1002/zamm.202100023.
  * Yang et al. (2019) Yang, X., Kweun, M., Kim, Y.Y., 2019.  Monolayer metamaterial for full mode-converting transmission of elastic waves.  Applied Physics Letters 115, 071901\.  URL: https://aip.scitation.org/doi/10.1063/1.5109758, doi:10.1063/1.5109758. publisher: American Institute of Physics.
  * Yuan et al. (2021) Yuan, Z., Cui, Z., Ju, J., 2021.  Micropolar homogenization of wavy tetra-chiral and tetra-achiral lattices to identify axial–shear coupling and directional negative Poisson’s ratio.  Materials & Design 201, 109483\.  URL: https://www.sciencedirect.com/science/article/pii/S0264127521000368, doi:10.1016/j.matdes.2021.109483.
  * Zadpoor (2016) Zadpoor, A.A., 2016.  Mechanical meta-materials.  Materials Horizons 3, 371–381.  URL: http://xlink.rsc.org/?DOI=C6MH00065G, doi:10.1039/C6MH00065G.
  * Zheng et al. (2019) Zheng, M., Liu, X., Chen, Y., Miao, H., Zhu, R., Hu, G., 2019.  Theory and Realization of Nonresonant Anisotropic Singly Polarized Solids Carrying Only Shear Waves.  Physical Review Applied 12, 014027\.  URL: https://link.aps.org/doi/10.1103/PhysRevApplied.12.014027, doi:10.1103/PhysRevApplied.12.014027. publisher: American Physical Society.
  * Zheng et al. (2020) Zheng, M., Park, C.I., Liu, X., Zhu, R., Hu, G., Kim, Y.Y., 2020.  Non-resonant metasurface for broadband elastic wave mode splitting.  Applied Physics Letters 116, 171903\.  URL: https://aip.scitation.org/doi/10.1063/5.0005408, doi:10.1063/5.0005408. publisher: American Institute of Physics.

Supplementary Information

##  Appendix S-I Difference in displacement fields between homogeneous and
heterogeneous materials

Figure S1: Comparison between the displacement fields obtained from finite
element simulations of a homogeneous specimen (left) and a heterogeneous
structure made of geometry #2 (center). For the homogeneous specimen a finite
element simulation using $10\times 10$ bilinear quadrilateral elements was
executed. The heterogeneous specimen was simulated using using $1000\times
1000$ bilinear quadrilateral elements. Afterwards, the displacement data at
the unit cell corners were extracted and interpolated with a bilinear
polynomial for each unit cell, to allow for a comparison with the homogeneous
specimen. The difference between the fields is shown on the right.. Figure S2:
Comparison between the displacement fields obtained from finite element
simulations of a homogeneous specimen (left) and a heterogeneous structure
made of geometry #3 (center). For the homogeneous specimen a finite element
simulation using $10\times 10$ bilinear quadrilateral elements was executed.
The heterogeneous specimen was simulated using using $1000\times 1000$
bilinear quadrilateral elements. Afterwards, the displacement data at the unit
cell corners were extracted and interpolated with a bilinear polynomial for
each unit cell, to allow for a comparison with the homogeneous specimen. The
difference between the fields is shown on the right. Figure S3: Comparison
between the displacement fields obtained from finite element simulations of a
homogeneous specimen (left) and a heterogeneous structure made of geometry #4
(center). For the homogeneous specimen a finite element simulation using
$10\times 10$ bilinear quadrilateral elements was executed. The heterogeneous
specimen was simulated using using $1000\times 1000$ bilinear quadrilateral
elements. Afterwards, the displacement data at the unit cell corners were
extracted and interpolated with a bilinear polynomial for each unit cell, to
allow for a comparison with the homogeneous specimen. The difference between
the fields is shown on the right.

##  Appendix S-II Comparison of simulated and experimentally measured full-
field displacement fields

Figure S4: Comparison between numerical (top) and experimentally
measured(bottom) full-field displacement and strain field data for the
$10\times 10$ tessellation of unit cell geometry #2 subjected to displacement-
controlled uniaxial tension test. Figure S5: Comparison between numerical
(top) and experimentally measured (bottom) full-field displacement and strain
field data for the $10\times 10$ tessellation of unit cell geometry #3
subjected to displacement-controlled uniaxial tension test. Figure S6:
Comparison between numerical (top) and experimentally measured (bottom) full-
field displacement and strain field data for the $10\times 10$ tessellation of
unit cell geometry #4 subjected to displacement-controlled uniaxial tension
test.

##  Appendix S-III Comparison of simulated and experimentally measured
displacement fields after postprocessing

Figure S7: Comparison between numerical (top) and experimentally measured
(bottom) full-field displacement and strain field data for the $10\times 10$
tessellation of unit cell geometry #1 subjected to displacement-controlled
uniaxial tension. Figure S8: Comparison between the synthetic (left) and
experimentally measured (right) displacement fields of the heterogeneous
structure made of geometry #2. Note that bilinear polynomials are used to
interpolate the displacement data at the unit cell corners. Figure S9:
Comparison between the synthetic (left) and experimentally measured (right)
displacement fields of the heterogeneous structure made of geometry 3. Note
that bilinear polynomials are used to interpolate the displacement data at the
unit cell corners. Figure S10: Comparison between the synthetic (left) and
experimentally measured (right) displacement fields of the heterogeneous
structure made of geometry 4. Note that bilinear polynomials are used to
interpolate the displacement data at the unit cell corners.