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

Journal Pre-proof

Single-test evaluation of directional elastic properties of anisotropic

structured materials

Jagannadh Boddapati, Moritz Flaschel, Siddhant Kumar, Laura De

Lorenzis, Chiara Daraio

PII:

S0022-5096(23)00275-2

DOI:

https://doi.org/10.1016/j.jmps.2023.105471

Reference:

MPS 105471

To appear in:

Journal of the Mechanics and Physics of Solids

Received date : 5 May 2023

Revised date :

29 July 2023

Accepted date : 14 October 2023

Please cite this article as: J. Boddapati, M. Flaschel, S. Kumar et al., Single-test evaluation of

directional elastic properties of anisotropic structured materials.

Journal of the Mechanics and

Physics of Solids

(2023), doi: https://doi.org/10.1016/j.jmps.2023.105471.

This is a PDF file of an article that has undergone enhancements after acceptance, such as the

addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive

version of record. This version will undergo additional copyediting, typesetting and review before it

is published in its final form, but we are providing this version to give early visibility of the article.

Please note that, during the production process, errors may be discovered which could affect the

content, and all legal disclaimers that apply to the journal pertain.

©

2023 Published by Elsevier Ltd.

Journal Pre-proof

Graphical Abstract

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

Jagannadh Boddapati, Moritz Flaschel, Siddhant Kumar, Laura De Lorenzis, Chiara Daraio

Revised Marked PDF

Journal Pre-proof

Highlights

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

Jagannadh Boddapati, Moritz Flaschel, Siddhant Kumar, Laura De Lorenzis, Chiara Daraio

•

Data-driven identification of anisotropic linear elastic constitutive material parameters for structured materials.

•

The technique is based on the virtual fields method and uses a single tension test.

•

Experimental validation performed on additively manufactured specimens.

•

Shear-normal coupling parameters are experimentally evaluated.

Journal Pre-proof

Single-test evaluation of directional elastic properties of anisotropic structured

materials

Jagannadh Boddapati

a,1

, Moritz Flaschel

b,1

, Siddhant Kumar

, Laura De Lorenzis

, Chiara Daraio

∗

Division of Engineering and Applied Science, California Institute of Technology, Pasadena, CA 91125, USA

Department of Mechanical and Process Engineering, ETH Z ̈urich, 8092 Z ̈urich, 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 in-

creases. For example, in two-dimensions, the anisotropic behavior of materials is described by up to 6 independent

elastic sti

ness 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 sti

ness tensor parameters of

two-dimensional anisotropic structured materials using just one tension test, thus eliminating the need for multiple ex-

periments, 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

∗

Corresponding author

Email address:

daraio@caltech.edu

(Chiara Daraio)

These authors contributed equally.

Preprint submitted to Journal of the Mechanics and Physics of Solids

July 28, 2023

Journal Pre-proof

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 e

ective 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 param-

eters 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 di

erent approaches to experimentally measure the elastic

parameters for di

erent anisotropy classes (Schittny et al., 2013; Considine et al., 2014; 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 sti

ness

tensor components when the o

-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 di

erent acoustic

wave speeds along di

erent material directions involve multiple tests and assume a certain material symmetry in pre-

dicting elastic parameters (Every and Sachse, 1990; Franc ̧ois et al., 1998). To date, there are no experimental methods

that can measure the sti

ness parameters of fully anisotropic structured materials 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 di

erent loading conditions are needed. Full-field identification methods allow extracting

Journal Pre-proof

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 (E

ffi

cient 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 ma-

terial 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

ediac (1989) (see also Gr

ediac et al. (2008); Pierron and Gr

ediac (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

ediac 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;

Journal Pre-proof

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

and the infinitesimal strain tensor

, is given by the generalized Hooke’s law (Rychlewski, 1984; Ting

and Chi-Tsai, 1996),

or (

i j

i jkl

)

(1)

where

is a fourth-order tensor, known as the elasticity tensor or the sti

ness 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



















1111

1122

1112

1122

2222

2212

1112

2212

1212



















(2)

where

1111

1122

2222

1112

2212

1212

are the elasticity tensor parameters in a given reference frame,

,ε

are the axial strains,

is the shear strain,

,σ

are the axial stresses, and

is the shear stress. For readability,

we combine the pair of indices as follows: ()

→

()

→

()

→

()

and write Eq. (2) as





































(3)

Our objective is to identify these six material parameters

from experimental measure-

ments while fulfilling certain constraints. From thermodynamic constraints, the elasticity tensor has to be positive

Journal Pre-proof

definite, which implies

(4a)

−

(4b)

The sti

ness parameter

represents the extension-to-extension deformation coupling. The sti

ness parameters

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

and

will be zero if the material has symmetry planes along the

and

axes.

Thus, the existence of shear-normal coupling and the maximum number of independent sti

ness tensor parameters

depend on the symmetries associated with the material microscopic topology (Ting and Chi-Tsai, 1996; Podest

a et al.,

2019). In plane elasticity, sti

ness tensors are categorized into four symmetry classes. They are denoted as

(2) for

Isotropic,

for Tetragonal,

for Orthotropic and

for Digonal (fully anisotropic) with 2, 3, 4 and 6 independent

parameters respectively. This categorization is based on the invariants of the sti

ness tensor (Forte and Vianello,

2014; Au

ray 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)

. We focus on identifying the e

ective 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

(

), as a linear superposition of cosine periodic functions:

(

)

∑

cos

(

)

∀

(

)

∈

[

−

∀

∈

[

−

(5)

where

are spatial frequencies, and

are the corresponding cosine function weights. The function is then

thresholded at a value

, 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

100 square mesh. In this pixelated representation, the gray phase

represents a sti

er material and the black phase represents a softer material (see Section 4.3.1).

Our methods are easily extendable to non-square unit cells.

Journal Pre-proof

The periodicity is ensured from the choice of the cosine functions directly. We randomly sample the weights

and the threshold value

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 sti

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

. A two-dimensional square anisotropic structured solid with 10

10 unit cell tessellation, with side length

, is subjected to a displacement-controlled tension test. The boundary

conditions are such that the bottom end is fixed, while a displacement of

]

is prescribed at the top end. The

reaction force components measured at the fixed end are denoted as

Figure 1: a) Design of an anisotropic unit cell geometry by thresholding a periodic function

(

). b) A two-phase unit cell geometry consisting

of a sti

er (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

- and

- 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

ediac, 1989; Gr

ediac et al., 2008; Pierron and Gr

ediac, 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

Journal Pre-proof

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 di

erent 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 sti

ness is set

to the homogenized sti

ness of the metamaterial. Under the same loading conditions, the two specimens exhibit

erent local displacements, which is likely caused by local size e

ects and the di

erent 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

ects 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).

3.2. Required data

To identify the unknown parameters, the VFM needs diverse local strain data, i.e., strain fields that are not ho-

mogeneous. Therefore, data that serve as input for the VFM are usually generated by testing complex specimen

geometries under complex loading conditions (Kim et al., 2014; Rossi et al., 2015). However, for the structured

materials considered in this study, generating complex specimen geometries (e.g., a plate with holes or notches) and

studying their behavior under complex loading conditions is not trivial. Local features like holes or notches may lead

We 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 di

ffi

cult.

Journal Pre-proof

to stress concentrations or singularities and hence more complex geometries may be prone to behave locally inelas-

tic or even fail. Further, we will show later that the VFM performs best if we do not feed the method with data at

the boundary of the specimen. Introducing more features like holes or notches to the geometry would increase the

boundary and thus reduce the amount of data suitable to be used by the VFM.

For our purposes we will show that, due to the anisotropy of the material, a clamped square plate under uniaxial

tension produces a su

ffi

ciently heterogeneous strain field. We hence consider a displacement-controlled uniaxial

tension experiment of a square-shaped specimen that consists of

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 (

(

unit cell corners and the net reaction forces. A quadrilateral finite element mesh is generated such that each of

the

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

(

) is hence approximated by

(

)

∑

(

)

(6)

where

(

denotes the number of nodes in the finite element mesh and

are the known nodal displacements,

while

(

) 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.,

(

)

(

∇

(

)

(

∇

(

)

3.3. Weak formulation of linear momentum balance

We denote the specimen domain and its boundary as

Ω

and

∂

Ω

, respectively, and the surface traction force acting

∂

Ω

. Assuming no inertia and body forces, the weak form of linear momentum balance reads

∫

Ω

(

) :

∇

(

) d

−

∫

∂

Ω

(

) d

(7)

which has to hold true for all admissible, i.e., su

ffi

ciently regular, test functions

(

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

(

) used to interpolate the displacement data

(

)

∑

(

)

(8)

Journal Pre-proof

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

∑







∫

Ω

∇

(

) d

︸

︷︷

︸

int

−

∫

∂

Ω

(

) d

︸

︷︷

︸

ext







(9)

where the first and second integral are the nodal internal forces

int

and nodal external forces

ext

, respectively. By

employing the constitutive relation Eq. (3), the nodal internal forces may be written as

int

∫

Ω

∇

∫

Ω













∫

Ω













∫

Ω













vec

(10)

where the elasticity tensor parameters

vec

[

]

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 di

erent choices of test

functions will yield di

erent 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 di

erent, a phenomenon that is best observed at the boundary and at the

corners of the specimen where local e

ects 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

{

,...,

}

as the set of all unit cells and

as the set of all

nodes corresponding to the unit cell

∈C

, and define a set of test functions as











(

)

∑

∈D

(

)

∈C

∈{

}











(11)

where

are the unit vectors in the corresponding

- and

-direction. Note that the test functions are normalized by

dividing by the number of nodes corresponding to the unit cell

(equal to 4 in our case).

Journal Pre-proof

Using the test functions in

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

ext

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

int











(

)

∑

∈D

(

)

∈C

int

∈{

}











(12)

where

int

⊂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

int

depend

, each field is zero in either

- or

-direction. The non-zero component of an exemplary virtual field in

int

shown in Fig. 2 (left).

Evaluating Eq. (9) for this set of functions leads to

∑

∈D

int

∀

∈C

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

int

(left) and non-zero component of a virtual field in

center

(right).

Equations (13) are not su

ffi

cient to identify the unknown material parameters, as the trivial solution

vec

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

Ω

∗

{

≤

}

and its boundary as

∂

Ω

∗

, the weak form of linear momentum balance for this domain reads

∫

Ω

∗

(

) :

∇

(

) d

−

∫

∂

Ω

∗

(

) d

(14)

Journal Pre-proof

Inserting the test function ansatz leads to

∑







∫

Ω

∗

∇

(

) d

︸

︷︷

︸

∗

int

−

∫

∂

Ω

∗

(

) d

︸

︷︷

︸

∗

ext







(15)

We define

center

{

}

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

, we consider instead

center

{

}

. We choose a set of virtual fields

center

that are constant along

center

and zero at all other

nodes

center











(

)

∑

∈D

center

(

)

∈{

}











(16)

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

∑

∈D

center

∗

int

∫

∂

Ω

center

(17)

where

∂

Ω

center

is the top boundary of

Ω

∗

. 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

, meaning that the sum of the internal forces at

∂

Ω

center

must equal the net reaction force.

3.6. Deterministic parameter identification

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

equations

int

vec

(18)

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

center

vec

(19)

where

int

and

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

opt

vec

arg min

vec

(

‖

int

vec

‖

center

vec

−

‖

)

(20)

where

‖·‖

is the Euclidean norm and

0 is a weighting parameter that scales the di

erent contributions to the

minimization problem. As there are less equations in the system (19) than in (18), the weighting parameter should be

Journal Pre-proof

chosen su

ffi

ciently larger than one (

1). Following previous works Flaschel et al. (2021, 2022, 2023), we choose

100 and keep it constant throughout this work. Based on our experience, the choice of

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

opt

vec

with

(

int

)

int

(

center

)

center

(

center

)

(21)

which leads to a determined system of equations that can be solved for

opt

vec

. The minimization problem in Eq. (20)

can alternatively be written as

opt

vec

arg min

vec

‖

vec

−

‖

(22)

where we have defined







int

√

center













√







(23)

The necessary condition for a minimum then reads

opt

vec

(24)

3.7. Bayesian inference

Besides the previously introduced deterministic approach, we further study the problem from a stochastic perspec-

tive. To this end, we construct a Bayesian linear regression model, for which we assume no intercept and a di

use

prior, as implemented in the Matlab

built-in function

bayeslm

We denote the number of rows in

and we define

with

∈ {

,...,

}

as the

row of

. For

each equation in the overdetermined system of equations

vec

, we assume the likelihood of obtaining

as a

Gaussian likelihood with mean

vec

and standard deviation

σ >

0, i.e.,

(

vec

,σ

)

√

πσ

exp

[

−

(

−

vec

)

]

(25)

where

vec

and

are treated as random variables. Assuming further that the likelihoods are conditionally indepen-

dent, we define the joint likelihood as

(

vec

,σ

)

∏

(

vec

,σ

)

(26)

Assuming here a di

use prior for the joint prior distribution of

vec

and

, i.e.,

(

vec

,σ

)

∝

(27)

the marginal posterior distributions of

vec

and

are analytically tractable and implemented in the Matlab

function

bayeslm

Journal Pre-proof

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, experi-

mental setup, testing and DIC.

4.1. Design and choice of unit cell geometries

We pick four unit cells with distinct

diverse e

ective sti

ness tensor parameters (all with six non-zero sti

ness

parameters). Table 1 shows the unit cells along with their symmetry class and homogenized sti

ness tensor. Geometry

#1 has

as the largest sti

ness parameter with

almost comparable to

and

. While geometry #2

has

as the largest sti

ness parameter with

, geometry #3 has negative values for all of the o

-diagonal

parameters. Geometry #4 has four independent sti

ness parameters with

as one of the largest values among other

sti

ness parameters, along with

and

. The fill fraction of the sti

phase for all the unit cells lies

between 60 and 70 %.

Unit Cell

Geometry

Name

Homogenized Sti

ff

ness Tensor

(

) [MPa]

Elastic

Symmetry Class

Geometry #1







131

98 63

98 198

38 83

87 95







Geometry #2







127

50 105

16 110







Geometry #3







−

107

−

105







Geometry #4







74 40

36 18

36 65

74 18

95 18

95 86







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.

Journal Pre-proof

For tessellation, we vary the number of unit cells

between 5 and 25.

Homogenization:

We compute the e

ective mechanical properties of the unit cells using the theory of homoge-

nization 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 di

ffi

cult to fabricate, we here pick 10

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

5 mm excluding the

portion that goes into the grips. We use Stratasys’ proprietary material DM8530 for the sti

phase and TangoBlack for

the soft phase. The material properties (DM8530: Young’s modulus

1000

90 MPa and Poisson’s ratio

0.35,

TangoBlack: Young’s modulus

0.7 MPa and Poisson’s ratio

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 test-

ing 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 ̃

02 and a global strain rate of 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

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.

4.3.3. Digital image correlation

Given a reference image

and and a deformed image

, the correlation algorithm aims at minimizing the sum of

squared di

erences over the considered domain

Ω

∫

Ω

(

))

−

(

)

(28)

Journal Pre-proof

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

where

is the position in the reference image and

(

) is the displacement field which is interpolated as

(

)

∑

(

)

(29)

where

are a set of shape functions and

the associated degrees of freedom. There are two approaches to determine

the unknowns

local DIC

and

global DIC

(Hild and Roux, 2012). In the

local

approach, the region of interest

(

Ω

) 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 displace-

ments (as in Agnelli et al. (2021)). We choose an edge length of 18 pixels (

∼

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. In the future, experimental errors may be further reduced using the global DIC technique

with quadratic interpolation (Blaysat et al., 2020). 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.

Journal Pre-proof

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 sti

ness is equal to the homogenized

sti

ness of the heterogeneous structure. To simulate the displacement of a homogeneous body, we assume a 10

bilinear quadrilateral finite element mesh. The displacement of the heterogeneous body is computed on a much finer

mesh with 1000

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 di

erences due to local e

ects in the heterogeneous structure,

which appear to be dominant at the boundary and corners of the specimen.

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 to S6 for the other geometries). We observe very

good agreement between the numerical and experimental data, especially for the variables

,ε

. However, the

experimentally measured

appears to be slightly higher than the numerical data, by about 0.1 mm, for all the

geometries. Also the two

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 ( ̃

) is 0.02.

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

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. B.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 e

ective material parameters, it is important to understand the continuum behavior of the heterogeneous structures