of 15
Mechanics of Materials 180 (2023) 104630
Available online 17 March 2023
0167-6636/© 2023 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (
http://creativecommons.org/licenses/by/4.0/
).
Contents lists available at
ScienceDirect
Mechanics of Materials
journal homepage:
www.elsevier.com/locate/mecmat
Mesh d-refinement: A data-based computational framework to account for
complex material response
Sacha Wattel
a
, Jean-François Molinari
a
, Michael Ortiz
b
,
c
, Joaquin Garcia-Suarez
a
,
∗
a
Institute of Civil Engineering, Institute of Materials, École Polytechnique Fédérale de Lausanne (EPFL), CH 1015 Lausanne, Switzerland
b
Division of Engineering and Applied Science, California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA
c
Hausdorff Center for Mathematics, Universität Bonn, Endenicher Allee 60, 53115 Bonn, Germany
A R T I C L E I N F O
Keywords:
Data-driven computational mechanics
Hybrid formulation
Non-linearity
FEM-DD coupling
Model-free
A B S T R A C T
Model-free data-driven computational mechanics (DDCM) is a new paradigm for simulations in solid mechanics.
The modeling step associated to the definition of a material constitutive law is circumvented through the
introduction of an abstract phase space in which, following a pre-defined rule, physically-admissible states are
matched to observed material response data (coming from either experiments or lower-scale simulations). In
terms of computational resources, the search procedure that performs these matches is the most onerous step
in the algorithm. One of the main advantages of DDCM is the fact that it avoids regression-based, bias-prone
constitutive modeling. However, many materials do display a simple linear response in the small-strain regime
while also presenting complex behavior after a certain deformation threshold. Motivated by this fact, we
present a novel refinement technique that turns regular elements (equipped with a linear-elastic constitutive
law) into data-driven ones if they are expected to surpass the threshold known to trigger material non-linear
response. We term this technique ‘‘data refinement’’, ‘‘d-refinement’’ for short. It works both with data-driven
elements based on either DDCM or strain–stress relations learned from data using neural networks. Starting
from an initially regular FEM mesh, the proposed algorithm detects where the refinement is needed and iterates
until all elements presumed to display non-linearity become data-driven ones. Insertion criteria are discussed.
The scheme is well-suited for simulations that feature non-linear response in relatively small portions of the
domain while the rest remains linear-elastic. The method is validated against a traditional incremental solver
(i.e., Newton–Raphson method) and we show that the d-refinement framework can outperform it in terms of
speed at no loss of accuracy. We provide an application that showcases the advantage of the new method:
bridging scales in architected metamaterials. For this application, we also succinctly outline how d-refinement
can be used in conjunction with a neural network trained on microscale data.
1. Introduction
Data-driven computational mechanics (DDCM) is a
rapidly-expanding research field. Data-driven solvers were originally
devised to deal with problems in small-strains statics (
Kirchdoerfer
and Ortiz
, 2016
; Conti et al.
, 2018
); extensions to dynamics (
Kirch-
doerfer and Ortiz
, 2018
; Salahshoor and Ortiz
, 2023
; Garcia-Suarez
et al.
, 2022
), large deformations (
Platzer et al.
, 2021
; Conti et al.
,
2020; Korzeniowski and Weinberg
, 2022
) and dissipative (inelastic)
material response (
Eggersmann et al.
, 2019
; Carrara et al.
, 2022
, 2021
;
Karapiperis et al.
, 2020
, 2021
) followed.
The data-driven paradigm relies on using ‘‘observed information’’
directly, avoiding the fitting/calibration exercise that turns these ob-
servations into a fully-defined mathematical function (‘‘material con-
stitutive model’’). Thus, part of the efforts in the field have been
∗
Corresponding author.
E-mail address:
joaquin.garciasuarez@epfl.ch
(J. Garcia-Suarez).
focused on generating those datasets either from experimental observa-
tions (
Leygue et al.
, 2018
) or from microstructural simulations (
Kara-
piperis et al.
, 2020
, 2021
; Korzeniowski and Weinberg
, 2022
). Further
efforts are currently underway to maximize its numerical efficiency (
Eg-
gersmann et al.
, 2021a
,b; Karapiperis et al.
, 2021
; Korzeniowski and
Weinberg
, 2021
). The mathematical foundations of the method have
been solidly laid out both in small strains (
Conti et al.
, 2018
) and large
deformations (
Conti et al.
, 2020
).
Arguably, the main benefit of the DDCM paradigm is that it enables
circumventing the definition of intricate phenomenological laws, which
are required to describe complex material response during numerical
simulations (
Kovachki et al.
, 2022
). However, most materials do display
a simple, consistent linear-elastic behavior as long as they remain in the
https://doi.org/10.1016/j.mechmat.2023.104630
Received 15 December 2022; Received in revised form 7 March 2023; Accepted 9 March 2023
Mechanics of Materials 180 (2023) 104630
2
S. Wattel et al.
small-strain regime, after which either non-linear or inelastic phenom-
ena are expected to kick in. Here, we consider that ‘‘complex response’’
is anything other than linear and elastic. As linear-elastic simulations
are both straightforward to implement and computationally efficient, it
is therefore desirable to stick to them to solve as much of the simulation
domain as possible. For this, one needs to devise a framework that can
detect when an element surpasses a deformation threshold, and, once
that trigger is indeed detected, proceed to enrich that element with a
data-driven formulation. This would also redound on simulation speed,
as it would reduce the number of phase-space searches, which remains
the most time-consuming step in the algorithm, despite the fact that
they are parallelizable.
A recent article of Yang et al. (2022) has been the first, to our
knowledge, that combined data-driven and traditional elements in
the same mesh. They used their coupling technique to analyze struc-
tures in which part of the main domain remains linear-elastic while
a portion of it weakens due to an external factor (corrosion). Tak-
ing advantage of this coupling scheme, we introduce a new kind
of mesh refinement technique that aims to improve both numerical
efficiency and accuracy while maintaining the advantages of DDCM:
‘‘data refinement’’, ‘‘d-refinement’’ for short, that automatically turns
regular elements into data-driven ones where and when necessary.
This eases the computational implementation, as it avoids intricacies
associated with numerically resolving explicitly either material non-
linearities, inelastic mechanisms or both. D-refinement also improves
accuracy: where constitutive modeling is more arduous, it replaces
assumption-ridden phenomenological descriptions by data, while keep-
ing the model in the region where it performs best (infinitesimal strain),
thus also partially avoiding errors associated with dealing with a dis-
crete material dataset (Kirchdoerfer and Ortiz, 2016). See Section 2.3.1
for further discussion on this point.
We will show that the d-refinement technique does not require
load stepping in the absence of inelasticity. Since it combines both
FEM and DDCM, the definition of the constants that metrize the space
in DDCM (Kirchdoerfer and Ortiz, 2016; Karapiperis et al., 2021) is
straightforward in d-refinement: we can make them equal to the elastic
constants of the FEM linear material.
The d-refinement framework forges a better union between model-
driven and data-driven ways of computing. It aspires to compete with
incremental solvers based on tangent operators as the tool-of-choice
when it comes to perform mechanics calculations that include either
non-linear behavior or inelasticity or both. Furthermore, if the dataset
comes from lower-scale simulations, d-refinement can also be regarded
as an optimized
퐹퐸
2
method (Feyel, 1999), in which (a) the microstruc-
tural response is obtained off-line before the simulation (Karapiperis
et al., 2021; Korzeniowski and Weinberg, 2022) and (b) microstructure
RVE response is directly used at the element level only when its
linear-elastic range is exhausted.
The paper is structured as follows: Section 2 reviews the basics of
DDCM and the ‘‘static’’ FEM-DD implementation, and introduces the
d-refinement solving procedure. Section 3 validates the methodology
via comparisons to conventional incremental solvers, we compare two
kinds of systems: 3D trusses and 2D plane-stress elements. Section 4
presents an application of d-refinement to multiscale analysis. Starting
from a material dataset, we use efficient data-driven simulations to
characterize the mechanical response of the unit cell of an architected
metamaterial, in terms of both linear-elastic constants and a dataset
that includes non-linear response. Then, this information is used to
study the K-field in a cracked macroscale sample (Shaikeea et al.,
2022). After conclusions in Section 5, we showcase in an Appendix A
how the d-refinement method can be modified to handle constitutive
models learned from data via trained neural networks.
2. Methods
2.1. Data-drivencomputationalmechanics
DDCM poses the solid mechanics boundary-value problem in an ab-
stract phase space of work-conjugate variables (e.g., strain and stress),
in which each element state is defined by a point in a ‘‘local’’ phase
space, and the overall system phase space is the Cartesian product
of that of its elements. The field equations (equilibrium and com-
patibility), along with the boundary conditions, represent a manifold
embedded in phase space. The infinite set that contains all possible
points that satisfy these constraints is termed
E
, and referred to as the
‘‘physically-admissible set’’. Traditionally, the constitutive law would
amount to a graph in this space (Conti et al., 2018), classical solutions
of the boundary-value problem would correspond to the intersection of
the aforementioned manifold and the graph. However, in DDCM, the
material behavior is represented by a ‘‘material discrete set’’, termed
D
.
As the intersection of admissible set
E
and
D
is unlikely, this concept is
replaced by the notion of ‘‘transversality’’ (Conti et al., 2018) between
E
and
D
.
The FEM-DD coupling scheme that we are to use and to introduce in
the next section converges to the original DD implementation (Kirch-
doerfer and Ortiz, 2016) when the whole mesh becomes data-driven;
hence, for the sake of brevity, we introduce the coupling scheme
algorithm directly in next section.
2.2. FEM-DDcoupling
2.2.1. Introduction
We begin with a brief digression concerning terminology. Model-
driven is a fair description of the current paradigm, which relies on
phenomenological models of material response to close the problem:
boundary conditions, force equilibrium and displacement compatibility
equations need to be supplemented with a force-deformation relation.
Nevertheless, using the acronym ‘‘MD’’ to refer to model-driven is
deemed confusing, as both in the mechanics and in the material science
communities it is associated with ‘‘molecular dynamics’’, a field where
data-driven approaches are also possible (Bulin et al., 2022). We would
rather use, at least in the context of solid and structural mechanics, the
name ‘‘DD-FEM coupling’’ in lieu of ‘‘DD-MD coupling’’. Small caveats
aside, the term is self-explanatory and, we believe, will be readily
understood by researchers in any of the fields mentioned before. In a
recent contribution by Korzeniowski and Weinberg (2022), the authors
refer to the conventional DDCM way of computing as ‘‘DD-FEM’’ as it is
‘‘data-driven in a finite-element mesh’’. However, since we tend to asso-
ciate the finite-element method not only with the discretization of the
domain, but also with a way to model materials (that is why we would
speak of ‘‘non-linear finite elements’’), we believe our terminology to
be appropriate.
The coupling implementation (Yang et al., 2022) divides the set of
all elements in the numerical model (
푆
) into two subsets:
푆
2
includes
the indices of DD elements, while
푆
1
contains those of the regular ones
(
푆
2
∩
푆
1
= ∅
and
푆
2
∪
푆
1
=
푆
). The size of each subset,
|
⋅
|
, is defined as
the number of elements in it, thus
|
푆
2
|
+
|
푆
1
|
=
푁
푒
, where
푁
푒
is the total
number of elements in the mesh. The concept of phase space distances
remain unchanged for DD elements: first, a point in the local phase
space
푍
푒
(of the
푒
th element,
푒
∈
푆
2
) is
풛
푒
= {
흐
푒
,
흈
푒
} ∈
R
2
푁
푐
, where
푁
푐
is the number of relevant components of stress/strain (one for uni-axial
elements, three for plane problems, six for 3D problems, and nine for
3D problems in generalized continua (Karapiperis et al., 2020)). This
local phase space is equipped with a metric defining a norm for each
point as
|
풛
푒
|
2
=
1
2
C
푒
흐
푒
⋅
흐
푒
+
1
2
C
−1
푒
흈
푒
⋅
흈
푒
,
(1)
where
C
푒
is a symmetric positive-definite matrix of constants not
necessarily related to any material property.
Mechanics of Materials 180 (2023) 104630
3
S. Wattel et al.
Given that the global phase space is
푍
=
푍
1
×
푍
2
...×
푍
|
푆
2
|
, the norm
associated to
풛
∈
푍
can be taken as
|
풛
|
=
(
∑
푒
∈
푆
2
푤
푒
|
풛
푒
|
)
1∕2
,
(2)
where
푤
푒
is the volume of the
푒
th element. The latter in turn provides
a notion of distance between two points,
풛
and
풚
, in phase space:
푑
(
풛
,
풚
) =
|
풛
−
풚
|
.
(3)
Note that the solution of traditional finite-element analysis also
represents a point in a phase space, and hence the global solution of
the system that contains both kinds of elements can be represented in
a larger phase space (i.e.,
푍
1
×
푍
1
×
⋯
×
푍
푒
).
By means of the constitutive law, assumed linear and elastic, we
can write
흈
푒
=
D
푒
흐
푒
for every element s.t.
푒
∈
푆
1
, with
D
푒
∈
R
푁
푐
×
푁
푐
containing the regular elastic constants of the
푒
th element. Using the
discretized compatibility relation, one can express the strain of an
element in terms of the global nodal displacement vector,
풖
∈
R
푁
dof
(
푁
dof
is the total number of degrees of freedom in the mesh),
흐
푒
=
푩
푒
풖
for
푒
= 1
,
...
,푁
푒
,
(4)
푩
푒
∈
R
푁
푐
×
푁
dof
encapsulates the nodal connectivity and the shape
functions used to perform nodal interpolation within each element. The
global equilibrium equation is expressed as
∑
푒
∈
푆
2
푤
푒
푩
⊤
푒
흈
푒
+
∑
푒
∈
푆
1
푤
푒
푩
⊤
푒
D
푩
푒
풖
=
풇
,
(5)
where
풇
∈
R
푁
dof
is the nodal external force vector.
Notice that the displacements of the regular elements enter the
nodal equilibrium equations along with the stress of the DD ones.
Using Lagrange multipliers,
휼
∈
R
푁
dof
, these equilibrium constraints
must be enforced concurrently with the distance minimization between
points in the admissible set,
풛
∈ E
, and the material set,
풛
∗
∈ D =
D
1
× D
2
×
⋯
× D
|
푆
2
|
, for the DD portion of the mesh (Kirchdoerfer and
Ortiz, 2016). Thus, we introduce a functional
훱
=
푑
2
(
풛
,
풛
∗
) +
휼
⋅
(
풇
−
∑
푒
∈
푆
2
푩
⊤
푒
흈
푒
−
∑
푒
∈
푆
1
푩
⊤
푒
D
푩
푒
풖
)
.
(6)
Compatibility is guaranteed by directly expressing, in the distance
function, the strains
흐
퐞
with displacements (Eq. (4)). The functional
obtained depends on
흈
퐞
,
풖
and
휼
. Enforcing the stationarity of
훱
renders a problem that can be solved iteratively, until the state of
the system that minimizes the distance between the solution and the
admissible points in the material dataset is attained.
Setting every variation of
훱
to zero, we obtain
훿
풖
⟹
∑
푒
∈
푆
2
푤
푒
푩
⊤
푒
C
푒
(
푩
푒
풖
−
흐
∗
푒
)
−
∑
푒
∈
푆
1
푩
⊤
푒
D
푩
푒
휼
= 0
,
(7a)
훿
흈
퐞
⟹
C
−1
푒
(
흈
퐞
−
흈
∗
퐞
) =
푩
푒
휼
,
(7b)
훿
휼
⟹
∑
푒
∈
푆
2
푤
푒
푩
⊤
푒
흈
푒
+
∑
푒
∈
푆
1
푤
푒
푩
⊤
푒
D
푩
푒
풖
=
풇
,
(7c)
which yields, after minor rearrangements of Eq. (7a) and introduc-
ing Eq. (7b) into Eq. (7c),
∑
푒
∈
푆
1
푤
푒
푩
⊤
푒
D
푒
푩
푒
풖
−
∑
푒
∈
푆
2
푤
푒
푩
⊤
푒
C
푒
푩
푒
휼
=
∑
푒
∈
푆
2
푤
푒
푩
⊤
푒
C
푒
흐
∗
푒
,
(8a)
∑
푒
∈
푆
1
푤
푒
푩
⊤
푒
D
푒
푩
푒
풖
+
∑
푒
∈
푆
2
푤
푒
푩
⊤
푒
C
푒
푩
푒
휼
=
풇
−
∑
푒
∈
푆
2
푤
푒
푩
⊤
푒
흈
∗
푒
,
(8b)
where we use the superindex ‘‘
∗
’’ to highlight elements in the material
dataset
D
. Given
풛
∗
푒
=
(
흐
∗
푒
,
흈
∗
푒
)
∈ D
푒
∀
푒
∈
푆
2
, Eq. (8) can be solved for
풖
and
휼
, and then physically-admissible strain and stresses for the DD
elements can be computed as Yang et al. (2022)
흈
푒
=
흈
∗
푒
+
D
푒
푩
푒
휼
,
(9a)
흐
푒
=
푩
푒
풖
.
(9b)
The subsequent solving of Eq. (8) and Eq. (9) is likened to a projec-
tion operation that takes the phase-space point
풛
∗
in
D
and projects it
onto the admissible set:
푃
퐸
풛
∗
∈ E
.
The next step of the algorithm involves projecting back onto the
material set, what is achieved through an element-wise search that
picks the point in
D
closest to the previously obtained
풛
∈ E
, and
it is symbolically expressed as
푃
퐷
풛
∈ D
. The necessary notion of
distance between points is provided by Eq. (3). In essence, this defines
the simplest data-driven solver, which allows fixed-point iterations
from the
푘
th admissible DD solution
풛
푘
to the next one, i.e.,
풛
(
푘
+1)
=
푃
퐸
푃
퐷
풛
(
푘
)
. The algorithm yields the final solution at the
푛
th iteration if
all the DD elements are assigned the same set of datapoints twice in a
row, i.e.,
풛
∗(
푛
)
=
푃
퐷
푃
퐸
풛
∗(
푛
−1)
=
풛
∗(
푛
−1)
.
The coupling is expected to converge with respect to the dataset
size since the linear-elastic finite elements can be thought as DD ones
in which
풛
=
풛
∗
always, and thus the theorems developed for pure DD
meshes (Conti et al., 2018, 2020) apply.
2.2.2. Implementation
Algorithm1
Fixed-point algorithm for DD-FEM coupling
Require:
∀
푒
= 1
,
...
,푁
푒
, compatibility matrices
푩
푒
;
∀
푒
∈
푆
1
,
D
푒
matrices containing elastic constants;
∀
푒
∈
푆
2
, local datasets
D
푒
and
C
푒
matrices containing distance constants;
∀
푖
= 1
,
...
,푁
dof
, external
forces
푓
푖
or boundary conditions.
(i) Set
푘
= 0
. Initial data assignation:
for
푒
∈
푆
2
do
Choose
풛
∗(0)
푒
= (
흈
∗(0)
푒
,
흐
∗(0)
푒
) ∈ D
푒
randomly
endfor
(ii) Project onto
E
:
Solve in eq. (8) for
풖
(
푘
)
and
휼
(
푘
)
for
푒
∈
푆
2
do
Compute
흐
푒
and
흈
푒
from eq. (9)
endfor
(iii) Project onto
D
:
for
푒
∈
푆
2
do
Choose
풛
∗(
푘
+1)
푒
∈ D
푒
closest to
풛
(
푘
)
푒
∈ E
푒
endfor
(iv) Convergence test:
if
풛
(
푘
+1)
푒
=
풛
(
푘
)
푒
∀
푒
∈
푆
2
then
Final displacement field
풖
=
풖
(
푘
)
exit
else
푘
←
푘
+ 1
,
goto
(ii)
endif
Algorithm 1 presents the pseudo-code implementation of a coupled
FEM-DD solver.
2.3. D-refinement
2.3.1. Motivation
As discussed in the introduction, we foresee data-driven simulations
in which some elements represent material behavior that does display a
meaningful, well-defined, linear-elastic range. Performing phase-space
searches for elements that do not abandon the small-strain regime
appears unnecessary, as no extra accuracy would be gained at the
cost of potential errors associated to dataset finite size. Since DDCM
requires no intersection between the sets but just transversality, there
is always some potential error that is introduced every time that a
material datum is assigned to an element (Kirchdoerfer and Ortiz,
2016). In summary, the phase space searches should be limited to
regions that clearly depart from linearity, what will be a small portion
of the total if the problem displays localization like, e.g., stark stress
Mechanics of Materials 180 (2023) 104630
4
S. Wattel et al.
Fig. 1.
Schematic representation of DD element initialization possibilities on a 2D phase space. The first circle (black) represents the initialization point and the second one (blue)
the point in the dataset that the algorithm converges to in the next iteration (see that each method may pick a different datum after the iteration). (For interpretation of the
references to color in this figure legend, the reader is referred to the web version of this article.)
concentrations or shear banding. Once again, let us state that ‘‘complex
material response’’ in the context of this text refers to anything beyond
linear elasticity.
Let us also remark that the searches become more involved and
time-consuming as the dimensionality of the space increases. From
1-dimensional problems requiring only one component of stress and
strain to define the per-element 2-dimensional phase space, to 3-
dimensional elements requiring six of each in a 12-dimensional phase
space, the resources invested in performing the projection onto
D
scale
rapidly. Of course, the size of
D
itself also plays a major role, but this
issue can be offset via smart searching procedures (Kirchdoerfer and
Ortiz, 2016; Eggersmann et al., 2021a), e.g., employing the k-d tree
algorithm.
2.3.2. Implementationdetails
For maximum efficiency, the raw material dataset should be sifted
to get rid of points that represent linear response. For the upcoming
simulations, we defined limit stresses
휎
lim
in terms of some combination
of the stress components (e.g., mean stress) and chose to include
each datum
풛
푒
whenever
푓
(
풛
푒
)
>
0
.
8
휎
lim
, where
푓
represents the
aforementioned stress combination. The reason for using a relaxed
inclusion is leaving enough wiggle room in the dataset to account also
for almost-non-linear response.
Likewise, to take into account the uncertainty associated with the
definition of the threshold itself, we choose to switch to DD whenever
푓
(
풛
)
>
0
.
9
휎
lim
.
When a regular element is swapped by a DD one, the strain–stress
state of the latter can be initialized in different ways (Fig. 1).
•
Method # 1: The last state of the regular element, call it
풛
퐹퐸
푒
, can
be swapped by its closest point in
D
, i.e., by
풛
∗
푒
=
argmin
푑
(
풛
퐹퐸
푒
,
D)
, or, using projection notation
풛
∗
푒
=
푃
퐷
풛
퐹퐸
푒
. This approach
initializes the element at a proper location in phase space, at the
expense of performing an extra search.
•
Method # 2: Simply initialize the new DD element state at the
origin, i.e.,
풛
∗
푒
= 0 ∈
R
2
푁
푐
. This method has the advantage of
avoiding the initial projection, but could cause the element to
follow a phase-space trajectory different from the one in method
#1 in subsequent steps. This approach appears more adequate
for history-dependent materials, but these are out of the scope of
this text. Herein, we will show that this initialization procedure
is faster (since it avoids extra phase-space searches) but can lead
to excessive refinement when used without adequate incremental
loading.
Incremental loading is unavoidable in some circumstances as, e.g.,
large deformation analysis (to recompute geometrical stiffness changes
associated to shape evolution), in the presence of dynamics (time
marching), or inelasticity (history-dependent material response). Nev-
ertheless, when dealing with quasi-static elastic yet non-linear re-
sponse, the load could be applied at once and the phase-space searches
would take care of finding the most convenient element states. We will
show that this is indeed the case, but that the final solution can depend
on the choice of aforementioned initialization methods.
We are assuming monotonously-increasing loading throughout this
text. If the material exhibits inelasticity, elements that have been re-
fined should not be swapped back to linear FE as the deformation they
suffered is irreversible. However, if the material is non-linear elastic,
one could define a threshold under which a DD-informed element is
converted back to a linear FE element.
The general version of the d-refinement solution procedure is devel-
oped in Algorithm 2.
3. Verification
Two verification exercises are presented in this section. For both,
an underlying constitutive law – from which a synthetic database
is generated – is assumed to exist, as this allows reaching solutions
with a traditional iterative solver, here specifically Newton–Raphson
(NR). This iterative solution is then considered as exact in order to
estimate the accuracy of the data-driven or d-refined solutions. The
two cases are used to study different behaviors of the method. In
the first one, constant-stress triangular (CST) elements are used under
the plane stress assumption. The main topics of interest are the wall
time speed-up, the spread of the refinement and the influence of the
initiation of newly refined elements. The second one is a 3D truss-
beam discretized with axial bar elements. This is used to study the
accuracy and convergence of the method with respect to the number
of datapoints, the influence of noise in the dataset and the impact of
the load stepping procedure.
3.1. 2Dplane-stresselements:platewithhole
We study first the capabilities of the one-load-step d-refinement
procedure. We create a dataset using the per-element phase-space
trajectories of the incremental NR solver, and then (a) we compare
to NR solver in terms of simulation speed assuming no incremental
loading in the d-refinement case. The effect of changing the number of
elements in the dataset is considered. (b) We also verify the accuracy
Mechanics of Materials 180 (2023) 104630
5
S. Wattel et al.
Algorithm2
D-refinement framework
Require:
∀
푒
= 1
,
...
,푁
푒
, compatibility matrices
푩
푒
,
D
푒
matrices
containing elastic constants, local datasets
D
푒
and
C
푒
matrices con-
taining distance constants;
∀
푖
= 1
,
...
,푁
dof
, final external forces
퐹
푖
or
boundary conditions.
(i) Set
푗
= 1
,
푆
(1)
2
= ∅
and
푆
(1)
1
= {1
,
2
,
3
,
...
,푁
푒
}
.
(ii) Incremental loading. Initialize
푙
= 1
for
푙
≤
푁
steps
do
Set loading level
풇
=
푙
∕
푁
steps
푭
if
|
푆
(
푗
)
2
|
= 0
then
Set
휼
= 0
and solve for
풖
in eq. (8b).
else
(ii.a) Initialize
푘
= 0
(ii.b) Project onto
E
:
Solve in eq. (8) for
풖
(
푘
)
and
휼
(
푘
)
for
푒
∈
푆
2
(
푗
)
do
Compute
흐
푒
and
흈
푒
from eq. (9)
endfor
(ii.c) Project onto
D
:
for
푒
∈
푆
(
푗
)
2
do
Choose
풛
∗(
푘
+1)
푒
∈ D
푒
closest to
풛
(
푘
)
푒
∈ E
푒
endfor
(ii.d) Convergence test:
if
풛
(
푘
+1)
푒
=
풛
(
푘
)
푒
∀
푒
∈
푆
(
푗
)
2
then
Final displacement field
풖
=
풖
(
푘
)
goto
(iii)
else
푘
←
푘
+ 1
,
goto
(ii.b)
endif
endif
(iii) Check for linear elements above threshold:
Set
푆
(
푗
+1)
1
=
푆
(
푗
)
1
,
푆
(
푗
+1)
2
=
푆
(
푗
)
2
for
푒
∈
푆
(
푗
)
1
do
Compute element strain and/or stress:
흐
푒
and
흈
푒
if
threshold criterion
then
Append
푒
to
푆
(
푗
+1)
2
, drop it from
푆
(
푗
+1)
1
Initial assignation: either choose
풛
∗(
푘
)
푒
∈ D
푒
closest to
(
흐
푒
,
흈
푒
)
or pre-defined
풛
∗
푒
∈ D
푒
endif
endfor
(vi) Global convergence condition
if
|
푆
(
푗
)
2
|
=
|
푆
(
푗
−1)
2
|
then
All DD elements converged (or no need for DD), no need of
further refinement
exit
else
푗
←
푗
+ 1
endif
푙
←
푙
+ 1
endfor
of d-refinement via comparison to NR, and (c) compare the two ini-
tialization methods (no incremental loading in the d-refinement case)
in terms of accuracy (comparison to NR) and computational time. For
each initialization method, the proportion of elements that become DD
and how many of them reach the same state as in the NR solution is
considered. Lastly, (d) we assess the impact of widely distributed non-
linearity (i.e., higher load, and hence less linear elements and more
refinement).
The phase-space searches are carried out in parallel using six pro-
cessors as long as the total number of DD elements is above twelve. This
problem is implemented in a Mathematica notebook (Wolfram, 2000),
which can be downloaded from on-line repositories (see supplementary
material section).
3.1.1. System
We resort to a classic 2D benchmark exercise: a thin square plate
(side length
퐿
= 0
.
2m
) with a circular hole (radius
푟
= 0
.
02m
) loaded
in tension on two opposite edges, see Fig. 3(a). The applied traction
푝
equals
100MPa
.
Material.
The material the plate is made of is assumed to change (de-
crease) stiffness when the mean stress (
휎
푚
= (
휎
I
+
휎
II
)∕2 = (
휎
xx
+
휎
yy
)∕2
)
surpasses a limit value
휎
lim
. No inelastic behavior is considered in these
simulations, deformations are considered recoverable; this does not
play out anyway since we study a quasi-static monotonously-increasing
loading scenario. The Young’s modulus of the material is assumed to
change according to an underlying model (Fig. 2):
퐸
(
휎
푚
) =
⎧
⎪
⎨
⎪
⎩
퐸
0
if
휎
푚
< 휎
lim
max
[(
휎
lim
휎
푚
)
퐸
0
,
0
.
5
퐸
0
]
if
휎
푚
≥
휎
lim
,
(10)
퐸
0
being the zero-strain modulus, while the Poisson’s ratio is assumed
to remain constant. For the simulations we are to show,
퐸
0
= 200GPa
,
휈
= 0
.
33
and
휎
lim
= 75MPa
. This model is ‘‘invisible’’ to the data-based
solvers, but not for the NR one. We use
C
푒
=
D
푒
=
D
∀
푒
, where
D
=
퐸
0
1 −
휈
2
⎡
⎢
⎢
⎣
1
휈
0
휈
1 0
0 0
1−
휈
2
⎤
⎥
⎥
⎦
.
(11)
We note that models of this type have been proposed as a means of
quantifying the effect of crack shielding by microcracking in brittle
materials such as ceramics (Hutchinson, 1987; Ortiz, 1987).
The mesh is generated with CST elements. The system is solved first
with linear-elasticity elements to check that at least some elements,
but not all, surpass the limit stress. A simulation where every element
remains linear or becomes non-linear would not be an interesting
case study for d-refinement. After this first verification, we move
to solve it also using pure DDCM, with d-refinement and with the
Newton–Raphson (NR) solver.
At each loading level, the incremental NR solver convergence condi-
tion is
|
풇
ext
−
풇
int
|
∕
|
풇
ext
|
<
tol
, where
풇
ext
is the nodal external force
vector and
풇
int
the internal one for unrestrained degrees of freedom.
The total external force is applied linearly over ten steps.
The datasets used for this first comparison proceed from the solution
of the incremental solver: at the end of each loading step, the strain
and stress state of each element is recorded as a point in phase space
that becomes part of the material set
D
according to the criterion in
Section 2.3.2.
3.1.2. Solutionanalysis
Qualitatively, Fig. 3(b), the method performs as expected: the top
frame displays the elements that overcome the mean-stress limit in
the linear-elastic simulation (vertically-hatched elements), while solid-
color elements in the same figure show the spread of non-linearity
predicted by the NR solution. The latter set is larger and contains most
of the former. The stress is more widely distributed over the domain
once the non-linear softer material is accounted for, what pulls more
elements past the threshold than those presumed by the linear-elastic
solver. The lower frame shows, in addition to the aforementioned NR
elements that go non-linear, the elements that become data-driven ac-
cording to d-refinement (horizontal hatching). See how these subsume
completely the others, meaning that this solver does predict where non-
linear behavior needs to be accounted for and proceeds in accordance
to implement the DD formulation therein.
Walltimespeed-up.
Regarding the results presented next, let us remark
that the time difference between realizations of the same simulation
differ by about than a
±10%
(maximum), in particular those solved with
pure DDCM, in which the initial assignation is random.
We begin by comparing the performance of the new approach
against NR solver for different values of
tol
while raising the load up to
Mechanics of Materials 180 (2023) 104630
6
S. Wattel et al.
Fig. 2.
Visualizing the material response of model Eq.
( 10). (a) Stiffness loss as mean stress increases. (b) Stress–strain relation for uniaxial loading (cf. Ref.
Ortiz
( 1987)).
Fig. 3.
System and simulation’s main outcome (results corresponding to tol
= 10
−3
and |D|
=
886 elements). (a) Scheme of the system, including loading and reduced mesh (taking
advantage of symmetry). (b) Results: elements above
휎
lim
in linear-elastic simulations (vertically-hatched gray), those above the same threshold in Newton–Raphson simulations
(solid pink) and elements that switched to DD in the course of the simulation (diagonally-hatched blue). (For interpretation of the references to color in this figure legend, the
reader is referred to the web version of this article.)
its final value,
푝
= 100
MPa, over ten increments at constant rate. Note
that this final value would yield a state of mean stress
휎
푚
≈ 50
MPa
<
휎
lim
away from the hole.
We take advantage of the most accurate NR simulation (
tol = 10
−5
)
to generate the dataset as described above and of the linear-elastic
solution to set the method constants:
C
푒
=
D
∀
푒
∈
푆
2
. The set
D
contains
at first 29240 data (including information about elastic loading), and
it reduces to just 886 (that only include information within the non-
linear regime) after the sifting (Section
2.3.2
). Let us re-iterate that
this convenient reduction is uncontroversial under the assumption of
monotonously-increasing load, but new considerations enter the picture
if one wants to consider non-monotonous loading (either elastic or
inelastic); for instance, one would have to take decisions as to how to
deal with DD elements that revert to the elastic regime.
The results concerning computational speed-up are consigned in
Table
1 . The d-refinement technique is faster than both NR solver (for
all the precision tolerances tried) and pure DD solution. For the latter,
a simulation with a reduced material set (
|
D
|
= 4598
, including linear
response) was sampled at random from the NR phase-space trajectories
(29240 points).
The evolution of this particular d-refinement simulation is repre-
sented in
Fig. 4 . All the elements that were presumed by the linear-
elastic solver to be above the mean-stress threshold were turned into
DD after the first iteration. Of course, all these elements were assigned
new dataset points, what is reflected in the other (red) curve first point.
As the simulation moved forth, the number of DD elements converged
quickly, while the material point assigned to each one of them did so
more slowly.
Mechanics of Materials 180 (2023) 104630
7
S. Wattel et al.
Fig. 4.
(a) Simulation progress in terms of portion of the mesh being refined and number of DD elements that changed material data by the end of each step. (b) Linear scaling
of d-refinement wall time (
훥푡
푑
−
푟푒푓
) with size of dataset (for a similar number of refined elements): the solid line (Method #1) represents initial assignation
풛
∗
푒
=
argmin
푑
(
풛
퐹퐸
푒
,
D)
(202 refined elements by the end of the simulation), the dashed line (Method #2)
풛
∗
푒
=
ퟎ
(270 refined elements on average). (For interpretation of the references to color in this
figure legend, the reader is referred to the web version of this article.)
Accuracyofd-refinement.
To assess the difference between the solution
using mesh d-refinement and the most accurate incremental solver
(tol
= 10
−5
), we choose to compare the phase-space distance between
the points corresponding to the Newton–Raphson solution (
풛
NR
∈
R
푁
푒
×2
푁
푐
,
푁
푐
= 3
for plane stress) and the ones of the d-refinement
solution (
풛
d-ref
∈
R
푁
푒
×2
푁
푐
):
|
풛
d-ref
−
풛
NR
|
|
풛
d-ref
|
=
(
∑
푁
푒
푒
=1
푤
푒
|
풛
d-ref
푒
−
풛
NR
푒
|
)
1∕2
(
∑
푁
푒
푒
=1
푤
푒
|
풛
NR
푒
|
)
1∕2
≈ 0
.
034
,
(12)
that is, the distance between solutions in phase space is less than 4% of
the one from the origin to the d-refinement solution. This is in spite of
the material set containing exactly the same points as the non-linear
portion of the Newton–Raphson trajectory. Actually, if we compare
where the refined elements land, only 23 out of 202 arrive at their final
configuration in the incremental solver. This does not translate in error
nonetheless (recall that the error we have estimated using distances is
less than a 4%), which necessarily means that the local density of the
material set in phase space is such that some elements can converge to
nearby points at virtually no cost.
For this kind of initialization (
풛
∗
푒
=
푃
퐷
풛
퐹퐸
, labeled ‘‘method #1’’), it
is interesting to note that the proportion of elements that become data-
driven is either completely independent of or very weakly-dependent
on the phase-space density. Increasing the precision of the incremental
solver by raising the number of load steps from 10 to 25 increases
the material dataset size to
|
D
|
= 1878
points, and this has no effect
either over the number of elements that become DD (202) or on the
proportion of these that reach the same destination (23). In terms of
computational time, Table 2, d-refinement again outperforms the most
accurate NR (
tol = 10
−5
).
Effect of dataset size on d-refinement performance.
This also motivates
a sensitivity study as to the influence of the size of the material set
D (number of material data): the incremental solver is kept at ten
load increments and
tol = 10
−3
, but we grow
D
by lowering the
threshold to facilitate points in the Newton–Raphson trajectory entering
D. We find the expected linear proportionality between total wall time
and dataset size, see solid line in Fig. 4. We discuss the influence of
the other initial data assignation (which happens to reduce time, see
dashed line in Fig. 4, at the expense of accuracy) in the next section,
but both methods display this linear scaling. To further decrease the
computational time, one could use some of the techniques discussed in
the literature (Eggersmann et al., 2021a,b).
Influence of initial data assignation.
Until now, reported results used
phase-space searches to assign an initial datum to each element that
evolves from FEM to DD. This involves additional searches that could
be avoided by always assigning a pre-defined initial datum, i.e., the
origin.
This creates a somewhat disjoint material set
D
, since there is a
single point away from the cluster that represents non-linear behavior.
The approach brings about a substantial speed improvement: based
on the results presented in Fig. 4, we find that this initial assignation
method (‘‘method #2’’) reduces the computation time by about 45%.
This comes at the expense of increasing the error when compared to
NR results: the distance ratio Eq. (12) doubled from 3.5% to 7%. So, as
often, there is here a trade-off between accuracy and time efficiency.
Depending on the priority, one would favor one method over the
other. However, this conclusion is contingent on localized non-linear
behavior, what leads to having similar final number of DD elements,
see next point.
Impactofspreadingnon-linearbehavior.
As the load increases, the actual
solution departs more meaningfully from the linear-elastic one. As
elements excursion deeper past the threshold, richer datasets are also
necessary to capture the response.
For instance, keeping the dataset inclusion criterion and increasing
by 20% the load we have been using leads to an increase in
|
D
|
from
886 to 3419 phase-space points.
This also increases both the NR computation time and the d-
refinement’s, as well as the relative error between them (Eq. (12)), see
Table 3. The spreading of the non-linearity induced by the intenser
load means more DD elements and hence more phase-space searches,
what slows d-refinement down to the point of NR being faster in this
case. The final number of refined elements has escalated from 202 to
620. Interestingly, initialization type 2 (
풛
∗
푒
=
ퟎ
) ends up defaulting to a
pure data-driven mesh due to accumulation of errors and the total time
skyrockets.
We remarked that the precision of either approach can be increased
by sharpening the refinement threshold criterion. In any case, such
an approach is a barren one, since in practice the actual threshold
will not be well-defined; the difference between d-refinement and the
incremental solver could be regarded not as a demerit of the former but
as modeling bias infused in the latter.
The main conclusion is that in quasi-static, monotonous-loading
cases where non-linearity is expected to be present well over the do-
main, d-refinement must be used with initialization #1, since this yields
Mechanics of Materials 180 (2023) 104630
8
S. Wattel et al.
Table 1
Solving the plane-stress problem with different methods: wall time comparison (
푝
=
100
MPa, initialization type 1, 10 NR increments).
훥푡
d-ref
(
|
D
|
=
886)
훥푡
NR
(
tol = 10
−3
)
훥푡
NR
(
tol = 10
−4
)
훥푡
NR
(
tol = 10
−5
)
훥푡
DD
(
|
D
|
=
4588)
[
푠
]
20.3
38.5
51.2
72.8
928.7
Increment –
× 1
.
9
× 2
.
5
× 3
.
6
× 45
.
8
better first guesses, thus boosting convergence and limiting the number
of extra elements to be refined. To gauge the extent of the non-linear
portion beforehand, comparing the results of a linear-elastic calculation
to the threshold may suffice. We have shown that initialization method
#2 is adequate and efficient when there is localized non-linearity and,
in the next section, we will also prove its suitability when used in
tandem with appropriate load increments.
3.2. 1Dbarelementsin3Dspace:octet-trussstructure
In this study, we create a dataset by sampling a constitutive law,
and compare d-refinement and NR solver (a) changing the density of
the dataset, and (b) changing the number of loading steps.
This application is implemented in a Jupyter notebook, which is
provided as supplementary material.
3.2.1. System
The second study features a truss made of octet unit cells (
Desh-
pande et al.
, 2001
). The model is made up by axially-loaded bar ele-
ments (no bending), with linearized kinematics, and a total of 6
×
2
×
2
unit octet truss cells submitted to a 3-point bending test as represented
in Fig. 5 . The parameters of the unit cell octet truss were inspired by
the work in Ref.
Shaikeea et al.
( 2022) and are summarized in
Table
4 .
It should be noted that the goal is not to reproduce their results: to do
so, either buckling in the case of slender member or joint stiffness for
thick member should have been accounted for. This model is intended
to test and explore the performance of the d-refinement solver. It has
two main advantages: bar elements are straightforward to discretize
and their local phase-space is two-dimensional:
풛
푒
= (
휖
푒
,휎
푒
) ∈
R
2
.
The beam displacements are fixed on both lower edges and a
displacement-controlled deflection is imposed on the top middle nodes.
While the load is applied in steps, the problem is considered to be quasi-
static and thus solved with the static solver presented above. See
Fig. 5
for a graphical representation.
Material.
The behavior of the bulk material is also inspired from
Ref. Shaikeea et al.
( 2022), in which trimethylolpropane triacrylate
(TMPTA) is used. The material displays similar Young’s modulus at
zero-strain
퐸
0
= 430
MPa and limit strength
휎
푓
= 11
MPa. However,
this artificial material is assumed to be non-linear elastic, as opposed to
brittle, and is represented by a hyperbolic-tangent stress–strain relation
(see Fig. 6 ):
퐸
(
휖
) =
퐸
0
[
1 − tanh
2
(
퐸
0
휎
푙
휖
)]
,
(13)
휎
(
휖
) =
휎
푓
tanh
(
퐸
0
휎
푙
휖
)
.
(14)
A synthetic database of 2700 points is generated by sampling uni-
formly the constitutive law between
휖
= −0
.
2
and
휖
= 0
.
2
. For the
purpose of d-refinement, the trigger to leave the linear elastic domain
is assumed to be
휎 > 휎
푙
= 7MPa
. Fig. 6 summarizes the material
behavior, with the underlying constitutive law, the dataset and the
assumed linear elastic domain.
Table 2
Solving the plane-stress problem: wall time comparison (
푝
= 100
MPa, initialization type
1, 25 NR increments).
훥푡
d-ref
(
|
D
|
=
1858)
훥푡
NR
(
tol = 10
−5
)
[
푠
]
55.3
122.3
Increment
–
× 2
.
2
Table 3
Solving the plane-stress problem: wall time comparison and accuracy (
푝
= 120
MPa,
different initialization types, 10 NR increments). The NR solution is used as reference
to compare in this case.
훥푡
d-ref
(
|
D
|
=
3419)
훥푡
NR
(
tol = 10
−5
)
Method #1
Method #2
[
푠
]
154.8
924.3
105.3
Increment
× 1
.
5
× 8
.
8
–
Error (Eq.
( 12))
0
.
05
0
.
14
–
Table 4
Summary of the different parameters used in the truss models with both the geometric
parameters and the bulk material properties.
Notation
Value
푙
530
μ
m
Strut length
푑
65
μ
m
Strut diameter
퐸
0
430
MPa
TMPTA Young’s modulus at zero strain
휎
푓
11
MPa
TMPTA tensile strength
Fig. 5.
Scheme of three-point bending configuration. Arrows represent the imposed
displacement at the top middle nodes.
Fig. 6.
Constitutive law, synthetic dataset and the limits of the linear elastic domain.
For the dataset, only 1% of the points are plotted to increase legibility.
3.2.2. Solutionanalysis
The following is a comparative study of the results obtained with
three different solvers: a Newton–Raphson informed by the constitutive
law (NR), a pure data-driven one informed by the dataset from the
constitutive law (DDCM) and a d-refinement one informed by the
dataset as well as the zero-strain Young’s modulus and linear elastic
limit (DR). Initialization method #2 is chosen in this case. Unlike the
Mechanics of Materials 180 (2023) 104630
9
S. Wattel et al.
Fig. 7.
Load–deflection curves (upper) obtained with the three solvers and a 20-step
loading protocol (left vertical axis). The deflection is normalized by the beam length.
The total force is computed through the internal forces in the bars that converge at
the nodes where the displacement is enforced. The lower curves (right vertical axis)
represent the fraction of rod elements that have been refined, i.e. switched to DD, for
both the 20-step loading and the 2-step loading. (For interpretation of the references
to color in this figure legend, the reader is referred to the web version of this article.)
previous plane-stress study, the incremental version of d-refinement is
used here.
Error analysis.
The load response of the beam for each solver is
recorded to be then compared. With a 2700-point dataset (
|
D
|
= 2700
)
and taking the Newton–Raphson solution as a reference, the relative
errors in load response are 2.8% and 1.8% for the DD and DR solvers,
respectively.
Convergencewithnoiselessandnoisydataset.
Then, the number of points
in the dataset is increased. The errors, as defined by Eq. (12), of
both the DD and DR solvers reduce at a rate given by a power law
with exponent one. This
∼
|
D
|
−1
scaling is coherent with the litera-
ture (Kirchdoerfer and Ortiz, 2016). However, for the DR solver, the
convergence to the NR solution is ultimately limited by the quality
of the linear assumption at small strains. The linear approximation
does not overlap perfectly with the underlying hyperbolic-tangent law,
hence some minor difference in the linear-elastic elements persists
independently of the number of datapoints. We stress that this is not a
flaw in the d-refinement method but a matter of modeling bias.
When random Gaussian noise is added to the datapoints (both in
strain and stress) with standard deviation inversely proportional to the
square root of the number of datapoints (Kirchdoerfer and Ortiz, 2016),
the same behavior is observed but the convergence rate slows down to
|
퐷
|
−1∕2
which is again in agreement with the literature results (Kirch-
doerfer and Ortiz, 2018, 2017). For data whose noise is independent
of the number of datapoints, max-ent data-driven solvers (Kirchdoerfer
and Ortiz, 2017) are more adequate than the fixed-point one presented
here. As for a more formal proof of convergence, the linear elastic
approximation can be seen as being part of the dataset (as a graph or
a sampling of a line with infinite density for example) and thus the
convergence proof derived in Ref. Conti et al. (2018) can be extended
to the d-refinement solver: if the linear approximation and the dataset
converge to the real underlying material behavior, so will do the
d-refinement solver.
Influence of load-stepping.
As the stepping procedure unfolds and the
load increases, more elements go over the threshold and are switched
to the DD solver. The size of the loading step influences which elements
are switched: a large step may cause the refinement to overshoot, i.e., to
refine elements that should have remained within the linear elastic
regime. The red lines (lower) in Fig. 7 highlight this phenomenon:
the fraction of refined elements at the end of the loading is 7% with
20 steps and 10% with 2 steps. In Fig. 8, with 20-step loading, the
refined elements correspond exactly to the ones that undergo large
straining according to the NR solution. In 2-step loading, more elements
are refined. This overshooting, associated to using initialization at the
origin and not enough load steps, does not influence the accuracy of
the final solution provided that (a) the dataset has points in the elastic
zone, at least in the proximity of the limit as suggested before, and
(b) the refined elements represent a fraction of the total (recall the
results obtained in the previous case study when increasing the load).
The only cost of the overshooting in refinement is computational, no
meaningful extra errors are induced. However, adding more loading
steps also have the cost of computing extra intermediate steps that
might not be of interest, so a balance should be met depending on the
type and purpose of the simulation. It should be noted that the standard
Newton–Raphson solver could diverge if the load steps were too large,
thus another advantage of the d-refinement method is its stability.
4. Application:bridgingscalestostudythefractureprocesszone
inarchitectedmetamaterials
While DDCM allows to use experimental data to make predictions
without resorting to phenomenological models, it can also be applied
as an efficient method to solve multiscale problems (Karapiperis et al.,
2021, 2020; Korzeniowski and Weinberg, 2022). It is particularly suited
in the case where a convincing microscale model exists, which, for
any reason, cannot be run at the larger scale of interest. In such a
case, a dataset consisting of RVE microscale simulations over different
loading paths can be computed offline and then passed to a macroscale
DDCM model. Since the overhead associated to running microscale
simulations is much lower than the one of performing experiments, the
development of the dataset can be done in several iterations: areas of
the phase space which have been found to not be dense enough can be
repopulated (Karapiperis et al., 2021); one can even foresee launching
the RVE microscale simulations on the fly (Karapiperis et al., 2021).
To showcase the merits of DDCM for multiscale modeling, the two
models presented in Section 3 are combined to devise an architected
metamaterial (Meza et al., 2015). Starting from the synthetic dataset
of TMPTA material behavior, microscale simulations of an RVE unit
cell of octet truss were run using the DDCM solver. Then, this dataset
was fed to a continuous macroscale model with the goal of resolving
the stress distribution around a crack tip with the d-refinement solver.
In summary, we will go from a dataset of the underlying material
to the response of the microstructure, to finally reach the macroscale
and perform mechanical analysis where the non-linear behavior of the
octet-truss architected metamaterial is taken into consideration with no
recourse to any constitutive modeling whatsoever.
This family of artificial materials has commanded the attention of
the research community in recent decades (Fleck et al., 2010). Due
to their manufacturability and remarkable properties, they are poised
to challenge conventional materials in a bevy of applications, from
thermal insulation (Dou et al., 2018) to impact absorption (Portela
et al., 2021).
4.1. System
Microscale.
The microscale model is the same used in Section 3.2. The
RVE is a single unit octet-truss cell (represented in the middle of Fig. 9).
The axes of the frame of reference (
풙
,
풚
,
풛
) are aligned with the lattice
planes ([100], [010], [001]), see Fig. 10(a). The cross-sections of the
bar elements that lie on the boundary are reduced as they are shared
with the neighboring cells.
To create the dataset, the stress response to specified strain states
(
휖
푥푥
휖
푦푦
휖
푥푦
) are simulated. For this, the strain state is converted into
imposed displacements at the ten boundary nodes. To compute the
homogenized stress states, the sum of the forces, on each facet, for each
direction, is divided by the apparent area of the cell facet. The material
behavior is again considered elastic. Since it was verified previously
that the accuracy of the final solution does not depend on the number
of loading steps when using initialization method #2, no incremental
loading is employed. The material response of the unit cell is illustrated
in Fig. 10 where strain–stress paths are plotted for uni-axial elongation,
pure shear and combined shear-stretching.