# Physically recurrent neural networks for path-dependent heterogeneous
materials: embedding constitutive models in a data-driven surrogate

M. A. Maia  I. B. C. M. Rocha  P. Kerfriden  F. P. van der Meer  Delft
University of Technology, Department of Civil Engineering and Geosciences,
P.O.Box 5048, 2600GA, Delft, The Netherlands  Mines ParisTech (PSL
University), Centre des matériaux, 63-65 Rue Henri-Auguste Desbruères BP87,
F-91003, Évry, France  Cardiff University, School of Engineering, Queen’s
Buildings, The Parade, Cardiff, CF24 3AA, United Kingdom

###### Abstract

Driven by the need to accelerate numerical simulations, the use of machine
learning techniques is rapidly growing in the field of computational solid
mechanics. Their application is especially advantageous in concurrent
multiscale finite element analysis (FE2) due to the exceedingly high
computational costs often associated with it and the high number of similar
micromechanical analyses involved. To tackle the issue, using surrogate models
to approximate the microscopic behavior and accelerate the simulations is a
promising and increasingly popular strategy. However, several challenges
related to their data-driven nature compromise the reliability of surrogate
models in material modeling. The alternative explored in this work is to
reintroduce some of the physics-based knowledge of classical constitutive
modeling into a neural network by employing the actual material models used in
the full-order micromodel to introduce non-linearity. Thus, path-dependency
arises naturally since every material model in the layer keeps track of its
own internal variables. For the numerical examples, a composite Representative
Volume Element with elastic fibers and elasto-plastic matrix material is used
as the microscopic model. The network is tested in a series of challenging
scenarios and its performance is compared to that of a state-of-the-art
Recurrent Neural Network (RNN). A remarkable outcome of the novel framework is
the ability to naturally predict unloading/reloading behavior without ever
seeing it during training, a stark contrast with popular but data-hungry
models such as RNNs. Finally, the proposed network is applied to FE2 examples
to assess its robustness for application in nonlinear finite element analysis.

###### keywords:

Artificial Neural networks (ANN) , Multiscale , Heterogeneous materials ,
Path-dependency

##  1 Introduction

Driven by the need to accelerate numerical simulations, machine learning (ML)
techniques are rapidly growing in the field of computational solid mechanics.
The use of surrogate models in particular is especially advantageous in
concurrent multiscale finite element analysis (FE2). In this approach, each
integration point at the macroscale is linked to a microscopic model, allowing
complex materials such as composite materials to be explicitly modeled using
relatively simple constitutive models. Although appealing, the computational
cost associated with the concurrent Finite Element (FE) simulations is often
prohibitive, hindering its widespread use in practical engineering-scale
applications. In that sense, another layer of approximation, for instance
based on ML techniques, that help accelerate FE2 simulations is key.

However, despite the successful use of this type of technique in several
fields (e.g., speech recognition, language processing, and biomedical
sciences), the use of surrogate models in material modeling is still filled
with challenges and open issues. One of them relates to the data scarcity that
exposes how data-driven models do not perform as well in extrapolation as they
do within their training range. In general, another shortcoming of black-box
models is their lack of interpretability. Although not a limiting aspect by
itself, the limited interpretability these models offer can make it difficult
to understand whether the fitted model is capable of yielding physically
consistent solutions for loading paths different than those seen during
training.

In this scenario, the definition of a Design of Experiments (DoE) strategy
plays an important role in the generalization capabilities of the surrogate
model. The designer needs to take into account the computational (or
experimental) budget available to collect the data and what type of loading
scenarios are expected to be experienced by the model. The former is known a
priori, but the latter is not known exactly until the finite element
simulation itself is run. As a consequence, large training and testing
datasets are usually employed in an attempt to cover all different loading
scenarios for a given strain/stress range.

When path-dependency is present, conceiving an efficient DoE becomes even more
complex since stresses now depend on the history of the material, leading to
potentially infinite parameter space. One way to bypass the computationally
intensive approach is to run and incorporate new load paths on-the-fly,
updating the surrogate model as necessary as in the works of Goury et al. [1]
and Rocha et al. [2]. In addition to the definition of the DoE, the physics-
based assumptions about the material should also be wisely taken into account
when choosing the modeling approach. For instance, conventional Artificial
Neural Networks (ANN) with strains and stresses taken as input and output,
respectively, will fail to capture phenomena such as elastic unloading. This
scenario is illustrated by Vlassis and Sun [3] and occurs due to the unique
mapping between inputs and outputs that does not reflect two different stress
states for the same strain. One way to differentiate the different paths is by
augmenting the feature space of the ANN with extra variables that carry
partial information about the history of stresses and/or strains [4, 5].

Another highly popular strategy to handle path-dependency is to use Recurrent
Neural Networks (RNNs), a variation of ANNs capable of learning from
sequential data. This is done by incorporating the previous state of the
network when making predictions for the current state. However, this approach
often suffers from short-term memory when dealing with long sequences. To
address that issue, more complex architectures with more model parameters and
fine control of the flow of information retained through a sequence of
information were proposed (e.g., GRU and LSTM). These networks gained
popularity and have been used to model a wide variety of materials that show
history-dependency [6, 7, 8, 9, 10, 11].

However, RNNs are still severely limited by the curse of dimensionality
associated with sampling arbitrarily long strain paths. The overview on the
potential of RNNs for modeling path-dependent plasticity presented by Gorji et
al. [12] illustrates the multiple ways in which one can define a DoE and argue
that GRUs, in particular, can be used to model the plastic response of
materials provided a rich enough training dataset is used. However, the
authors do not investigate further how these networks would perform in
scenarios slightly different than those seen in training or their efficiency
in an FE2 framework.

Recently, a new class of ANNs called Physics-Informed Neural Networks (PINNs)
was proposed by [13] to solve any given laws of physics described by general
nonlinear partial differential equations. For training PINNs, in addition to
the data used to model the governing equation, the loss function is augmented
with information on the imposed physical constraints (e.g. initial and
boundary conditions). In this approach, automatic differentation plays a key
role in allowing the outputs of the network to be differentiated with respect
to the input.

Applications of PINNs as surrogate constitutive models are still relatively
new, but a few works showcase their potential in modeling elastoplastic [14,
15] and elastic-viscoplastic behaviors [16]. In [14] and [16], the loss
function is also incremented with terms related to the physics constraints
assumed by the material models considered (e.g. the Karush-Kuhn-Tucker
conditions on the yield function). While both works approximate the
displacements and the stress based on the position of the material point at
the macroscale and on the material properties, Haghighat et al. [14] further
investigates the applicability of their surrogate model in a sensitivity
analysis study. Finally, Eghbalian et al. [15] proposed an architecture where
the nonlinear incremental elasticity and the strain decomposition are
hardwired into the networks’s architecture and the loss function. Thus leading
to improved capability in predicting unseen loading scenarios over
conventional ANNs.

Inspired by PINNs, Masi et al. [17] tailored a network architecture capable of
ensuring thermodynamically consistent predictions by evaluating the numerical
derivatives of the network with respect to its inputs. Later, Masi and
Stefanou [18] proposed a key extension to the approach, in which the
identification of the internal variables is no longer user-dependent. For that
purpose, an encoder-decoder architecture is incorporated to identify the
minimum number of internal variables from the set of internal coordinates of
the system (e.g., displacement fields, internal forces, etc.) in an
unsupervised manner. The feature allows the recovery of the full-field
information of the microstructure by the decoder. Another strategy that takes
into account the derivatives of the approximated functions to impose
thermodynamically consistency is proposed by Vlassis and Sun [3]. The authors
treat the yield function as a signed distance function level set and formulate
a supervised learning task that deduces the learned evolving yield function
against a monotonically increasing accumulated plastic strain. As a result,
cyclic loading paths can be predicted based only on monotonic data. For a
recent and comprehensive review on the state-of-the-art literature of ANNs in
the constitutive modeling of composite materials, the reader is referred to
Liu et al. [19].

Different from regular ANNs or RNNs, the Deep Material Network (DMN) proposed
by Liu et al. [20] assigns physical meaning to the hyperparameters. The DMN
learns the hidden topology representation of the micromodel (RVE) based on the
stiffness matrices and residual stresses of the different material phases. A
major advantage of DMNs is the ability to extrapolate to nonlinear behavior
based only in elastic snapshots. Several follow-up publications and
improvements on DMN can be found in the literature. The author’s latest work
[21] is dedicated to tackling multiscale failure analysis, a far less explored
avenue by the community. Other succesful applications in this field can be
found in the works of Kerfriden et al. [22], Bessa et al. [23] and Oliver et
al. [24], where different reduced-order modeling strategies were employed.

In an alternative approach that incorporates some level of physics in the
surrogate model, Fuhg et al. [25] combines a physical model that accounts
partially for the constitutive behavior of the material and a Gaussian Process
(GP) correction to locally improve the baseline physical model. Although this
boosts accuracy inside the training space, in practice, it does not translate
to benefits in terms of extrapolation. Another alternative that retains some
of the physics embedded in the full-order solution is model order reduction.
These methods are frequently employed to alleviate the computational cost of
FE2 [26, 27] and usually offer better generalization properties to unseen
points than common surrogate models. The drawback is that they are inherently
slower [28]. Such techniques can also be coupled with RNNs [29], ANNs [5] and
GPs [30] to quickly infer the coefficients of the reduced basis approximation
for arbitrary parameter values.

Although a large body of literature has been devoted to applying ML techniques
in the solid mechanics field, the gaps and issues left by data-driven
surrogate models are still an open issue, compromising their reliable and
widespread use in practical applications. In this work, a new network design
is proposed specifically for the modeling of path-dependent materials and to
accelerate concurrent finite element simulations. In Section 2 the FE2 method
is presented, while in Section 3 one of the most popular approaches to tackle
the computational bottleneck originated from it is briefly discussed. In
Section 4, the main features of the novel neural network are described. In
Section 5, the Design of Experiments and methodology adopted for the
comparative study shown in Section 6 is described. In this study, the
performance of the proposed network is compared to a RNN for a single-scale
problem. In Section 7 the novel approach is integrated into an FE2 framework
and tested in two applications for robustness and accuracy. Finally,
conclusions are presented in Section 8.

##  2 Concurrent multiscale analysis

Let $\Omega$ define the macroscopic domain being modeled. To find the internal
stresses and displacement field of such body in absence of body forces, a
boundary value problem that satisfies the following equilibrium equations is
defined:

$\displaystyle\nabla\bm{\sigma}^{\Omega}=\mathbf{0}$ (1)

where $\nabla$ is the divergence operator and $\bm{\sigma}^{\Omega}$ is the
macroscopic stress, which depends on the macroscopic displacement field
$\mathbf{u}^{\Omega}$ (for simplicity, this dependence is omitted). The
governing equations are subjected to the boundary conditions:

$\displaystyle\bm{\sigma}^{\Omega}\ \mathbf{n}=\mathbf{t}^{\Gamma_{f}}\ \ \
\text{on}\ \Gamma_{f}\qquad\mathbf{u}^{\Omega}=\mathbf{u}^{\Gamma_{u}}\ \ \
\text{on}\ \Gamma_{u}$ (2)

where $\mathbf{n}$ is the normal to the surface $\Gamma_{f}$ and
$\mathbf{u}^{\Gamma_{u}}$ and $\mathbf{t}^{\Gamma_{f}}$ represent a set of
Neumann and Dirichlet boundary conditions acting on the body surface such that
$\Gamma_{u}\cap\Gamma_{f}=\emptyset$ as illustrated in Fig. 1a. To relate
strains and stresses, a constitutive model $\mathcal{D}$ is required:

$\bm{\sigma}^{\Omega}=\mathcal{D}^{\Omega}\
(\bm{\varepsilon}^{\Omega},\bm{\alpha}^{\Omega})$ (3)

where $\bm{\alpha}^{\Omega}$ are history variables that account for path-
dependency and $\bm{\varepsilon}^{\Omega}$ is the macroscopic strain defined
under small displacement assumptions as:

$\bm{\varepsilon}^{\Omega}=\frac{1}{2}\Big{(}\nabla\mathbf{u}^{\Omega}+(\nabla\mathbf{u}^{\Omega})^{\text{T}}\Big{)}$
(4)

In the concurrent multiscale approach, the model $\mathcal{D}^{\Omega}$ is not
directly formulated but is instead obtained by nesting a lower scale finite
element model to each integration point. In that scale, the microscopic
stucture of complex materials can be explicitly modeled using simpler
constitutive models for each of the components. Further discussion on how to
solve the microscopic problem and link both scales is shown in Section 2.1 and
2.2.

(a) Schematic representation of FE2

(b) Controlling nodes of RVE

Figure 1: Scheme of FE2 framework and definition of the boundary value problem
on RVE

To solve the boundary value problem at the macroscale, the Finite Element (FE)
method is employed to discretize the domain $\Omega$ into a number of elements
connected by nodes with $N$ degrees of freedom. The global equilibrium is
solved iteratively in its discretised weak form:

$\mathbf{r}=\mathbf{f}^{\Gamma}-\mathbf{f}^{\Omega}(\mathbf{u}^{\Omega})=\mathbf{0}$
(5)

where $\mathbf{r}\in\mathbb{R}^{N}$ is a residual vector that goes to zero
when equilibrium is reached, $\mathbf{f}^{\Gamma}\in\mathbb{R}^{N}$ is the
global external vector that represents the Neumann boundary conditions and
$\mathbf{f}^{\Omega}\in\mathbb{R}^{N}$ is the global internal force vector
given by a volume integral:

$\mathbf{f}^{\Omega}=\operatorname*{A}_{e=1}^{n_{e}}\int_{\Omega_{e}}\mathbf{B}_{e}^{\text{T}}\
\bm{\sigma}^{\Omega}(\mathbf{u}^{\Omega}_{e})\ \text{d}\Omega$ (6)

where $\operatorname*{A}$ is an assembly operator that takes into account the
connectivities between the elements and the global system and $\mathbf{B}$ is
a matrix with the spatial derivatives of the shape functions used to
interpolate nodal displacements. Finally, an iterative procedure is adopted to
solve Eq. (5) for the macroscopic displacement field:

$\Delta\mathbf{u}^{\Omega}=\mathbf{u}^{\Omega}_{n}-\mathbf{u}^{\Omega}_{o}=-\mathbf{K}_{o}^{-1}\mathbf{r}_{o}$
(7)

where the subscripts $o$ and $n$ refer to old and new analysis increments,
respectively and $\mathbf{K}\in\mathbb{R}^{N\times N}$ is the global tangent
stiffness matrix given by:

$\mathbf{K}=\operatorname*{A}_{e=1}^{n_{e}}\int_{\Omega_{e}}\mathbf{B}_{e}^{\text{T}}\
\mathbf{D}_{e}^{\Omega}(\mathbf{u}^{\Omega}_{e})\ \mathbf{B}_{e}\
\text{d}\Omega$ (8)

and $\mathbf{D}^{\Omega}$ is the constitutive tangent matrix, discussed in
Section 2.2.

The key difference to a classical FE simulation lies in the embedding of
another FE model in the macroscopic integration points. Here, to obtain the
internal forces in Eq. (6) and the tangent stiffness matrix in Eq. (8) for a
single integration point of the macroscale, one needs to run an entire FE
model instead of a single evaluation of a homogenous material model. This is
the most computationally expensive part of the framework and is where the
proposed network aims to tackle. The approach here is to replace the solution
to the microscopic problem (discussed in Section 2.1) with a surrogate model,
specifically a neural network. The homogenization procedure required to
upscale the responses to the macroscale is discussed in Section 2.2.

###  2.1 Microscopic scale

Let $\omega$ be a Representative Volume Element (RVE) of the microscopic
material features whose behavior is to be upscaled. Assuming that the
principle of separation of scales (i.e. $\Omega\gg\omega$) holds, the two
scales can be linked by enforcing:

$\mathbf{u}^{\omega}=\bm{\varepsilon}^{\Omega}\mathbf{x}^{\omega}+\tilde{\mathbf{u}}$
(9)

where the linear displacement field is the result of the imposed macroscopic
strains $\bm{\varepsilon}^{\Omega}$ and the fluctuation field
$\mathbf{\tilde{u}}$ is the result of microscopic inhomogeneities. The
principle of separation of scales implies the stain averaging theorem that
states that the macroscopic strains are considered uniform over the RVE
domain:

$\bm{\varepsilon}^{\Omega}(\mathbf{x}^{\Omega})=\frac{1}{\omega}\int_{\omega}\bm{\varepsilon}^{\omega}(\mathbf{x}_{\omega})\
\text{d}\omega$ (10)

where $\bm{\varepsilon}^{\omega}$ is the microscopic strain tensor. Therefore,
the microscopic displacement field in Eq. (9) can only satisfy (10) if the
fluctuation displacement field vanishes at the RVE boundary when upscaling
quantities. An additional requirement on the fluctation field having zero
resultant work at the boundaries arises from the Hill-Mandel principle. Both
requirements are met using Periodic Boundary Conditions (PBC) to represent the
behavior of a macroscopic bulk material point. Fig. 1b illustrates the node
groups and boundary edges needed to implement the PBC. In Section 5, the
generation of $\bm{\varepsilon}\text{-}\bm{\sigma}$ paths for the training of
the surrogate models is obtained by setting a user-defined function to set the
prescribed displacements $\mathbf{u}_{1}$ and $\mathbf{u}_{2}$.

Finally, keeping the hypothesis of small strains, the stress equilibrium
problem is described as:

$\displaystyle\nabla\bm{\sigma}^{\omega}=\mathbf{0}\qquad\bm{\sigma}^{\omega}=\mathcal{D}^{\omega}\
(\bm{\varepsilon}^{\omega},\bm{\alpha}^{\omega})\qquad\bm{\varepsilon}^{\omega}=\frac{1}{2}\Big{(}\nabla\mathbf{u}^{\omega}+(\nabla\mathbf{u}^{\omega})^{\text{T}}\Big{)}$
(11)

where $\mathbf{u}^{\omega}$ is the microscopic displacement field and
$\bm{\sigma}^{\omega}$ and $\bm{\varepsilon}^{\omega}$ are the microscopic
stress and strain tensors, respectively. An analoguous procedure to the one
detailed in Section 2 is used to find the microscopic displacement field
(subjected to the periodic boundary conditions). Note that at this scale,
regular physics-based material models $\mathcal{D}^{\omega}$ (e.g.
elastoplasticity, viscoelasticity, etc.) are employed to represent the
constitutive behavior of the homogeneous material of the discretized elements.

###  2.2 Homogenization procedure

After convergence of the microscopic displacement field $\mathbf{u}^{\omega}$,
the upscaling procedure is performed based on the Hill-Mandel principle. The
principle postulates that the macroscopic stress power must equal the volume
average of the microscopic power over the RVE. Considering the definition in
Eq. (9), an expression similar to Eq. (10) is obtained for the homogenized
stresses:

$\displaystyle\bm{\sigma}^{\Omega}=\frac{1}{\omega}\int_{\omega}\bm{\sigma}^{\omega}\text{d}\omega$
(12)

As for the macroscopic constitutive tangent stiffness $\mathbf{D}^{\Omega}$, a
probing operator $\mathcal{P}$ is applied on the global microscopic tangent
stiffness matrix $\mathbf{K}^{\omega}$ without the need to invert it as
proposed by Nguyen et al. [31].

##  3 Recurrent Neural Networks

In this section, a brief overview of the working mechanisms of Recurrent
Neural Networks is presented. Although part of the paper is dedicated to
comparing them with the novel approach, here the idea is to use well-known
concepts from ANNs and RNNs to illustrate features of the proposed network in
the following sections.

As a starting point, consider a conventional feed-forward neural network to
surrogate the nonlinear constitutive relationship of a path-independent
material given by the following parametric regression model:

$\widehat{\bm{\sigma}}^{\Omega}=\mathcal{NN}\
(\bm{\varepsilon}^{\Omega},\mathbf{W},\mathbf{b})$ (13)

where $\mathbf{W}$ and $\mathbf{b}$ are weights and biases calibrated through
a fitting procedure based on observations of the actual microscopic model.
During training, the strains are fed to the first neural layer (input layer)
and values are propagated until the final layer (output layer) to give the
predicted stresses $\widehat{\bm{\sigma}}^{\Omega}$. These are in turn
compared to the ground truth value according to a loss function. Based on
that, the model parameters are adjusted so that the error between the
predicted stresses and the actual stresses is minimized:

$\mathbf{W},\mathbf{b}=\textrm{argmin}\ \overline{\mathbf{W}}\
\overline{\mathbf{b}}\sum_{i\ \in\ \mathbf{X}}\|\
\bm{\sigma}_{i}(\bm{\varepsilon}^{\Omega}_{i})-\widehat{\bm{\sigma}}_{i}(\bm{\varepsilon}^{\Omega}_{i},\overline{\mathbf{W}},\overline{\mathbf{b}})\
\|^{2}$ (14)

where $\mathbf{X}\in\mathbb{R}^{n_{\varepsilon}\times N}$ is a snapshot matrix
with $N$ pairs of $\bm{\varepsilon}^{\Omega}$-$\bm{\sigma}^{\Omega}$ obtained
from microscopic simulations. This setting is the most straight-forward way to
map pairs of macroscopic strains and stresses but does not offer good
generalization properties once path-dependency is introduced. In that case,
one way to overcome the lack of history information is to extend their feature
space with e.g. previous (incremental) strains and/or stresses [4, 5].

As an alternative, RNNs offer additional parameters (i.e. the hidden state)
and mechanisms (i.e. the gates that control the flow of information being
propagated) to learn history information from sequential data in an implicit
way. These parameters describe the evolution of the so-called hidden state and
can encapsulate information from previous iterations without the need to
introduce history variables in the feature space. In a regular RNN, the
outputs and hidden state are given by:

$\displaystyle\mathbf{h}^{t}=\upphi(\mathbf{W}_{1}\mathbf{v}^{t}+\mathbf{W}_{s}\mathbf{h}^{t-1}+\mathbf{b}_{s})$
(15)
$\displaystyle\widehat{\bm{\sigma}}^{t}=\upphi(\mathbf{W}_{2}\mathbf{h}^{t}+\mathbf{b}_{2})$

where $\upphi(\cdot)$ is an activation function, $\mathbf{W}_{s}$ and
$\mathbf{b}_{s}$ are the additional model parameters (compared to conventional
feed-forward neural networks), $\mathbf{v}^{t}$ are the current neuron values
coming from the last layer and $\mathbf{h}^{t}$ and $\mathbf{h}^{t-1}$ are
current and previous states, respectively. This arrangement allows the network
to learn how stress evolves for a sequence of strains instead of bulding a
regression model from independent stress-strain pairs and is illustrated in
Fig. 2a. However, in practice, the efficiency of RNNs are impeded by vanishing
gradient problems and are not suitable for long-term history dependent
problems.

To overcome that, more sophisticated architectures (more popularly known as
cells) have been proposed. Among the most popular ones are the Gated Recurrent
Unit (GRU) and the Long-Short Term Memory (LSTM), illustrated in Fig. 2b and
2c, respectively. The internal mechanisms, also known as gates, used to
control the flow of information passing from one state to another are
represented by the colored circles. For each gate, additional parameters need
to be learned by the network in a way that the element-wise application of the
sigmoid (red circles) or tanh (purple circles) functions can wisely retain
what should be preserved and what can be forgotten in a long sequence.

Despite best efforts, these architectures are still vulnerable to overfitting,
compromising their ability to generalize well to new data. Potential solutions
to prevent this phenomenon include regularization techniques such as L2
penalty, early stopping and dropout. In this work, a special type of dropout
proposed by Kingma et al. [32] is used in combination with a GRU architecture
is considered to perform the comparison with the proposed network. In this
Bayesian GRU, the regular dropout with continuous noise (i.e. Gaussian
dropout) is reintepreted as a variational method that allows optimal dropout
rates to be inferred from the data as opposed to it being fixed and defined in
advance as usual. This circumvents the need for a validation set during model
selection. For more details, the interested reader is referred to [32].

(a) Classical Recurrent Unit

(b) Gated Recurrent Unit (GRU)

(c) Long-Short Term Memory Unit (LSTM)

Figure 2: Different architectures for recurrent units/cells: red circles
correspond to the element-wise application of the sigmoid function, while
purple circles correspond to the tanh function.

##  4 Physically Recurrent Neural Network

This section presents the neural network proposed to capture path-dependent
behavior of heterogeneous microscopic models. The core task of the network is
to learn how the macroscopic strain $\bm{\varepsilon}^{\Omega}$ can be
dehomogenized into a small set of representative material points and how their
responses can be combined to obtain the homogenized macroscopic stresses
$\bm{\sigma}^{\Omega}$. For this, the parametric regression model
$\mathcal{R}$ illustrated in Fig. 3 is proposed: a combination of a data-
driven encoder, a material layer with embedded physics-based material models
and a data-driven decoder. Each of these components are discussed in detail in
Sections 4.1, 4.2 and 4.3, respectively. Finally, the training process is
described in Section 4.4 and the use of this network as constitutive model in
FE2 frameworks is discussed in Section 4.5.

Figure 3: The proposed architecture

###  4.1 Encoder

The encoder consists of all parameters that convert the macroscopic strain
from the input layer to the values used as input of the material layer, which
corresponds to the grey lines in Fig. 3. Since these values are the inputs of
actual material models that are later combined by the decoder into the
prediction of the macroscopic stresses, we interpret the role of the encoder
as being the microscopic periodic boundary-value problem (BVP) solved with FE,
only at a much lower computational cost. On the other hand, no information on
the displacement field of the micromodel is retrieved by the network as the
encoder is learned based exclusively on snapshots of macroscopic stresses.
This understanding is depicted in Fig. 4 by the grey curved line linking the
strains from the macroscopic scale and the fictitious strains seen by the
material points in the network.

Figure 4: Interpretation of the proposed network with respect to a full-order
solution

As for the architecture, an arbitrary number of layers and units (with
conventional activation functions as illustrated in Fig. 3) can be used. In
case regular dense layers are employed, the neuron states ($\mathbf{a}_{i-1}$)
from the previous layer $i$-1 are propagated to the following layer $i$
according to:

$\displaystyle\mathbf{v}_{i}=\mathbf{W}_{i}\mathbf{a}_{i-1}+\mathbf{b}_{i}\qquad\mathbf{a}_{i}=\upphi(\mathbf{v}_{i})$
(16)

where $\mathbf{W}_{i}\in\mathbb{R}^{n_{i}\times n_{i-1}}$ is a weight matrix
and $\mathbf{b}$ is a bias term with $n_{i}$ being the number of neurons of
layer $i$, and $\upphi$ is an activation function applied in an element-wise
manner to the neuron values of $i$ to introduce nonlinearity into the network.
In the particular case where the dense layer is either the input or the output
layers, no activation function is applied and $\mathbf{a}_{i}=\mathbf{v}_{i}$.
Popular activation functions include the sigmoid, tanh and ReLU. In the
present investigation, a single dense layer is considered before the material
layer, which results in a linear relationship between macroscopic strains and
local strains in the material layer.

###  4.2 Material layer

The material layer is responsible for introducing explicitly the same physics-
based material models used in the RVE that the network will be a surrogate
for. To properly incorporate them and take full advantage of its outputs,
important changes on how neurons are evaluated compared to regular dense
layers are proposed. First, instead of introducing nonlinearity in a element-
wise manner with a scalar-to-scalar activation function, neurons are grouped
in $m$ sets of the size of the input layer (see colored boxes in Fig. 3) and
then evaluated as a subgroup. Each subgroup is referred to as a fictitious
material point and its size is equal to the length of the strain vector (i.e.
length 3 for the present investigation in two dimensions). In this
arrangement, each neuron of the subgroup $j$ represents one component of the
strain vector $\bm{\varepsilon}_{j}$, as illustrated in Fig. 5a.

(a) Fictitious material point

(b) Material cell

Figure 5: Schemes of (a) fictitious material point and (b) material layer as a
cell

Supposing the micromodel contains $n$ material models
$\mathcal{D}^{\omega}_{1}$, $\dots\mathcal{D}^{\omega}_{n}$, several
combinations of them can be employed in the material layer. The choice on
which material models should be used to evaluate the fictitious material
points depends on the types of non-linearity embedded in each of these models.
For instance, in this work, the micromodel consists of a composite
microstructure with two material models: a linear elastic model
$\mathcal{D}^{\omega}_{1}$ to describe the fibers and an elastoplastic model
with isotropic hardening $\mathcal{D}^{\omega}_{2}$ to describe the matrix.
However, in the material layer, only model $\mathcal{D}^{\omega}_{2}$ is
employed in all fictitious material points. This is possible because the
network still can make any of the subgroups to behave linear elastically by
passing small strains and scaling the stresses to give a significant elastic
contribution.

However, if instead of a linear elastic model, a nonlinear hyperelastic model
was used to describe the fibers the network would not perform optimally. In
that case, although the elastoplastic model does introduce nonlinearity to the
network, the (nonlinear) contribution from the fibers is no longer embedded in
that model. Furthermore, if the nonlinear hyperelastic model was the one
chosen to evaluate all subgroups, the network would essentially become a feed-
forward one with no history information taken into account (explicitly or
implicitly), losing the ability to predict elastic unloading. For such a
micromodel, both material models would need to be considered.

Another interesting hypothetical case is that of a micromodel with two
elastoplastic phases with different material properties. This time, depending
on the contrast of the material properties, a single material model with a
fixed set of properties coming from one of the two phases might be enough to
reproduce the homogenized response of the micromodel. In that case, making the
material properties of each fictitious material point a trainable feature
would also be advantageous. This could also be a valuable feature when dealing
with experimental data or with a micromodel with continuously varying material
properties. Naturally, far more complex arrangements than those discussed here
are found in practice. Addressing these extensions is object of ongoing
research. A general guideline is to employ all different sources of
nonlinearity with their respective known material properties in the material
layer.

For simplicity, we choose to illustrate the network with a general material
model $\mathcal{D}^{\omega}$ for all fictitious material points. Such model
can take the form of any of the material models with known material properties
used in the micromodel. Here, $\mathcal{D}^{\omega}$ takes as input the
current strain $\bm{\varepsilon}^{t}\in\mathbb{R}^{n_{\varepsilon}}$ and the
internal variables from previous time step
$\bm{\alpha}^{t-1}\in\mathbb{R}^{n_{IntVar}}$, where $n_{IntVar}$ is the
number of internal variables of the material model. With that, the model is
used to evaluate current stress state
$\bm{\sigma}^{t}\in\mathbb{R}^{n_{\epsilon}}$ and updated internal variables
$\bm{\alpha}^{t}$. These quantities motivated the tailor-made architecture of
the proposed layer.

To store the internal variables used as input/output of the material model, an
auxiliary vector $\mathbf{h}_{j}\in\mathbb{R}^{n_{IntVar}}$ referred as
history vector is defined. In the particular case of a subgroup with a
material model with no internal variables (e.g. linear elastic model),
$\mathbf{h}_{j}$ does not exist since $n_{IntVar}=0$. For the first time step,
the history vector is initialized as zero for all $m$ subgroups. As
information reaches the material layer and the material model is called, three
outputs are made available: the stresses $\bm{\sigma}^{t}_{j}$, the updated
internal variables $\bm{\alpha}^{t}_{j}$ and the tangent stiffness matrix
$\mathbf{D}^{t}_{j}\in\mathbb{R}^{n_{\varepsilon}\times\varepsilon}$. In this
layer, only the stresses are propagated forward. To do this, each stress
component is associated to a unit of the subgroup, as illustrated in Fig. 5a.
Then, the updated internal variables $\bm{\alpha}^{t}$ are stored in
$\mathbf{h}^{t}_{j}$ so that when new strains $\bm{\varepsilon}^{t+1}_{j}$ are
fed to the fictitious material point, the material model is aware of its own
history so far, making the $\bm{\varepsilon}\text{-}\bm{\sigma}$ path of each
subgroup unique. This architecture is illustrated in Fig. 5b.

Note that $\mathbf{h}_{j}$ has no weighted connections with the neurons in the
material layer (or any other layer). It is merely an object to store
information about the history of the fictitious material point. It is also
worth mentioning that since no data from the microscale has been collected and
imposed in the network, the paths seen by the fictitious material points do
not need to hold any similarity with actual integration points of the
microscopic model.

Using standard machine learning notation, the material layer propagates
previous neuron states ($\mathbf{a}_{k-1}$) and applies the material model
$\mathcal{D}^{\omega}$ as follows:

$\mathbf{v}_{k}=\mathbf{W}_{k}\mathbf{a}_{k-1}+\mathbf{b}_{k}\Rightarrow\mathbf{a}_{k},\mathbf{h}=\mathcal{D}^{\omega}(\mathbf{v}_{k},\mathbf{h}^{t-1})$
(17)

where $\mathbf{W}_{k}\in\mathbb{R}^{n_{k}\times n_{k-1}}$ is the weight matrix
connecting layers $k$-1 and $k$, $\mathbf{b}_{k}\in\mathbb{R}^{n_{k}}$ is a
bias term. In addition, $\mathbf{v}_{k}$ are the neuron values (correspond to
the concatenated vector of all microscopic strains $\bm{\varepsilon}_{j}$),
$\mathbf{h}^{t-1}$ and $\mathbf{h}^{t}$ are history-related term (correspond
to concatenated vector of all internal variables $\bm{\alpha}_{j}^{t-1}$)
resulting from the material models with path-dependent behavior from past and
current time step and $\mathbf{a}_{k}$ are the current neuron states
(correspond to the concatenated vector of all microscopic stresses
$\bm{\sigma}_{j}$).

###  4.3 Decoder

The decoder consists of all parameters that convert the outputs from the
material layer to the predicted macroscopic stress
$\widehat{\mathbf{\sigma}}^{\Omega}$ in the output layer, which corresponds to
the brown lines in Fig. 3. Similar to the encoder, an arbitrary number of
conventional layers and units can be employed. In the full-order solution,
after convergence of the microscopic BVP, the macroscopic stresses are
obtained by the volume average of microscopic stresses over the entire RVE. In
the network, since the solution of the microscopic BVP is replaced by the
encoder and the microscopic material points in the RVE are replaced by the few
fictitious material points, the decoder is then analogous to the
homogenization operator that transforms local stresses to macroscopic
stresses, as illustrated by the brown curved line in Fig. 4.

In the present work, we use a single dense layer (output) with linear
activation and physics-motivated modifications to perform the task. With this,
all the nonlinearity of the network arises from the models in the material
layer. As discussed previously, the decoder can be understood as the averaging
operator in a multiscale approach and with the chosen architecture (dense-
material-dense), the weights of the output layer can be seen as the relative
contribution of each fictitious material point to the macroscopic stress.
Based on that, a constraint on the positivity of the weights of the output
layer is considered. For that, a softplus function $\rho(\cdot)$ is applied
element-wise on the weights matrix before computing the neuron values of the
last layer:

$\mathbf{v}_{l}=\rho(\mathbf{W}_{l})\ \mathbf{a}_{l-1}+\mathbf{b}_{l}$ (18)

where $\mathbf{b}_{l}$ is set to zero and $\mathbf{a}_{l-1}$ corresponds to
the stresses coming from the material layer. This procedure guarantees that,
after the tranformation, weights will always be positive.

###  4.4 Training

The goal of the training phase is to minimize a loss function given by:

$L=\frac{1}{N}\ \sum_{i=1}^{N}\frac{1}{2}\ \|\
\bm{\sigma}^{\Omega}(\bm{\varepsilon}^{\Omega}_{i})-\widehat{\bm{\sigma}}^{\Omega}(\bm{\varepsilon}^{\Omega}_{i})\
\|^{2}$ (19)

where $N$ is the number of snapshots. Based on it, a Stochastic Gradient
Descent (SGD) optimization algorithm is used to update the trainable
parameters $\mathbf{W}$ and $\mathbf{b}$:

$\displaystyle\mathbf{W}^{n}=\mathbf{W}^{o}-\mathcal{A}\Big{(}\frac{1}{B}\sum_{i=1}^{B}\frac{\partial{L_{i}}}{\partial{\mathbf{W}}}\Big{)}$
(20)
$\displaystyle\qquad\mathbf{b}^{n}=\mathbf{b}^{o}-\mathcal{A}\Big{(}\frac{1}{B}\sum_{i=1}^{B}\frac{\partial{L_{i}}}{\partial{\mathbf{b}}}\Big{)}$

where $L_{i}$ is the loss of the $i$-th sample, $o$ indicates current values,
$n$ indicates updated values and $B$ is the size of the sample mini-batch used
in the update. Finally the $\mathcal{A}$ operator depends on the solver. In
this work, the Adam optimizer [33] is used.

To compute the gradients appearing in Eq. (20), backpropagation in time is
employed in a similar fashion as done to RNNs: based on the network state
($\mathbf{v}$ and $\mathbf{a}$) after computing each training curve with $n$
pairs of $\bm{\sigma}-\bm{\varepsilon}$, the chain rule is used to propagate
the derivative of the loss functions starting from the output layer and
progressively moving back through the network and through time. Commonly, this
process is dealt with by automatic differentiation, but we present the
expressions to allow for integrating the network directly into existing FE
software. For this, two auxiliary quantities are defined. The first is defined
for each layer and helps propagating the error through the network
$\mathbf{d}_{i}\in\mathbb{R}^{n_{i}}$. Starting from the output layer $l$, it
is defined as:

$\mathbf{d}_{l}=\frac{\partial{L}}{\partial{\mathbf{a}_{l}}}=\widehat{\bm{\sigma}}^{\Omega}-\bm{\sigma}^{\Omega}$
(21)

Next, the effect of the activation function is taken into account as:

$\bar{\mathbf{d}_{i}}=\mathbf{d}_{i}\odot\frac{\partial{\upphi(\mathbf{v}_{i})}}{\partial{\mathbf{v}_{i}}}$
(22)

after which it is possible to compute the gradients of the trainable
parameters:

$\displaystyle\frac{\partial{L}}{\partial{\mathbf{W}_{i}}}=\bar{\mathbf{d}}_{i}\mathbf{a}_{i}^{T}\qquad\frac{\partial{L}}{\partial{\mathbf{b}_{i}}}=\bar{\mathbf{d}}_{i}$
(23)

Finally, the values $\mathbf{d}$ of the previous layer (the next layer to be
backpropagated) can be computed as:

$\mathbf{d}_{i-1}=\mathbf{W}_{i}^{T}\bar{\mathbf{d}}_{i}$ (24)

and the algorithm moves to Eq. (22) for layer $i$-1.

When reaching the material layer, recall that the internal variables are
stored in $\mathbf{h}$ and used for keeping track of the evolution of the
internal variable through time. For that reason, a second auxiliary quantity
is introduced and Eq. (22) is replaced by:

$\bar{\mathbf{d}_{i}}=\mathbf{d}_{i}\odot\frac{\partial{\mathbf{a}_{i}}}{\partial{\mathbf{v}_{i}}}+\mathbf{d}_{h}^{t+1}\odot\frac{\partial{\mathbf{h}}}{\partial{\mathbf{v}_{i}}}$
(25)

where the first term concerns the derivatives of stresses with respect to
strains, the second term concerns the derivatives of the current internal
variables with respect to strains and $\mathbf{d}_{h}\in\mathbb{R}^{n_{i}}$ is
given by:

$\mathbf{d}_{h}=\mathbf{d}_{i}\odot\frac{\partial{\mathbf{a}_{i}}}{\partial{\mathbf{h}}}+\mathbf{d}_{h}^{t+1}\odot\frac{\partial{\mathbf{h}}}{\partial{\mathbf{h}^{t-1}}}$
(26)

Note that the derivatives of the stresses with respect to the strains of
material point $j$ are an output of the material model: the tangent stiffness
matrix $\mathbf{D}_{j}$. The remaining derivatives in Eqs. (25) and (26) are
evaluated using central finite differences. Naturally, computing gradients
with other methods would also be possible. For instance, if the material model
used in the network supports automatic differentation, storing the internal
variables in $\mathbf{h}$ for backpropagation can be bypassed as the
derivatives are automatically obtained in this approach.

Finally, to obtain the gradients of the trainable parameters including the
history-dependence coming from the material layer and compute the values
$\mathbf{d}$ of the previous layer, we consider Eq. (25) instead of (22) in
the expressions shown in Eqs. (23) and (24), respectively.

###  4.5 Use as constitutive model

To make new stress predictions, the macroscopic strain
$\bm{\varepsilon}^{\Omega}$ is fed to the input layer and a complete forward
pass is performed. The final activated neuron values of the output layer give
the predicted stress. To obtain the macroscopic consistent tangent stiffness
matrix $\mathbf{D}^{\Omega}$, a complete backward pass is required:

$\mathbf{D}^{\Omega}=\frac{\partial{\widehat{\bm{\sigma}}^{\Omega}}}{\partial{\bm{\varepsilon}^{\Omega}}}=\frac{\partial{\mathbf{a}_{l}}}{\partial{\mathbf{v}_{0}}}=\mathbf{J}$
(27)

which is obtained with a backward pass through the network:

$\mathbf{J}_{i}=\mathbf{J}_{i+1}\mathbf{I}_{i}^{\upphi}\mathbf{W}_{i}\ \ \
\textrm{with}\ \ \ \mathbf{J}_{l+1}=\mathbf{I}$ (28)

where $\mathbf{I}_{i}^{\upphi}$ is a matrix whose diagonal contains the
derivatives of the activation function with respect to the neuron values
$\mathbf{v}$:

$\mathbf{I}_{i}^{\upphi}=\textrm{diag}(\frac{\partial{\upphi(\mathbf{v}_{i})}}{\partial{v}})$
(29)

except for the material layer. In that case, such matrix is full and consists
of the concatenation of the tangent stiffness matrix of all fictitious
material points. It is worth mentioning that despite the linear dependency on
the tangent stiffness matrices of the material models, the Jacobian matrix of
the network does not inherit their spectral properties.

###  4.6 Analogies to other methods

In this section, the parallels between features of the proposed network and
related works in the literature are briefly discussed. One possible analogy
comes from hyper-reduced-order models [34]. With the architecture chosen for
the present investigation, both methods work on a reduced number of material
points with modified (integration) weights. However, in the network, these
points are only fictitious and learned by the encoder based on snapshots of
the homogenized stresses. Moreover, each stress component is associated with a
different weight. By contrast, the material points in the hyper-reduction
approach exist in the microscopic model and a single modified integration
weight of each material point selected is used to compute all its
stress/internal force components.

Following the discussion on the encoder, it is worth highlighting how this
feature would be framed with respect to asymptotic homogenization schemes such
as Mori-Tanaka [35]. In this type of solution, the microscopic problem is also
not solved explictly and only average fields are calculated. Relying on the
equivalent inclusion idea and on Eshelby’s solution [36], the strain
concentration tensor is obtained analytically and yields the full solution of
the microcopisc model as it correlates the average field of the phases in the
micromodel with its average field. In our network, although the macroscopic
stresses are also obtained by relating macroscopic and (fictitious)
microscopic strains through an encoder, here no average field is calculated
for each of the phases. Indeed, not every phase needs to be included in the
material layer and multiple strain paths for the same phase are considered.
Furthermore, while Mori-Tanaka is accurate for moderate volume fractions of
the inclusions, such restriction is not present in our method.

Compared to PINNs, in which physical constraints are explicitly included in
the loss function, here, most physical constraints are naturally taken care of
by the physics-based material models directly embedded in the material layer.
The proposed approach is also more general as it can be directly used for
arbitrary material models and is not particularly tailored to a single type of
model (e.g. elastoplastic behavior [14, 15]). As an added benefit, our model
selection procedure makes physical sense: we add more material points or
material models to the network.

Another noteworthy strategy with relevant analogies to our method is the DMN
[20]. In this approach, the contribution of a few material points evaluated
using the classical constitutive models in the RVE is also employed to make
predictions in the online phase. On the same reasoning as discussed in Section
4.2, since the inputs come directly from actual material models, path-
dependency is captured naturally. However, the main concept and architecture
of DMNs are different from the ones explored here. In the offline phase, the
goal of the DMN is to find a topological representation of the RVE with fewer
degrees of freedom (i.e. material points) based only on the elastic stiffness
matrices of the different material phases that compose the original
micromodel. For the online phase, the feature space is increased to include
residual stresses of the micromodel components, and an iterative procedure is
implemented. The authors compare the incremental strains of the material
points at the beginning of the iteration with the one obtained by a de-
homogenization process that backpropagates the macroscopic incremental strain
from the output layer to the bottom layer (i.e. input layer). Upon
convergence, the set of internal variables of each material point at the
bottom layer is therefore updated.

In the present work, the feature space is the same in both phases and no
iterative procedure is employed in the forward pass, which simplifies
implementation and reduces even further the number of material model calls.
Here, the strain path each fictitious material point follows is simply
described by the encoder and not all phases need to be included in the
network. The homogenized stresses and tangent stiffness are obtained in a
single forward and backward pass, respectively. Furthermore, the
backpropagation in our approach is considerably simpler than the DMN. Although
the use of homogenization (and de-homogenization) operations in the DMN
assigns physical interpretation to the model, it also makes training a rather
intricate process.

Finally, to draw a parallel with LSTMs, one might understand $\mathbf{h}$ as
the cell state $\mathbf{c}$, but instead of using bijective and smooth
functions such as the sigmoid and tanh functions to describe the evolution of
the material response, the material model itself is directly employed. This
bypasses the need to learn new parameters to regulate the flow of information
kept or forgotten throughout time (see Fig. 2b) and has important implications
for the training process. The most important one is the ability to mirror
physical behaviors such as elastic unloading/reloading without ever seeing the
pattern during training, a stark contrast with LSTMs and GRUs that usually
require extensive training sets with multiple cycles of loading and reloading
at different strain levels with different step sizes. The physical
interpretation of the nonlinearity is directly embedded in the network. In the
numerical examples of this work, the nonlinearity is due to plasticity, but
other effects such as hyperelasticity, visco-plasticity, stiffness
degradation, or any combination thereof, could be embedded by adapting the
constitutive model that is used in the material layer.

##  5 Design of Experiments

One critical aspect of the training and testing of surrogate models is the
formulation of a sampling plan. Typically, a uniform distribution of the
sampling points is desirable, but that task becomes more complex when path-
dependent behavior is present. In this case, pairs of strains and stresses are
collected and processed as sequences, which leads to potentially infinite-
dimensional parameter spaces.

In this work, two strategies are considered. In the first approach,
proportional loading paths are generated, which means that the stress ratio
between the components is constant. Here, the sequence of strains is created
based on two features: the loading function $\lambda(\Delta\varepsilon,t)$ and
the loading direction given by the unit vector $\mathbf{n}$, where
$\Delta\varepsilon$ is the step size and $t$ is the current time step. For
each time step, $\mathbf{n}$ is multiplied by the scalar-valued loading
function $\lambda$ creating a new set of strains, which is in turn applied at
the controlling nodes of the microscopic model.

For monotonic loading, the loading function is as depicted in Fig. 7a. The
values in the unit vector can come from prior knowledge of the material as
illustrated in Fig. 6b, in which only fundamental cases such as uniaxial
strain, pure shear, and biaxial cases are considered, or from random
distributions as represented by the purple line Fig. 6a. In the present work,
the random directions are obtained by sampling values from $n_{\varepsilon}$
independent Gaussian distributions (X $\sim\mathcal{N}(0,1)$) and subsequently
normalizing the vectors.

(a) Types of loading used in this study

(b) Proportional loading with a priori known directions

Figure 6: Different design of experiments strategies

Despite the simplicity in creating such paths, RNNs trained exclusively on
monotonic cannot predict cylic responses. Thus, to create non-monotonic
sequences, a linear piecewise function as the one depicted in Fig 7b is used.
Note that even though the loading function is changed, the unit vector is kept
constant for the entire strain sequence, yielding proportional loading.
However, to cover the entire space of possible cyclic responses, a large (and
a priori unknown) number of curves comprehending different unloading points
with different duration of unloading/reloading and step sizes is necessary. In
this work, that matter is first handled in a simplified way by only sampling
two different cycles of unloading/reloading.

Finally, in a more general approach, a second strategy to create the
$\bm{\varepsilon}\text{-}\bm{\sigma}$ paths is considered: the random walks.
These are typically defined by sampling random strain increments with random
loading directions for each time step, resulting in non-proportional loading.
In this work, random walks are created by associating the prescribed strains
to independent Gaussian Processes (GPs) with X
$\sim\mathcal{N}(\mu,\,\sigma^{2})$ and covariance function given by:

$k(\mathbf{x}_{p},\mathbf{x}_{q})=\sigma_{f}^{2}\exp\Big{(}-\frac{1}{2\ell^{2}}\|\mathbf{x}_{p}-\mathbf{x}_{q}\|^{2}\Big{)}$
(30)

where $\sigma_{f}^{2}$ is the variance and $\ell$ is a length scale. In this
setting, the lengthscale controls the smoothness of the strain path and the
variance controls how large the step size can be for each prescribed degree of
freedom of the controlling nodes. Similar approaches were employed by Mozaffar
et al. [8] and Logarzo et al. [11].

Input : lengthscale $\ell$, variance $\sigma_{f}^{2}$, number of strain
components $N_{\textrm{components}}$, number of time steps
$N_{\textrm{steps}}$

Output : macroscopic strains $\mathfrak{D}_{\varepsilon}$ and macroscopic
stresses $\mathfrak{D}_{\sigma}$

1 initialize datasets: $\mathfrak{D}_{\varepsilon}\leftarrow\varnothing$,
$\mathfrak{D}_{\sigma}\leftarrow\varnothing$

2 for  _$i\in[1,2,...,N_{\text{components}}]$_ do

3 initialize input and output datasets for $\textrm{GP}_{i}$:
$\mathbf{X}_{\textrm{GP}_{i}}\leftarrow 0$,
$\mathbf{Y}_{\textrm{GP}_{i}}\leftarrow 0$

4 initialize $\textrm{GP}_{i}$: $\textrm{GP}_{i}\leftarrow$ initGP(
$\mathbf{X}_{\textrm{GP}_{i}},\mathbf{Y}_{\textrm{GP}_{i}},\ell,\sigma_{f}^{2}$)

5

6for  _$t\in[1,2,...,N_{\text{steps}}]$_ do

7 initialize current macroscopic strain:
$\bm{\varepsilon}_{\textrm{current}}\leftarrow\varnothing$

8 for  _$i\in[1,2,...,N_{\text{components}}]$_ do

9 sample from posterior distribution: $\varepsilon_{i}\leftarrow$
GPi::samplePosterior ( $t$ )

10 add value to strain vector:
$\bm{\varepsilon}_{\textrm{current}}\leftarrow\bm{\varepsilon}_{\textrm{current}}\cup\varepsilon_{i}$

11 solve micromechanical BVP: $\bm{\sigma}_{\textrm{current}}\leftarrow$
fullModel::materialUpdate ( $\bm{\varepsilon}_{\textrm{current}}$ )

12 if  _convergence_ then

13 store equilibrium solution of micromodel: fullModel::storeSolution( )

14 add macroscopic strains and stresses to dataset:
$\mathfrak{D}_{\varepsilon}\leftarrow\mathfrak{D}_{\varepsilon}\cup\bm{\varepsilon_{\textrm{current}}}$,
$\mathfrak{D}_{\sigma}\leftarrow\mathfrak{D}_{\sigma}\cup\bm{\sigma}_{\textrm{current}}$

15 for  _$i\in[1,2,...,N_{\text{components}}]$_ do

16 add time step and current strain to dataset of GPi:
$\mathbf{X}_{\textrm{GP}_{i}}\leftarrow\mathbf{X}_{\textrm{GP}_{i}}\cup t$,
$\mathbf{Y}_{\mathrm{GP}_{i}}\leftarrow\mathbf{Y}_{\textrm{GP}_{i}}\cup\varepsilon_{i}$

17 update GPi with new data: GPi::update ( $\mathbf{X}_{\textrm{GP}_{i}}$,
$\mathbf{Y}_{\textrm{GP}_{i}}$ )

return ($\mathfrak{D}_{\varepsilon}$, $\mathfrak{D}_{\sigma}$)

Algorithm 1 Generation of random loading path using GPs

The details of the present implementation are given in Algorithm 1. Note that
instead of drawing the entire strain sequence for a given component, we sample
it step by step and update the GP dataset before sampling again. This strategy
results in the same strain sequence given a fixed random seed throughout the
steps, but in this way the GPs can also be used in applications where the
number of loading steps is changed on-the-fly. Following the work of Logarzo
et al. [11], we define the mean of all GPs to be zero and include $t=0$ and
$\varepsilon_{i}=0$ as a prior. In addition to that, references to fullModel
(i.e. the full-order microscopic model) in Algorithm 1 are kept as minimal and
general as possible. One example of loading path resulting from Algorithm 1 is
illustrated in Fig. 6a.

In this paper, both strategies generate $\bm{\varepsilon}\text{-}\bm{\sigma}$
curves containing 60 time steps, unless stated otherwise. To summarize the
types of loading studied in the following sections:

  * 1.

Type I: monotononic and proportional loading paths with a priori known
directions. The 18 directions used to train the proposed network are
illustrated in Fig. 6b and include uniaxial strains, pure shear, biaxial cases
and biaxial with shear cases.

  * 2.

Type II: monotonic and proportional loading paths randomly spread across the
design space. The loading directions are generated randomly and the loading
function is as shown in Fig. 7a.

  * 3.

Type III: non-monotonic and proportional loading paths randomly spread across
the design space. Again, the loading directions are random, but the loading
function is now given by Fig. 7b and includes one cycle of unloading;

  * 4.

Type IV: Variations to Type III:

    * (a)

Type IVa: same loading directions as the test set of Type III, but
unloading/reloading takes place at a different point in time as shown in Fig.
7c;

    * (b)

Type IVb: same loading directions as the test set of Type III, but time step
is 10 $\times$ smaller. Thus, to reach the same norm as the original curve in
Type II, 600 time steps are evaluated, as depicted in Fig. 7d;

  * 5.

Type V: non-monotonic and non-proportional loading paths randomly spread
across the design space. A GP-based path described by Eq. (30) is illustrated
in Fig. 6a. Fig. 7e ilustrates the strain paths of each component using this
approach with lengthscale $\ell=20$ and $\sigma_{f}=1.0\times 10^{-3}$.

(a) Monotonic loading (Types I and II)

(b) Non-monotonic loading (Type III)

(c) Different unloading/reloading (Type IVa)

(d) Non-monotonic with smaller step size (Type IVb)

(e) Random non-proportional and non-monotonic with $\ell=20$ (Type V)

Figure 7: Proportional and non-proportional loading functions

##  6 Assessing network performance

In this section, the performance of the proposed network is compared to a
state-of-the-art RNN trained on different training dataset sizes and methods
to sample the design space. The comparison is done for a single micromodel.
Specifically, four scenarios are investigated: (i) predicting
unloading/reloading behavior from monotonic data, (ii) predicting
unloading/reloading behavior from non-monotonic data, (iii) predicting unseen
patterns from non-monotonic data, and (iv) training with non-monotonic and
non-proportional loading paths. In the first three scenarios, our network is
trained exclusively on the fundamental loading cases of Type I (18 curves),
while the training of the RNN is an open question to be addressed in the
following sections.

From now on, the network presented in this work will be referred to as
Physically Recurrent Neural Network, or simply PRNN. The PRNN was trained for
80000 epochs, while the RNN was trained for 60000 epochs with an early
stopping criterion, which consists of interrupting training if the best
training loss so far is not improved over a given period (in this work, 5000
epochs). The Adam optimizer is used in all cases with batch size of 9 and
default parameters suggested by Kingma and Ba [33], with the exception of the
learning rate of 0.01 for the RNNs. The layer sizes are chosen through model
selection to provide optimal results and fair comparison with our approach to
the best of our knowledge. The methodology adopted is briefly described in
Section 6.1.

The microscopic model consists of an RVE with 36 elastic fibers (volume
fraction = 0.6) with properties E = 74000 MPa and $\nu$ = 0.2 embedded in an
elastoplastic matrix with isotropic hardening. The geometry and the mesh with
7048 elements are depicted in Fig. 8. The elastoplastic matrix is modeled
using the von Mises yield criterion with properties E = 3130 MPa, $\nu$ = 0.3
and yield stress given by:

$\sigma_{y}=64.8-33.6\exp^{-\varepsilon^{p}_{eq}/0.0003407}$ (31)

where $\varepsilon^{p}_{eq}$ is the equivalent plastic strain defined as:

$\varepsilon^{p}_{eq}=\sqrt{\frac{2}{3}\
\bm{\varepsilon}^{p}:\bm{\varepsilon}^{p}}$ (32)

and $\bm{\varepsilon}^{p}$ is the plastic strain. Plane stress conditions are
assumed.

(a)  Figure 8: Geometry and mesh discretization of microscopic model adopted
in this work.

###  6.1 Model selection

In this section, a grid-search strategy is employed to choose the best
architecture for the networks. The goal is to find the optimum layer size
before heading to the testing sections. Here, four different layer sizes and
10 different weight initializations are considered to mitigate the effect of
randomness.

For the PRNN, a single material layer without biases and four different sizes
are considered. In this case, only the elastoplastic model is used so that the
size of the material layer is the only variable in the model selection. Recall
that the network still can make a subgroup to behave elastically by passing
small strains to the material model. The training and the validation sets
consist of 18 Type I curves and 54 Type II curves, respectively. Fig. 9a shows
the boxplot with the average validation error of each run alongside the mean
error value. The networks with 6 units (or two fictitious material points)
performed best.

For the Bayesian RNN, a single GRU cell is employed and in each
initialization, the training set consists of the fixed set of 18 Type I
curves, 90 Type II curves randomly chosen from a pool of 1800 curves, and 90
Type III curves also randomly chosen from a pool of 1800 curves, amounting to
198 loading paths. That way, all types of curves used for training in the
following sections are covered. Fig. 9b shows the boxplot with the average
training error of each run alongside the mean error value of all runs
represented by the x marker. In this study, it is found that the GRU with 128
units performs best. In this case, no validation set is needed to determine
the best dropout rate as the type of RNN used in this investigation infers it
from the training data by default (see Section 3).

(a) Validation error for PRNNs trained on $\mathfrak{D}_{\mathrm{PRNN}}=\\{$18
Type I curves$\\}$ and $\mathcal{V}_{\mathrm{PRNN}}=\\{$54 Type II curves$\\}$

(b) Training error for RNNs trained on $\mathfrak{D}_{\mathrm{RNN}}=\\{$18
Type I curves, 90 Type II and 90 Type III curves$\\}$

Figure 9: Model selection of PRNN and RNN

###  6.2 Predicting unloading/reloading from monotonic data

In this section, the training process of the RNNs on Type II curves (i.e.
without unloading) is reported. The first test set consists of 100 Type II
curves. Fig. 10 shows the average error of the RNNs over the test set compared
to the best (blue triangle), worst (upside-down blue triangle), and average
error (blue circle) found by the PRNNs. Note that the secondary axis starts
with 18 curves, this is because the known directions used for training the
PRNNs are also a fixed set in the training of the RNNs. The training of the
RNNs is stopped with 288 curves. At that stage, a similar level of accuracy
between the PRNNs and the RNNs is obtained (although with a training set 16
times larger) and the addition of new curves only yields a marginal gain in
accuracy.

Figure 10: Absolute error over monotonic test set
$\mathcal{T}_{\mathrm{II}}=\\{$100 Type II curves$\\}$ for RNNs trained on
Types I and II and PRNNs on Type I only

Next, a new test set with 100 Type III curves is evaluated by the same
networks trained on the 288 monotonic curves. This time, the RNN fails to
capture unloading and the addition of more monotonic data is ineffective, as
shown in Fig. 11a. This outcome is not new to the literature and it is not
surprising that RNNs need to see unloading behavior during training in order
to be able to describe it. However, in contrast to the RNNs, the PRNNs provide
the same level of accuracy for the test sets with and without unloading, even
when not exposed to unloading data during training. In Fig. 11b a single
representative case from test set $\mathcal{T}_{\mathrm{III}}$ is plotted
using the best RNN and PRNN. Both networks show good agreement with the
reference solution until unloading starts (a feature not covered during
training), but only the PRNN is capable of capturing the elastic
unloading/reloading.

(a) Absolute error over test set $\mathcal{T}_{\mathrm{III}}=\\{$100 Type III
curves$\\}$ for RNNs trained on Types I and II and PRNNs on Type I only

(b) Stress-time view of representative case from test set
$\mathcal{T}_{\mathrm{III}}$ using the best RNN and PRNN

Figure 11: Absolute error over non-monotonic test set
$\mathcal{T}_{\mathrm{III}}$ for RNNs and PRNNs trained only on monotonic data
and representative case

###  6.3 Predicting unloading/reloading behavior from non-monotonic data

Following the conclusions of Section 6.2, the training set of the RNN is
expanded to include curves with the same unloading behavior as the one
observed in the test set $\mathcal{T}_{\mathrm{III}}$. The 288 monotonic
curves of Types I and II from the previous section are combined with an
increasing number of non-monotonic curves of Type III. This time, with the
right features included in the training set, Fig. 12 shows a monotonic
decrease of the average error for the RNN on $\mathcal{T}_{\mathrm{III}}$
curves. However, the performance of the RNN only meets the one obtained by the
PRNN with around 32 times more data.

Figure 12: Absolute error over non-monotonic test set
$\mathcal{T}_{\mathrm{III}}=\\{$100 Type III curves$\\}$ for RNNs trained on
Types I, II and II and PRNNs on Type I only

###  6.4 Predicting unseen patterns from non-monotonic data

In this section, three additional test sets are considered for the RNN trained
on Types I, II and III and the PRNN trained on Type I only. The goal is to
test the ability of the networks to predict the macroscopic stress with
patterns different from those seen during training.

First, we consider a test set with 100 unseen curves of Type IVa, which
consists of proportional curves in random directions with a different
predefined unloading/reloading behavior than that of Type III. The average
error for that set is shown in Fig. 13a. Here, the 576 curves from the
previous section are no longer enough to provide good accuracy when predicting
a different unloading/reloading. By adding more Type III curves to the
training set of the RNN, the average error decreases from 6.4 MPa to around
4.0 MPa but no significant gain in the accuracy is observed when the total
number of curves is larger than 864 curves. Based on that, a representative
case from test set $\mathcal{T}_{\mathrm{IVa}}$ is shown in Fig. 13b. Despite
the relative low error from both networks, note that the RNN looses
performance once unloading starts while the PRNN continues to show good
agreement throughout the entire loading path.

(a) Absolute error over test set $\mathcal{T}_{\mathrm{IVa}}=\\{$100 Type IVa
curves$\\}$ for RNNs trained on Types I, II and III and PRNN on Type I only

(b) Stress-time view of representative case from test set
$\mathcal{T}_{\mathrm{IVa}}$ using the best RNN and PRNN

Figure 13: Absolute error over non-monotonic test set
$\mathcal{T}_{\mathrm{IVa}}=\\{$100 Type IVa curves$\\}$ for RNNs trained on
Types I, II and II and PRNNs on Type I only and representative case

For the next scenario, a test set with 100 Type IVb curves are considered.
These curves have the same unloading/reloading behavior as Type III, but with
a 10$\times$ smaller time step. Fig. 14a illustrates the average error of 10
networks over that test set and again, it is clear that the addition of new
curves with patterns different from the exact one being tested is not
beneficial to the RNN. Again, the PRNN provides good accuracy. Essentially,
the PRNN is only as sensitive to step size as the material models embedded in
it. Fig. 14b illustrates the networks’ predictions for a curve in the test set
$\mathcal{T}_{\mathrm{IVb}}$.

(a) Absolute error over test set $\mathcal{T}_{\mathrm{IVb}}=\\{$100 Type IVb
curves$\\}$ for RNNs trained on Types I, II and III and PRNNs on Type I only

(b) Representative $\bm{\varepsilon}\text{-}\bm{\sigma}$ curve from test set
$\mathcal{T}_{\mathrm{IVb}}$ using the best RNN and PRNN

Figure 14: Absolute error over different step size test set
$\mathcal{T}_{\mathrm{IVb}}=\\{$100 Type IVb curves$\\}$ for RNNs trained on
Types I, II and III and PRNNs on Type I only and representative case

As a final test, a set of 100 Type V curves, which corresponds to non-
proportional and non-monotonic paths, is considered. This type of curve
combines the two previous features: different step sizes and different
unloading/reloading locations. Fig. 15a shows the average error for additional
non-monotonic curves in the training of the RNNs. It is clear that the RNN
completely fails to capture non-proportional paths (lowest error around 32
MPa) and that the addition of more data with features different than those
being tested is a waste of resources. Although in different levels, a similar
trend of loss of accuracy is also observed in the PRNNs, where the best,
average and worse performances result in errors around 8.9 MPa, 17.0 MPa and
11.2 MPa respectively. Fig. 15b illustrates the different order of error
between the PRNN and the RNN on a representative case from test set
$\mathcal{T}_{\mathrm{V}}$.

(a) Absolute error over test set $\mathcal{T}_{\mathrm{V}}=\\{$100 Type V
curves$\\}$ for RNN trained on Types I, II and III and PRNN on Type I only

(b) Representative $\bm{\varepsilon}\text{-}\bm{\sigma}$ curve from test set
$\mathcal{T}_{\mathrm{V}}$ using the best RNN and PRNN

Figure 15: Absolute error over non-proportional and non-monotonic test set
$\mathcal{T}_{\mathrm{V}}=\\{$100 Type V curves$\\}$ for RNNs trained on Types
I, II and III and PRNNs on Type I only and representative case

###  6.5 Training both networks on non-monotonic and non-proportional loading

In this section, both networks are trained on the most generic set of curves,
i.e. random non-monotonic and non-proportional curves of Type V. In addition
to that, we trained the PRNNs on the known and proportional loading cases for
comparison purposes. In that case, the size of the material layer is kept at 6
units and three training dataset sizes are considered. First, only the pure
uniaxial cases are included, which yields 6 loading cases. Then, the 4 biaxial
cases are added to the previous training dataset, resulting in 10 loading
cases. And finally, we add the 8 cases with biaxial and shear loading, which
amounts to the 18 fundamental paths shown in Fig. 6b.

When training both networks on non-proportional paths a new model selection
procedure was carried out to determine the optimum size of the material layer
and the GRU cell, respectively. In this preliminary study, 10 different
weights initialization are considered again. For training the RNNs and the
PRNNs, 2304 and 198 Type V curves are used, respectively. Figs. 16a and 16b
show the boxplot with the average error of each run alongside the mean error
value. In this case, the networks with 18 units (which corresponds to six
fictitious material points) performed better. For the Bayesian RNN, the GRUs
with 128 units continue to provide the best performance. Therefore, this
architecture is the one used in the comparison presented below.

(a) Validation error for PRNNs trained on
$\mathfrak{D}_{\mathrm{PRNN}}=\\{$198 Type V curves$\\}$ and
$\mathcal{V}_{\mathrm{PRNN}}=\\{$54 Type V curves$\\}$

(b) Training error for RNNs trained on $\mathfrak{D}_{\mathrm{RNN}}=\\{$2304
Type V curves$\\}$

Figure 16: Model selection of PRNN and RNN for non-monotonic and non-
proportional loading

This time, all test sets discussed in Section 6.2, 6.3 and 6.4 are used again
to assess the accuracy of the networks with the new sampling strategy. Figures
17a-17e show the best, worst and average error for the cases studied so far in
order. Based on this study, a few important insights are worth mentioning:
after a certain point (around 576 curves), the RNNs reach an optimum level of
accuracy and the addition of new curves no longer boosts predictions for
proportional loading cases ($\mathcal{T}_{\mathrm{II}}$,
$\mathcal{T}_{\mathrm{III}}$, $\mathcal{T}_{\mathrm{IVa}}$ and
$\mathcal{T}_{\mathrm{IVb}}$). This is in line with the behavior observed in
the previous sections, in which the RNNs would only perform well when trained
with the same features as in the test set. And more importantly, changing the
sampling strategy also showed to have limited effect on improving their
performance. Granted, increasing even further the number of curves used for
training as well as the complexity of the RNN might help in that task.
However, the point stands that with the PRNN, this is not necessary. Note that
for the same training set sizes, the PRNN with either the known or the random
curves perfoms better than using the RNNs.

Finally, when choosing between known and random curves for training the PRNN,
the latter shows comparable errors with the first when predicting proportional
loading, but is significantly more accurate (see detail in Fig. 17e) for non-
proportional loading. For that reason, the PRNN trained on Type V is chosen to
illustrate the network’s capacity in the following FE2 examples. On average,
the accuracy of the PRNN reaches a plateau around 36 curves. From that point
on, the benefit of adding new data is limited.

(a) Test set $\mathcal{T}_{\mathrm{II}}$ = {100 Type II curves}

(b) Test set $\mathcal{T}_{\mathrm{III}}$ = {100 Type III curves}

(c) Test set $\mathcal{T}_{\mathrm{IV}_{\mathrm{a}}}$ = {100 Type IVa curves}

(d) Test set $\mathcal{T}_{\mathrm{IV}_{\mathrm{b}}}$ = {100 Type IVb curves}

(e) Test set $\mathcal{T}_{\mathrm{V}}$ = {100 Type V curves}

Figure 17: Absolute error of networks trained on different sampling strategies
and different test sets

##  7 FE2 applications

In this section, the PRNN trained with 36 non-proportional and non-monotonic
Type V curves and lowest test error for test set $\mathcal{T}_{\mathrm{V}}$ in
Section 6.5 is employed as the constitutive model in two numerical examples.
Results obtained with the PRNN as constitutive model are compared against
results obtained with full FE2 with the same micromodel that the network was
trained to be a surrogate for. Both types of analysis are performed with an
in-house Finite Element code using the open-source Jem/Jive C++ numerical
analysis library [37].

At the macroscale, an arc-length method with an adaptive-stepping scheme [38]
is adopted to tackle potential convergence issues. This way, if a loading step
does not converge with the given (full) step size, a reduction factor is
applied to it until the loading step converges or until a maximum number of
reductions in the initial step has been reached, terminating the analysis.
Upon convergence, in the following step, the analysis is resumed with the full
step size. In this work, each load step can be reduced by a factor of 0.4 for
a maximum number of 5 times.

All simulations, including the PRNN training, were executed on a single core
of a Xeon E5-2630V4 processor on a cluster node with 128 GB RAM running CentOS
7.

###  7.1 Tapered bar

The first example consists in a tapered composite specimen with length of 128
mm and height of 8 mm loaded in transverse tension. In this setting, the
36-fiber RVE model used to train the networks in the previous sections is
embedded at each integration point of the macroscale. The geometry and the
boundary conditions are shown in Fig. 18a. The FE2 problem is solved for 110
load steps with unloading according to the function shown in Fig. 18b. At this
point, the macroscopic response is already in the plastic regime.

(a) Geometry and boundary conditions

(b) Loading function

Figure 18: Tapered bar FE2 example with (a) geometry and boundary conditions
and (b) loading function

(a) Strain field at the end of the analysis using the PRNN

(b) Detailed view of PRNN for integration point 1

Figure 19: Strain field using the PRNN on the left and detailed view of PRNN
for a single macroscopic integration point on the right.

The strain field at the end of the analysis is shown in Fig. 19a along with
the location of one macroscopic integration point. This point is used to
illustrate the state of the PRNN throughout the time steps. Recall that the
network used in this section consists of 18 units (i.e. 6 fictitious material
points). Thus, each row in Fig. 19b corresponds to a fictitious material
point, each with its own stress path and internal variables (even though only
one of them is plotted).

In this case, each component of the macroscopic response (i.e. homogenized
stresses) is simply the linear combination of the local stresses of the 6
material models. It can be observed that the two macroscopic responses are in
excellent agreement, with minor deviations in stress components with low
magnitude. This is only visualized for a single point, but it is emphasized
that in this multiscale problem, agreement in the evolution of the stresses in
a single integration point indicates that the whole problem is solved
accurately. Moreover, the equivalent plastic deformation of each material
point is plotted in the last column. Note that despite the plastic response of
the RVE after time step 20, three of the fictitious material points of the
network ($m_{1}$, $m_{2}$ and $m_{6}$) remain in the elastic regime.

(a) Load-displacement curve obtained with the full-order FE2 approach and with
PRNN for different macroscopic mesh discretizations

(b) Average RMSE at each time step of analysis with
$\Delta_{\mathrm{elem}}^{\Omega}=$ 8 mm

Figure 20: Tapered bar FE2 example with (a) load-displacement curve with the
full-order solution and best PRNN and (b) average error of PRNN’s predictions
at each time step of the analysis with $\Delta_{\mathrm{elem}}^{\Omega}=$ 8
mm.

The accuracy of the PRNN is further assessed by inspecting the load-
displacement curve at the top edge of the bar. Fig. 20a shows the load-
displacement curve using the full-order solution and the network’s response
for different macroscopic mesh discretizations. For the refinement studied so
far ($\Delta_{\mathrm{elem}}^{\Omega}=\ $ 8 mm), it is clear that the proposed
network can capture accurately the entire nonlinear response, albeit with
minor deviations when the tapered bar changes from tension to compression and
then back again to tension.

Next, the global accuracy of the method is verified. Since the surrogate model
is not as accurate as the full-order model, a different equilibrium solution
at a certain time step affects the equilibrium in the following loading steps,
leading to accumulated error and diverging
$\bm{\varepsilon}\text{-}\bm{\sigma}$ paths. In this case, since no reduction
in the step size was observed at any moment, the simple average error between
the PRNN prediction and the full-order solution is calculated for each time
step averaged over all integrations points at the macroscale, as illustrated
in Fig. 20b. For most of the simulation, the average absolute error in the
predictions remains below 1 MPa with two peaks around 4 MPa and 6 MPa when the
loading is reversed, following the trends observed in Fig. 20a. For reference,
the lowest error of the network for test set $\mathcal{T}_{\mathrm{V}}$ is
plotted.

Table 1:  Computational cost for different mesh discretizations and efficiency of network in FE2 approach. | Macroscale element size ($\Delta_{\mathrm{elem}}^{\Omega}$) [mm] | 8 | 4 | 2
---|---|---|---|---
| Number of elements at the macroscale | 64 | 134 | 454
Online | FE2 wall-clock time [s] | 21574 | 47644 | 178800
PRNNFE2 wall-clock time [s] | 0.81 | 1.55 | 8.41
Speed-upa [-] | 26560 | 30746 | 21526
Offline | Av. wall-clock time per curve (dataset gen.) [s] | 265b | N/A | N/A
Av. training time (excl. dataset gen.) [s] | 38045b | N/A | N/A
aEvaluated as FE2 wall-clock time/PRNNFE2 wall-clock time and averaged over 5 runs |
bOne-off cost regardless of macroscopic mesh discretization |

For the purpose of assessing the efficiency of the network in accelerating the
FE2 simulations, three different levels of mesh refinement of the tapered bar
are taken into account, being the coarses the one shown in Fig. 18. Table 1
summarizes the wall-clock time spent in the analysis of the different
discretizations, as well as the speed-up in comparison to the full-order
solution. For the mesh used to illustrate this section
($\Delta_{\mathrm{elem}}^{\Omega}=$ 8 mm), replacing the solution of the BVP
of the micromodel with the network led to a speed-up over 26000 with the
accuracy reported in Fig. 20b. Considering the offline costs, the training
time is still lower than that of using the full-order solution with 134
macroscopic elements, which is a very modest number of elements for a
multiscale problem.

Since the network is trained to replace the solution of the microscopic model,
no additional training is required for the analysis of more complex cases
where the macroscale problems require more elements and time steps. Hence, in
general, higher speed-ups should be achieved with denser meshes. However, in
this particular problem, this is not always the case. An increase from the
coarsest to the intermediary mesh is observed, but no gain is achieved when
refining even further. In that case, the reduction in performance due to the
higher number of iterations caused by the necessity of adaptively reducing the
step size in order to ensure convergence. In constrast, the full-order
solution was more numerically stable for this mesh density and the adaptive-
stepping scheme was not triggered.

(a) Macroscopic strains in $x$ and $y$

(b) Macroscopic strains in $x$ and $xy$

(c) Macroscopic strains in $xy$ and $y$

Figure 21: Joint distribution of strains from the training set of the best
PRNN and strain distribution obtained by PRNNFE2 for the tapered bar problem
with different macroscopic mesh discretizations.

As the mesh is refined and strain localization takes place (see the red region
in Fig. 19a), even higher strain levels are achieved, pushing the network to
make predictions in unexplored regions during training, as illustrated in Fig.
21. Note that the network is already making far-reaching predictions in the
coarsest discretization, although in a less extensive way. In the mesh with
$\Delta_{\mathrm{elem}}^{\Omega}=\ $8 mm, the maximum strain does not exceed
0.2, while $\Delta_{\mathrm{elem}}^{\Omega}=\ $ 2 mm leads to strains higher
than 0.3. In spite of these complicating aspects, it is worth mentioning this
is still a significant speed-up. Moreover, a far less severe effect on the
global accuracy is observed, as illustrated by the almost overlapping load-
displacement curves in Fig. 20a. In that sense, the adaptive-stepping scheme
plays an important role to help overcome convergence issues.

###  7.2 Plate with multiple holes

As a final example, a composite plate with multiple cutouts with geometry and
boundary and loading conditions as illustrated in Fig. 22 is studied. Again,
an FE2 approach is employed to solve the problem for the same microscopic
model with which all the networks in Section 6 were trained for. This time, no
unloading is imposed and 150 load steps with $\Delta$s $=5.0\times 10^{-3}$
are considered. The load-displacement curve at the right edge of the plate is
plotted in Fig. 23a using both the full-order solution and the network. Again,
good agreement is observed between the macroscopic responses. The slight
inaccuracy between those are quantified in Fig. 23b, in which the average
absolute error of the component with the highest magnitudes ($\sigma_{x}$) is
around 1 MPa for almost the entire simulation.

Figure 22: Plate with cutouts: geometry and boundary and loading conditions.

(a) Load-displacement curve

(b) Average RMSE

Figure 23: Plate with cutouts: (a) load-displacement curve and (b) average
error of PRNN’s predictions at each time step of the analysis.

The displacement field at the end of the analysis is shown in Fig. 24a, where
the location of five macroscopic integration points are marked for further
inspection. The stress paths for each of these points are illustrated in Fig.
24b, where the full-order solution and the network prediction are plotted in
black and gray, respectively. Note how the stress paths are non-proportional
even for the relatively simple loading condition observed in the macroscale.
The integration point on the edge of one of the cutouts, namely point 4, is
also the one with the highest stress magnitude and the closest to a uniaxial
state in the $x$ direction while the other points experience multiaxial
loading more strongly.

(a) Displacement field at the end of the analysis with PRNN and selection of
integration points for inspection

(b) Stress-time view of macroscopic integration points. Full-order solution in
black and network’s prediction in gray

Figure 24: Plate with cutouts: (a) displacement field at the end of the
analysis and (b) selected integration points shown in stress-time view.

In terms of efficiency, the solution using the network is around 26705 times
faster than the full-order solution, which took approximately 258780 s (around
72 hours). The order of magnitude in the speed-up is similar to that obtained
in the tapered bar problem. Although no additional offline costs are incurred
because the network has been trained before for the same microscopic model and
it does not depend on the macroscopic problem at hand, it is worth stressing
that the runtime of the full-order solution exceeds the sum of the online and
offline costs of the PRNN.

Finally, this example shows that the network can capture multiaxial stress
states and non-proportional loading as obtained in FE2 simulations accurately.
No convergence issues were encountered in the PRNNFE2 simulation which points
to the smoothness of the predictions that is not always guaranteed with
surrogate models (see e.g. RNN curves in Fig. 15b).

##  8 Concluding remarks

In this paper, a novel network with embedded physics-based constitutive models
is proposed as surrogate model for the behavior of path-dependent materials in
FE2 simulations. The central idea is to address common problems in modeling
path-dependent materials using black-box models (e.g. unique mapping between
input and output and limited extrapolation abilities) by taking a step back
and reintroducing physics into the network in a way that requires very little
extra coding effort with respect to existing FE2 frameworks. This is done by
employing the same material models used for the microscopic level as part of
one of the layers of the network.

To accommodate this non-standard neural layer the following changes with
respect to standard neural network architectures are proposed. First, neurons
are assembled in groups of the size of the number of strain/stress components
of the problem. These are referred to as fictitious material points. Secondly,
to take advantage of all the information coming from the physics-based
material model, we store the updated internal variables used to fully describe
the state of the fictitious material point in an auxiliary vector. With that,
when new strain values are fed, the material point will start from the last
stored internal variables. As a consequence, each subgroup follows a unique
path without the need to increase the feature space with extra history
variables.

The properties and assumptions made by the physics-based constitutive model
are inherited by the network and play a major role in reducing the amount of
data required to mirror physical and complex behaviors such as elastic
unloading/reloading. Here, the decomposition of the strain in elastic and
plastic parts is an assumption built in the material model used to describe
the nonlinear microscopic material phase and is also observed in the network
when the local stresses of the fictitious material points are evaluated. This
simple but highly-flexible arrangement allows the network to capture arbitrary
unloading behaviors with only monotonic data, a stark contrast with other
popular models such as RNNs. The PRNN inherits from FE2 the idea that complex
behavior of heterogeneous materials can be accurately described by letting
simpler constitutive models that represent the microscopic constituents
interact. The difference is that the interaction between the constituents is
not based on micromechanics diretly but learned from data obtained from
micromodel simulations.

Based on that, an extensive numerical comparison involving a state-of-the-art
black-box model, namely a Bayesian Recurrent Neural Network (referred to in
this work as RNN), was carried out in order to elucidate the abilities of the
proposed network (referred to as PRNN). First, we trained both networks only
on 18 monotonic curves with known directions and proportional loading in a
similar fashion as done to calibrate classical mesomodels. Such strategy led
to poor performance when trying to predict other random directions from the
RNN, but good accuracy from our method (Fig. 10). Following that conclusion,
the size of the RNN’s training dataset was sequentially increased until the
addition of more data did not result in significant gains in accuracy. At that
stage, the PRNN performed with the same level of accuracy but with a factor of
16 times less data.

Next, both networks were used to predict non-monotonic loading. For that
scenario, the PRNN showed the same level of accuracy as before with the same
minimal training dataset (Fig. 11a). Such outstanding result is not observed
in the RNNs, which again required a larger training set. This time, non-
monotonic loading curves were added until the RNN’s accuracy could no longer
be significantly improved. As a result, a 32 times larger training dataset in
comparison to the one used to train our network was necessary. Furthermore,
while our network continues to perform well in all the scenarios tested so
far, two other situations exposed the pitfalls of RNNs: (i) when trying to
predict unloading in a different location than the one seen in training (Fig.
13a) and (ii) when the step size was modified (Fig. 14a). This is typically
tackled by sampling different unloading behaviors with different step sizes,
leading to the choice of arbitrarily long sequences. However, we showed that
this is a trivial scenario for the PRNN. The network is only as sensitive to
step size as the material models it includes.

In the last test, both networks were used to predict non-proportional and non-
monotonic paths and neither succeeded, although they failed at very different
levels. While the lowest error of the RNN was around 32 MPa, the best PRNN led
to an error around 9 MPa error (Fig. 15a). Based on that, a second approach to
generate the dataset was considered. Random strain paths were generated from
Gaussian Processes priors, which produces non-proportional and non-monotonic
loading as opposed to the proportional loading previously considered for
training. This time, the size of the training dataset and the type of loading
was also a variable for the PRNN. It was found that training the PRNN on
random non-proportional and non-monotonic curves yields higher accuracy than
training with known, proportional, and monotonic curves for all loading
scenarios (Fig. 17). Although training with known directions is appealing,
having a network that provides lower errors and consistent performance with
random directions is also interesting. Ultimately, the PRNN consistently
outperformed the RNN with 64 times less data.

Finally, after ensuring the PRNN capacity in several challenging scenarios for
black-box models, one of the networks trained on non-proportional and non-
monotonic curves was chosen to surrogate the microscopic model in a set of two
FE2 examples. The first example concerned a tapered bar in transverse tension
and was used to illustrate how the different material models in the PRNN
behave for a single macroscopic integration point (Fig. 18). For different
discretizations, speed-ups between 21000 and 31000 were obtained for the
online phase (Table 1). Such substantial efficiency gain is explained from the
dramatically reduced number of material model calls and the bypassing of
solving the nonlinear microscopic system of equations for macroscopic stress
evaluation for a single load step of each macroscopic integration point.

In the last example, a similar order of magnitude of speed-up was observed and
the accuracy of the PRNN was illustrated by comparing the
$\bm{\varepsilon}\text{-}\bm{\sigma}$ paths of different macroscopic
integration points of a plate with multiple cutouts subjected to tension (Fig.
24b). For the analyzed cases, the time needed for a single FE2 analysis
exceeded the total offline and online time for the PRNN analysis, even though
the selected problems had a very modest number of macroscopic elements.
Moreover, performing subsequent macroscale analysis with the same material
would not require any additional offline training and would therefore leverage
the complete speed-up of four orders of magnitude.

## Acknowledgements

The authors acknowledge the TU Delft AI Initiative for their support through
the SLIMM AI Lab. FM acknowledges financial support from the Dutch Research
Council (NWO) under Vidi grant 16464.

## References

  * Goury et al. [2016] O. Goury, D. Amsallem, S. Bordas, W. Liu, P. Kerfriden,  Automatised selection of load paths to construct reduced-order models in computational damage micromechanics: from dissipation-driven random selection to bayesian optimization,  Computational Mechanics 58 (2016) 213–234. doi:10.1007/s00466-016-1290-2.
  * Rocha et al. [2021] I. B. C. M. Rocha, P. Kerfriden, F. P. van der Meer,  On-the-fly construction of surrogate constitutive models for concurrent multiscale mechanical analysis through probabilistic machine learning,  Journal of Computational Physics: X 9 (2021) 1–30. doi:10.1016/j.jcpx.2020.100083. arXiv:2007.07749.
  * Vlassis and Sun [2021] N. N. Vlassis, W. Sun,  Sobolev training of thermodynamic-informed neural networks for interpretable elasto-plasticity models with level set hardening,  Computer Methods in Applied Mechanics and Engineering 377 (2021) 113695\. doi:https://doi.org/10.1016/j.cma.2021.113695.
  * Lefik et al. [2009] M. Lefik, D. Boso, B. Schrefler,  Artificial neural networks in numerical modelling of composites,  Computer Methods in Applied Mechanics and Engineering 198 (2009) 1785–1804. doi:https://doi.org/10.1016/j.cma.2008.12.036, advances in Simulation-Based Engineering Sciences – Honoring J. Tinsley Oden.
  * Huang et al. [2020] D. Huang, J. N. Fuhg, C. Weißenfels, P. Wriggers,  A machine learning based plasticity model using proper orthogonal decomposition,  Computer Methods in Applied Mechanics and Engineering 365 (2020) 113008\. doi:10.1016/j.cma.2020.113008. arXiv:2001.03438.
  * Ghavamian and Simone [2019] F. Ghavamian, A. Simone,  Accelerating multiscale finite element simulations of history-dependent materials using a recurrent neural network,  Computer Methods in Applied Mechanics and Engineering 357 (2019) 112594\. doi:10.1016/j.cma.2019.112594.
  * Wu et al. [2020] L. Wu, V. D. Nguyen, N. G. Kilingar, L. Noels,  A recurrent neural network-accelerated multi-scale model for elasto-plastic heterogeneous materials subjected to random cyclic and non-proportional loading paths,  Computer Methods in Applied Mechanics and Engineering 369 (2020) 113234\. doi:https://doi.org/10.1016/j.cma.2020.113234.
  * Mozaffar et al. [2019] M. Mozaffar, R. Bostanabad, W. Chen, K. Ehmann, J. Cao, M. A. Bessa,  Deep learning predicts path-dependent plasticity,  Proceedings of the National Academy of Sciences 116 (2019) 26414–26420. doi:10.1073/pnas.1911815116. arXiv:https://www.pnas.org/content/116/52/26414.full.pdf.
  * Chen [2021] G. Chen,  Recurrent neural networks (RNNs) learn the constitutive law of viscoelasticity,  Computational Mechanics 67 (2021) 1009–1019. doi:10.1007/s00466-021-01981-y.
  * Koeppe et al. [2021] A. Koeppe, F. Bamer, M. Selzer, B. Nestler, B. Markert, Explainable artificial intelligence for mechanics: physics-informing neural networks for constitutive models, 2021. arXiv:2104.10683.
  * Logarzo et al. [2021] H. J. Logarzo, G. Capuano, J. J. Rimoli,  Smart constitutive laws: Inelastic homogenization through machine learning,  Computer Methods in Applied Mechanics and Engineering 373 (2021) 113482\. doi:https://doi.org/10.1016/j.cma.2020.113482.
  * Gorji et al. [2020] M. B. Gorji, M. Mozaffar, J. N. Heidenreich, J. Cao, D. Mohr,  On the potential of recurrent neural networks for modeling path dependent plasticity,  Journal of the Mechanics and Physics of Solids 143 (2020) 103972. doi:https://doi.org/10.1016/j.jmps.2020.103972.
  * Raissi et al. [2019] M. Raissi, P. Perdikaris, G. Karniadakis,  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,  Journal of Computational Physics 378 (2019) 686–707. doi:https://doi.org/10.1016/j.jcp.2018.10.045.
  * Haghighat et al. [2021] E. Haghighat, M. Raissi, A. Moure, H. Gomez, R. Juanes,  A physics-informed deep learning framework for inversion and surrogate modeling in solid mechanics,  Computer Methods in Applied Mechanics and Engineering 379 (2021) 113741\. doi:https://doi.org/10.1016/j.cma.2021.113741.
  * Eghbalian et al. [2022] M. Eghbalian, M. Pouragha, R. Wan, A physics-informed deep neural network for surrogate modeling in classical elasto-plasticity, 2022\. doi:10.48550/ARXIV.2204.12088.
  * Arora et al. [2022] R. Arora, P. Kakkar, B. Dey, A. Chakraborty, Physics-informed neural networks for modeling rate- and temperature-dependent plasticity, 2022. doi:10.48550/ARXIV.2201.08363.
  * Masi et al. [2021] F. Masi, I. Stefanou, P. Vannucci, V. Maffi-Berthier,  Thermodynamics-based Artificial Neural Networks for constitutive modeling,  Journal of the Mechanics and Physics of Solids 147 (2021). arXiv:2005.12183.
  * Masi and Stefanou [2022] F. Masi, I. Stefanou,  Multiscale modeling of inelastic materials with thermodynamics-based artificial neural networks (tann),  Computer Methods in Applied Mechanics and Engineering 398 (2022) 115190\. doi:https://doi.org/10.1016/j.cma.2022.115190.
  * Liu et al. [2021] X. Liu, S. Tian, F. Tao, W. Yu,  A review of artificial neural networks in the constitutive modeling of composite materials,  Composites Part B: Engineering 224 (2021) 109152. doi:https://doi.org/10.1016/j.compositesb.2021.109152.
  * Liu et al. [2019] Z. Liu, C. T. Wu, M. Koishi,  A deep material network for multiscale topology learning and accelerated nonlinear modeling of heterogeneous materials,  Computer Methods in Applied Mechanics and Engineering 345 (2019) 1138–1168. doi:10.1016/j.cma.2018.09.020. arXiv:1807.09829.
  * Liu [2021] Z. Liu,  Cell division in deep material networks applied to multiscale strain localization modeling,  CoRR abs/2101.07226 (2021). arXiv:2101.07226.
  * Kerfriden et al. [2013] P. Kerfriden, O. Goury, T. Rabczuk, S. Bordas,  A partitioned model order reduction approach to rationalise computational expenses in nonlinear fracture mechanics,  Computer Methods in Applied Mechanics and Engineering 256 (2013) 169–188. doi:https://doi.org/10.1016/j.cma.2012.12.004.
  * Bessa et al. [2017] M. Bessa, R. Bostanabad, Z. Liu, A. Hu, D. W. Apley, C. Brinson, W. Chen, W. Liu,  A framework for data-driven analysis of materials under uncertainty: Countering the curse of dimensionality,  Computer Methods in Applied Mechanics and Engineering 320 (2017) 633–667. doi:https://doi.org/10.1016/j.cma.2017.03.037.
  * Oliver et al. [2017] J. Oliver, M. Caicedo, A. Huespe, J. Hernández, E. Roubin,  Reduced order modeling strategies for computational multiscale fracture,  Computer Methods in Applied Mechanics and Engineering 313 (2017) 560–595. doi:https://doi.org/10.1016/j.cma.2016.09.039.
  * Fuhg et al. [2021] J. N. Fuhg, C. Böhm, N. Bouklas, A. Fau, P. Wriggers, M. Marino,  Model-data-driven constitutive responses: Application to a multiscale computational framework,  International Journal of Engineering Science 167 (2021) 103522. doi:https://doi.org/10.1016/j.ijengsci.2021.103522.
  * Ghavamian et al. [2017] F. Ghavamian, P. Tiso, A. Simone,  Pod–deim model order reduction for strain-softening viscoplasticity,  Computer Methods in Applied Mechanics and Engineering 317 (2017) 458–479. doi:https://doi.org/10.1016/j.cma.2016.11.025.
  * Rocha et al. [2018] I. B. C. M. Rocha, F. P. van der Meer, L. Sluys,  Efficient micromechanical analysis of fiber-reinforced composites subjected to cyclic loading through time homogenization and reduced-order modeling,  Computer Methods in Applied Mechanics and Engineering 345 (2018). doi:10.1016/j.cma.2018.11.014.
  * Rocha et al. [2020] I. B. C. M. Rocha, P. Kerfriden, F. P. van der Meer,  Micromechanics-based surrogate models for the response of composites: A critical comparison between a classical mesoscale constitutive model, hyper-reduction and neural networks,  European Journal of Mechanics, A/Solids 82 (2020) 103995. doi:10.1016/j.euromechsol.2020.103995.
  * Vijayaraghavan et al. [2021] S. Vijayaraghavan, L. Wu, L. Noels, S. P. A. Bordas, S. Natarajan, L. A. A. Beex, Neural-network acceleration of projection-based model-order-reduction for finite plasticity: Application to rves, 2021\. arXiv:2109.07747.
  * Guo et al. [2021] T. Guo, O. Rokoš, K. Veroy,  Learning constitutive models from microstructural simulations via a non-intrusive reduced basis method,  Computer Methods in Applied Mechanics and Engineering 384 (2021) 113924\. doi:https://doi.org/10.1016/j.cma.2021.113924.
  * Nguyen et al. [2012] V. P. Nguyen, O. Lloberas-Valls, M. Stroeven, L. J. Sluys,  Computational homogenization for multiscale crack modeling. Implementational and computational aspects,  International Journal for Numerical Methods in Engineering 89 (2012) 192–226. doi:10.1002/nme.3237.
  * Kingma et al. [2015] D. P. Kingma, T. Salimans, M. Welling, Variational dropout and the local reparameterization trick, 2015. doi:10.48550/ARXIV.1506.02557.
  * Kingma and Ba [2014] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization, 2014. doi:10.48550/ARXIV.1412.6980.
  * Hernández et al. [2017] J. Hernández, M. Caicedo, A. Ferrer,  Dimensional hyper-reduction of nonlinear finite element models via empirical cubature,  Computer Methods in Applied Mechanics and Engineering 313 (2017) 687–722. doi:https://doi.org/10.1016/j.cma.2016.10.022.
  * Mori and Tanaka [1973] T. Mori, K. Tanaka,  Average stress in matrix and average elastic energy of materials with misfitting inclusions,  Acta Metallurgica 21 (1973) 571–574. doi:https://doi.org/10.1016/0001-6160(73)90064-3.
  * Eshelby and Peierls [1957] J. D. Eshelby, R. E. Peierls,  The determination of the elastic field of an ellipsoidal inclusion, and related problems,  Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 241 (1957) 376–396. doi:10.1098/rspa.1957.0133. arXiv:https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1957.0133.
  * Nguyen-Thanh et al. [2020] C. Nguyen-Thanh, V. P. Nguyen, A. de Vaucorbeil, T. Kanti Mandal, J.-Y. Wu,  Jive: An open source, research-oriented c++ library for solving partial differential equations,  Advances in Engineering Software 150 (2020) 102925. doi:https://doi.org/10.1016/j.advengsoft.2020.102925.
  * van der Meer [2012] F. P. van der Meer,  Mesolevel modeling of failure in composite laminates: Constitutive, kinematic and algorithmic aspects,  Archives of Computational Methods in Engineering 19 (2012) 381–425. doi:10.1007/s11831-012-9076-y.