# Thermodynamics-based Artificial Neural Networks for constitutive modeling

Filippo Masi<EMAIL_ADDRESS>Ioannis Stefanou  ioannis.stefanou@ec-
nantes.fr Paolo Vannucci<EMAIL_ADDRESS>Victor Maffi-Berthier
<EMAIL_ADDRESS>Institut de Recherche en Génie Civil et
Mécanique,
UMR 6183, CNRS, Ecole Centrale de Nantes, Université de Nantes,
1 rue de la Nöe, F-44300, Nantes, France.
Ingérop Conseil et Ingénierie,
18 rue des Deux Gares, F-92500, Rueil-Malmaison, France.
LMV, UMR 8100, Université de Versailles et Saint-Quentin,
55 avenue de Paris, F-78035, Versailles, France.

###### Abstract

Machine Learning methods and, in particular, Artificial Neural Networks (ANNs)
have demonstrated promising capabilities in material constitutive modeling.
One of the main drawbacks of such approaches is the lack of a rigorous frame
based on the laws of physics. This may render physically inconsistent the
predictions of a trained network, which can be even dangerous for real
applications.

Here we propose a new class of data-driven, physics-based, neural networks for
constitutive modeling of strain rate independent processes at the material
point level, which we define as Thermodynamics-based Artificial Neural
Networks (TANNs). The two basic principles of thermodynamics are encoded in
the network’s architecture by taking advantage of automatic differentiation to
compute the numerical derivatives of a network with respect to its inputs. In
this way, derivatives of the free-energy, the dissipation rate and their
relation with the stress and internal state variables are hardwired in the
network. Consequently, our network does not have to identify the underlying
pattern of thermodynamic laws during training, reducing the need of large
data-sets. Moreover the training is more efficient and robust, and the
predictions more accurate. Finally and more important, the predictions remain
thermodynamically consistent, even for unseen data. Based on these features,
TANNs are a starting point for data-driven, physics-based constitutive
modeling with neural networks.

We demonstrate the wide applicability of TANNs for modeling elasto-plastic
materials, with strain hardening and strain softening. Detailed comparisons
show that the predictions of TANNs outperform those of standard ANNs. TANNs ’
architecture is general, enabling applications to materials with different or
more complex behavior, without any modification.

###### keywords:

Data-driven modeling; Machine learning; Artificial neural network;
Thermodynamics; Constitutive model.

††journal: Journal

##  1 Introduction

A large spectrum of constitutive models have been proposed in the literature,
based on observations and experimental testing. Existing constitutive laws can
account for phenomena taking place in various length scales. This is achieved
either through heuristic approaches and assumptions or through asymptotic
approximations and averaging (e.g. [38, 48, 15, 3]). The history and the state
of a material is commonly taken into account through ad hoc enrichment of
simpler constitutive laws and extensive calibration. For this purpose, the
laws of thermodynamics offer a useful framework for deriving more
sophisticated laws, by intrinsically respecting the energy balance and the
entropy production requirements (see e.g. [25, 14, 24] among others).

An important limitation in constitutive modeling is the availability of data
at different time- and length-scales. However, with the increase of
computational power, it is nowadays possible to foresee micromechanical
simulations that can account for realistic physics and explore stress paths
and non-linear phenomena, which are experimentally inaccessible with the
current methods. Of course, some constitutive assumptions will be always
necessary, but these might be at a smaller scale, where the material
properties are measurable and easier to identify. This scale is for instance
the scale of the microstructure of a material (e.g. the scale of sand grains,
crystals, alloys’ grains, composites’ fibers, masonry bricks’ etc. including
their topological configuration).

However, it is likely that the existing constitutive models might not be
sufficient for describing complex material behaviors emerging from the
microstructure. Therefore, calibration (parameter fitting) of known
constitutive descriptors might be insufficient for representing the full space
of material response, provided by sophisticated micromechanical simulations.
Moreover, micromechanical simulations have currently a tremendous calculation
cost, which is impossible to afford in large-scale, non-linear, incremental
simulations (e.g. Finite Elements) that are usually needed in applications
(cf. [41, 42, 50, 51, 8, 38, 48, 12, 15]).

A promising solution to this issue seems to be Machine Learning. According to
[18], “Machine Learning is the science (and art) of programming computers so
they can learn from data”. In the context of computer programming, learning is
defined by E. Tom Mitchell [46] as follows: “A computer program is said to
learn from experience $E$ with respect to some task $T$ and some performance
measure $P$, if its performance on $T$, as measured by $P$, improves with
experience”. In the frame of constitutive modeling, a Machine Learning program
can learn the stress-strain behavior of a material, given examples of stress-
strain increments, which are either determined experimentally or through
detailed micromechanical simulations. The data that the system uses to learn
are called the training data-set and each training example is called a
training instance (or sample). In our case, the task $T$, for instance, can be
the prediction of the stress for a given increment and internal state of the
material. The experience $E$ is the training data-set and the performance
measure $P$ can be the prediction error. Machine Learning is a general term to
describe a large spectrum of numerical methods. Some of them offer very rich
interpolation spaces, which, in theory, could be used for approximating
complicated functions belonging to uncommon spaces. Here we focus on the
method of Artificial Neural Networks (ANNs), which is considered to be a sub-
class of Machine Learning methods. According to Chen and Chen [7] and Cybenko
[9], ANNs have proved to be universal approximators, due to their rich
interpolation space. Therefore, they seem to be a useful and promising tool
for data-driven constitutive modeling of many materials (e.g. sand, masonry,
alloys, ceramics, composites etc.).

Recognizing this potential, there is an increasing amount of new literature
employing ANNs in constitutive modeling. Starting form the seminal work of
Ghaboussi et al. [20] and without being exhaustive, we refer to Lefik and
Schrefler [36], Ghaboussi and Sidarta [21], Jung and Ghaboussi [30], Heider et
al. [23], Settgast et al. [52], Ghavamian and Simone [22], Liu and Wu [37], Lu
et al. [39], Xu et al. [54], Huang et al. [27], Mozaffar et al. [47], Frankel
et al. [16], Liu and Wu [37], Gajek et al. [17] and references therein. The
main idea in these works is to appropriately train ANNs, feeding them with
material data, and predict the material response at the material point level.
In this sense ANNs can be seen as rich interpolation spaces, able to represent
complex material behavior. The Boundary Value Problem (BVP), set to determine
the behavior of a solid under mechanical and/or multiphysics couplings, is
then solved by replacing the standard constitutive equations or algorithms by
the trained ANN. This replacement is straightforward and non-intrusive in
Finite Element codes. It is worth emphasizing that the aforementioned data-
driven approaches are different from another promising data-driven method
(i.e., data driven computing [32]) in which the BVP is solved directly from
experimental material data (measurements), bypassing the empirical material
modeling step, involving the calibration of constitutive parameters [32, 29,
33, 34, 28, 11]. While data-driven computing can be extremely powerful in many
applications [11], the first class of methods above-mentioned (based on the
constitutive behavior at the material point level) can be advantageous when
modeling complex and abstract constitutive behaviors, which are not a priori
known. Moreover, they can be used even if the BVP does not have a unique
solution due to important non-linearities and bifurcation phenomena (e.g. loss
of uniqueness, strain localization at the length of interest, multiphysics,
runway instabilities etc.).

Nevertheless, until now ANNs for constitutive modeling are mainly used as a
‘black-box’ mathematical operator, which once trained on available data-sets,
does not embody the basic laws of thermodynamics. As a result, vast amount of
high quality data (e.g. with reduced noise and free of outliers) are needed to
enable ANNs to identify and learn the underlying thermodynamic laws. Moreover,
nothing guarantees that the predictions of trained ANNs will be
thermodynamically consistent, especially for unseen data.

In this paper, we encode in the ANN architecture the two basic laws of
thermodynamics. This assures thermodynamically consistent predictions, even
for unseen data (which can exceed the range of training data-sets). Therefore,
we assure thermodynamically consistent network’s predictions, both for seen
and unseen data (which can exceed the range of the training data-sets).
Moreover, our network does not have to identify/learn the underlying pattern
of thermodynamical laws. Consequently, smaller data-sets are needed in
principle, the training is more efficient and the accuracy of the predictions
higher. The price to pay, in comparison with existing approaches, is the need
of two additional scalar functions (outputs) in the training data-set. These
are the free-energy and the dissipation rate. However, these quantities are
easily accessible in micromechanical simulations (e.g. [48, 12, 15]) and can
also be obtained experimentally in some cases. Then, based on classical
derivations in thermodynamics (e.g. [25, 13]) specific interconnections are
programmed inside our ANN architecture to impose the necessary thermodynamic
restrictions. These thermodynamic restrictions concern the stresses and
internal state variables and their relation with the free-energy and the
dissipation rate. Our approach is inspired by the so-called Physics-Informed
Neural Networks (PINNs) [49], in which reverse-mode autodiff [4] is used,
allowing the numerical calculation of the derivative of an ANN with respect to
its inputs.

The calculation of these derivatives, imposes some numerical requirements
regarding the mathematical class of the activation functions to be used. More
specifically, the internal ANN restrictions, derived from the first law of
thermodynamics, require activation functions whose second gradient does not
vanish. Otherwise, the problem of second-order vanishing gradients, as it is
called here (cf. classical vanishing gradients problem in ANNs [e.g. 18]), can
inhibit back-propagation and make training to fail. This new problem and its
remedy is extensively explored and discussed herein.

For the sake of simplicity and for distinguishing our approach from existing
ones, we call the proposed ANN architecture Thermodynamics-based Artificial
Neural Networks (TANNs). In our opinion TANN should be the starting point for
data-driven and physics-based constitutive modeling at the material point
level. For the implementation of TANNs, we leverage Tensorflow v2.0, an open-
source symbolic tensor manipulation software library [1]. Other
libraries/packages can be used as well.

The paper is structured as follows. Section 2 presents a brief summary of the
theoretical background of thermodynamics. In Section 3 an overview of the
methodology proposed and architecture of TANN is given. The main differences
with classical, standard ANNs for material constitutive modeling are also
discussed. Particular attention is given to the choice of activation functions
and the issue of second-order vanishing gradient is investigated in detail.
Generation of material data-sets, with which the training of ANNs is
performed, is presented in Section 4. Some applications of TANN for uni-
dimensional and three-dimensional elasto-plastic material models are presented
in Section 5. Extensive comparisons with standard ANNs which are not based on
thermodynamics are also presented. Supplementary figures and data are
available in Supplementary Material (SM) file. All code accompanying this
manuscript is available on request.

##  2 Thermodynamics principles: energy conservation and dissipation
inequality

###  2.1 Energy conservation

A convenient way to express the (local) energy conservation is

$\rho\dot{\text{{e}}}=\sigma\cdot\text{D}^{\text{Sym}}\text{v}-\text{div}q+\rho
h,$ (1)

with $\rho$ being the material density; e the specific internal energy (per
unit mass); $\sigma$ the Cauchy stress tensor; Dv the spatial velocity
gradient tensor; $q$ the rate of heat flux per unit area; $h$ the specific
energy source (supply) per unit mass, and ”$\cdot$” denotes contraction of
adjacent indices.

###  2.2 Second principle

The second law of thermodynamics can be formulated in terms of the local
Clausius-Duhem inequality

$\rho\dot{\text{{s}}}\geq\frac{\rho h}{\theta}-\text{div}\,\left(\frac{q\cdot
n}{\theta}\right),$ (2)

with s being the specific (per unit mass) entropy; $h/\theta$ and $-(q\cdot
n)/\theta$ the rate of entropy supply and flux, respectively. By removing the
heat supply $h$ between the energy equation (1) and the entropy inequality (2)
leads to

$\rho\left(\theta\dot{\text{{s}}}-\dot{\text{{e}}}\right)+\sigma\cdot\text{D}^{\text{Sym}}\text{v}-\dfrac{q\cdot\text{D}\theta}{\theta}\geq
0,$ (3)

where the first two terms represent the rate of mechanical dissipation
$\text{{D}}=\rho\left(\theta\dot{\text{{s}}}-\dot{\text{{e}}}\right)+\sigma\cdot\text{D}^{\text{Sym}}\text{v}$
and the latter the thermal dissipation rate, i.e.,
$\text{{D}}^{th}=-\frac{q\cdot\text{D}\theta}{\theta}$. The thermal
dissipation is non-negative because heat only flows from regions of higher
temperature to lower temperature$-$that is, the heat flux $q$ is always in the
direction of the negative thermal gradient. As it follows we argue that the
mechanical dissipation rate must itself be non-negative (point-wise), i.e.,
$\text{{D}}\geq 0$.

###  2.3 Dissipation function

The definition of the (mechanical) dissipation rate D leads to

$\rho\dot{\text{{e}}}=\rho\theta\dot{\text{{s}}}+\sigma\cdot\text{D}^{\text{Sym}}\text{v}-\text{{D}}.$
(4)

Let define the specific (per unit volume) internal energy
$\text{{E}}=\rho\text{{e}}$ and entropy $\text{{S}}=\rho\text{{s}}$ and
further assume constant material density, i.e.,
$\frac{\text{d}}{\text{d}t}\rho=0-$that is,
$\dot{\text{{E}}}=\rho\dot{\text{{e}}}$ and
$\dot{\text{{S}}}=\rho\dot{\text{{s}}}$. We shall assume a small strain
regime, i.e., $\text{D}u\ll 1$, with $\varepsilon:=\text{D}^{\text{Sym}}u$ the
small strain tensor, where $u$ is the displacement vector field, and
$\dot{\varepsilon}:=\text{D}^{\text{Sym}}\text{v}$ its rate of change.
Equation (4) hence becomes

$\dot{\text{{E}}}=\theta\dot{\text{{S}}}+\sigma\cdot\dot{\varepsilon}-\text{{D}}.$
(5)

Let assume a strain-rate independent material such that

$\text{{E}}:=\widetilde{\text{{E}}}\left(\text{{S}},\varepsilon,\mathcal{Z}\right),$
(6)

and

$\text{{D}}:=\widetilde{\text{{D}}}\left(\text{{S}},\varepsilon,\mathcal{Z},\dot{\mathcal{Z}}\right),$
(7)

where $\mathcal{Z}=(\zeta_{i},\dots,\zeta_{N})$ denotes a set of $N$
(additional) internal state variables, $\zeta_{i}$, $i=1,\dots,N$. We define
here (thermodynamic) state variables those macroscopic quantities
characterizing the state of a system, see e.g. [43]. The physical
representation of $\zeta_{i}$ is not a priori prescribed. For instance, in the
case of isotropic damage, $\zeta$ is a scalar; in anistotropic damage, a
tensor; in the case of elasto-plasticity, a second order tensor, etc. The
generalization to a finite-strain formulation can be achieved by considering
the deformation gradient, $F$, and the first Piola-Kirchhoff tensor, $P$, as
strain and stress measures, respectively (see e.g. [40] and [2]).
Nevertheless, as it would presented in Section 3, an incremental formulation
of the material response is herein adopted. Therefore, the hypothesis of a
small strain regime is usually realistic, at least for a large class of
materials and an updated Lagrangian scheme.
Time differentiation of the internal energy gives

$\dot{\text{{E}}}=\dfrac{\partial\text{{E}}}{\partial\text{{S}}}\cdot\dot{\text{{S}}}+\dfrac{\partial\text{{E}}}{\partial\varepsilon}\cdot\dot{\varepsilon}+\sum_{i=1}^{N}\dfrac{\partial\text{{E}}}{\partial\zeta}_{i}\cdot\dot{\zeta}_{i},$
(8)

which is equal to (5) and, grouping terms, it leads to

$\left(\dfrac{\partial\text{{E}}}{\partial\text{{S}}}-\theta\right)\dot{\text{{S}}}+\left(\dfrac{\partial\text{{E}}}{\partial\varepsilon}-\sigma\right)\cdot\dot{\varepsilon}-\left(\sum_{i=1}^{N}\dfrac{\partial\text{{E}}}{\partial\zeta}_{i}\cdot\dot{\zeta}_{i}+\text{{D}}\right)=0.$
(9)

The arbitrariness of $\dot{\text{{S}}}$, $\dot{\varepsilon}$, and
$\dot{\zeta}$ leads to the following relations

$\displaystyle\theta=\dfrac{\partial\text{{E}}}{\partial\text{{S}}},$ (10a)
$\displaystyle\sigma=\dfrac{\partial\text{{E}}}{\partial\varepsilon},$ (10b)
$\displaystyle\sum_{i=1}^{N}\dfrac{\partial\text{{E}}}{\partial\zeta}_{i}\cdot\dot{\zeta}_{i}+\text{{D}}=0.$
(10c)

Further introducing the thermodynamic stress, conjugate to $\zeta_{i}$,
$\mathcal{X}=(\chi_{1},\ldots,\chi_{N})$, with

$\chi_{i}:=-\frac{\partial\text{{E}}}{\partial\zeta_{i}}\qquad\forall\>i\in\left[1,N\right],$
(11)

we obtain the following, alternative definition of the dissipation

$\text{{D}}=\sum_{i=1}^{N}\chi_{i}\cdot\dot{\zeta}_{i}$ (12)

###  2.4 Isothermal processes

In the case of isothermal process, the (specific) Helmholtz free-energy,
$\text{{F}}:=\text{{E}}-\text{{S}}\theta=\widetilde{\text{{F}}}(\theta,\varepsilon,\mathcal{Z})$,
which is the Legendre transform conjugate of e, is preferable. In this case,
the dissipation rate is such that
$\text{{D}}:=\widetilde{\text{{D}}}(\theta,\varepsilon,\mathcal{Z},\dot{\mathcal{Z}})$.
The equations presented above (9-31) still hold (by replacing E with F)

$\text{{S}}=-\dfrac{\partial\text{{F}}}{\partial\theta},\quad\sigma=\dfrac{\partial\text{{F}}}{\partial\varepsilon},\quad\text{{D}}=-\sum_{i}\dfrac{\partial\text{{F}}}{\partial\zeta_{i}}\cdot\dot{\zeta}_{i}=\sum_{i}\chi_{i}\cdot\dot{\zeta}_{i}.$
(13)

##  3 Thermodynamics-based Artificial Neural Networks

Within the framework of ANN material models, we can distinguish two main
classes. The first consists of direct, so-called “black-box”, approaches,
where the information flow passes through the machine learning tool (usually
feed-forward or recurrent artificial neural networks) which operates as a mere
regression operator, see e.g. [20, 36]. The second class coincides with ANN
models incorporating some knowledge in an informed, guided graph, see e.g.
[23]. Both classes, however, are affected by the lack of physics, being the
predictions not always compatible with thermodynamic principles (at least).

(a) black-box (BB) network.

(b) informed neural network (i-NN1).

(c) informed neural network (i-NN1).

Figure 1: Examples of direct, black-box (BB) (left) and informed (right)
neural networks for material laws modeling. Inputs are highlighted in gray (
), outputs in black ( ).

Figure 1(a) depicts the direct approach (BB), in which ANNs, either forward
fed or recurrent, are used to predict the stress increment, (output,
$\mathcal{O}$) $\mathcal{O}=\Delta\sigma=\sigma^{t+\Delta t}-\sigma^{t}$, from
the input $\mathcal{I}=(\varepsilon^{t},\Delta\varepsilon)$, being
$\varepsilon^{t}$ the precedent strain state and $\Delta\varepsilon$ its
increment. In concise form, we write $\mathcal{O}=\textsf{BB}@\,\mathcal{I}$.
In this scheme, $\varepsilon^{t}$ and $\Delta\varepsilon$ can be regarded as
the state variables, namely the ANN state variables (not necessarily
coinciding with those introducted in Sect. 2), on which the updated material
stress depends on. Two examples of guided, informed ANNs are illustrated in
Figures 1(b) and 1(c). In both cases, the ANN intrinsically accounts for path-
dependency, see e.g. [23], making sequence of predictions of the main output.
The network i-NN1 makes use of the last predicted output, i.e., $\sigma^{t}$,
to make predictions of the next output, $\mathcal{O}=\Delta\sigma$. The inputs
are hence $\mathcal{I}=(\varepsilon^{t},\Delta\varepsilon,\sigma^{t})$. We
shall notice that, differently from BB, the stress at the precedent state,
$\sigma^{t}$, is also considered to be an ANN state variable. Other
alternatives exist in the selection of the ANN variables of state. One may
chose, as we shall see in Section 5, (thermodynamic) state variables to be ANN
state variables.
In the case of temperature-dependent material response, the second case
(i-NN2) allows to make predictions that depend on the precedent temperature
state, $\theta^{t}$, namely $\mathcal{O}=\textsf{i-NN2}@\,\mathcal{I}$, with
$\mathcal{I}=(\varepsilon^{t},\Delta\varepsilon,\sigma^{t},\theta^{t})$ and
$\mathcal{O}=(\Delta\sigma,\Delta\theta)$.

The main aim of this work is to change the classical paradigm of data-driven
ANN material modeling into physics-based ANN material modeling. We propose a
new class of ANN based on thermodynamics, which are thermodynamics-based
artificial neural networks (TANN). By exploiting the theoretical background
presented in Section 2, we propose neural networks which, by definition,
respect the thermodynamic principles, holding true for any class of material.
In this framework, TANN posses the special feature that the entire
constitutive response of a material can be derived from definition of only two
potential functions: an energy potential and a dissipation (pseudo-) potential
[25]. TANNs are fed with thermodynamics ”information” by relying on the
automatic differentiation technique [4] to differentiate neural networks
outputs with respect to their inputs. This strategy allows to construct a
general framework of neural networks material models which, in principle, can
be exploited to predict the behavior of any material and assure that the
predictions of TANN will be thermodynamically consistent even for inputs that
exceed the training range of data. In this paper, we only focus on strain-rate
independent processes. Moreover, our approach can be extended, following the
developments in [24], to materials showing viscosity and strain-rate
dependency.
The model relies on an incremental formulation and can be used in existing
Finite Element formulations (among others), see e.g. [36]. Figure 2
illustrates the scheme of TANN. The model inputs are the strain increment, the
previous material state at time $t$, which is identified herein through the
material stress, $\sigma^{t}$, temperature, $\theta^{t}$, and the internal
state variables, $\zeta_{i}^{t}$, as well as the time increment $\Delta t$,
namely
$\mathcal{I}=(\varepsilon^{t},\Delta\varepsilon,\sigma^{t},\theta^{t},\zeta_{i}^{t},\Delta
t)$. The primary outputs, $\mathcal{O}_{1}$, are internal variables increment,
$\Delta\zeta_{i}$, the temperature increment, $\Delta\theta$, and the energy
potential at time $t+\Delta t$, $\textsf{F}^{t+\Delta t}$, i.e.
$\mathcal{O}_{1}=(\Delta\zeta_{i},\Delta\theta,\textsf{F}^{t+\Delta t})$.
Secondary outputs, $\mathcal{O}_{2}-$that is, outputs computed by
differentiation of the neural network with respect to the inputs$-$are the
stress increment, $\Delta\sigma$, and the dissipation rate,
$\textsf{D}^{t+\Delta t}$, which we denote as
$\mathcal{O}_{2}=\nabla_{\mathcal{I}}\mathcal{O}_{1}=(\Delta\sigma,\textsf{D}^{t+\Delta
t})$.
The class of neural network we propose differs from the previous ones by the
fact that the quantity of main interest, i.e., the stress increment, is
obtained as a derived one, which intrinsically satisfies the first principle
of thermodynamics (and, as we shall see, the second principle, as well). In
the following, we briefly recall the basic concepts of artificial neural
networks (paragraph 3.1), we then focus on the issue of the second-order
vanishing gradients that may afflict the training and the performance of an
ANN model (paragraph 3.2). In particular, it is shown that, in the framework
of thermodynamics-based artificial neural networks, particular attention has
to be paid to the selection of activation functions. Finally, we present in
detail the architecture of our model (paragraph 3.3).

Figure 2: Schematic architecture of TANN. Inputs are highlighted in gray ( );
outputs in black, ( ) and ( ); and intermediate quantities in white ( ).
Dashed lines represent definitions, while arrows are used to denote ANN.

###  3.1 Artificial neural networks overview

We give herein a brief overview of the basic concepts of artificial neural
networks (ANNs). For more details, we refer to [26, 19]. ANNs can be regarded
as non-linear operators, composed of an assembly of mutually connected
processing units$-$nodes$-$, which take an input signal $\mathcal{I}$ and
return the output $\mathcal{O}$, namely

<EMAIL_ADDRESS>(14)

ANN consist of at least three types of layers: input, output and hidden
layers, with equal or different number of nodes. Figure 3 depicts a network
composed of one hidden layer, with 3 nodes, an input layer with 2 inputs, and
an output layer with 1 node. When an ANN has two or more hidden layers, it is
called a deep neural network [19]. Denoting the input array with
$\mathcal{I}=\left(i_{t}\right)$, with $t=1,2\ldots,n_{\mathcal{I}}$
($n_{\mathcal{I}}$ is the number of inputs), and the outputs with
$\mathcal{O}=\left(o_{j}\right)$, with $j=1,2\ldots,n_{\mathcal{O}}$
($n_{\mathcal{O}}$ is the number of outputs), the signal flows from layer
$(l-1)$ to layer $(l)$ according to

$p^{(l)}_{k}=\mathcal{A}^{(l)}\left(z_{k}^{(l)}\right),\quad\text{with}\quad
z_{k}^{(l)}=\sum_{s}^{n^{(l-1)}_{\mathcal{N}}}\left(w^{(l)}_{ks}p^{(l-1)}_{s}\right)+b^{(l)}_{k},$
(15)

where $p^{(l)}_{k}$ is the output of node $k$, at layer $(l)$;
$\mathcal{A}^{(l)}$ is the activation function of layer $(l)$;
$n^{(l-1)}_{\mathcal{N}}$ is the number of neurons in layer $(l-1)$;
$w^{(l)}_{ks}$ is the weight between the $s$-th node in layer $(l-1)$ and the
$k$-th node in layer $(l)$; and $b^{(l)}_{k}$ are the biases of layer $(l)$.
With reference to Figure 3, the output is given by

$\displaystyle\mathcal{O}=$
$\displaystyle\mathcal{A}^{(\text{o})}\left(z^{(\text{o})}\right)\quad\text{with}\quad
z^{(\text{o})}=\sum_{r}w_{r}^{(2)}p_{r}^{(1)}+b^{(2)}$ $\displaystyle
p_{r}^{(1)}=$
$\displaystyle\mathcal{A}^{(1)}\left(z^{(1)}_{r}\right)\quad\text{with}\quad
z^{(1)}_{r}=\sum_{t}w_{rt}^{(1)}i_{t}+b_{r}^{(1)},$

where the activation function of the output layer,
$\mathcal{A}^{(\text{out})}$, in a regression problem, is a linear function,
in the most part of applications. The weights and biases of interconnections
are adjusted, in an iterative procedure (gradient descent algorithm [19]), to
minimize the error between the benchmark, $\mkern
1.5mu\overline{\mkern-1.5mu\mathcal{O}\mkern-1.5mu}\mkern 1.5mu$, and
prediction, $\mathcal{O}$, that is measured by a loss function, $\mathcal{L}$.
In the following, the Mean (over a set of $N$ samples) Absolute Error is used
as loss function, i.e.,

$\mathcal{L}=\dfrac{\sum_{i=1}^{N}|\mkern
1.5mu\overline{\mkern-1.5mu\mathcal{O}\mkern-1.5mu}\mkern
1.5mu_{i}-\mathcal{O}_{i}|}{N},$ (16)

where $i=1,2,\ldots N$. The errors related to each node of the output layer
are hence back-propagated to the nodes in the hidden layers and used to
calculate the gradient of the loss function, namely

$\dfrac{\partial\mathcal{L}}{\partial
w_{ks}^{(l)}}=\dfrac{\partial\mathcal{L}}{\partial o_{j}}\dfrac{\partial
o_{j}}{\partial z_{r}^{(l+m)}}\dfrac{\partial z_{r}^{(l+m)}}{\partial
p_{r}^{(l+m-1)}}\cdots\dfrac{\partial z_{k}^{(l+1)}}{\partial
p_{k}^{(l)}}\dfrac{\partial p_{k}^{(l)}}{\partial w_{ks}^{(l)}},$ (17)

which is also used to update weights and biases, and force the minimization of
the loss function values, i.e.

$w_{ks}^{(l)-{\text{new}}}:=w_{ks}^{(l)}-\epsilon\dfrac{\partial\mathcal{L}}{\partial
w_{ks}^{(l)}},$ (18)

where $\epsilon$ is the so-called learning rate.

Figure 3: Graph illustration of an ANN structure with two inputs, one output,
and one hidden layer with three nodes.

The weights and biases updating, the so-called training process, is performed
on a subset of the input-output data-set, defined as training set, known from
experimental tests or numerical simulations of the phenomenon investigated.
The ANN is trained. The training process is stopped as the loss function is
below a specific tolerance. Then a test set, a subset of the input-output
data-set different to the training set, is used to check the error of the
network predictions. Once the ANN is trained, it is used in recall mode to
obtain the output of the problem at hand.
Due to their rich interpolation space, ANNs have proved to be universal
approximators, see e.g. [7, 9], although the choice of hyper-parameters, such
as the number of neurons, the network topology, the weights, etc. are problem-
dependent. The same stands for the activation functions, which may be chosen
to have some desirable properties of non-linearity, differentiation,
monotonicity, etc. Most of these properties mostly stem from issues related to
the gradient descent algorithm and the so-called (first-order) vanishing
gradient problem. As it follows, we briefly present this well-known issue and
we further give insights in a variation of it: the second-order vanishing
gradient.

###  3.2 First- and second-order vanishing gradients

During the training process, if the gradient of the loss function with respect
to a certain weight tends to zero$-$that is, see Eq. (18), when
$\mathcal{A^{\prime}}^{(l)}=\partial p_{j}^{(l)}/\partial z_{j}^{(l)}\approx
0$ (with $\mathcal{A}^{\prime}$ the first-derivative of the activation
function with respect to its arguments)$-$the update operation can fail, and
the weight and biases values are not updated. In this case, we have the so-
called first-order vanishing gradient [19].

Figure 4: Some of the most common activation functions and their first-order
gradient. From left to right: the logistic (sigmoid) function, the hyperbolic
tangent, the Rectified Linear Unit (ReLU), and the Exponential Linear Unit
(ELU).

Figure 4 displays some of the most common activation functions and their
derivatives$-$that is, the logistic (sigmoid) function, the hyperbolic
tangent, the Rectified Linear Unit (ReLU), and the Exponential Linear Unit
(ELU). The sigmoid function is S-shaped, continuous, differentiable, its
output values range from 0 to 1, and its first-order gradient (derivative)
assumes values much smaller than 1. When inputs become large (negative or
positive), the function saturates at 0 or 1, with a derivative extremely close
to 0. Thus when backpropagation kicks in, it has virtually no gradient to
propagate back through the network, which is problematic for training. The
hyperbolic tangent activation function is very similar to the sigmoid, but it
is centered at zero allowing to maintain the output values within a normalized
range (between -1 and 1). Nevertheless, it suffers from saturated gradients
(at $z=0$, for $z<<-1$ and $z>>1$). ReLU is continuous but not differentiable
at $z=0$. Nevertheless it is an unsaturated activation function for positive
values of $z$ (its gradient has no maximum) and, therefore, it allows to avoid
vanishing gradient issues for $z>0$. Nevertheless, it suffers from a problem
known as the dying ReLUs: during training, some neurons are effectively
deactivated, meaning they stop outputting anything other than 0 (for $z<0$).
To this purpose many variants exist. The ELU activation, for instance, takes
on negative values when $z<0$, which allows the unit to have an average output
closer to 0. This helps alleviate the vanishing gradient problem, as discussed
earlier. Second, it has a nonzero gradient for $z<0$, which avoids the dying
units issue. Finally, the function is smooth everywhere, including z = 0,
which helps speed up gradient descent.

When dealing with TANN, second-order vanishing gradients can appear. This is a
new concept and, in order to illustrate it, we will use a simle example.
Assume an ANN which takes as input some $\mathcal{I}=x$ and returns (a)
$\mathcal{O}_{1}=x^{2}$ and (b) its derivative with respect to the input,
i.e., $\mathcal{O}_{2}=\nabla_{\mathcal{I}}\mathcal{O}_{1}=2x$ (see Figure 5).
Let consider one hidden layer, with activation function $\mathcal{A}$ and
$N_{n}$ nodes. The activation function of the single output layer, which
returns $x^{2}$, is assumed to be linear. In this case, the output (a) is
given by

$\begin{split}\mathcal{O}_{1}=&\>p^{(\text{o})}=\mathcal{A}^{(\text{o})}\left(z_{k}^{(\text{o})}\right)\\\
\mathcal{O}_{1}=&\>\sum_{j}w_{j}^{(\text{o})}p_{j}^{(1)}+b^{(\text{o})}\\\
\mathcal{O}_{1}=&\>\sum_{j}w_{j}^{(\text{o})}\mathcal{A}\left(w_{j}^{(1)}i+b_{j}^{(1)}\right)+b^{(\text{o})}.\end{split}$
(19)

The derivatives of the outputs with respect to the inputs can be easily
computed, in this simple example, by taking advantage of the automatic
(numerical) differentiation [4]. Output (b) is hence computed by the ANN as

$\begin{split}\mathcal{O}_{2}=\nabla_{\mathcal{I}}\mathcal{O}_{1}=\frac{\partial\mathcal{O}_{1}}{\partial\mathcal{I}}=&\sum_{j}\dfrac{\partial
p^{(\text{o})}}{\partial z^{(\text{o})}}\dfrac{\partial
z^{(\text{o})}}{\partial p_{j}^{(1)}}\dfrac{\partial p_{j}^{(1)}}{\partial
z_{j}^{(1)}}\dfrac{\partial z_{j}^{(1)}}{\partial\mathcal{I}}\\\
\frac{\partial\mathcal{O}_{1}}{\partial\mathcal{I}}=&\sum_{j}w_{j}^{(\text{o})}w_{j}^{(1)}\mathcal{A}^{\prime}\left(z^{(1)}_{j}\right).\end{split}$
(20)

Figure 5: ANN which takes as input $x$ and returns (a) $\mathcal{O}=x^{2}$ and
(b) its derivative with respect to the input, i.e.,
$\nabla_{\mathcal{I}}\mathcal{O}=2x$, with one hidden layer whose activation
function is $\mathcal{A}$.

Consider the following loss function

$\mathcal{L}=\text{w}_{o}\mathcal{L}_{o}+\text{w}_{\nabla_{\mathcal{I}}\mathcal{O}}\mathcal{L}_{\nabla_{\mathcal{I}}\mathcal{O}},$

where $\mathcal{L}_{o}$ and $\mathcal{L}_{\nabla_{\mathcal{I}}\mathcal{O}}$
are the loss functions corresponding to output $\mathcal{O}$ and
$\nabla_{\mathcal{I}}\mathcal{O}$, respectively. Regularized weights,
$\text{w}_{o}$ and $\text{w}_{\nabla_{\mathcal{I}}}$, can be used to obtain
comparable order of magnitude of the two loss functions. During training,
weights and biases are updated according to Eq. (18) where the computed
gradients are

$\displaystyle\dfrac{\partial\mathcal{L}}{\partial w_{j}^{(\text{o})}}=$
$\displaystyle\>\mathcal{A}\mathcal{L}^{\prime}_{o}+w_{j}^{(1)}\mathcal{A}^{\prime}\mathcal{L}^{\prime}_{\nabla_{\mathcal{I}}\mathcal{O}}$
(21a) $\displaystyle\dfrac{\partial\mathcal{L}}{\partial w_{j}^{(\text{1})}}=$
$\displaystyle\>iw_{j}^{(\text{o})}\mathcal{A}^{\prime}\mathcal{L}^{\prime}_{o}+\left(w_{j}^{(\text{o})}\mathcal{A}^{\prime}+iw_{j}^{(1)}\mathcal{A}^{\prime\prime}\right)\mathcal{L}^{\prime}_{\nabla_{\mathcal{I}}\mathcal{O}}$
(21b) $\displaystyle\dfrac{\partial\mathcal{L}}{\partial b^{(\text{o})}}=$
$\displaystyle\>\mathcal{L}^{\prime}_{o}$ (21c)
$\displaystyle\dfrac{\partial\mathcal{L}}{\partial b_{j}^{(1)}}=$
$\displaystyle\>w_{j}^{(\text{o})}\mathcal{A}^{\prime}\mathcal{L}^{\prime}_{o}+w_{j}^{(\text{o})}w_{j}^{(1)}\mathcal{A}^{\prime\prime}\mathcal{L}^{\prime}_{\nabla_{\mathcal{I}}\mathcal{O}}.$
(21d)

It follows, from relations (21b) and (21d), that the gradient descent
algorithm needs the computation of both first- and second-order gradients of
the activation function $\mathcal{A}$. This particular result is a direct
consequence of the minimization of the error between the gradient of the
outputs with respect to the inputs, i.e.
$\mathcal{O}_{2}=\nabla_{\mathcal{I}}\mathcal{O}_{1}$, and the corresponding
benchmark values, $2x$. This is what we call second-order vanishing gradient
problem. It is tantamount to the first-order variant, but it involves the
second derivatives (and not only the first) of the activation functions in an
ANN. With reference to Figure 4, none of the depicted, classical activation
functions is suitable for such class of problems. Consequently, care must be
taken in selecting activation functions that do not have second-order
vanishing gradients.

####  3.2.1 Understanding second-order vanishing gradient

In the following, we investigate the performance and influence of different
activation functions on the computational time to train an ANN with input
$\mathcal{I}$, primary output $\mathcal{O}_{1}$, and secondary output
$\mathcal{O}_{2}=\nabla_{\mathcal{I}}\mathcal{O}_{1}$. Consider the above
discussed example with $\mathcal{I}=x$, $\mathcal{O}_{1}=x^{2}$, and
$\mathcal{O}_{2}=2x$. The ANN has one hidden layer, with $N_{n}=6$ nodes, and
activation functions as reported in Table 1. The output layer has linear
activation and null bias. The absolute error is selected as loss function for
both $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$. Training is performed on 1000
samples, normalized between -1 and 1. A very small value for the learning rate
is selected, i.e., $\epsilon=10^{-5}$ in order to facilitate the gradient
descent algorithm in reaching small values of the loss function. We use early-
stopping. In other words, training is stopped as the error of a validation set
(500 samples) starts to increase while the learning error still decreases
[19]. The validation set is used to avoid over-fitting of the training data.

Table 1: Set of activation functions considered to investigate the performance
of the network with outputs $\mathcal{O}=x^{2}$ and
$\nabla_{\mathcal{I}}\mathcal{O}=2x$, with $\mathcal{I}=x$, in the framework
of first- and second-order vanishing gradients.

Function $z$ range $\mathcal{A}(z)$ $\mathcal{A}^{\prime}(z)$
$\mathcal{A}^{\prime\prime}(z)$ $\text{ReLU}_{z}$ $z<0$ $0$ $0$ $0$ $z\geq 0$
$z$ $1$ $0$ $\text{ReLU}_{0.5z^{2}+z}$ $z<0$ $0$ $0$ $0$ $z\geq 0$
$0.5z^{2}+z$ $z+1$ $1$ $\text{ReLU}_{z^{2}}$ $z<0$ $0$ $0$ $0$ $z\geq 0$
$z^{2}$ $2z$ $2$ $\text{ELU}_{\text{e}}$ $\forall z$ $\text{e}^{z}-1$
$\text{e}^{z}$ $\text{e}^{z}$ $\text{ELU}_{z}$ $z<0$ $\text{e}^{z}-1$
$\text{e}^{z}$ $\text{e}^{z}$ $z\geq 0$ $z$ $1$ $0$ $\text{ELU}_{0.5z^{2}+z}$
$z<0$ $\text{e}^{z}-1$ $\text{e}^{z}$ $\text{e}^{z}$ $z\geq 0$ $0.5z^{2}+z$
$z+1$ $1$ $\text{ELU}_{z^{2}}$ $z<0$ $\text{e}^{z}-1$ $\text{e}^{z}$
$\text{e}^{z}$ $z\geq 0$ $z^{2}$ $2z$ $2$ $\text{ELU}_{z^{4}}$ $z<0$
$\text{e}^{z}-1$ $\text{e}^{z}$ $\text{e}^{z}$ $z\geq 0$ $z^{4}$ $4z^{3}$
$12z^{2}$ $\text{ELU}_{z^{4}+0.5z^{2}+z}$ $z<0$ $\text{e}^{z}-1$
$\text{e}^{z}$ $\text{e}^{z}$ $z\geq 0$ $z^{4}+0.5z^{2}+z$ $4z^{3}+z+1$
$12z^{2}+1$

For each tested activation function, Table 2 shows the Mean Absolute Error
(MAE) calculated using a set of new, unseen data (500 samples) of input-output
predictions for $x^{2}$ and $2x$. The advancement of training is quantified
herein as the number of epochs, i.e., the number with which the training
algorithm works with the training data-set [19]. Activation functions with
quadratic terms, or of higher degree, perform very well, compared to their
linear equivalents. RELU${}_{z^{2}}$, ELU${}_{z^{2}}$ outperform as their
shape is very similar to the input-output regression they are trained to
learn. Nevertheless, it is worth noticing that training fails when activation
functions with vanishing second gradient are used (e.g. RELUz and ELUz).
Figure 6 compares the ANN predictions for a selection of activation functions
with the analytical (exact) results. Whilst RELUz is clearly inadequate, ELUz
predictions overall agree with the analytical values. This is due to the fact
that the ANN takes advantage of the exponential term, for negative $z$ and
thus successfully manage to satisfy both $\mathcal{O}$ and
$\nabla_{\mathcal{I}}\mathcal{O}$. Additional hidden layers may improve the
performance of the network. It can be further noticed that activation function
of high degree, e.g. $\text{ELU}_{\text{e}}$, $\text{ELU}_{z^{4}}$, and
$\text{ELU}_{z^{4}+0.5z^{2}+z}$, even if successful, require a large number of
epochs.

Table 2: Activation functions and performance with unseen data.

Activation function $\mathcal{A}$ $\mathcal{L}$ $\mathcal{L}_{\mathcal{O}}$
$\mathcal{L}_{\nabla_{\mathcal{I}}\mathcal{O}}$ no. epochs $(10^{-4})$
$(10^{-4})$ $(10^{-4})$ (-) $\text{ReLU}_{z}$ 1521.2 205.98 1315.18 920
$\text{ReLU}_{0.5z^{2}+z}$ 762.4 93.58 668.85 8054 $\text{ReLU}_{z^{2}}$ 0.061
0.0241 0.0371 148 $\text{ELU}_{\text{e}}$ 127.2 26.83 100.38 19477
$\text{ELU}_{z}$ 108.56 12.12 96.44 17280 $\text{ELU}_{0.5z^{2}+z}$ 65.5 10.91
54.63 12178 $\text{ELU}_{z^{2}}$ 0.13 0.067 0.067 88 $\text{ELU}_{z^{4}}$
65.36 33.75 31.61 20051 $\text{ELU}_{z^{4}+0.5z^{2}+z}$ 12.94 1.81 11.13 9683

(a) $x^{2}$ predictions, $\mathcal{O}_{1}$, using ReLU and $\text{ELU}_{z}$.

(b) $2x$ predictions, $\mathcal{O}_{2}$, using ReLU and $\text{ELU}_{z}$.

(c) $x^{2}$ predictions, $\mathcal{O}_{1}$, using $\text{ELU}_{z}$,
$\text{ELU}_{0.5z^{2}+z}$, $\text{ELU}_{z^{2}}$, $\text{ELU}_{\text{e}}$, and
$\text{ELU}_{z^{4}}$.

(d) $2x$ predictions, $\mathcal{O}_{2}$, using $\text{ELU}_{z}$,
$\text{ELU}_{0.5z^{2}+z}$, $\text{ELU}_{z^{2}}$, $\text{ELU}_{\text{e}}$, and
$\text{ELU}_{z^{4}}$.

Figure 6: Comparison of different activation functions for the prediction of
the primary output, $x^{2}$ (a), and secondary output, $2x$ (b). From top to
bottom the range of $z$ decreases from larger to smaller values, to observe
the behavior at $z\approx 0$.

###  3.3 Architecture of Thermodynamics-based ANN

Herein we detail the architecture and the internal steps/definitions TANN is
relying on. The architecture is detailed in Figure 7. The input vector is
$\mathcal{I}=(\varepsilon^{t},\Delta\varepsilon,\sigma^{t},\theta^{t},\zeta_{i}^{t},\Delta
t)$, the primary and secondary outputs are
$\mathcal{O}=(\Delta\zeta_{i},\Delta\theta,\textsf{F}^{t+\Delta t})$ and
$\nabla_{\mathcal{I}}\mathcal{O}=(\Delta\sigma,\textsf{D}^{t+\Delta t})$,
respectively. TANN involves the following steps:

  * 1.

computation of the updated strain (definition): $\varepsilon^{t+\Delta
t}:=\varepsilon^{t}+\Delta\varepsilon$

  * 2.

prediction of the kinematic variables and temperature increments with two sub-
ANNs:

$\Delta\zeta=\textsf{sNN}_{\zeta}@\left(\varepsilon^{t+\Delta
t},\Delta\varepsilon^{t},\sigma^{t},\theta^{t},\zeta^{t}\right)$

and

$\Delta\theta=\textsf{sNN}_{\theta}@\left(\varepsilon^{t+\Delta
t},\Delta\varepsilon,\sigma^{t},\theta^{t},\zeta^{t}\right)$

  * 3.

computation of

    * $(a)$

the updated kinematic variables rates (backward finite difference
approximation):

$\dot{\zeta}^{t+1}\approx\frac{\Delta\zeta}{\Delta t}$

    * $(b)$

the updated kinematic variables (definition):
$\zeta^{t+1}:=\zeta^{t}+\Delta\zeta^{t}$

    * $(c)$

the updated temperature (definition): $\theta^{t+1}:=\theta^{t}+\Delta\theta$

  * 4.

prediction of the updated energy potential:

$\text{{F}}^{t+\Delta t}=\textsf{sNN}_{\textsf{F}}@\\{\varepsilon^{t+\Delta
t}\quad\zeta^{t+\Delta t}\quad\theta^{t+\Delta t}\\}$

  * 5.

computation of the updated dissipation (definition, Eq. (13)):
$\text{{D}}^{t+\Delta t}:=-\frac{\partial\text{{F}}^{t+\Delta
t}}{\partial\zeta^{t+\Delta t}}\cdot\dot{\zeta}^{t+\Delta t}$

  * 6.

computation of

    * $(a)$

the updated stress (definition, Eq. (13)): $\sigma^{t+\Delta
t}:=\frac{\partial\text{{F}}^{t+\Delta t}}{\partial\varepsilon^{t+\Delta t}}$

    * $(b)$

the stress increment (definition): $\Delta\sigma:=\sigma^{t+\Delta
t}-\sigma^{t}$

(a) non-isothermal processes.

(b) isothermal processes.

Figure 7: Architecture of TANN: general case (a) and for isothermal processes
(b). Inputs are highlighted in gray ( ); outputs in black, ( ) for direct ANN
predictions and ( ) for derived outputs; and intermediate quantities
(definitions) are in white ( ) and ( ). Relationships obtained from
definitions are represented with dashed lines, while arrows denote ANNs.

TANN is thus composed of three sub-ANNs; $\textsf{sNN}_{\zeta}$ predicts the
internal variables increment, $\textsf{sNN}_{\theta}$ predicts the temperature
increment (note that in case of the isothermal conditions, this component can
be removed from the architecture, see Fig. 7(b)), and
$\textsf{sNN}_{\textsf{F}}$ predicts the Helmholtz free-energy. The main
output, the increment in stress, is computed according to expression (13),
which stems from thermodynamic requirements. By virtue of the fact that the
entire constitutive response of a material can be derived from definition of
only two potential functions, the model is able to predict the stress
increment from the knowledge of the energy potential (and the internal
variables $\zeta_{i}$). It is worth noticing that, differently from common
approaches (cf. Sect. 3), the sub-network $\textsf{sNN}_{\textsf{F}}$ is
required to learn a scalar quantity$-$that is, the Helmholtz free-energy
potential. This offer compelling advantages. When dealing with ANNs, the curse
of dimensions (increasing effort in training and large amount of training data
required) is an important issue when the studied problem passes to higher
dimensions, see e.g. [5]. Passing from 1D to 3D, for instance, increases the
number of variables the ANN need to learn. For stresses, from one single
scalar value, in 1D, we pass to a vector with six-components, in 3D. The
computational effort is thus not trivial. Nevertheless, TANN is, in principle,
less affected by these issues as the two potentials, on which the entire set
of predictions relies on, are scalar functions.
The computation of dissipation, from expression (13), plays a double role.
First, it assures thermodynamic consistency of the predictions of TANN (first
law). Second, it brings the information to distinguish between reversible and
irreversible processes, e.g. elasticity from plasticity/damage, etc., and it
is trained to be positive or zero (second law).
It is worth noticing that further improvements of the performance of TANNs may
be obtained, as suggested in the work of Karpatne et al. [31], by adding a
physical inconsistency term to the loss functions (e.g., with respect to
dissipation).

##  4 Generation of data

We present the procedure used to generate material data TANN is trained with.
Herein, data are obtained by numerical integration of an incremental form of
the constitutive relations. To this purpose, we assume the Ziegler’s
orthogonality condition (see paragraph 4.1 and [56, 24, 25]), which, in
general, it is not a strict requirement. Nevertheless, it is worth noticing
that this restriction applies only on the generated data, and not on the ANN
class here proposed. More precisely, TANN architecture still holds even for
materials for which the Ziegler’s normality condition does not apply. We shall
recall that the aim is to demonstrate the advantages of thermodynamics-based
neural networks with respect to classical approaches. Hence the restrictions,
imposed by the orthogonality hypothesis for the generation of data, are
expected not to affect the comparisons presented in Section 5.

###  4.1 Incremental formulation

Following the hyperplasticity framework proposed in [14], the thermo-
mechanical, non-linear, incremental constitutive relation for strain-rate
independent materials, undergoing infinitesimal strains, is here derived in
the framework of isothermal processes ($\theta=\text{cost}$). By
differentiating the energy expressions (13) and rearranging the terms, we
obtain the following non-linear incremental relations

$\displaystyle\dot{\sigma}=\partial_{\varepsilon\varepsilon}\text{{F}}\cdot\varepsilon+\sum_{k}\partial_{\varepsilon\zeta_{k}}\text{{F}}\cdot\dot{\zeta}_{k},$
(22a)
$\displaystyle-\dot{\chi}_{i}=\partial_{\zeta_{i}\varepsilon}\text{{F}}\cdot\varepsilon+\sum_{k}\partial_{\zeta_{i}\zeta_{k}}\text{{F}}\cdot\dot{\zeta}_{k}.$
(22b)

where the following notation is adopted

$\partial_{\varepsilon\varepsilon}\text{{F}}=\dfrac{\partial^{2}\text{{F}}}{\partial\varepsilon_{ij}\partial\varepsilon_{kl}},\quad\partial_{\varepsilon\zeta_{k}}\text{{F}}=\dfrac{\partial^{2}\text{{F}}}{\partial\varepsilon_{ij}\partial\zeta_{k}},\quad\partial_{\zeta_{i}\zeta_{k}}\text{{F}}=\dfrac{\partial^{2}\text{{F}}}{\partial\zeta_{i}\partial\zeta_{k}}.$

Further, introducing the thermodynamic dissipative stresses
$\mathcal{X}^{\dagger}=(X_{1},\ldots,X_{N})$ and assuming the Ziegler’s
orthogonality condition ([56]), the following non-linear, incremental
constitutive relation can be found

$\dot{\Xi}=\begin{cases}\mathcal{M}|_{y=0}\;\cdot\dot{\varepsilon}&\text{if
}y=0\\\ \mathcal{M}|_{y<0}\;\cdot\dot{\varepsilon}&\text{else }\end{cases}$
(23)

with

$\dot{\Xi}=\begin{bmatrix}\dot{\sigma}\\\ -\dot{X}_{i}\\\ \dot{\zeta_{i}}\\\
\lambda\end{bmatrix},\quad\mathcal{M}|_{y=0}=\begin{bmatrix}\partial_{\varepsilon\varepsilon}\text{{F}}-\sum_{k}\partial_{\varepsilon\zeta_{k}}\text{{F}}\cdot\left(\frac{\mathcal{C}_{\varepsilon}}{B}\cdot\frac{\partial
y}{\partial X_{k}}\right)\\\
\partial_{\zeta_{i}\varepsilon}\text{{F}}-\sum_{k}\partial_{\zeta_{i}\zeta_{k}}\text{{F}}\cdot\left(\frac{\mathcal{C}_{\varepsilon}}{B}\cdot\frac{\partial
y}{\partial X_{k}}\right)\\\
-\dfrac{\mathcal{C}_{\varepsilon}}{B}\cdot\frac{\partial y}{\partial X_{i}}\\\
-\dfrac{\mathcal{C}_{\varepsilon}}{B}\end{bmatrix},\quad\text{and}\quad\mathcal{M}|_{y<0}=\begin{bmatrix}\partial_{\varepsilon\varepsilon}\text{{F}}\\\
\partial_{\zeta_{i}\varepsilon}\text{{F}}\\\ \emptyset\\\
\emptyset\end{bmatrix},$ (24)

and $\cdot$ denotes the contraction of adjacent indices. In the above
relations (23-24), whose derivation is presented in Appendix A,
$y=\tilde{y}(\varepsilon,\mathcal{Z},\mathcal{X}^{\dagger})$ is the yield
function, $\emptyset$ denotes a quantity (scalar or tensorial, depending on
the dimensionality of the internal variable set) equal to zero,

$\mathcal{C}_{\varepsilon}=\frac{\partial
y}{\partial\varepsilon}-\sum_{i=1}^{N}\frac{\partial y}{\partial
X_{i}}\cdot\partial_{\zeta_{i}\varepsilon}\text{{F}},$
$B=\sum_{i=1}^{N}\frac{\partial y}{\partial\zeta_{i}}\cdot\frac{\partial
y}{\partial X_{i}}-\sum_{i=1}^{N}\frac{\partial y}{\partial
X_{i}}\left(\sum_{k=1}^{N}\partial_{\zeta_{k}\varepsilon}\text{{F}}\cdot\frac{\partial
y}{\partial X_{k}}\right).$

###  4.2 Data generation

Data are generated in a Python environment [44], where SymPy [45] and SciPy
[53] libraries are used for symbolic calculations and numerical integration.
Data are generated by identifying an initial state for the material at time
$t$,

$\text{state at time }t:\qquad\Xi^{t}=\begin{bmatrix}\sigma^{t}\\\
-X_{i}^{t}\\\ \zeta_{i}^{t}\\\
0\end{bmatrix}\quad\text{and}\quad\varepsilon^{t},$

and a given strain increment $\dot{\varepsilon}^{t}$, assuming constant and
unitary time increment $\Delta t=1$
($\dot{\varepsilon}^{t}=\Delta\varepsilon^{t}$). Numerical integration of the
ordinary differential equations (23) is performed with an explicit solver [6]
to obtain the state at the new time $t+\Delta t$, i.e.,

$\text{state at time }t+\Delta t:\qquad\Xi^{t+\Delta
t}=\begin{bmatrix}\sigma^{t+\Delta t}\\\ -X_{i}^{t+\Delta t}\\\
\zeta_{i}^{t+\Delta t}\\\ \lambda^{t+\Delta t}\end{bmatrix}$

The training data play a crucial role for both the accuracy of the predictions
and the generalization with respect to the ANN state variables, e.g., strain
increments. The generalization capability of a network is here defined as the
ability to make predictions for loading paths different from those used in the
training operation. Nevertheless, a significant dependency on the ANN state
variables is usually observed. This may result in a poor network
generalization. In Lefik and Schrefler [36], an improvement of the
generalization capability of the ANN is proposed. Artificial sub-sets of data,
with zero strain increments, are added in the set of training data to force
the network in learning that to zero input increments correspond zero output
increments.
In the available literature, strain-stress loading paths are commonly used in
training. If recursive neural networks are used, feeding them with history
variables (loading paths) is the only possible solution (see e.g. [47]).
Nevertheless, ANNs do not necessary need the data-sets to be (historical)
paths.

Herein, we generate data randomly, and not following prescribed loading paths.
Conversely, this allows us to (1) improve the representativeness of the
material data and (2) improve the generalization of the network on the strain
increments. The initial state, $\Xi^{t}$ and $\varepsilon^{t}$, and the strain
increment, $\Delta\varepsilon^{t}$, are randomly generated from standard
distributions with mean value equal to zero and standard deviation equal to
$\Xi_{\text{max}}$, $\varepsilon^{t}_{\text{max}}$, and
$\Delta\varepsilon^{t}_{\text{max}}$, respectively. The Cauchy and
thermodynamic stresses, $\sigma^{t}$ and $X_{i}^{t}$, as well as the internal
variables $\zeta_{i}^{t}$ are then calculated to satisfy the constraint
$y^{t}\leq 0$. This incremental procedure is repeated for
$N_{\text{samples}}$, resulting in a set of $N_{\text{samples}}$ ordered pairs
$\\{\Xi^{t},\varepsilon^{t},\Delta\varepsilon^{t};\Xi^{t+\Delta t}\\}$, from
which the corresponding energy potential and dissipation rate at time
$t+\Delta t$ are evaluated.
The choice of the standard deviations $\Xi^{t}_{\text{max}}$,
$\varepsilon^{t}_{\text{max}}$, and $\Delta\varepsilon^{t}_{\text{max}}$
depends on the investigated problem. As it follows, they are selected in a way
such that $50\div 60\%$ of data samples, $N_{\text{samples}}$, lie on the
yield surface, i.e., $y^{t}=0$. Figure 8 depicts the sampling for one of the
investigated applications (see paragraph 5.1).

(a) Cauchy stress, $\sigma^{t}$, $\sigma^{t+\Delta t}$, and increment
$\Delta\sigma$.

(b) thermodynamic stress, $X^{t}$, $X^{t+\Delta t}$, and increment $\Delta X$.

(c) total strain, $\varepsilon^{t}$, $\varepsilon^{t+\Delta t}$, and increment
$\Delta\varepsilon$.

(d) internal variable (inelastic strain), $\zeta^{t}$, $\zeta^{t+\Delta t}$,
and increment $\Delta\zeta$.

(e) Helmholtz free-energy $F^{t+\Delta t}$.

(f) dissipation rate $D^{t+\Delta t}$.

Figure 8: Sampling for one of the studied applications (paragraph 5.1). 1D
elasto-plastic material with Young’s modulus $E=200$ GPa and yield strength
$\sigma_{y}=200$ MPa. Standard deviations
$\sigma^{t}_{\text{max}}=X^{t}_{\text{max}}=2\sigma_{y}$,
$\varepsilon^{t}_{\text{max}}=\zeta^{t}_{\text{max}}=10^{-2}$, and
$\Delta\varepsilon^{t}_{\text{max}}=10^{-3}$.

##  5 Applications

Herein we use TANNs to the modeling of multi-dimensional elasto-plastic
materials and demonstrate the wide applicability and effectiveness of TANNs.
It is worth noticing that, even though the applications here investigated
consist of elasto-plastic materials, the proposed class of ANN can be
successfully applied (without any modification) to materials with different or
more complex behavior, accounting e.g. for damage and/or other non-linearities
(in the framework of strain-rate independent processes). In paragraph 5.1, a
one-dimensional motivating example, a spring-slider model, is studied. Then,
the more general cases of three-dimensional elasto-plasticity (paragraph 5.2),
accounting for perfect-plasticity, hardening and softening behaviors, are
presented.
As it follows, the hyper-parameters (i.e., number of hidden layers, neurons,
activation functions, etc.) of the networks are selected to give the best
predictions, while requiring minimum number of hidden layers and nodes per
layer. This is accomplished by comparing the learning error on the set of test
patterns, per each trial choice of the hyper-parameters. In each training
process, we use the commonly used technique of early-stopping (see par.
3.2.1): the iterative update of weights and biases is stopped as the test
error starts to increase while the learning error still decreases (indicator
of an over-fitting of the training data).
Throughout this Section relatively simple deep feed-forward neural networks
architectures are used (with, at maximum, two hidden layers) and no additional
regularization techniques are employed (e.g., L1/L2 penalties, dropout, etc.).
Each numerical example is accompanied with a detailed discussion about the
network architecture.

###  5.1 1D elasto-plasticity with kinematic hardening

The Helmholtz free-energy potential, dissipation rate, and yield function that
define the elasto-plastic 1D model with kinematic hardening (1D spring-slider,
[25]) are:

$\begin{split}\textsf{F}=&\frac{E}{2}\left(\varepsilon-\zeta\right)^{2}+\frac{H}{2}\zeta^{2},\quad\textsf{D}=k|\dot{\zeta}|,\\\
\text{and}\quad y=&\frac{|\sigma-H\zeta|}{k}-1\leq 0,\end{split}$ (25)

with $H$ being the kinematic hardening/softening-parameter and $k$ being the
yield strength (slider threshold). The internal variable, $\zeta$, represents
herein the plastic deformation. Table 3 displays the choice of the material
parameters, selected to represent a steel-like material with either (1)
perfect-plastic, (2) hardening, or (3) softening behavior.

Table 3: Material parameters for 1D elasto-plastic materials.

case $E$ $k$ $H$ (GPa) (MPa) (GPa) 1D-1 200 200 0 1D-2 200 200 10 1D-3 200 200
-10

####  5.1.1 Training

According to the procedure detailed in Section 4, 2000 data (random increments
at random states) are generated, for each (material) case, with the procedure
detailed in Section 4. Training is performed with 50 % of them (i.e., 1000). A
validation set of 500 samples (validation data) is used to avoid over-fitting
together with early stopping rules [7]. The performance of the predictions, at
the end of the training, is evaluated on a set of 500 samples (test data). The
sampling for material case 1D-1 is shown in Figure 8, while the samples
distribution for cases 1D-2 and 1D-3 are presented in SM. Adam optimizer with
Nesterov’s acceleration gradient [10] is selected and a batch size of 10
samples is used. We use the Mean Absolute Error (MAE) as loss functions for
each output in order to assure the same precision between data of low and high
numerical values (cf. Mean Square Error). Regularized weights are used to have
consistent order of magnitude of different quantities involved in the loss
functions. The architecture of TANN for all 1D cases consists of one hidden
layer with 6 neurons (and leaky ReLU activation function) for the predictions
of $\Delta\zeta$ and one hidden layer with 9 neurons (activation
$\text{ELU}_{z^{2}}$) to predict $\textsf{F}^{t+\Delta t}$. The output layers
for both sub-networks have linear activation function and biases set to zero.
The corresponding number of degrees of freedom, i.e., the number of the hyper-
parameters, is 72. Higher number of hidden layers could be used as well, but
this is out of the scope of our investigations. Figure 9 displays the loss
functions of each output as the training is performed, i.e., in number of
epochs. The early stopping rule assures convergence with MAEs of the same
order of magnitude for the 4 outputs, $\Delta\zeta$, $\textsf{F}^{t+\Delta
t}$, $\Delta\sigma$, and $\textsf{D}^{t+\Delta t}$. Similar behaviors in
training and performance with respect to the set of test data are also found
for cases 1D-2 and 1D-3. For this reason, they are not presented herein.

Figure 9: Errors of the predictions of TANN (loss functions), as the training
is being performed, evaluated with respect to the training (train) and
validation (val) sets. Weights and biases update are computed only on the
training set.

####  5.1.2 Predictions in recall mode

Once the neural networks have been trained, we use them in recall mode to
predict the stress increment for a given strain, strain increment, and
possibly other variables, and we compare the predictions with the
corresponding targets. The results of the numerical integration scheme
presented in Section 4 are here considered as the exact solution of the
material response. In particular, starting from an initial configuration, we
make cyclic (or random) increments of the strain, $\Delta\varepsilon$. TANN
hence predicts the corresponding increments, $\\{\Delta\zeta,\Delta\sigma\\}$,
which will be transformed into the inputs in the successive call, as well as
the energy and dissipation rates, $\\{\textsf{F}^{t+\Delta
t},\textsf{D}^{t+\Delta t}\\}$. This procedure is applied recursively. The
neural network is so self-fed. Figures 10(a), 10(b), and 10(c)
illustrate$-$for cases 1D-1, 1D-2, and 1D-3, respectively$-$the predictions of
TANN for cyclic paths with strain increments
$\Delta\varepsilon^{n}=\Delta\varepsilon\text{
sgn}\left(\cos\frac{n\pi}{2N}\right)$, where $n=1,2,\ldots$,
$N=\varepsilon_{\text{max}}/\Delta{\varepsilon}$,
$\varepsilon_{\text{max}}=2\times 10^{-3}$, and $\Delta{\varepsilon}=10^{-4}$.

(a) material case 1D-1.

(b) material case 1D-2.

(c) material case 1D-3.

Figure 10: Predictions of TANN due to cyclic loading, compared with the target
constitutive model for cases (a) 1D-1, (b) 1D-2, and (c) 1D-3 (see Tab. 3),
with strain increments $\Delta{\varepsilon}=10^{-4}$.

TANN is found to successfully predict all quantities of interest. Moreover,
and most important, the architecture and the training of the network allows to
obtain thermodynamically consistent results. The first law of thermodynamics
is automatically satisfied as a result of the structure of TANN and the
predicted dissipation rate is always positive. Indeed, even if the second
principle of thermodynamics is not explicitly assured by the TANN
architecture, the fact that the training has been performed with consistent
material data (i.e., positive dissipation rate) results automatically in the
fulfillment of the second principle. Moreover, the linear dependency of
dissipation with respect to $\dot{\zeta}$, property that stems from the strain
rate independent material formulation, is also automatically recovered by
TANN.

####  5.1.3 Generalization of the network

Herein we investigate the generalization capability of TANN (i.e., the ability
to make predictions for loading paths different from those used in the
training operation). This is achieved by feeding the trained network with
input values that not necessarily belong to the training range. With reference
to Figures 8 and S1 and S2 (SM file), the training ranges of the inputs are
represented in Table 4. Figure 11 displays the predictions for a cycling
loading path $\Delta\varepsilon^{n}=\Delta\varepsilon\text{
sgn}\left(\cos\frac{n\pi}{2N}\right)-$with $\Delta\varepsilon\in(10^{-5},1)$.
We clearly see that for input variables outside the training range, the
predictions of the network become less accurate. Nevertheless, the predictions
are always thermodynamically consistent. Moreover, the quantities of primary
interest, such as the stress, the internal state variable, and the energy are
in extremely good agreement with the reference model. The same stands also for
the dissipation rate. We notice, once more, that its values are always
positive, even when the network is used for predictions beyond the training
range. In the SM file we present the predictions for the plastic hardening and
softening behaviors, respectively.

Table 4: Range of the value of inputs used for training, for material cases
1D-1, 1D-2, and 1D-3.

case $\max|\varepsilon^{t}|$, $\min|\varepsilon^{t}|$
$\max|\Delta\varepsilon|$, $\min|\Delta\varepsilon|$ $\max|\zeta^{t}|$,
$\min|\zeta^{t}|$ $\max|\sigma^{t}|$, $\min|\sigma^{t}|$ (-) (-) (-) (MPa)
1D-1 $4\times 10^{-2}$, $1\times 10^{-7}$ $1.4\times 10^{-2}$, $1\times
10^{-7}$ $4\times 10^{-2}$, $6\times 10^{-6}$ $200$, $0.168$ 1D-2 $4\times
10^{-2}$, $4.\times 10^{-6}$ $1.2\times 10^{-2}$, $1\times 10^{-7}$ $3.6\times
10^{-2}$, $6\times 10^{-7}$ $495$, $0.300$ 1D-3 $4\times 10^{-2}$, $1\times
10^{-7}$ $1.4\times 10^{-2}$, $1\times 10^{-7}$ $3.5\times 10^{-2}$, $6\times
10^{-6}$ $505$, $0.060$

It is worth noticing that in all the cases, even for very large strain
increments$-$for which the predictions of the network in terms of dissipation
rate, energy potential, and (for some values) stress and internal variable
differ from the target values$-$, TANN successfully predicts the Jacobian,
i.e., $\frac{\partial\sigma}{\partial\epsilon}$, in very good agreement with
the reference model. This is true even when the error in the stress prediction
is not negligible. This is of particular importance for numerical simulations
with implicit algorithms. Therefore, TANN can successfully replace complicated
constitutive models or multiscale approaches, but considerably and safely
decreasing the calculation cost, even when the requested increments are
outside the training range.

(a) $\sigma$ prediction.

(b) $\zeta$ prediction.

(c) F prediction.

(d) D prediction.

Figure 11: Sensitivity on the inputs for material case 1D-1, with strain
increments $\Delta\varepsilon=\Delta{\varepsilon}\text{
sgn}\left(\cos\frac{n\pi}{2N}\right)-$with
$N=\varepsilon_{\text{max}}/\Delta{\varepsilon}$,
$\varepsilon_{\text{max}}=2\times 10^{-3}\div 2$, and $\Delta{\varepsilon}$
varying from $\times 10^{-5}$ (top) to $1$ (bottom). Each column displays the
response (from left to right) in term of $\Delta\sigma$ (a), $\Delta\zeta$
(b), $\textsf{F}^{t+\Delta t}$ (c), and $\textsf{D}^{t+\Delta t}$ (d). Each
row represents the prediction at different $\Delta\varepsilon$.

####  5.1.4 TANN vs ANN

We compare herein the performance, in recall mode, of TANN with respect to the
classical approach of ANN for constitutive modeling [20, 36]. Figure 12
displays the architecture of the network, ANN, with inputs
$\mathcal{I}=(\varepsilon^{t},\Delta\varepsilon,\sigma^{t},\zeta^{t})$ and
output $\mathcal{O}=(\Delta\zeta,\Delta\sigma)$. As for TANN, the stress
increment is derived by assuming ANN state variables (see Sect. 3) such that
they coincide with the thermodynamic state variables, i.e., $\varepsilon^{t}$
and $\zeta^{t}$, and their increments, i.e., $\Delta\varepsilon$ and
$\Delta\zeta$, as in [35, 55]. ANN is thus composed of two sub-ANNs;
$\textsf{aNN}_{\zeta}$ predicts the internal variables increment and
$\textsf{sNN}_{\sigma}$ predicts the stress increment, i.e.,
$\Delta\sigma=\textsf{aNN}_{\sigma}(\varepsilon^{t+\Delta},\Delta\varepsilon,\zeta^{t+\Delta
t},\Delta\zeta)$. The architecture of the network is selected to give the best
performance while assuring the same amount of degrees of freedom, hyper-
parameters, of TANN. Both sub-networks, $\textsf{aNN}_{\zeta}$ and
$\textsf{aNN}_{\sigma}$, consist of one hidden layer, with 6 neurons each and
leaky ReLU activation function. Same with TANN, the output layers have linear
activation function and zero bias. Training is performed on the same set of
samples that are used for the thermodynamics-based network. Figure 13 displays
the error of the predictions of ANN, as training is performed, and compares it
with TANN.

(a) ANN scheme.

(b) ANN architecture.

Figure 12: Schematic (a) and full architecture (b) of the network, not based
on thermodynamics, ANN. Inputs are highlighted in gray ( ), outputs in black (
).

We present in Figure 14 the comparisons between the predictions of TANN and
ANN, in terms of stress, for a cycling loading path
$\Delta\varepsilon^{n}=\Delta\varepsilon\text{
sgn}\left(\cos\frac{n\pi}{2N}\right)$. The cases with strain hardening and
strain softening are presented in the SM file. TANN is clearly superior in
terms of (a) accuracy of the prediction and (b) generalization with respect to
the inputs. Moreover, ANN predictions do not fulfill the principles of
thermodynamics, even though the training of the network has been performed on
consistent material data. This is clearly shown by computing from the
predictions of ANN the increment of the Helmholtz free-energy and dissipation
rate using the corresponding definitions, Eq. (25). Figures 14(c) and 14(d)
display the computed quantities, F and D, for material case 1D-1. The
predictions of the standard ANN clearly do not respect the thermodynamics
principles (both the first and second laws).

(a) mean absolute error of $\Delta\zeta$ prediction

(b) mean absolute error of $\Delta\sigma_{i}$ prediction

Figure 13: Training of ANN compared with TANN evaluated with respect to the
training (train) and validation (val) sets.

(a) $\sigma$ prediction.

(b) $\zeta$ prediction.

(c) F prediction (TANN), Eq. (25) (ANN).

(d) D prediction (TANN), Eq. (25) (ANN).

Figure 14: Comparison of the predictions of TANN and those of standard ANN,
for material case 1D-1 (perfect plasticity). Each row represents the
prediction at different $\Delta\varepsilon$.

(a) material case 1D-1 (perfect plasticity).

(b) material case 1D-2 (hardening).

(c) material case 1D-3 (softening).

Figure 15: Comparison of the predictions of TANN and standard ANN for a random
loading path, for material cases 1D-1, perfect plasticity (a), 1D-2, hardening
(b), and 1D-3, softening (c). The loading path is displayed at the first row,
in each case.

Figure 15 displays the predictions of both TANN and ANN for a random loading
path, for a perfectly plastic behavior (Fig. 15(a)), hardening (Fig. 15(b)),
and softening (Fig. 15(c)). Once more, the performance of the thermodynamics-
based network, as well as its generalization capabilities, are significantly
better than those of standard ANN. One could, of course, increase the number
of layers and neurons in order to assure better predictions, but still there
will not be guarantee that the predictions of standard ANN would be
thermodynamically consistent. An increase of the number of samples used in the
training operation may as well improve the predictions. This is not the case
for TANN.

###  5.2 3D elasto-plasticity

In order to illustrate the performance of TANN in three dimensions we use the
simple von Mises elasto-plastic model with kinematic hardening (and
softening). The model can be derived from the following expressions of the
energy potential and dissipation rate

$\begin{split}\textsf{F}=&\frac{9K}{2}\left(\varepsilon_{p}-\zeta_{p}\right)\cdot\left(\varepsilon_{p}-\zeta_{p}\right)+\\\
&+G\left(e-z\right)\cdot\left(e-z\right)+\frac{H}{2}z\cdot z,\\\
\textsf{D}=&k\sqrt{2}\sqrt{\dot{z}\cdot\dot{z}},\end{split}$ (26)

where $k$ represents the elastic limit in simple shear; $K$ and $G$ are the
bulk and shear moduli; $\varepsilon_{p}$ and $\zeta_{p}$ are, respectively,
the mean total and plastic deformation; and $e$ and $z$ are, respectively, the
total and plastic deviatoric strain tensors. The yield surface can be derived
as shown in Appendix A [25] and is defined as

$y=D-X^{\prime}\cdot z=\sqrt{X^{\prime}\cdot X^{\prime}}-\sqrt{2}k\leq 0,$
(27)

with $X^{\prime}_{ij}=2G\left(e_{ij}-z_{ij}\right)+Hz_{ij}$.

Table 5: Material parameters for 3D elasto-plastic von Mises material.

case $K$ $G$ $k$ $H$ (GPa) (GPa) (MPa) (GPa) 3D-1 167 77 140 0 3D-2 167 77 140
10 3D-3 167 77 140 -10

####  5.2.1 Training

Data are generated as detailed in Section 4. A total of 6000 data with random
increments of deformation are generated. In order to improve the performance
of the network in recall mode, additional sampling with random uni-axial and
bi-axial loading paths are also used. The samples are split into training
(50%), validation (25%), and test (25%) sets. The sampling in terms of the
mean and deviatoric stresses, $p$ and $q$, and deformations, $\varepsilon_{p}$
and $e$, is presented in the SM file. For the sake of simplicity, stress and
deformation are converted in the principal axes frame of reference.
The network architecture is adapted to the size of the inputs and outputs,
with respect to the mono-dimensional case. In particular, the sub-network
$\textsf{sNN}_{\zeta}$ consists of two hidden layers, with 48 neurons (leaky
ReLU activation function), and three output layers, one per each (principal)
component of (increment of) $\zeta$. The sub-network
$\textsf{s-NN}_{\textsf{F}}$ has one hidden layer with 36 neurons (activation
$\text{ELU}_{z^{2}}$). The output layers for both sub-networks have linear
activation functions and biases set to zero. The resulting number of hyper-
parameters is $\approx 3000$. The loss functions of each output as the
training is performed, for material case 3D-1 (perfect plasticity), is
presented in the SM file. Similar behaviors are also recovered for cases 3D-2
(hardening), 3D-3 (softening).

####  5.2.2 Predictions in recall mode

Once the network has been trained, it is used, in recall mode, to make
predictions. We briefly present the performance of TANN in predicting the
material response for a cyclic loading path. Figure 16 depicts the comparison
with the target material model for material case 3D-1. The predictions for the
same loading path with material cases 3D-2 and 3D-3 are presented in the SM
file. In all the cases, the network shows good performance.
Similarly with the 1D case, the generalization capabilities of the network are
presented together with the comparison of the thermodynamics-based network
with the standard ANN approach.

(a) loading path.

(b) $\sigma_{i}$ prediction.

(c) $\zeta_{i}$ prediction.

(d) energy (left) and dissipation rate (right) predictions.

Figure 16: Predictions of TANN due to cyclic loading, compared with the target
constitutive model, case 3D-1, perfect plasticity: loading path (16(a)), in
terms of principal deformations; principal stress predictions
($\sigma_{1},\sigma_{2},\sigma_{3}$) (16(b)); plastic deformation predictions
($\zeta_{1},\zeta_{2},\zeta_{3}$) (16(c)); energy (F) and dissipation rate
predictions (D) (16(d)).

####  5.2.3 TANN vs standard ANN. Generalization of the network

Herein we investigate the performance of TANN with respect to the classical
approach of ANN, as well as the sensitivity with respect to the input
variables range. Similarly to the comparisons in paragraph 5.1.4, we select a
network with inputs
$\mathcal{I}=(\varepsilon^{t}_{i},\Delta\varepsilon_{i},\sigma_{i}^{t},\zeta_{i}^{t})$
and output $\mathcal{O}=(\Delta\zeta_{i},\Delta\sigma_{i})$, with $i=1,2,3$
denoting the principal components. The architecture is selected to give the
best performance, preserving the same number of hyper-parameters between TANN
and standard ANN. The network, ANN, consists of the two sub-networks,
$\textsf{aNN}_{\zeta}$ and $\textsf{aNN}_{\sigma}$, with two hidden layers,
each one, leaky ReLU activation functions, and number of neurons per layer
equal to 48. As for $\textsf{sNN}_{\zeta}$ and $\textsf{sNN}_{\sigma}$, in
$\textsf{aNN}_{\zeta}$ and $\textsf{aNN}_{\sigma}$ three output layers (1
neuron each) are used, with linear activation functions and zero biases. In
the SM file we present the comparison of the MAE of the network predictions
with respect to the target values (training and validation data-sets).

We first compare the performance of both networks, TANN and standard ANN, in
predicting the material response for cyclic isotropic loading paths (material
case 3D-1, cf. Table 5). A linear elastic material response is expected and
retrieved. Figure 17 displays the stress predictions of TANN and ANN, compared
with the target values, for different strain increments. It is worth
mentioning that the standard approach of ANN does not succeed in accurately
predicting the elastic deformation range. Moreover, contrary to TANN, the
stress predictions of standard ANN, depend strongly on the cyclic loading. As
the network is used recursively, in recall mode, the stress predictions
rapidly become less and less precise, due to error accumulation.

(a) strain increment $\Delta{\varepsilon}=1\times 10^{-3}$.

(b) strain increment $\Delta{\varepsilon}=1\times 10^{-2}$.

(c) strain increment $\Delta{\varepsilon}=1\times 10^{-1}$.

Figure 17: Comparison of the stress predictions of TANN and standard ANN with
respect to the target values, for the cyclic, isotropic loading path
$\Delta\varepsilon_{1}^{n}=\Delta\varepsilon_{2}^{n}=\Delta\varepsilon_{3}^{n}=\Delta{\varepsilon}\text{
sgn}\left(\cos\frac{n\pi}{2N}\right)-$with
$N=\varepsilon_{\text{max}}/\Delta{\varepsilon}$,
$\varepsilon_{\text{max}}=2\times 10^{-3}$ (a),
$\varepsilon_{\text{max}}=\times 10^{-1}$ (b), and
$\varepsilon_{\text{max}}=1$ (c), for material case 3D-1 (perfect plasticity).
Each row represents the prediction at different $\Delta\varepsilon$
increments.

The performance of both networks is further compared for the following three
loading paths

$\begin{split}\text{uni-
axial:}\qquad&\Delta\varepsilon_{1}^{n}=\Delta{\varepsilon}\text{
sgn}\left(\cos\frac{n\pi}{2N}\right),\quad\Delta\varepsilon_{2}^{n}=\Delta\varepsilon_{3}^{n}=0;\\\
\text{bi-axial:}\qquad&\Delta\varepsilon_{1}^{n}=\Delta{\varepsilon}\text{
sgn}\left(\cos\frac{n\pi}{2N}\right),\quad\Delta\varepsilon_{2}^{n}=-\Delta{\varepsilon}\text{
sgn}\left(\cos\frac{n\pi}{4N}\right),\quad\Delta\varepsilon_{3}^{n}=0;\\\
\text{tri-axial:}\qquad&\Delta\varepsilon_{1}^{n}=\Delta{\varepsilon}\text{
sgn}\left(\cos\frac{n\pi}{2N}\right),\quad\Delta\varepsilon_{2}^{n}=\Delta\varepsilon_{3}^{n}=\Delta{\varepsilon}\text{
sgn}\left(\sin\frac{n\pi}{2N}\right),\end{split}$ (28)

with $N=\frac{\varepsilon_{\text{max}}}{\Delta{\varepsilon}}$,
$\varepsilon_{\text{max}}=2\times 10^{-3}\div 2$, and
$\Delta\varepsilon=10^{-4}\div 1$.
In Figures 18 and 19 the results obtained for a uni-axial loading scheme are
presented, for different values of the strain increment, and material case
3D-1 (cf. Table 5). Figures 18(a) and 18(b) display the material response in
terms of the principal stress, $\sigma_{1}$, and inelastic strain,
$\zeta_{1}$, over the principal strain, $\varepsilon_{1}$. Figures 18(c) and
18(d) compare the energy and dissipation rate predicted by TANN with those
computed, with standard ANN, directly using the corresponding definitions for
the free-energy and dissipation rate, Eq. (26). The predictions of TANN are in
good agreement with the constitutive model, independently from the strain
increment, which exceeds considerably the training range. Nevertheless, the
performance of ANN is found to be strongly affected by the values of
$\Delta\varepsilon$. Standard ANN performs poorly for strain increments
smaller and larger than the ones at which it was trained
($\Delta\varepsilon=1\times 10^{3}$). Furthermore, standard ANN predicts
thermodynamically inconsistent outputs.
We emphasize that, even though for relatively large strain increments TANN
predictions are less accurate (see e.g. Fig. 18(d)), predictions remain
thermodynamically consistent. Moreover, the network successfully learns the
Jacobian, $\frac{\partial\sigma}{\partial\epsilon}$, contrary to standard ANN.
As discussed for the 1D case, TANN hence offers compelling capabilities in
replacing complicated constitutive models in numerical simulations, with
smaller computational cost. Figure 19 illustrates the stress predictions, but
in terms of the mean and deviatoric stress (computed from the principal stress
predictions).

(a) $\sigma$ prediction

(b) $\zeta$ prediction

(c) F prediction (TANN), Eq. (26) (ANN).

(d) D prediction (TANN), Eq. (26) (ANN).

Figure 18: Comparison of the predictions of TANN and standard ANN with respect
to the target values, for the uni-axial cyclic loading path, Eq. (28), for
material case 3D-1 (perfect plasticity). Each row represents the prediction at
different $\Delta\varepsilon$ increments.

(a) $p$ computation.

(b) $q$ computation.

Figure 19: Comparison of the stress predictions of TANN and ANN in terms of
mean and deviatoric stress, $p$ (left) and $q$ (right), for the uni-axial
loading path in Figure 18, Eq. (28). Each row represents the prediction at
different $\Delta\varepsilon$ increments.

Figures 20 and 21 present the predictions of TANN and standard ANN for the bi-
axial cyclic path (Eq. (28)), in terms of the principal stresses, for the
material case 3D-1 (cf. Table 5). Similarly, Figures 22 to 24 show the
predictions for the tri-axial loading path. The predictions obtained for the
same loading path but for material cases 3D-2 (hardening) and 3D-3 (softening)
are presented in the SM file. As mentioned above, TANN capabilities of
generalizing the predictions with respect to the values assumed by the input
variables are remarkably good. Standard ANN succeeds in correctly predicting
the stress increments only in a reduced range of the strain increment values,
close to the range of its training data. It is worth noticing that in the
computation of the mean and deviatoric stress, from the principal stress
components,TANN gives relatively higher errors, but still much lower when
compared with those of ANN). The performance of TANN in predicting $p$ and $q$
can be further improved by small modifications to include these invariants as
outputs. However, this kind of optimization exceeds the scope of the current
work.

Similar to the 1D case, it should be stressed that the performance of TANN and
standard ANN can be improved by increasing the dimension of the training data-
set, the number of the hyperparameters (e.g. numbers of hidden layers, etc.).
Nevertheless, the fundamental gap between the two approaches in assuring
thermodynamically consistent quantities still persist.

(a) $\sigma$ prediction.

(b) $\zeta$ prediction.

Figure 20: Comparison of the stress and internal variable predictions of TANN
and standard ANN with respect to the target values, for the bi-axial cyclic
loading path, Eq. (28), for material case 3D-1 (perfect plasticity). Each row
represents the prediction at different $\Delta\varepsilon$ increments.

(a) F prediction (TANN), Eq. (26) (ANN).

(b) D prediction (TANN), Eq. (26) (ANN).

Figure 21: Comparison of the energy and dissipation rate predictions of TANN
and computation according to Eq. (26) for standard ANN with respect to the
target values, for the bi-axial cyclic loading path, Eq. (28), for material
case 3D-1 (perfect plasticity). Each row represents the prediction at
different $\Delta\varepsilon$ increments. Figure 22: Comparison of the stress
predictions of TANN and standard ANN with respect to the target values, for
the tri-axial cyclic loading path, Eq. (28), for material case 3D-1 (perfect
plasticity). Each row represents the prediction at different
$\Delta\varepsilon$ increments. Figure 23: Comparison of the internal variable
predictions of TANN and standard ANN with respect to the target values, for
the tri-axial cyclic loading path, Eq. (28), for material case 3D-1 (perfect
plasticity). Each row represents the prediction at different
$\Delta\varepsilon$ increments.

(a) F prediction (TANN), Eq. (26) (ANN).

(b) D prediction (TANN), Eq. (26) (ANN).

Figure 24: Comparison of the energy and dissipation rate predictions of TANN
and computation according to Eq. (26) for standard ANN with respect to the
target values, for the tri-axial cyclic loading path, Eq. (28), for material
case 3D-1 (perfect plasticity). Each row represents the prediction at
different $\Delta\varepsilon$ increments.

##  6 Concluding remarks

A new class of artificial neural networks models to replace constitutive laws
and predict the material response at the material point level was proposed.
The two basic laws of thermodynamics were directly encoded in the architecture
of the model, which we refer to as Thermodynamics-based Neural Network (TANN).
Our approach was inspired by the so-called Physics-Informed Neural Networks
(PINNs) [49], where the automatic differentiation was used to perform the
numerical calculation of the derivative of a neural network with respect to
its inputs.

The numerical requirements regarding the mathematical class of appropriate
activation functions to be used together with automatic differentiation were
investigated. More specifically, the internal restrictions, derived from the
first law of thermodynamics, require activation functions whose second
gradient does not vanish. This new problem and its remedy was extensively
explored and discussed in the manuscript.

TANN, relying on an incremental formulation and on the theoretical
developments in [25], posses the special feature that the entire constitutive
response of a material can be derived from definition of only two scalar
functions: the free-energy and the dissipation rate. This assures
thermodynamically consistent predictions both for seen and unseen data.
Differently from the standard ANN approaches, TANN does not have to identify,
through learning, the underlying thermodynamic laws. Indeed, predictions of
standard ANNs may be thermodynamically inconsistent, even though the training
of the network has been performed on consistent material data. Being aware of
physics, TANNs are found to be a robust approach. If the training data-sets
are not thermodynamically consistent, the training operation of the network
will be unsuccessful, contrary to the standard ANN approach (see SM).

For the cases here investigated, we showed that TANNs are characterized by
high accuracy of the predictions, higher than those of standard approaches.
The integration of thermodynamic principles inside the network renders TANN’s
ability of generalization (i.e., make predictions for loading paths different
from those used in the training operation) remarkably good. Consequently, TANN
is an excellent candidate for replacing constitutive calculations at Finite
Element incremental formulations. Moreover, thanks to the implementation of
the free-energy in the network predictions and its thermodynamical relation
with the stresses, the Jacobian
$\frac{\partial\Delta\sigma}{\Delta\varepsilon}$ at the material point level
is better predicted even for increments far outside the training data-set
range. As a result quadratic convergence in implicit formulations can be
preserved, reducing the calculation cost.

Further extensions of TANN in a wide range of applications, for complex
materials, are straightforwards, as the thermodynamics principles hold true
for any known class of material, at any length (micro- and macro-scale).

## Acknowledgments

The author I.S. would like to acknowledge the support of the European Research
Council (ERC) under the European Union Horizon 2020 research and innovation
program (Grant agreement ID 757848 CoQuake).

## Appendix A. Derivation of the incremental material formulation

By differentiating the energy expressions (13) and rearranging the terms, we
obtain the following non-linear incremental relations

$\displaystyle\dot{\sigma}=\partial_{\varepsilon\varepsilon}\text{{F}}\cdot\varepsilon+\sum_{k}\partial_{\varepsilon\zeta_{k}}\text{{F}}\cdot\dot{\zeta}_{k}+\partial_{\varepsilon\theta}\text{{F}}\,\dot{\theta}$
(29a)
$\displaystyle-\dot{\chi}_{i}=\partial_{\zeta_{i}\varepsilon}\text{{F}}\cdot\varepsilon+\sum_{k}\partial_{\zeta_{i}\zeta_{k}}\text{{F}}\cdot\dot{\zeta}_{k}+\partial_{\zeta_{i}\theta}\text{{F}}\,\dot{\theta}$
(29b)
$\displaystyle-\dot{\text{{S}}}=\partial_{\theta\varepsilon}\text{{F}}\cdot\varepsilon+\sum_{k}\partial_{\theta\zeta_{k}}\text{{F}}\cdot\dot{\zeta}_{k}+\partial_{\theta\theta}\text{{F}}\,\dot{\theta},$
(29c)

where the following notation is adopted

$\displaystyle\partial_{\varepsilon\varepsilon}\text{{F}}=$
$\displaystyle\dfrac{\partial^{2}\text{{F}}}{\partial\varepsilon_{ij}\partial\varepsilon_{kl}},\quad\partial_{\varepsilon\zeta_{k}}\text{{F}}=$
$\displaystyle\dfrac{\partial^{2}\text{{F}}}{\partial\varepsilon_{ij}\partial\zeta_{k}},$
$\displaystyle\partial_{\varepsilon\theta}\text{{F}}=$
$\displaystyle\dfrac{\partial^{2}\text{{F}}}{\partial\varepsilon_{ij}\partial\theta},\quad\partial_{\theta\theta}\text{{F}}=$
$\displaystyle\dfrac{\partial^{2}\text{{F}}}{\partial\theta^{2}}.$

We introduce the thermodynamic dissipative stresses
$\mathcal{X}^{\dagger}=(X_{1},\ldots,X_{N})$ with

$X_{i}:=\dfrac{\partial\text{{D}}}{\partial\dot{\zeta}_{i}}\qquad\forall\>i\in[1,N].$
(30)

For a rate-independent material, the dissipation is a homogeneous first-order
function in the internal variable rates $\dot{\zeta}_{i}$ [25]. This
homogeneity can be expressed by the Euler’s relation

$\text{{D}}=\sum_{i=1}^{N}\dfrac{\partial\text{{D}}}{\partial\dot{\zeta}_{i}}\cdot\dot{\zeta}_{i}=\sum_{i}X_{i}\cdot\dot{\zeta}_{i},$
(31)

which, together with (11), implies

$\sum_{i=1}^{N}\left(X_{i}-\chi_{i}\right)\cdot\dot{\zeta}_{i}=0$ (32)

Ziegler’s orthogonality condition [56] is further assumed, i.e.,
$X_{i}=\chi_{i}$ $\forall\>i\in[1,N]$.
Being D homogeneous first-order function in $\dot{\zeta}_{i}$, the Legendre
transform, conjugate to $X_{i}$, is degenerate, that is equal to zero, and
represents the yield function
$y=\tilde{y}(\theta,\varepsilon,\mathcal{Z},\mathcal{X}^{\dagger})$, i.e.

$\lambda y=\sum_{i}X_{i}\cdot\dot{\zeta}_{i}-\text{{D}}=0,$ (33)

where $\lambda$ is a non-negative multiplier. From the properties of Legendre
transform, the following flow rules must hold

$\dot{\zeta}_{i}=\lambda\dfrac{\partial y}{\partial
X_{i}}\quad\forall\>i\in[1,N].$ (34)

Since $\lambda\geq 0$ and $\lambda y=0$, $y\leq 0$. If $y=0$, the following
consistency equation is met

$\dot{y}=\frac{\partial
y}{\partial\varepsilon}\cdot\dot{\varepsilon}+\sum_{i=1}^{N}\frac{\partial
y}{\partial\zeta_{i}}\cdot\dot{\zeta}_{i}+\sum_{i=1}^{N}\frac{\partial
y}{\partial X_{i}}\cdot\dot{X}_{i}+\frac{\partial
y}{\partial\theta}\;\dot{\theta}=0.$ (35)

By further using the flow rules (34) and Ziegler’s normality condition, we
obtain

$\lambda=-\dfrac{\mathcal{C}_{\varepsilon}}{B}\cdot\dot{\varepsilon}-\dfrac{\mathcal{C}_{\theta}}{B}\cdot\dot{\theta},$
(36)

with

$\mathcal{C}_{\varepsilon}=\frac{\partial
y}{\partial\varepsilon}-\sum_{i=1}^{N}\frac{\partial y}{\partial
X_{i}}\cdot\partial_{\zeta_{i}\varepsilon}\text{{F}},$
$\mathcal{C}_{\theta}=\frac{\partial
y}{\partial\theta}-\sum_{i=1}^{N}\frac{\partial y}{\partial
X_{i}}\cdot\partial_{\zeta_{i}\theta}\text{{F}},$

and

$B=\sum_{i=1}^{N}\frac{\partial y}{\partial\zeta_{i}}\cdot\frac{\partial
y}{\partial X_{i}}-\sum_{i=1}^{N}\frac{\partial y}{\partial
X_{i}}\left(\sum_{k=1}^{N}\partial_{\zeta_{k}\varepsilon}\text{{F}}\cdot\frac{\partial
y}{\partial X_{k}}\right).$

Finally, we arrive to the following, incremental non-linear formulation, for
$y=0$,

$\dot{\Xi}=\mathcal{M}|_{y=0}\;\dot{\xi},\quad\text{ with
}\quad\dot{\Xi}=\begin{bmatrix}\dot{\sigma}\\\ -\dot{X}_{i}\\\
-\dot{\text{{S}}}\\\ \dot{\zeta_{i}}\\\
\lambda\end{bmatrix},\quad\dot{\xi}=\begin{bmatrix}\dot{\varepsilon}\\\
\dot{\theta}\end{bmatrix},\quad\mathcal{M}|_{y=0}=\begin{bmatrix}\textsf{M}_{\varepsilon\varepsilon}&\textsf{M}_{\varepsilon\theta}\\\
\textsf{M}_{\zeta\varepsilon}&\textsf{M}_{\zeta\theta}\\\
\textsf{M}_{\theta\varepsilon}&\textsf{M}_{\theta\theta}\\\
-\dfrac{\mathcal{C}_{\varepsilon}}{B}\cdot\frac{\partial y}{\partial
X_{i}}&-\dfrac{\mathcal{C}_{\theta}}{B}\cdot\frac{\partial y}{\partial
X_{i}}\\\
-\dfrac{\mathcal{C}_{\varepsilon}}{B}\cdot&-\dfrac{\mathcal{C}_{\theta}}{B}\end{bmatrix},$
(37)

and

$\displaystyle\textsf{M}_{\varepsilon\varepsilon}=$
$\displaystyle\partial_{\varepsilon\varepsilon}\text{{F}}-\sum_{k}\partial_{\varepsilon\zeta_{k}}\text{{F}}\cdot\left(\frac{\mathcal{C}_{\varepsilon}}{B}\cdot\frac{\partial
y}{\partial X_{k}}\right),$ $\displaystyle\textsf{M}_{\varepsilon\theta}=$
$\displaystyle\partial_{\varepsilon\theta}\text{{F}}-\sum_{k}\partial_{\varepsilon\zeta_{k}}\text{{F}}\cdot\left(\frac{\mathcal{C}_{\theta}}{B}\cdot\frac{\partial
y}{\partial X_{k}}\right),$ $\displaystyle\textsf{M}_{\zeta\varepsilon}=$
$\displaystyle\partial_{\zeta_{i}\varepsilon}\text{{F}}-\sum_{k}\partial_{\zeta_{i}\zeta_{k}}\text{{F}}\cdot\left(\frac{\mathcal{C}_{\varepsilon}}{B}\cdot\frac{\partial
y}{\partial X_{k}}\right),$ $\displaystyle\textsf{M}_{\zeta\theta}=$
$\displaystyle\partial_{\zeta_{i}\theta}\text{{F}}-\sum_{k}\partial_{\zeta_{i}\zeta_{k}}\text{{F}}\cdot\left(\frac{\mathcal{C}_{\theta}}{B}\cdot\frac{\partial
y}{\partial X_{k}}\right),$ $\displaystyle\textsf{M}_{\theta\varepsilon}=$
$\displaystyle\partial_{\theta\varepsilon}\text{{F}}-\sum_{k}\partial_{\theta\zeta_{k}}\text{{F}}\cdot\left(\frac{\mathcal{C}_{\varepsilon}}{B}\cdot\frac{\partial
y}{\partial X_{k}}\right),$ $\displaystyle\textsf{M}_{\theta\theta}=$
$\displaystyle\partial_{\theta\theta}\text{{F}}-\sum_{k}\partial_{\theta\zeta_{k}}\text{{F}}\cdot\left(\frac{\mathcal{C}_{\theta}}{B}\cdot\frac{\partial
y}{\partial X_{k}}\right).$

In case of $y<0$, relation (37) becomes

$\dot{\Xi}=\mathcal{M}|_{y<0}\;\dot{\xi},\quad\text{ with
}\quad\mathcal{M}|_{y<0}=\begin{bmatrix}\partial_{\varepsilon\varepsilon}\text{{F}}&\partial_{\varepsilon\theta}\text{{F}}\\\
\partial_{\zeta_{i}\varepsilon}\text{{F}}&\partial_{\zeta_{i}\theta}\text{{F}}\\\
\partial_{\theta\varepsilon}\text{{F}}&\partial_{\theta\theta}\text{{F}}\\\
0&0\\\ 0&0\end{bmatrix}.$ (38)

## References

## References

  * Abadi et al. [2016] Abadi, M., Barham, P., Chen, J., Chen, Z., Davis, A., Dean, J., Devin, M., Ghemawat, S., Irving, G., Isard, M., et al., 2016. Tensorflow: A system for large-scale machine learning. In: 12th USENIX Symposium on Operating Systems Design and Implementation. pp. 265–283.
  * Anand et al. [2012] Anand, L., Aslan, O., Chester, S. A., 2012. A large-deformation gradient theory for elastic–plastic materials: Strain softening and regularization of shear bands. International Journal of Plasticity 30-31, 116 – 143.
  * Bakhvalov and Panasenko [1989] Bakhvalov, N., Panasenko, G., 1989. Homogenisation: Averaging Processes in Periodic Media: Mathematical Problems in the Mechanics of Composite Materials.
  * Baydin et al. [2017] Baydin, A. G., Pearlmutter, B. A., Radul, A. A., Siskind, J. M., 2017. Automatic differentiation in machine learning: a survey. The Journal of Machine Learning Research 18 (1), 5595–5637.
  * Bessa et al. [2017] Bessa, M., Bostanabad, R., Liu, Z., Hu, A., Apley, D. W., Brinson, C., Chen, W., Liu, W. K., 2017. A framework for data-driven analysis of materials under uncertainty: Countering the curse of dimensionality. Computer Methods in Applied Mechanics and Engineering 320, 633–667.
  * Bogacki and Shampine [1989] Bogacki, P., Shampine, L. F., 1989. A 3 (2) pair of Runge-Kutta formulas. Applied Mathematics Letters 2 (4), 321–325.
  * Chen and Chen [1995] Chen, T., Chen, H., 1995. Universal approximation to nonlinear operators by neural networks with arbitrary activation functions and its application to dynamical systems. IEEE Transactions on Neural Networks 6 (4), 911–917.
  * Collins-Craft et al. [2020] Collins-Craft, N. A., Stefanou, I., Sulem, J., Einav, I., 2020. A cosserat breakage mechanics model for brittle granular media. Journal of the Mechanics and Physics of Solids, 103975.
  * Cybenko [1989] Cybenko, G., 1989. Approximation by superpositions of a sigmoidal function. Mathematics of control, signals and systems 2 (4), 303–314.
  * Dozat [2016] Dozat, T., 2016. Incorporating Nesterov momentum into Adam.
  * Eggersmann et al. [2019] Eggersmann, R., Kirchdoerfer, T., Reese, S., Stainier, L., Ortiz, M., 2019. Model-free data-driven inelasticity. Computer Methods in Applied Mechanics and Engineering 350, 81–99.
  * Eijnden et al. [2016] Eijnden, A. P. V. D., Bésuelle, P., Collin, F., Chambon, R., Desrues, J., 2016\. Modeling the strain localization around an underground gallery with a hydro-mechanical double scale model ; effect of anisotropy. Computers and Geotechnics.
  * Einav [2012] Einav, I., 2012. The unification of hypo-plastic and elasto-plastic theories. International Journal of Solids and Structures 49 (11-12), 1305–1315.
  * Einav et al. [2007] Einav, I., Houlsby, G., Nguyen, G., 2007. Coupled damage and plasticity models derived from energy and dissipation potentials. International Journal of Solids and Structures 44 (7-8), 2487–2508.
  * Feyel [2003] Feyel, F., 2003. A multilevel finite element method (FE2) to describe the response of highly non-linear structures using generalized continua. Computer Methods in Applied Mechanics and Engineering 192 (28-30), 3233–3244.
  * Frankel et al. [2019] Frankel, A. L., Jones, R. E., Alleman, C., Templeton, J. A., 2019. Predicting the mechanical response of oligocrystals with deep learning. Computational Materials Science 169, 109099.
  * Gajek et al. [2020] Gajek, S., Schneider, M., Böhlke, T., 2020. On the micromechanics of deep material networks. Journal of the Mechanics and Physics of Solids, 103984.
  * Geron [2015] Geron, A., 2015. Hands-on MachineLearning with Scikit-Learn & Tensorflow. Vol. 1. O’Reilly Media.
  * Géron [2019] Géron, A., 2019. Hands-On Machine Learning with Scikit-Learn, Keras, and TensorFlow: Concepts, Tools, and Techniques to Build Intelligent Systems. O’Reilly Media.
  * Ghaboussi et al. [1991] Ghaboussi, J., Garrett, J. H., Wu, X., 1991. Knowledge‐based modeling of material behavior with neural networks. Journal of Engineering Mechanics 117 (1), 132–153.
  * Ghaboussi and Sidarta [1998] Ghaboussi, J., Sidarta, D., 1998. New nested adaptive neural networks (nann) for constitutive modeling. Computers and Geotechnics 22 (1), 29–52.
  * Ghavamian and Simone [2019] Ghavamian, F., Simone, A., 2019. Accelerating multiscale finite element simulations of history-dependent materials using a recurrent neural network. Computer Methods in Applied Mechanics and Engineering 357, 112594.
  * Heider et al. [2020] Heider, Y., Wang, K., Sun, W., 2020. SO(3)-invariance of informed-graph-based deep neural network for anisotropic elastoplastic materials. Computer Methods in Applied Mechanics and Engineering 363, 112875.
  * Houlsby and Puzrin [2000] Houlsby, G., Puzrin, A., 2000. A thermomechanical framework for constitutive models for rate-independent dissipative materials. International journal of Plasticity 16 (9), 1017–1047.
  * Houlsby and Puzrin [2007] Houlsby, G. T., Puzrin, A. M., 2007. Principles of hyperplasticity: an approach to plasticity theory based on thermodynamic principles. Springer Science & Business Media.
  * Hu and Hwang [2002] Hu, Y. H., Hwang, J.-N., 2002. Handbook of neural network signal processing.
  * Huang et al. [2020] Huang, D. Z., Xu, K., Farhat, C., Darve, E., 2020. Learning constitutive relations from indirect observations using deep neural networks. Journal of Computational Physics, 109491.
  * Ibanez et al. [2018] Ibanez, R., Abisset-Chavanne, E., Aguado, J. V., Gonzalez, D., Cueto, E., Chinesta, F., 2018. A manifold learning approach to data-driven computational elasticity and inelasticity. Archives of Computational Methods in Engineering 25 (1), 47–57.
  * Ibañez et al. [2017] Ibañez, R., Borzacchiello, D., Aguado, J. V., Abisset-Chavanne, E., Cueto, E., Ladevèze, P., Chinesta, F., 2017. Data-driven non-linear elasticity: constitutive manifold construction and problem discretization. Computational Mechanics 60 (5), 813–826.
  * Jung and Ghaboussi [2006] Jung, S., Ghaboussi, J., 2006. Neural network constitutive model for rate-dependent materials. Computers & Structures 84 (15-16), 955–963.
  * Karpatne et al. [2017] Karpatne, A., Watkins, W., Read, J., Kumar, V., 2017. Physics-guided neural networks (pgnn): An application in lake temperature modeling. arXiv preprint arXiv:1710.11431.
  * Kirchdoerfer and Ortiz [2016] Kirchdoerfer, T., Ortiz, M., 2016. Data-driven computational mechanics. Computer Methods in Applied Mechanics and Engineering 304, 81–101.
  * Kirchdoerfer and Ortiz [2018a] Kirchdoerfer, T., Ortiz, M., 2018a. Data-driven computing in dynamics. International Journal for Numerical Methods in Engineering 113 (11), 1697–1710.
  * Kirchdoerfer and Ortiz [2018b] Kirchdoerfer, T., Ortiz, M., 2018b. Data-driven computing in dynamics. International Journal for Numerical Methods in Engineering 113 (11), 1697–1710.
  * Lee [2018] Lee, J., 2018. Analysis methodology of inelastic constitutive parameter using state space method and neural network. International Journal of Engineering & Technology 7 (3.34).
  * Lefik and Schrefler [2003] Lefik, M., Schrefler, B. A., 2003. Artificial neural network as an incremental non-linear constitutive model for a finite element code. Computer methods in applied mechanics and engineering 192 (28-30), 3265–3283.
  * Liu and Wu [2019] Liu, Z., Wu, C., 2019. Exploring the 3d architectures of deep material network in data-driven multiscale mechanics. Journal of the Mechanics and Physics of Solids 127, 20–46.
  * Lloberas Valls et al. [2019] Lloberas Valls, O., Raschi Schaw, M., Huespe, A. E., Oliver Olivella, X., 2019. Reduced finite element square techniques (rfe2): towards industrial multiscale fe software. In: COMPLAS 2019: XV International Conference on Computational Plasticity: Fundamentals and Applications. International Centre for Numerical Methods in Engineering (CIMNE), pp. 157–169.
  * Lu et al. [2019] Lu, X., Giovanis, D. G., Yvonnet, J., Papadopoulos, V., Detrez, F., Bai, J., 2019\. A data-driven computational homogenization method based on neural networks for the nonlinear anisotropic electrical response of graphene/polymer nanocomposites. Computational Mechanics 64 (2), 307–321.
  * Mariano and Galano [2015] Mariano, P. M., Galano, L., 2015. Fundamentals of the Mechanics of Solids. Springer.
  * Masi et al. [2020] Masi, F., Stefanou, I., Maffi-Berthier, V., Vannucci, P., 2020. A discrete element method based-approach for arched masonry structures under blast loads. Engineering Structures 216, 110721.
  * Masi et al. [2018] Masi, F., Stefanou, I., Vannucci, P., 2018. A study on the effects of an explosion in the Pantheon of Rome. Engineering Structures 164, 259–273.
  * Maugin and Muschik [1994] Maugin, G. A., Muschik, W., 1994. Thermodynamics with internal variables. Part I. General concepts.
  * McKinney et al. [2011] McKinney, W., et al., 2011. Pandas: a foundational Python library for data analysis and statistics. Python for High Performance and Scientific Computing 14 (9).
  * Meurer et al. [2017] Meurer, A., Smith, C. P., Paprocki, M., Čertík, O., Kirpichev, S. B., Rocklin, M., Kumar, A., Ivanov, S., Moore, J. K., Singh, S., et al., 2017\. SymPy: symbolic computing in Python. PeerJ Computer Science 3, e103.
  * Mitchell et al. [1997] Mitchell, T. M., et al., 1997. Machine learning. 1997. Burr Ridge, IL: McGraw Hill 45 (37), 870–877.
  * Mozaffar et al. [2019] Mozaffar, M., Bostanabad, R., Chen, W., Ehmann, K., Cao, J., Bessa, M., 2019. Deep learning predicts path-dependent plasticity. Proceedings of the National Academy of Sciences 116 (52), 26414–26420.
  * Nitka et al. [2011] Nitka, M., Combe, G., Dascalu, C., Desrues, J., mar 2011. Two-scale modeling of granular materials: a DEM-FEM approach. Granular Matter 13 (3), 277–281.
  * Raissi et al. [2019] Raissi, M., Perdikaris, P., Karniadakis, G. E., 2019. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics 378, 686–707.
  * Rattez et al. [2018a] Rattez, H., Stefanou, I., Sulem, J., 2018a. The importance of thermo-hydro-mechanical couplings and microstructure to strain localization in 3d continua with application to seismic faults. part i: Theory and linear stability analysis. Journal of the Mechanics and Physics of Solids 115, 54 – 76\.
  * Rattez et al. [2018b] Rattez, H., Stefanou, I., Sulem, J., Veveakis, M., Poulet, T., 2018b. The importance of thermo-hydro-mechanical couplings and microstructure to strain localization in 3d continua with application to seismic faults. part ii: Numerical implementation and post-bifurcation analysis. Journal of the Mechanics and Physics of Solids 115, 1 – 29.
  * Settgast et al. [2019] Settgast, C., Abendroth, M., Kuna, M., 2019. Constitutive modeling of plastic deformation behavior of open-cell foam structures using neural networks. Mechanics of Materials 131, 1–10.
  * Virtanen et al. [2020] Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., et al., 2020. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature methods 17 (3), 261–272.
  * Xu et al. [2020] Xu, K., Huang, D. Z., Darve, E., 2020. Learning constitutive relations using symmetric positive definite neural networks. arXiv preprint arXiv:2004.00265.
  * Yu and Carrillo [2019] Yu, W., Carrillo, D., 2019. Chapter 1 - hierarchical dynamic neural networks for cascade system modeling with application to wastewater treatment. In: Alanis, A. Y., Arana-Daniel, N., López-Franco, C. (Eds.), Artificial Neural Networks for Engineering Applications. Academic Press, pp. 1 – 8.
  * Ziegler [2012] Ziegler, H., 2012. An introduction to thermomechanics. Elsevier.