A comparative study on different neural network architectures to model
  inelasticity
The mathematical formulation of constitutive models to describe the
path-dependent, i.e., inelastic, behavior of materials is a challenging task
and has been a focus in mechanics research for several decades. There have been
increased efforts to facilitate or automate this task through data-driven
techniques, impelled in particular by the recent revival of neural networks
(NNs) in computational mechanics. However, it seems questionable to simply not
consider fundamental findings of constitutive modeling originating from the
last decades research within NN-based approaches. Herein, we propose a
comparative study on different feedforward and recurrent neural network
architectures to model inelasticity. Within this study, we divide the models
into three basic classes: black box NNs, NNs enforcing physics in a weak form,
and NNs enforcing physics in a strong form. Thereby, the first class of
networks can learn constitutive relations from data while the underlying
physics are completely ignored, whereas the latter two are constructed such
that they can account for fundamental physics, where special attention is paid
to the second law of thermodynamics in this work. Conventional linear and
nonlinear viscoelastic as well as elastoplastic models are used for training
data generation and, later on, as reference. After training with random walk
time sequences containing information on stress, strain, and, for some models,
internal variables, the NN-based models are compared to the reference solution,
whereby interpolation and extrapolation are considered. Besides the quality of
the stress prediction, the related free energy and dissipation rate are
analyzed to evaluate the models. Overall, the presented study enables a clear
recording of the advantages and disadvantages of different NN architectures to
model inelasticity and gives guidance on how to train and apply these models.
Max Rosenkranz, Karl A. Kalina, Jörg Brummund, Markus Kästner
# A comparative study on different neural network architectures to model inelasticity

###### Abstract

The mathematical formulation of constitutive models to describe the path-dependent, i. e., inelastic, behavior of materials is a challenging task and has been a focus in mechanics research for several decades. There have been increased efforts to facilitate or automate this task through data-driven techniques, impelled in particular by the recent revival of neural networks (NNs) in computational mechanics. However, it seems questionable to simply not consider fundamental findings of constitutive modeling originating from the last decades research within NN-based approaches. Herein, we propose a comparative study on different feedforward and recurrent neural network architectures to model inelasticity. Within this study, we divide the models into three basic classes: black box NNs, NNs enforcing physics in a weak form, and NNs enforcing physics in a strong form. Thereby, the first class of networks can learn constitutive relations from data while the underlying physics are completely ignored, whereas the latter two are constructed such that they can account for fundamental physics, where special attention is paid to the second law of thermodynamics in this work. Conventional linear and nonlinear viscoelastic as well as elastoplastic models are used for training data generation and, later on, as reference. After training with random walk time sequences containing information on stress, strain, and - for some models - internal variables, the NN-based models are compared to the reference solution, whereby interpolation and extrapolation are considered. Besides the quality of the stress prediction, the related free energy and dissipation rate are analyzed to evaluate the models. Overall, the presented study enables a clear recording of the advantages and disadvantages of different NN architectures to model inelasticity and gives guidance on how to train and apply these models.

neural networks recurrent neural networks enforcing physics constitutive modeling thermodynamic consistency viscoelasticity plasticity

## 1 Introduction

Accurately describing the behavior of materials under mechanical loading by constitutive models has been a focus in mechanics research for several decades now. The formulation and parametrization of constitutive models is however still a challenging task especially for materials showing path-dependent, i. e., inelastic, behavior. As an alternative to traditional models, _data-based_ or _data-driven_ techniques are very promising and have the potential to improve or replace conventional models. These techniques have become increasingly popular in the computational mechanics community
during the last years [1, 2], where the application of neural networks (NNs) is probably the most common technique. In the following, a brief overview on NNs in constitutive modeling is given.

The concept of using NNs in constitutive modeling was initially put out by Ghaboussi et al. [3] in the early 1990s. However, in this early stage, generally pure black-box techniques were employed, i.e., networks that do not account for any physical principles and can therefore only accurately recreate the training data, in this case composed of stress-strain pairs, but perform badly when extrapolating. To address this issue, a relatively new approach in _NN-based constitutive modeling_, and scientific _machine learning_ (ML) in general, is to integrate crucial underlying physics in either a strong or weak form. These methods, known as _physics-informed_[4, 5], _mechanics-informed_[6, 7], _physics-augmented_[8, 9], _physics-constrained_[10], or _thermodynamics-based_[11], improve extrapolation capability and allow for the use of sparse training data.

The easiest material behavior to model is elasticity, since here a suitable model only needs to predict the stresses for specific deformation states. In the context of ML, the works [12, 13] seek to approximate the elastic potential by using a _feedforward neural network_ (FNN) with three deformation-type invariants as input. Thus, a number of constitutive requirements is fulfilled by construction, e.g., _thermodynamic consistency_, objectivity, or material symmetry. However, training of these models directly requires the elastic potential. Meanwhile, FNNs using invariants as input and the hyperelastic potential as output are a very well established approach [14, 15, 16, 17, 18, 8, 9]. Thereby, an improved training is applied that allows calibration of the network directly by tuples of stress and strain, i.e., the derivative of energy with respect to the deformation is included into the loss, which is also called Sobolev training [19, 20]. Alternatively, a network previously trained to predict stress coefficients can be used to construct a pseudopotential, thus ensuring thermodynamic consistency of NN-based elastic models a posteriori [21].

Compared to elasticity, the modeling of path-dependent, i.e., inelastic, constitutive behavior by NN-based approaches is more complex. Some early proposals [22, 23, 24, 25, 26] can already achieve quite good predictions by, for example, adding stress and strain states from previous time steps into the input layer of an FNN [25]. This allows the network to indirectly learn a kind of evolution equation. The model [22] uses internal variables of the material to reliably reproduce stress-strain curves of a viscoplastic material. Alternatively, load history dependent behavior can also be represented without the availability of the internal variables by so-called _recurrent neural networks_ (RNNs), which have been shown to be universal and accurate, particularly for more sophisticated recurrent cells, e.g., according to Hochreiter and Schmidhuber [27]. RNNs, especially _long short-term memory_ (LSTM) cells, have been intensively used to model inelasticity, e.g., in the works [28, 29, 30, 31], and are very promising regarding their prediction quality. Very recently, spiking LSTMs, which enable a massive reduction in memory and energy consumption over conventional neural networks, have been applied to model isotropic hardening plasticity [32]. In addition, a new type of RNN named linearized minimal state cell [33] prevents that its response depends on path-sampling and is therefore advantageous in modeling elastoplasticity. This approach has been used for both 2D and 3D datasets matching to the real mechanical behavior of an aluminum alloy as determined by simulations of crystal plasticity [34]. A further promising approach to represent anisotropic elastoplasticity combines Lie algebra with RNNs [35]. Finally, although being trained entirely on monotonic data, a hybrid model [36] combining a data-driven encoder and a physics-based decoder allows for accurate predictions of elastoplastic unloading/reloading paths.

Despite the great progress in ML-based constitutive modeling, the NN approaches to describing inelasticity mentioned so far are united by their lack of knowledge of the _second law of thermodynamics_. However, following Masi et al. [11], the incorporation of such fundamental physical principles offers decisive advantages as a more targeted and therefore faster training, that also requires only a small amount of data, and a significantly improved extrapolation capability. A data-driven framework called _deep material network_ is shown in the works [37, 38, 39, 40]. Thereby, the response of a representative volume element is reproduced by a network including a collection of connected mechanistic building blocks with analytical homogenization solutions, which enables to describe a complex effective response without the loss of essential physics. Within the works [41, 42, 43, 44, 45] the idea to replace parts of classical models with NNs is pursued to achieve this. E. g., the yield function or the evolution equations are described by FNNs instead of using a particular model. A coupling of NNs to the so-called micro-sphere approach is shown in the work [46]. A more freely formulated approach for the consideration of rate-independent inelasticity based on an adapted network architecture consisting of two FNNs is presented in the works [47, 11]. Thereby, the first network is used to learn the internal variables' evolution and the second for the approximation of the free energy, where the training procedure requires internal state variables. To account for the thermodynamic consistency, i.e., that the rate of dissipation \(\mathbf{D}\) is always greater equal to zero, this term, which follows from the free energy, is added into the loss function. A similar approach tailored for the modeling of inelasticity is shown in the work [48]. In contrast to the former model, internal state variables capturing the path-dependency are inferred automatically from the hidden state of an RNN. Thus, this method has the advantage of requiring only stresses and strains for training. The two mentioned models [11, 48] nevertheless have the weakness that the requirement \(\mathbf{D}\geq 0\) is not satisfied by design for arbitrary load cases, but is merely enforced by adding a penalty term to the loss function. Within the works [49, 7], on the other hand, thermodynamic consistency is fulfilled by design of the network architecture. This is achieved by combining the concept of _generalized standard materials_[50] with _input convex neural networks_ (ICNNs) [51]. Within the mentioned works, an application to viscoelasticity is shown. In a similar approach, finite viscoelasticity is
modeled by replacing the Helmholtz free energy function and dissipation potential with data-driven functions that a priori satisfy the second law of thermodynamics, using neural ordinary differential equations (NODEs) [52].

After the brief overview given above, it can be summarized that there are a variety of NN-based approaches to modeling inelasticity, with very different levels of incorporated physics. Most approaches were applied exclusively for describing one specific material class, elastoplasticity or viscoelasticity. Thus, this work aims on bringing the different approaches into a uniform framework and comparing them by applying to both elastoplastic as well as viscoelastic data in the 1D case. Thereby, special attention is paid to the fulfillment of the second law of thermodynamics. With regard to this, a division of the models into three basic classes is done: _black box NNs, NNs enforcing physics in a weak form_, and _NNs enforcing physics in a strong form_. Networks belonging to the first class learn constitutive relations from data while the underlying physics are completely ignored, whereas the latter two are constructed such that they can account for fundamental physical principles. However, NNs enforcing physics in a weak form do not necessarily satisfy the second law for arbitrary load cases, which is due to the fact that \(D\geq 0\) is only integrated into the loss function by a penalty term. In contrast, the network architecture of NNs enforcing physics in a strong form is designed in such a way that this condition is fulfilled in every case, i.e., by construction. In this paper, conventional linear and nonlinear viscoelastic as well as elastoplastic models are used for training data generation and, later on, as reference. After training with random walk time sequences containing information on stress, strain, and - for some models - internal variables, the NN-based models are compared to the reference solution. Besides the quality of the stress prediction, the predicted free energy and dissipation rate are analyzed to evaluate the models for both interpolation and extrapolation. In addition to the provided comparison, some of the NN-based models are extended and/or modified at several points, in particular the approaches belonging to NNs enforcing physics in a strong form.

The organization of the paper is as follows: In Sects. 2 and 3, the basics of constitutive modeling in continuum solid mechanics as well as artificial neural networks are given, respectively. After this, the considered NN-based constitutive models are introduced in Sect. 4. The generation of the database for training is given in Sect. 5. This is followed by a study of the prediction quality of the various NN-based models in Sect. 6. After a discussion of the results, the paper is closed by concluding remarks and an outlook to necessary future work in Sect. 7.

## 2 Classical constitutive models

The description of the behavior of materials requires constitutive equations. A framework for the formulation of these equations for different types of material behavior is shown in this section. According to Haupt [53], in terms of their behavior, materials can be categorized into the four classes (i) _elasticity_, (ii) _elastplasticity_, (iii) _viscoelasticity_, and (iv) _viscoelastoplasticity_. Classes (ii) to (iv) are referred to as _inelastic_ and exhibit dissipative behavior. In order to ensure the irreversibility of such processes, _thermodynamic consistency_ must be taken into account during the formulation of the material laws.

### General framework

#### 2.1.1 Dissipation inequality

Using entropy balance, energy balance and the second law of thermodynamics, an expression known as Clausius-Duhem inequality can be formulated, which, assuming isothermal 1D processes, takes the form

\[\sigma\dot{\epsilon}-\psi\geq 0\quad, \tag{1}\]

with \(\sigma\) being the stress, \(\epsilon\) the strain and \(\psi\) the Helmholtz free energy density, which for the sake of brevity will simply be referred to as free energy in the following. Starting from Eq. (1), depending on the particular choice of \(\psi\), different constitutive models which strictly satisfy the second law can be derived. Generally, the free energy \(\psi:=\psi\left(\epsilon,\xi^{1},\xi^{2},\ldots,\xi^{N}\right)\) is a function of \(\epsilon\) and the internal variables \(\xi^{a},\alpha\in\left\{1,\ldots,N\right\}\). This set of internal variables is required to describe the load history dependent internal state of a material point and does not necessarily represent measurable physical quantities. To shorten notation, the internal variables are summarized in the generalized vector \(\mathbf{\xi}\in\mathbb{R}^{N}\) where appropriate in the following. Applying the principle of equipresence [53], the stress \(\sigma:=\sigma\left(\epsilon,\mathbf{\xi}\right)\) is assumed to be a function of the same set of variables. Evaluating Eq. (1) yields

\[D:=\left(\sigma-\frac{\partial\psi}{\partial\epsilon}\right)\dot{\epsilon}- \underbrace{\frac{\partial\psi}{\partial\mathbf{\xi}}}_{=:-\tau}\cdot\dot{\mathbf{\xi}} \geq 0\quad, \tag{2}\]

where \(\mathbf{a}\cdot\mathbf{b}\) is the scalar product of two vectors \(\mathbf{a},\mathbf{b}\in\mathbb{R}^{N}\). The quantities \(D\) and \(\tau\in\mathbb{R}^{N}\) denote the dissipation rate and the vector of _thermodynamic conjugate forces_, also called _internal forces_[50], with respect to the internal variables \(\mathbf{\xi}\)
In order to comply with inequality (2), the necessary and sufficient conditions

\[\mathcal{D}\geq 0\quad\Longleftrightarrow\quad\sigma=\frac{\partial\psi}{ \partial\epsilon}\quad\wedge\quad\mathbf{\tau}\cdot\mathbf{\dot{\xi}}=\mathcal{D}\geq 0 \tag{3}\]

arise. Thus, the evolution equations for the internal variable are yet to be defined such that \(\mathcal{D}\geq 0\) is ensured at all times.

#### 2.1.2 Generalized standard materials

A common way to formulate the necessary evolution equations is to use the concept of _generalized standard materials_[50, 54, 55] which is briefly explained in the following. Within this concept, in addition to \(\psi\), a _dissipation potential_\(\phi:=\phi(\mathbf{\xi},\mathbf{\xi},\varepsilon)\) which is defined to be (i) _convex_ with respect to its first argument \(\mathbf{\dot{\xi}}\) and is additionally normalized with respect to \(\mathbf{\xi}\), i.e., it fulfills the conditions (ii) \(\phi(\mathbf{0},\mathbf{\zeta},\varepsilon)=0\) and (iii) \(\phi(\mathbf{0},\mathbf{\zeta},\varepsilon)\geq 0\) is a minimum in \(\mathbf{\dot{\xi}}\), is introduced. The dissipation potential may be non-smooth for rate independent, i. e. elastoplastic, materials. The internal forces can now be determined from \(\mathbf{\phi}\) according to \(\mathbf{\tau}\in\partial_{\xi}\phi\), where the operator \(\partial_{\xi}(\cdot)\) denotes the subdifferential of a non smooth convex function. On the other hand, if \(\phi\) is smooth, the relation changes to \(\mathbf{\tau}=\partial_{\xi}\phi\), where the introduced operator now represents the standard partial derivative. Thus, using Eq. (2), it follows the Biot equation which describes the internal variables' evolution and is given by

\[\mathbf{0}\in\partial_{\xi}\psi+\partial_{\xi}\phi\quad\text{or}\quad\mathbf{0}= \partial_{\xi}\psi+\partial_{\xi}\phi\quad\text{with}\quad\mathbf{\xi}(t=0)=\mathbf{ \xi}_{0} \tag{4}\]

for rate independent and rate dependent constitutive behavior, respectively. Note that the inequality given in Eq. (3), i. e., \(\mathcal{D}\geq 0\), is automatically fulfilled by Eq. (4) due to the stated requirements on \(\phi\), i. e., convexity and normalization. An alternative formulation follows with the _dual dissipation potential_\(\phi^{*}\) obtained by the Legendre-Fenchel transformation

\[\phi^{*}(\mathbf{\tau},\mathbf{\xi},\varepsilon):=\sup_{\mathbf{\dot{\xi}}}\left[\mathbf{\tau }\cdot\mathbf{\dot{\xi}}-\phi(\mathbf{\dot{\xi}},\mathbf{\xi},\varepsilon)\right]. \tag{5}\]

By using the dual dissipation one gets

\[\mathbf{\dot{\xi}}\in\partial_{\mathbf{\tau}}\phi^{*}\quad\text{or}\quad\mathbf{\dot{\xi }}=\partial_{\mathbf{\tau}}\phi^{*}\quad\text{with}\quad\mathbf{\xi}(t=0)=\mathbf{\xi}_{0}\quad, \tag{6}\]

respectively, instead of Eq. (4).

**Remark 1** It should be noted, that Eqs. (4) and (6) in this stated form only hold if the rates of the internal variables or the internal forces, respectively, are independent from each other. If there are constraints between the individual quantities, these must be explicitly taken into account during the evaluation of Eqs. (4) and (6). Otherwise, the dissipation potential might for example be expressed in terms of a reduced set of internal variable rates and the partial derivatives with respect to the omitted rates yield zero. This is the case, e.g., for the elastoplastic model as formulated in Tab. 1.

Within the special case of _rate-independent_ constitutive behavior, the dissipation function \(\phi\) is now obtained by using the _concept of maximum dissipation_[50, 55]. Thus, it is defined by the constrained maximization problem

\[\phi(\mathbf{\dot{\xi}},\mathbf{\xi},\varepsilon):=\sup_{\mathbf{\tau}\in\mathcal{E}}\left( \mathbf{\tau}\cdot\mathbf{\dot{\xi}}\right)\quad\text{with}\quad\mathcal{E}:=\left\{ \mathbf{\tau}\in\mathbb{R}^{N}\ |\ f(\mathbf{\tau},\mathbf{\xi})\leq 0\right\}. \tag{7}\]

In the equation above, \(\mathcal{E}\) denotes the admissible domain of internal forces with the yield function \(f(\mathbf{\tau},\mathbf{\xi})\) which is assumed to be convex with respect to \(\mathbf{\tau}\), normalized and homogeneous of degree one. The solution of Eq. (7) yields the evolution equations of internal variables together with the Karush-Kuhn-Tucker conditions:

\[\dot{\mathbf{\xi}}=\lambda\partial_{\mathbf{\tau}}f\ \wedge\ \lambda\geq 0\ \wedge\ f\leq 0\ \wedge\ \lambda f=0. \tag{8}\]

Therein, the scalar \(\lambda\in\mathbb{R}_{+}\) denotes the plastic multiplier.

### Specific constitutive models

The outlined framework given in Sect. 2.1 can be applied to describe, e.g., viscoelasticity or elastoplasticity, for which the corresponding models are briefly summarized in the following.

#### 2.2.1 Viscoelasticity

Viscoelastic behavior of solids is characterized by a strain rate dependent stress response with elastic equilibrium curve, see Fig. 1c. These properties can be described by a _generalized Maxwell model_[53] as shown in Fig. 1a, where \(E\) and \(E_{a}\) are the Young's moduli of the springs within the rheological model, \(\eta_{a}\) the respective viscosities, \(\epsilon_{a}^{\text{el}}\) and \(\epsilon_{a}^{\text{vi}}\) the elastic and viscous strains, and \(\sigma_{a}^{\text{ov}}\) the non-equilibrium stresses also denoted as overstresses. The viscosity \(\eta_{a}\) may be a function of the overstress \(\sigma_{a}^{\text{ov}}\), see Fig. 1b. A model with overstress dependent viscosities is referred to as nonlinear in the following and linear otherwise. The behaviour of such a generalized Maxwell model is described by the governing equations in Tab. 1.
#### 2.2.2 Elastoplasticity

In contrast to viscoelasticity, elastoplastic material behavior is characterized by strain rate independence but the presence of an equilibrium hysteresis [53]. Within an initial region \(\sigma\in\left(-\sigma_{y_{0}},\,\sigma_{y_{0}}\right)\), the material acts purely elastic. If the limits of this initial region are exceeded, i.e., \(|\sigma|>\sigma_{y_{0}}\), elastic-plastic deformation occurs and the elastic region can be both, shifted in the respective direction (_kinematic hardening_) as well as expanded (_isotropic hardening_) [56]. In order to model these properties, a so-called yield function \(f:=f\left(\tau,\xi\right)\leq 0\) is defined. It determines the boundaries of the elastic region such that the deformation is purely elastic if \(f<0\) and elastic-plastic if \(f=0\). The case \(f>0\) is not admissible. The evolution equations can be obtained from the maximum dissipation principle, see Sect. 2.1. Note, that the rates of the internal variables are not independent from each other, which has to be taken into account when evaluating Eq. (4), as explained in Remark 1. With that, the governing equations of elastoplastic behavior with mixed kinematic and isotropic hardening are derived, cf. Tab. 1. The corresponding rheological model is shown in 2a. Therein, \(E\) is the Young's modulus, \(H\) and \(\hat{H}\) are the kinematic and isotropic hardening modulus, respectively, and \(\sigma_{y_{0}}\) is the initial yield stress. Moreover, \(p\) denotes the back stress and \(\epsilon^{\text{el}}\) and \(\epsilon^{\text{pl}}\) the elastic and plastic strain.

Setting the hardening modules \(H\) or \(\hat{H}\) to zero results in the special cases with ideal plasticity (\(H=\hat{H}=0\)), only kinematic hardening (\(\hat{H}=0\)) or only isotropic hardening (\(H=0\)) see Figs. 2b-2e. The corresponding internal variables of the hardening modules set to zero then still take on values unequal to zero, but no longer have any influence on \(\sigma\), \(\psi\) or \(\mathcal{D}\) and can therefore be neglected. For instance, in the case of ideal plasticity \(H=\hat{H}=0\), the set of internal variables can be reduced to the plastic strain \(\epsilon^{\text{pl}}\).

\begin{table}
\begin{tabular}{l l l} \hline  & **Viscoelasticity** & **Elastoplasticity** \\
**Free energy**: & \(\psi:=\frac{1}{2}E\epsilon^{2}+\sum_{s=1}^{N}\frac{1}{2}E_{a}\left(\epsilon- \epsilon_{a}^{\text{si}}\right)^{2}\) & \(\psi:=\frac{1}{2}E\left(\epsilon-\epsilon^{\text{pl}}\right)^{2}+\frac{1}{2}Ha ^{2}+\frac{1}{2}\hat{H}\dot{a}^{2}\) \\
**Dissipation potential**: & \(\phi=\sum_{s=1}^{N}\frac{\gamma_{e}}{2}\left(\dot{a}_{a}^{\text{si}}\right)^{2}\) & \(\phi=\sigma_{y_{0}}\left|\dot{a}^{\text{pl}}\right|\) with \(\dot{a}=\dot{a}^{\text{pl}}\) and \(\dot{\dot{a}}=\left|\dot{a}^{\text{pl}}\right|\) \\
**Dual dissipation potential**: & \(\phi^{*}=\sum_{s=1}^{N}\frac{1}{2\gamma_{e}}\left(\sigma_{s}^{\text{ox}}\right)^ {2}\) & \(\phi^{*}=\left\{0,\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad \text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad \text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad \text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad \text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if }\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad \text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if }\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad \text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\quad \text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad \text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\quad\text{if}\quad\text{if}\quad\quad\text{if}\quad \text{if}\quad\text{if}\quad\text{if}\quad\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad \text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad \text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad \text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if} \quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad \text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\quad \text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\quad\text{if}\quad\text{if}\quad\text{if}\quad\text{if}\quad\quad\text
## 3 Basics of artificial neural networks

NNs can be divided into different classes, each adapted to different types of tasks [57, 58]. For constitutive modeling, the use of FNNs and RNNs is particularly suitable.

### Feedforward neural networks

FNNs are composed of artificial neurons arranged in several layers as shown in Fig. 2(a). The information of the input variables (features) \(i_{n}\), summarised in the input \(\mathbf{I}:=(i_{1},...,i_{N_{1}})\), is passed on to each neuron of the subsequent layer 2 and is used in these to calculate the outputs of this layer. The entire output of layer 2 in turn serves as input for each neuron in the next layer and so forth. The forward flow of information from the input layer to the final output \(\mathbf{O}:=(\,o_{1},\dots,o_{N_{L}})\) of layer \(L\) without feedback of intermediate results is called _feedforward_. The \(n\)th neuron of layer \(l\) is shown in Fig. 2(b). The output \(\mathcal{N}_{n}^{l}\) of this neuron is determined via

\[\mathcal{N}_{n}^{l}=\mathcal{A}_{n}^{l}\left(\chi_{n}^{l}\right)=\mathcal{A}_ {n}^{l}\left(\sum_{m=0}^{N_{l-1}}\mathcal{N}_{m}^{l-1}\mathbf{\omega}_{nm}^{l} \right)\quad. \tag{9}\]

Therein, \(\mathcal{A}_{n}^{l}\) denotes the activation function of this neuron, \(N_{l-1}\) the number of neurons in layer \(l-1\) and \(\mathbf{\omega}_{nm}^{l}\) the weight of the output \(\mathcal{N}_{m}^{l-1}\) in the argument of the activation function \(\mathcal{A}_{n}^{l}\). The weight \(\mathbf{\omega}_{n0}^{l}\) of the output of the imaginary neuron \(\mathcal{N}_{0}^{l-1}\equiv 1\) is referred to as the bias of the neuron. Common activation functions are the hyperbolic tangent \(\tanh(x)\), the rectifier \(\operatorname{ReLU}(x):=\max\left(0,x\right)\) or the softplus activation \(\operatorname{SP}(x):=\ln(1+\exp(x))\)[11, 57, 58].

Particularly, the outputs \(\mathcal{N}_{n}^{L}\) of the output layer \(L\) can also be expressed with (9). Recursively, the \(\mathcal{N}_{n}^{l}\) can be replaced until layer 1 with its inputs \(\mathcal{N}_{n}^{1}=i_{n}\) is reached. This clarifies that such a network may be arbitrarily nested, but for fixed \(\mathbf{\omega}_{nm}^{l}\) is a well-defined function of the input \(\mathbf{I}\). In this sense, the network as a function

\[F:\mathbb{R}^{N_{1}}\rightarrow\mathbb{R}^{N_{L}}\quad,\quad\mathbf{I}\mapsto \mathcal{F}(\mathbf{I})=\mathbf{O} \tag{10}\]

Figure 3: Structure of an FNN: **(a)** General representation of a full network with arbitrary number of layers and neurons and **(b)** functionality of the \(n\)th neuron of layer \(l\).

Figure 2: Elastoplastic constitutive model: **(a)** Rheological model with isotropic and kinematic hardening and typical hystereses for **(b)** ideal plasticity (\(\mathbf{H}=\hat{\mathbf{H}}=0\)), **(c)** kinematic hardening (\(\hat{\mathbf{H}}=0\)), **(d)** isotropic hardening (\(\mathbf{H}=0\)) and **(e)** kinematic and isotropic hardening.
maps the input \(\mathbf{I}\) to the output \(\mathbf{O}\). In order to adapt the weights \(w_{nm}^{l}\) such that the predictions of the NN match the expected values, a training data set \(S:=\left\{\mathcal{T}_{1},\ldots,\mathcal{T}_{N^{\text{th}}}\right\}\) consisting of \(N^{\text{ds}}\in\mathbb{N}\) data tuples \(\mathcal{T}_{a}\) is required. In each of these data tuples \(\mathcal{T}_{a}:=\left(\mathbf{I}_{a}\,,\,\mathbf{O}_{a}\right)\), an Input \(\mathbf{I}_{a}\) is assigned its expected output \(\mathbf{\tilde{O}}_{a}\). The expected output usually contains the desired values of the neurons in the last layer, but may also contain additional information, such as desired derivatives with respect to a certain input. The error of the predictions of an NN with respect to this set of data tuples is summarized in the loss function \(\mathcal{L}\left(\mathbf{w},\mathcal{S}\right)\), where \(\mathbf{w}\) denotes the vector of all weights \(w_{nm}^{l}\) of the NN. For the applications shown here, the loss function is mostly composed of mean absolute error terms, which for FNN read

\[\text{MAE}\left(q\right):=\frac{1}{N^{\text{ds}}}\sum_{i=1}^{N^{\text{ds}}} \left|q(\mathbf{I}_{i},\mathbf{w})-\bar{q}_{i}\right|\quad. \tag{11}\]

Therein, \(q\) is the quantity, whose prediction is evaluated. The minimization of this error function with respect to the weights \(\mathbf{w}\)

\[\mathbf{w}=\arg\min_{\mathbf{k}}\mathcal{L}\left(\mathbf{x},\mathcal{S}\right) \tag{12}\]

is called the training process of the NN. Various methods can be used to solve this minimization problem, e.g. _Stochastic Gradient Descent_ (SGD) or _Adam_, to name two of the most common optimizers. Within this work, _Sequential Least Squares Programming_ (SLSQP) is used for the optimization.

A special class of FNNs are ICNNs initially proposed by Amos et al. [51], i.e., networks that have the property to be convex with respect to their input arguments. This is achieved by using a convex and non-decreasing activation function and non-negative weights.

### Recurrent neural networks

In contrast to FNNs, RNNs allow intermediate results to be fed back, see Fig. 3(a). This feature enables the network to take into account not only one, but several time steps for the output and thus to learn to interpret the input variables in their (pseudo) temporal context. This can be particularly useful, for example, in text translation, speech recognition or prediction of the evolution of physical quantities. Thus, in RNNs, entire sequences of successive data are evaluated, which in the applications shown within this work are sets of time steps. RNNs are thus very different from FNNs in the way they operate. Nevertheless, their internal structure is also composed of a single or multiple FNNs.

These FNNs receive the current input vector \({}^{n}\mathbf{I}\) and the so called hidden state \({}^{n-1}\mathbf{h}\) of the previous time step and output the hidden state \({}^{n}\mathbf{h}\) of the current time step. Outside of the recurrent cell, the new hidden state is subsequently processed further to obtain the output \({}^{n}\mathbf{O}\). Fig. 3(b) illustrates this process for a simple RNN cell unfolded in time.

The depicted cell consists of only a single FNN and is not capable of incorporating many preceding time steps into the output of the current hidden state. Therefore, more complex internal structures have been developed using multiple FNNs in order to enable the cell to detect long term dependencies. The RNNs presented in this work use the so-called _long short-term memory_ (LSTM) [27] cell to overcome this problem. Besides the hidden state, this LSTM cell passes another set of information, the so-called cell state \({}^{n}\mathbf{c}\), from time step to time step. The number \(N^{\text{c}}\) of entries in \({}^{n}\mathbf{c}\) can be interpreted as the memory capacity of the cell and is a crucial parameter for the network performance. In contrast to FNNs, the training data set \(S:=\left\{\mathcal{T}_{1},\ldots,\mathcal{T}_{N^{\text{sqg}}}\right\}\) contains \(N^{\text{sqg}}\in\mathbb{N}\) tuples \(\mathcal{T}_{s}:=(\mathcal{T}_{s}^{1},\tilde{\mathcal{T}}_{s}^{0})\) of input sequences \(\mathcal{T}_{s}^{1}\) and assigned expected output sequences \(\tilde{\mathcal{T}}_{s}^{0}\). Each input sequence \(\mathcal{T}_{s}^{1}:=(^{1}\mathbf{I}_{s},\ldots,^{N^{\text{us}}}\mathbf{I}_{s})\) is an ordered series of \(N^{\text{ts}}_{s}\in\mathbb{N}\) input vectors, each representing one of the \(N^{\text{ts}}_{s}\) time steps of the sequence. The expected output sequence \(\tilde{\mathcal{T}}_{s}^{0}:=(^{1}\mathbf{\tilde{O}}_{s},\ldots,^{N^{\text{ts}}} \mathbf{\tilde{O}}_{s})\)

Figure 4: Functionality and internal structure of a standard RNN cell in two equivalent illustrations: the output is **(a)** fed back into the RNN cell and **(b)** passed from one time step to the next. I.e., (b) can be understood as (a) unfolded in time.
contains the same number of expected output vectors.

The error of the RNNs predictions is measured in the loss function \(\mathcal{L}\left(\mathbf{w},\mathcal{S}\right)\), where \(\mathbf{w}\) now contains the weights of all FNNs inside and outside the recurrent cell. For RNN, the mean absolute error measure of a quantity \(q\) takes the form

\[\text{MAE}\left(q\right)\,:=\frac{1}{N^{\text{seq}}}\sum_{s=1}^{N^{\text{seq}}} \left(\frac{1}{N_{s}^{\text{ts}}}\sum_{n=1}^{N_{s}^{\text{ts}}}\left|q(T_{s}^{ \text{t}},\mathbf{w})\,-\,^{n}\bar{q}_{s}\right|\right)\quad. \tag{13}\]

The training process, similar to the FNNs, is performed as minimization of the loss function with respect to all weights present in the network.

## 4 NN-based constitutive models

Although the models presented in Sect. 2 are able to reproduce essential aspects of real material behavior in a physically consistent manner, the applied approach to find a set of governing equations has some disadvantages. On the one hand, this set of equations has to be formulated manually, which becomes increasingly difficult for more complex materials in three dimensions and, on the other hand, it can only reproduce the phenomena taken into account during modeling. In this section, some possible methods are presented that, based on NNs, are intended to find arbitrary correlations automatically, i.e., without any additional manual modeling effort. These methods are divided into three categories, which describe the amount of incorporated physics: (i) the _black box models_ with no physical information, (ii) _NNs enforcing physics in a weak form_, which are trained to satisfy the second law and (iii) _NN-based models enforcing physics in a strong form_, in which the second law is satisfied a priori. Each method requires a training data set, whose generation is outlined in Sect. 5. The implementation is done in Python using the Tensorflow library, with the results being presented in Sect. 6.

### Black box models

Black box models are characterized by their lack of incorporated physical knowledge, i.e., they take a set of input values (depending on the concrete model) and output the desired quantity (the stress) without respecting relationships known from continuum mechanics. Two simple architectures are examined, one using a feedforward network, the other using a recurrent cell.

#### 4.1.1 Feedforward architecture with stress as output (FNN\({}^{\sigma}\))

This architecture is the simplest of the six architectures presented in this work and consists of only a single FNN, which outputs the stress directly. It receives the new strain \({}^{n+1}\epsilon\) as well as strains \({}^{m}\epsilon\), stresses \({}^{m}\sigma\) and time increments \({}^{m}\Delta t\,:=\,^{m+1}t\,-\,^{m}t\) of an arbitrary number of \(N^{\text{pt}}\geq 1\) preceding time steps and outputs the new stress, that is

\[\mathbf{I}\,:=\,(^{n+1}\epsilon,^{n}\epsilon,^{n}\sigma,^{n}\Delta t,\ldots,\,^{n +1-N^{\text{pt}}}\epsilon,^{n+1-N^{\text{pt}}}\sigma,^{n+1-N^{\text{pt}}} \Delta t)\mapsto\mathcal{O}\,:=\,^{n+1}\sigma\quad. \tag{14}\]

The time steps \({}^{m}\Delta t\) in the input are only necessary for rate-dependent behavior, i.e., for viscoelasticity. The loss function contains only the stress prediction,

\[\mathcal{L}\,:=\,\mathcal{L}^{\sigma}\,:=\,\text{MAE}(\sigma)\quad. \tag{15}\]

Similar approaches, e.g., the nested adaptive neural networks for constitutive modeling, have been described earlier [23, 24, 26]. However, as the results in Sect. 6 show, this method can only be applied to specific classes of material behavior. An improved architecture based on a recurrent cell can be used to overcome this restriction.

#### 4.1.2 Recurrent architecture with stress as output (RNN\({}^{\sigma}\))

The second black box architecture consists of a _recurrent cell_ followed by an FNN as shown in Fig. 6. Such models are used in the works [29]. For each time step in every sequence, the RNN cell takes as input the new strain and if necessary the time increment and outputs the new hidden state, which is subsequently fed into the FNN to receive the new stress,

\[{}^{n}\mathbf{I}^{\text{RNN}}\,:=\,(^{n}\epsilon,^{n-1}\Delta t)\mapsto{}^{n}\mathbf{ h}=:\,^{n}\mathbf{I}^{\text{FNN}}\mapsto{}^{0}\mathbf{O}^{\text{FNN}}\,:=\,^{n} \sigma\quad. \tag{16}\]

The initial state before the first time step of hidden state as well as cell state are set to \({}^{0}\mathbf{h}=\mathbf{0}\) and \({}^{0}\mathbf{c}=\mathbf{0}\). The loss function

\[\mathcal{L}\,:=\,\mathcal{L}^{\sigma}\,:=\,\text{MAE}(\sigma) \tag{17}\]

averages the error in the stress prediction over all time steps and sequences.

This architecture is capable of modeling a broad variety of materials. Nevertheless, neither is any physical knowledge provided to the network nor can any additional physical information be obtained besides the stress response. Therefore, two physically informed models are presented below.
### Neural networks enforcing physics in a weak form

In contrast to black box models, _NN-based models enforcing physics in a weak form_ take into account known relations, namely Eq. (3)\({}_{2}\) and Eq. (3)\({}_{3}\). This aims at making predictions that are actually physically consistent and gaining additional information, particularly the free energy and the dissipation rate. The networks are trained to respect these relations, but the architecture itself cannot guarantee \(D\geq 0\) a priori. To incorporate the relations into the training process, the network must be differentiated with respect to its inputs. The differentiation of a network results in restrictions regarding the choice of activation functions: since gradients of the loss function are required during the optimization, the activation functions are differentiated a second time. Consequently, the activations must be differentiable twice to enable an optimization [10, 11, 47]. Therefore, the hyperbolic tangent and the softplus activation are used for the corresponding networks in the following.

2.1 Feedforward architecture with internal variables and free energy as output (FNN\({}^{\sharp+\psi}\))

The architecture described below is based on Masi et al. [11] and is studied in detail therein. It uses two FNNs with different tasks, as shown in Fig. 7. The first subnetwork sFNN\({}^{\sharp}\) receives the new and old strain, old stress as well as the time increment and internal variables of the previous time step and is trained to predict new internal variables:

\[\mathbf{I}^{\text{FNN}^{\sharp}}:=(^{n+1}\epsilon,{}^{n}\epsilon,{}^{n}\sigma,{}^ {n}\Delta t,{}^{n}\xi^{1},\dots,{}^{n}\xi^{N})\mapsto\mathbf{O}^{\text{FNN}^{ \sharp}}:={}^{n+1}\mathbf{\xi}\quad. \tag{18}\]

The following second subnetwork sFNN\({}^{\varphi}\) uses the output of the first network alongside the new strain to determine the free energy, i.e.,

\[\mathbf{I}^{\text{FNN}^{\varphi}}:=(^{n+1}\epsilon,{}^{n+1}\xi^{1},\dots,{}^{n+1} \xi^{N})\mapsto\mathbf{O}^{\text{sFNN}^{\varphi}}:={}^{n+1}\psi\quad. \tag{19}\]

Since sFNN\({}^{\sharp}\) takes \({}^{n}\mathbf{\xi}\) as input, the availability of internal variables in the training data set is essential. The free energy is not required but improves the network performance. Within this study, the goal is to use as few information as possible, thus considering the free energy as not available during training. This applies to all of the following architectures as well.

In contrast to the original work [11], the training process is not performed as a single optimization considering both networks, but is split into two smaller steps, speeding up the optimization: First, sFNN\({}^{\sharp}\) is trained independently, i.e., detached from sFNN\({}^{\varphi}\), using

\[\mathcal{L}^{\text{sFNN}^{\sharp}}:=\mathcal{L}^{\xi}\quad\text{with}\quad \mathcal{L}^{\xi}:=\frac{1}{N}\sum_{\alpha=1}^{N}\text{MAE}\left(\xi^{\alpha}\right) \tag{20}\]

to learn to predict the internal variables. Subsequently, both networks are used together in order to train the second network via the loss function

\[\mathcal{L}^{\text{sFNN}^{\varphi}}:=w^{\sigma}\mathcal{L}^{\sigma}+w^{\psi} \mathcal{L}^{\psi}+w^{D\geq 0}\mathcal{L}^{D\geq 0}\quad\text{with}\quad \mathcal{L}^{\sigma}:=\text{MAE}\left(\sigma\right)\quad, \tag{21}\]

\[\mathcal{L}^{\varphi}:=\text{MAE}\left(\psi\right)\quad\text{and}\quad \mathcal{L}^{D\geq 0}:=\sum_{k=1}^{N^{\text{old}}}\text{ReLU}\left(-{}^{n+1}D_{k} \right)\quad\text{where}\quad{}^{n+1}D=-\sum_{\alpha=1}^{N}\frac{\partial^{n+ 1}\psi}{\partial^{n+1}\xi^{\alpha}}\frac{{}^{n+1}\xi^{\alpha}-{}^{n}\xi^{ \alpha}}{{}^{n}\Delta t}\quad. \tag{22}\]

The non-trainable parameters \(w^{\sigma}\), \(w^{\psi}\) and \(w^{D\geq 0}\) regulate the influence of each term on the value of the loss function. During this optimization, the loss is only differentiated with respect to the weights of sFNN\({}^{\varphi}\), so that only these weights
are adjusted and the weights of sFNN\({}^{\xi}\) retain their values as obtained from the first training step.

However, the provision of internal variables for the training data for more complex material behavior is not a trivial task, which is addressed in Masi and Stefanou [47]. In order to avoid providing internal variables, a similar architecture seems appropriate, in which sFNN\({}^{\xi}\) is replaced by a recurrent cell.

.2.2 Recurrent architecture with internal variables and free energy as output (RNN\({}^{\xi+\psi}\))

Since RNNs are capable of storing information from time steps far in the past, it seems reasonable to exploit this property to mimic internal variables. Fig. 8 shows a possible implementation of this idea, which is similar to He and Chen [48]. Therein, at each time step, the RNN cell takes the new strain as argument as well as the time increment in case of viscous behavior and produces a new vector-valued output \({}^{n}\mathbf{h}\), which carries the history information and is initialized with \({}^{0}\mathbf{h}={}^{0}\mathbf{c}=\mathbf{0}\). Subsequently, these values are fed into an FNN and are reduced to freely selectable number of internal variables \(\mathbf{N}^{\xi}\):

\[{}^{n}\mathbf{I}^{\text{RNN}}\ :=({}^{n}\varepsilon,^{-1}\Delta t)\mapsto{}^{n}\mathbf{h}= \cdot{}^{n}\mathbf{I}^{\text{FF}^{\xi}}\mapsto{}^{n}\mathbf{O}^{\text{FF}^{\xi}}:={}^{ n}\mathbf{\xi}\in\mathbb{R}^{N^{\xi}}\quad. \tag{23}\]

Within another FNN, these values along with the new strain are projected onto the free energy, i.e.,

\[\mathbf{I}^{\text{FF}^{\psi}}\ :=({}^{n}\varepsilon,^{n}\underline{\varepsilon}^{1}, \ldots,{}^{n}\underline{\varepsilon}^{N^{\xi}})\mapsto O^{\text{FF}^{\psi}}\ :={}^{n}\psi\quad. \tag{24}\]

The loss function

\[\mathcal{L}^{\text{FF}^{\psi}}\ :=w^{\sigma}\mathcal{L}^{\sigma}+w^{\xi}\mathcal{L}^{\xi}+w^{ \psi}\mathcal{L}^{\psi}+w^{D\geq 0}\mathcal{L}^{D\geq 0}\quad\text{with} \tag{25}\]

\[\mathcal{L}^{\sigma}\ :=\text{MAE}\left(\sigma\right)\quad,\quad\mathcal{L}^{\xi}\ :=\frac{1}{N_{\text{av}}}\sum_{a=1}^{N_{\text{av}}}\text{MAE}\left( \underline{\varepsilon}_{\text{av}}^{\sigma}\right)\quad,\quad\mathcal{L}^{ \psi}\ :=\text{MAE}\left(\psi\right)\quad\text{and} \tag{26}\]

\[\mathcal{L}^{D\geq 0}\ :=\frac{1}{N^{\text{seq}}}\sum_{s=1}^{N^{\text{seq}}} \left(\frac{1}{N^{\text{s}}_{s}}\sum_{n=1}^{N^{\text{s}}_{s}}\text{ReLU} \left(-{}^{n}D_{s}\right)\right)\quad\text{where}\quad{}^{n}\mathcal{D}=-\sum _{a=1}^{N}\frac{\partial^{n}\psi}{\partial^{n}\xi^{a}}\frac{{}^{n}\xi^{a}-{}^{ n-1}\underline{\varepsilon}^{a}}{{}^{n}\Delta t}\quad, \tag{27}\]

is composed of four terms that control the prediction of stress, internal variables and free energy as well as the compliance with thermodynamic consistency, i.e., \(D\geq 0\). In the error term for the internal variables \(\mathcal{L}^{\xi}\), only those internal variables

Figure 8: Functionality of the RNN\({}^{\xi+\psi}\) architecture: For each time step in a sequence, new strain and time increment are fed into the RNN cell, whose output is passed on to the feedforward network sFNN\({}^{\xi}\) to obtain the new set of internal variables. Together with the new strain, these serve as input in a final feedforward network sFNN\({}^{\psi}\) which predicts the new free energy. Differentiating sFNN\({}^{\psi}\) yields the new stress.

Figure 7: Functionality of the FNN\({}^{\xi+\psi}\) architecture: Based on the input values, the first subnetwork sFNN\({}^{\xi}\) predicts a new set of internal variables, which alongside the new strain is passed on to the second subnetwork sFNN\({}^{\psi}\) to predict the new free energy. Differentiating with respect to the new strain yields the new stress.
can be considered that are actually available in the training data set, in Eq. (26) denoted as \(\xi_{\text{av}}\) for \(N_{\text{av}}\) available internal variables. This might be the measurable plastic strain, for example. Consequently, if no internal variables are known, this term is omitted and the network has to find reasonable representations of internal variables on its own. The free energy term \(w^{\prime\prime}\mathcal{L}^{\psi}\) is also optional and can be omitted if the free energy is unknown.

### Neural networks enforcing physics in a strong form

The last category of models, denoted as _NNs enforcing physics in a strong form_ here, in contrast to the previous model class, satisfies \(D\geq 0\) a priori by construction of the model. To achieve this, the concept of generalized standard materials summarized in Sect. 2.1.2 is adopted into the data-driven paradigm by applying ICNNs [51]. This has been proposed for viscoelasticity by Huang et al. [49] as well as As'ad and Farhat [7]. In the following, three models, FNN\({}^{w+\phi}\), FNN\({}^{w+\phi^{*}}\), and FNN\({}^{w+\phi+\xi}\), are introduced in detail.

**Remark 2**.: It should be noted that an application of the introduced approaches FNN\({}^{w+\phi}\) and FNN\({}^{w+\phi^{*}}\) to data belonging to a rate-independent material automatically leads to a regularization [59], i.e., an approximation of the data by a rate-dependent model. This is due to the fact that the chosen activation function cannot represent the non-smooth dissipation potentials typical for plasticity, cf. Tab. 1.

#### 4.3.1 FNN-architecture with free energy and dissipation potential as output (FNN\({}^{w+\phi}\))

The first architecture of this category, FNN\({}^{w+\phi}\), is composed of two FNNs, modeling the free energy \(\psi\) and the dissipation potential \(\phi\). With the requirements on \(\phi\) from Sect. 2.1, i.e., (i) convexity in \(\xi\), (ii) \(\phi(\xi=0,\xi,\epsilon)=0\), and (iii) \(\phi(\xi=0,\xi,\epsilon)\) is a minimum in \(\xi\), this model can be used to make predictions that are _a priori thermodynamically consistent_. Convexity of the dissipation potential with respect to only \(\xi\), but not \(\xi\) and \(\epsilon\), is achieved through the multiplicative split

\[\phi^{\text{NN}}(\xi,\xi,\epsilon):=\phi^{\text{con}}(\xi)\,\phi^{+}(\epsilon, \xi) \tag{28}\]

into a convex part \(\phi^{\text{con}}(\xi)\) depending on only the rate of the internal variables and a positive part \(\phi^{+}(\epsilon,\xi)\) depending on the strain and the internal variables themselves. Each part is modeled by a single FNN. For \(\phi^{\text{con}}\) to be convex, it requires convex and non-decreasing activation functions, here the Softplus activation \(\text{SP}(x):=\ln(1+\exp(x))\), and non-negative weights across the whole network [51, 16, 5]. Positivity of \(\phi^{+}\) in turn requires only the weights and bias of the output layer to be non-negative and positive activation functions, here also the Softplus activation. As another, more general method to construct a network, that is convex in only some of its inputs, a _partially input convex neural network_ (PICNN) as presented by Amos et al. [51] could be used as well. To meet normalization conditions (ii) and (iii), the output \(\phi^{\text{NN}}\) of the combined network given in Eq. (28) is modified via

\[\phi(\xi,\xi,\epsilon):=\phi^{\text{NN}}(\xi,\xi,\epsilon)-\phi^{\text{NN}}( \xi=0,\xi,\epsilon)-\frac{\partial\phi^{\text{NN}}(\xi=0,\xi,\epsilon)}{ \partial\xi}\cdot\xi\geq 0\ \forall\xi,\xi,\epsilon\quad. \tag{29}\]

Likewise, the ansatz

\[\psi(\epsilon,\xi):=\psi^{\text{NN}}(\epsilon,\xi)-\psi^{\text{NN}}(\epsilon= 0,\xi=0)-\frac{\partial\psi^{\text{NN}}(\epsilon=0,\xi=0)}{\partial\epsilon} \epsilon-\frac{\partial\psi^{\text{NN}}(\epsilon=0,\xi=0)}{\partial\xi}\cdot\xi \tag{30}\]

guarantees that free energy, stress and internal forces equal zero in the initial, unloaded state [49].

Now, after the model formulation given above, it is explained how the training algorithm and the prediction of new time steps with a calibrated FNN\({}^{w+\phi}\) architecture work. Thereby, tuples of stress \(\sigma\), strain \(\epsilon\), and internal variables \(\xi\) are needed for the training of \(\phi\) and \(\psi\) with Eqs. (29) and (30).2

Footnote 2: To be able to use only stress and strain for training, it is also possible to integrate the full path of \(\xi\) in each training epoch. However, since this requires backpropagation over all time steps, the training becomes much more time consuming [52]. Alternatively, following As’ad and Farhat [7], a further network for the approximation of internal variables \(\xi\) during training can be applied. This also enables to end up with a training data set containing only tuples of \(\sigma\) and \(\epsilon\), whereby the training time is reduced compared to the first method. Another possibility is to determine the internal variables in advance within a preprocessing step. How this can be done is described in Ladeveze et al. [60] or Gerbaud et al. [61]. Here, the second method is used, see Sect. 4.3.3.

Training of the three networks \(\phi^{\text{con}}(\xi)\), \(\phi^{+}(\epsilon,\xi)\), and \(\psi^{\text{NN}}(\epsilon,\xi)\) is performed as follows: The training starts by predicting \({}^{w}\psi\) with Eq. (30) and obtaining \({}^{n}\sigma\) and \({}^{n}\sigma^{4}\) via differentiation. Following Eq. (4), the internal forces shall equal the dissipation potential differentiated with respect to the rate of the corresponding internal variable. This rate is approximated with \({}^{n}\xi^{\dot{\alpha}}:=(^{n}\xi^{\alpha}-^{n-1}\xi^{\alpha})/^{n-1}\Delta t\). The difference of the predicted internal forces \({}^{n}\tau^{\alpha}=-\frac{\partial\psi}{\partial^{2}\xi^{\alpha}}\) and
\({}^{n}\hat{\mathbf{\varepsilon}}^{a}:=\frac{\partial\phi}{\partial^{\frac{n}{2}}\phi^{a}}\) as well as the stress prediction are included in the loss function

\[\mathcal{L}:=\omega^{\sigma}\mathcal{L}^{\sigma}+\omega^{\text{Bitt}}\mathcal{L}^ {\text{Biot}}\quad\text{with}\quad\mathcal{L}^{\sigma}:=\text{MAE}\left(\sigma \right)\quad\text{and}\quad\mathcal{L}^{\text{Biot}}:=\frac{1}{N}\sum_{a=1}^{N }\text{MAE}\left(\tau^{a}\right)\quad, \tag{31}\]

which is to minimize in the process.

Once training is finished and good representations of \(\psi\) and \(\phi\) are found, these potentials can be used to predict the material response for a given strain path. Therefore, in each time step \(n\), the rates of the internal variables \({}^{n}\hat{\mathbf{\varepsilon}}^{a}\) are adapted iteratively using a _Newton-Raphson scheme_, such that \(\max_{\mathbf{\varepsilon}}|{}^{n}\tau^{a}-{}^{n}\hat{\mathbf{\varepsilon}}^{a}|<e\) with a given tolerance \(\mathbf{\varepsilon}\), where the new internal variables are obtained via \({}^{n}\hat{\mathbf{\varepsilon}}^{a}={}^{n-1}\hat{\mathbf{\varepsilon}}^{a}+{}^{n-1} \Delta{}^{n}\hat{\mathbf{\varepsilon}}^{a}\), see the scheme given in Fig. 9.

3.2 FNN-architecture with free energy and dual dissipation potential as output (FNN\({}^{\psi+\phi^{*}}\))

The second architecture presented here is the FNN\({}^{\psi+\phi^{*}}\) model which makes use of the _dual dissipation potential_\(\phi^{*}(\mathbf{\tau},\mathbf{\xi},\mathbf{\varepsilon})\) according to Eq. (5). The model equations are thus similar to Eqs. (28)-(30) with the difference that \(\phi\) has to be replaced by \(\phi^{*}\). The dual dissipation potential \(\phi^{*}\) has to guarantee (i) convexity in \(\mathbf{\tau}\), (ii) \(\phi^{*}(\mathbf{\tau}=\mathbf{0},\mathbf{\xi},\mathbf{\varepsilon})=0\), and (iii) \(\phi^{*}(\mathbf{\tau}=\mathbf{0},\mathbf{\xi},\mathbf{\varepsilon})\) is a minimum in \(\mathbf{\tau}\), to get an a priori thermodynamically consistent model. Again, convexity of \(\phi^{*}\) with respect to only \(\mathbf{\tau}\) is achieved through a multiplicative split similar to Eq. (28).

Besides the difference in the model equations, the FNN\({}^{\psi+\phi^{*}}\) architecture differs in its training algorithm and the method to predict the new time step. In FNN\({}^{\psi+\phi^{*}}\), training with respect to a data set consisting of tuples of \(\sigma\), \(\mathbf{\varepsilon}\), and \(\mathbf{\xi}\) is performed as follows: First, the free energy \({}^{n}\psi\) is calculated and \({}^{n}\sigma\) and \({}^{n}\tau^{a}\) are obtained via differentiation. Subsequently, the internal forces \({}^{n}\tau^{a}\) are fed into the combined network of the dual dissipation potential and \({}^{n}\phi^{*}\) is received. Differentiating \({}^{n}\phi^{*}\) with respect to the internal forces yields the rate of the internal variables \({}^{n}\hat{\mathbf{\varepsilon}}^{a}\). The loss function

\[\mathcal{L}:=w^{\sigma}\mathcal{L}^{\sigma}+w^{\hat{\mathbf{\varepsilon}}} \mathcal{L}^{\hat{\xi}}\quad\text{with}\quad\mathcal{L}^{\sigma}:=\text{MAE} \left(\sigma\right)\quad\text{and}\quad\mathcal{L}^{\hat{\xi}}:=\frac{1}{N} \sum_{a=1}^{N}\text{MAE}\left(\hat{\xi}^{a}\right) \tag{32}\]

now compares the predictions of the stress and the rates of the internal variables with their expected values \({}^{n}\bar{\sigma}\) and \({}^{n}\hat{\mathbf{\varepsilon}}^{a}\), respectively, where \({}^{n}\hat{\mathbf{\varepsilon}}^{a}=({}^{n}\hat{\mathbf{\varepsilon}}^{a}-{}^{n-1}\hat {\mathbf{\varepsilon}}^{a})/{}^{n-1}\Delta t\). Once the three networks are trained, the material state for a time step can be obtained using the pattern in Fig. 10, where due to the implicit description of \({}^{n}\hat{\mathbf{\varepsilon}}^{a}={}^{n-1}\hat{\mathbf{\varepsilon}}^{a}+{}^{n-1} \Delta t^{n}\hat{\mathbf{\varepsilon}}^{a}\), the evolution of the internal variables is determined iteratively using a Newton-Raphson scheme, such that \(\max_{\mathbf{\varepsilon}}|{}^{n}\hat{\mathbf{\varepsilon}}^{a}={}^{n}\hat{\mathbf{ \varepsilon}}^{a}|<\mathbf{\varepsilon}\) with a given tolerance \(\mathbf{\varepsilon}\).

3.3 Adapted training for FNN-architecture with free energy and dissipation potential as output (FNN\({}^{\psi+\phi+\xi}\))

The training of the former two architectures, i.e., FNN\({}^{\psi+\phi}\) and FNN\({}^{\psi+\phi^{*}}\), requires the internal variables to be present in the training data set. However, since these quantities are usually unknown, As'ad and Farhat [7] proposed an efficient training method for a model based on the dual dissipation potential, in which the internal variables are no longer required. This method is adopted analogously here for the architecture FNN\({}^{\psi+\phi}\).

To no longer have to specify the internal variables directly while keeping the training fast, another FNN, which takes the time \(t\) as input and outputs the internal variables \(\xi^{\text{NN}}(t)\), is integrated only for the training process. The

Figure 9: Functionality of the FNN\({}^{\psi+\phi}\) architecture: With the rate of internal variables of the current iteration, the internal forces are calculated using the network of the free energy (\(\mathbf{\tau}\), left) and the networks of the dissipation potential (\(\hat{\mathbf{\tau}}\), right). The rate of internal variables is now adapted iteratively, such that \(\mathbf{\tau}\approx\hat{\mathbf{\tau}}\).
values for the internal variables and their rates entering the networks for the free energy and the dissipation potential are taken from this auxiliary network, where the rate is approximated consistently to the prediction process by \(\sfrac{\hat{\xi}^{\text{NN}}}{\hat{\xi}}=\frac{1}{\Delta t}\left(\xi^{\text{NN}}( \sfrac{n}{t})-\xi^{\text{NN}}(\sfrac{n-1}{t})\right)\). The condition \(\xi^{\text{NN}}(0)=\mathbf{0}\) is enforced by subtracting \(\xi^{\text{NN}}(0)\) from the actual output. During the training process, the weights of this network are adapted such that the predicted temporal courses of \(\xi^{\text{NN}}(t)\) allow the errors for stress (\(\mathcal{L}^{\sigma}\)) and Biot term (\(\mathcal{L}^{\text{Biot}}\)) to become as small as possible. This leads to a reasonable representation of the internal variables without explicitly specifying their values. Once the training is finished, the auxiliary network is no longer necessary and the prediction process can be carried out in the same manner as described in Fig. 9. To make the training process easier for the optimizer, the courses of \(\xi^{\text{NN}}(t)\) that are to be represented with the auxiliary network should to be smooth and rather simple. This is done by choosing a smooth strain path as explained in Sect. 5.1. As'ad and Farhat [7] further propose an additive split of the free energy into an equilibrium and a dissipative part. In order to retain comparability and consistency with \(\text{FNN}^{\psi+\phi}\), this split is not performed herein. However, it should be noted that this additional assumption significantly facilitates training.

In addition to the loss function Eq. (31), another term \(w^{[\xi]}\mathcal{L}^{[\xi]}\) is added. This term has the purpose to keep the internal variables small. Otherwise, the chosen internal variables become unboundedly large. To illustrate this, consider a constitutive model with free energy function \(\psi=\psi(\epsilon,\xi)=\frac{1}{2}E(\epsilon-\xi)^{2}\) with corresponding stress and internal force \(\sigma=\tau=E(\epsilon-\xi)\). An equivalent expression with another internal variable \(\tilde{\xi}\) is \(\tilde{\psi}(\epsilon,\tilde{\xi})=\frac{1}{2}E(\epsilon-\frac{1}{A}\tilde{\xi})\) where \(\tilde{\xi}=A\xi\) and \(\tilde{\sigma}=E(\epsilon-\frac{1}{A}\tilde{\xi})=\sigma\) and \(\tilde{\tau}=\frac{1}{A}E(\epsilon-\frac{1}{A}\tilde{\xi})\). The loss term \(\mathcal{L}^{\text{Biot}}\) rewards small absolute values of the internal forces, i.e., a representation with large scaling factor \(|A|\), resulting in \(|\tilde{\xi}|\gg|\xi|\). This leads to difficulties in the prediction process later on. The new loss thus reads

\[\mathcal{L}\,:=w^{\sigma}\mathcal{L}^{\sigma}+w^{\text{Biot}}\mathcal{L}^{ \text{Biot}}+w^{[\xi]}\mathcal{L}^{[\xi]}\quad\text{with}\quad\mathcal{L}^{[ \xi]}\,:=\frac{1}{N^{\text{ds}}}\sum_{n=1}^{N^{\text{ds}}}\sum_{\alpha=1}^{N} |\xi_{\alpha}^{\text{NN}}(\sfrac{n}{t})|\quad, \tag{33}\]

where \(\mathcal{L}^{[\text{NN}]}\) penalizes large absolute values of the chosen internal variable. Note that this does not imply, that \(\text{sign}(\tilde{\xi})=\text{sign}(\xi)\), since \(A<0\) is still a valid option. Within the scope of this study, only a single internal variable is predicted. Another network for another internal variable or a network with two outputs could be used as well, but this is not shown herein.

## 5 Generation of the database for training and validation

Prior to studying the usability of the presented material models based on NNs, a database must be generated for the training process and, later on, as a reference for validation. Thereby, four different models according to Sects. 2.2.1 and 2.2.2 are used to generate material states belonging to a prescribed strain path: a _nonlinear viscoelastic model with one Maxwell element_ (V1), a _linear viscoelastic model with two Maxwell elements_ (V2), an _elastoplastic model with kinematic hardening_ (P1), as well as an _elastoplastic model with mixed kinematic-isotropic hardening_ (P2). The models' governing equations given in Tab. 1 are solved by applying an implicit Euler scheme. The chosen parameters of these reference constitutive models can be found in Tab. 2.

Figure 10: Functionality of the \(\text{FNN}^{\psi+\phi^{*}}\) architecture: Using the rate of internal variables of the current iteration, the internal variables alongside the strain are fed into the network of the free energy. Stress and internal forces are obtained via differentiation. The internal forces are passed on to the convex network for \(\phi^{*\text{con}}\), strain and internal variables to the positive network for \(\phi^{*+}\). Multiplying \(\phi^{*\text{con}}\) and \(\phi^{*+}\) yields the dual dissipation potential, which is differentiated with respect to the internal forces to receive \(\dot{\tilde{\xi}}\). The rate of internal variables is now adapted iteratively, such that \(\dot{\tilde{\xi}}\approx\dot{\tilde{\xi}}\).
### Generation of training data

Two different methods to generate training data are used. The first method can be applied to all architectures except for FNN\({}^{\psi+\phi+\xi}\) and uses random walk sequences. The FNN\({}^{\psi+\phi+\xi}\) requires a smooth strain path, which is given by a cubic spline combining a chosen set of knots. Both methods are described in detail below.

#### 5.1.1 Random walk strain paths

Four different data bases are generated, one for each examined material V1, V2, P1 and P2. Since the recurrent architectures require several independent sequences in their training data, the data points are generated in multiple sequences. For the feedforward architectures, these sequences are decomposed back into their separate time steps to obtain the required input-output pairs or data tuples. Each sequence is created from a random walk regarding the strain path, starting from the initial material state with \({}^{0}\sigma={}^{0}\epsilon={}^{0}\epsilon{}^{\alpha}=0\). A strain increment \({}^{0}\Delta\epsilon\) as a sample of a normal distribution with standard deviation \(s^{\Delta\epsilon}\) around mean 0 is applied to the initial state in a time increment \({}^{0}\Delta t\) sampled from a uniform distribution \({}^{0}\Delta t\in(\Delta t_{\text{min}}\,,\,\Delta t_{\text{max}})\). The strain rate is constant within this interval. The next material state is obtained by applying another strain increment \({}^{1}\Delta\epsilon\) within another time increment \({}^{1}\Delta t\) and so on, until the sequence contains \(N^{\text{ts}}\) timesteps. Here, sequences of length \(N^{\text{ts}}=100\) have shown to perform well for the recurrent architectures, independent from the examined material. In order to limit the strain to a reasonable range, the absolute value of the strain may not exceed \(|\epsilon|_{\text{max}}\). That is, if \(|{}^{n}\epsilon+{}^{n}\Delta\epsilon|>|\epsilon|_{\text{max}}\), the strain increment \({}^{n}\Delta\epsilon\) is sampled again until \(|{}^{n}\epsilon+{}^{n}\Delta\epsilon|\leq|\epsilon|_{\text{max}}\). The parameters of the random walk are chosen to be \(s^{\Delta\epsilon}=0.25\,\%\), \(|\epsilon|_{\text{max}}=2\,\%\), \(\Delta t_{\text{min}}=0.02\,\text{s}\) and \(\Delta t_{\text{max}}=0.1\,\text{s}\). Exemplarily, the strain path and respective stress response for the first sequence of the data base for V1 is shown in Fig. 11a.

The actual training data sets, i.e., the data that is effectively used during training, comprise a number of \(N^{\text{seq}}\) sequences or \(N^{\text{ds}}\) data tuples taken from this database. It was found that \(N^{\text{seq}}=100\) sequences for the recurrent or \(N^{\text{ds}}=1000\) data tuples for the feedforward architectures are sufficient without a loss of prediction quality, independent from the material behavior to be modeled. The selection of data from each database starts with the first sequence or the first time step of the first sequence and is continued chronologically. This means that the training data set of all feedforward architectures consists of information from the identical material states taken from the first 10 sequences and the recurrent architectures receive the same set of 100 sequences.

#### 5.1.2 Smooth strain path

In order to enable the NN representation of the internal variable with the auxiliary network, the FNN\({}^{\psi+\phi+\xi}\) architecture requires the training data to be generated from a simpler, more smooth strain path. This smooth strain path is described by a cubic spline connecting a set of manually chosen knots in the \(\epsilon\)-\(t\) space [7]. Once the spline is defined, the data points are generated by applying this strain path to the analytical models with uniformly distributed time increments \(\Delta t\in(0.01\,\text{s},0.02\,\text{s})\) until the data set contains 900 time steps. The same strain path is used for all materials V1, V2, P1 and P2 and is shown in Fig. 11b with the corresponding stress response for V1. The knots were chosen such that the data generated with this path are comparable to the data from the random walk sequences.

### Generation of validation data

In order to evaluate and compare the network performances for unseen data, the networks have to predict the constitutive response for a set of consecutive time steps taken from the strain path shown in Fig. 12a, where a constant time increment of 0.05 s is used. The path is chosen such that strain, strain rate, and time increments do not exceed the training range.

\begin{table}
\begin{tabular}{c c c} \hline \hline Description & Label & Parameters \\ \hline Nonlinear viscoelastic model with & V1 & \(E=1\,\text{GPa}\), \(E_{1}=10\,\text{GPa}\), \(\hat{\eta}_{1}=200\,\text{GPa}\,\text{s}\), \(a_{1}=-0.7\), \(b_{1}=0.1\) \\ one Maxwell element & & \\ Linear viscoelastic model with two & V2 & \(E=1\,\text{GPa}\), \(E_{1}=10\,\text{GPa}\), \(\eta_{1}=10\,\text{GPa}\,\text{s}\), \(E_{2}=20\,\text{GPa}\), \(\eta_{2}=5\,\text{GPa}\,\text{s}\) \\ Maxwell elements & & \\ \hline Elastoplastic model with kinematic & P1 & \(E=20\,\text{GPa}\), \(\sigma_{y}=100\,\text{MPa}\), \(H=10\,\text{GPa}\) \\ Elastoplastic model with mixed & P2 & \(E=20\,\text{GPa}\), \(\sigma_{y}=100\,\text{MPa}\), \(H=3\,\text{GPa}\), \(\hat{H}=3\,\text{GPa}\) \\ kinematic-isotropic hardening & & \\ \hline \hline \end{tabular}
\end{table}
Table 2: Chosen parameters of the four testing material models used as for training data generation and validation.
The networks thus do not have to extrapolate for this benchmark test. Reference data points are, analogously to the training data, obtained using an implicit Euler scheme.

Secondly, data for a further validation path with strain, strain rate and time increments exceeding the training range are generated, see Fig. 11(b). This will be used to evaluate the extrapolation capabilities of the different models later on. This path includes a small hysteresis beyond the training limits of the strain (\(2\,\mathrm{s}<t\leq 4\,\mathrm{s}\)) and a short interval (\(5\,\mathrm{s}<t\leq 5.08\,\mathrm{s}\)) of 16 time steps with \(\Delta t=0.005\,\mathrm{s}\) and \(\epsilon=62.5\,\mathrm{\char 37}\,\mathrm{s}^{-1}\), undercutting the limits of \(\Delta t\) and exceeding the limits of \(\epsilon\) in another hysteresis. Afterwards, a short relaxation with constant strain \(\epsilon=0\,\mathrm{\char 37}\) is performed until \(t=6\,\mathrm{s}\). Subsequently, large time increments outside of the training range are applied (\(\Delta t=0.125\,\mathrm{s}\) and \(\Delta t=0.2\,\mathrm{s}\) in \(6\,\mathrm{s}<t\leq 7\,\mathrm{s}\) and \(7\,\mathrm{s}<t\leq 8\,\mathrm{s}\), respectively). For all other time steps, the time increment equals \(\Delta t=0.05\,\mathrm{s}\). Finally, the strain is increased linearly until \(t=10\,\mathrm{s}\) and \(\epsilon=6\,\mathrm{\char 37}\).

### Normalization

The first and second derivatives of most common activation functions vanish for large absolute values of their arguments, rendering training impossible. In order to ensure nonzero gradients, the training data set, which for example contains stresses in the order of several MPa, has to be scaled down to values of magnitude 1. Since in general not all relevant features of the training data set are independent from each other, the scaling factors must be chosen appropriately, which is shown in the following.

Generally, an independent quantity \(f\) can be normalized to the range \((-1,1)\ni\tilde{f}\) via

\[\tilde{f}=\frac{f-m_{f}}{s_{f}}\quad\Longleftrightarrow\quad f=s_{f}\tilde{f} +m_{f}\quad\text{with}\quad m_{f}=\frac{1}{2}(f_{\text{max}}+f_{\text{min}}) \quad\text{and}\quad s_{f}=\frac{1}{2}(f_{\text{max}}-f_{\text{min}})\quad, \tag{34}\]

wherein \(\tilde{f}\) is the scaled value of \(f\) and \(f_{\text{max}}\) and \(f_{\text{min}}\) denote the maximum and minimum value of \(f\) across the whole data set [21].

Figure 11: Stresses and strains over time in the training data for the nonlinear viscoelastic material (V1): **(a)** the first sequence taken from the non-smooth random walk data base and **(b)** the smooth path generated with cubic splines used for \(\text{FNN}^{\text{w}+\phi+\frac{\pi}{2}}\).

Figure 12: Validation paths to evaluate the network performances for **(a)** interpolation and **(b)** extrapolation. The path in (a) is chosen such that neither strain nor time increments exceed the training range. Path (b) in turn exceeds those limits separately and includes sections of large strain rates. Black dashed lines indicate the training data limits.
For the blackbox models FNN\({}^{\sigma}\) and RNN\({}^{\sigma}\), the necessary quantities \(\varepsilon\), \(\sigma\) and \(\Delta t\) can be treated as independent and the scaling factors are thus obtained using Eq. (34). In all other architectures, not every feature is independent from the others. This is due to the occurring differential operators. Furthermore, it is reasonable to choose the scaling factors in such a way, that all relations, that hold for the unscaled quantities also hold for their scaled values to avoid back transformations during training. This simplifies and speeds up the calculation of the loss function. For the sake of simplicity, all \(m_{f}\) are set to zero in the following. Considering all remaining architectures, the relevant quantities are \(\varepsilon\), \(\sigma\), \(\Delta t\), \(\xi\), \(\tau\), \(\psi\), \(\phi\) and \(\phi^{*}\). As chosen independent quantities, the scaling factors of \(\varepsilon\), \(\Delta t\) and \(\psi\) are obtained using Eq. (34)\({}_{4}\). All other scaling factors follow from the evaluation of the relevant equations under the stated condition, resulting in the following relations:

\[s_{\sigma}=\frac{s_{\psi}}{s_{\varepsilon}}\quad,\quad s_{\xi}=s_{\varepsilon} \quad,\quad s_{\tau}=\frac{s_{\psi}}{s_{\varepsilon}}\quad,\quad s_{\phi}= \frac{s_{\psi}}{s_{\Delta t}}\quad\text{and}\quad s_{\phi^{*}}=\frac{s_{\psi} }{s_{\Delta t}}\quad. \tag{35}\]

The scaling is applied to all features in the respective data set and the training is performed using only these scaled values. The network itself consequently predicts only scaled values. The back transformation Eq. (34)\({}_{2}\) yields the values in familiar physical units. The trained network can now be applied to predict the material response for a given strain path.

## 6 Applications

Within this section, building on the data generated according to Sect. 5, each of the seven architectures presented in Sect. 4, i.e., FNN\({}^{\sigma}\), RNN\({}^{\sigma}\), FNN\({}^{\xi+\psi}\), RNN\({}^{\xi+\psi}\), FNN\({}^{\nu+\phi}\), FNN\({}^{\nu+\phi^{*}}\) and FNN\({}^{\nu+\phi+\xi}\), is tested on its capability to rebuild the constitutive response of the materials V1, V2, P1, and P2 after training. This is evaluated by the NN-based models' predicted stress responses for an unknown strain path given in Fig. 11(a), where the path is chosen such that strain, strain rate, and time increments do not exceed the training range. In addition, to analyze the extrapolation capability, the models' predicted responses are considered for a second validation path with strain, strain rate, and time increments exceeding the training range, cf. Fig. 11(b). In order to reduce the scope of the presented study to a reasonable level, the second validation path is investigated for V1 only. All of the used FNNs consist of an input layer, a single hidden layer and an output layer. The number of neurons in the hidden layer is denoted as \(N_{2}\).

### Training results and validation path without extrapolation

#### 6.1.1 Black box models

Fnn\({}^{\sigma}\)For each test material model, an FNN\({}^{\sigma}\) as described in Sect. 4.1.1 is created with the hyperparameters and training results given in Tab. 3. The model predictions are presented in Fig. 13.

The results for viscoelasticity show that the prediction of the stress behavior with FNN\({}^{\sigma}\) for the materials V1 is possible without further difficulties. For V2, in contrast, a precise prediction is only possible with two preceding time steps in the input, i.e., \(N^{\text{pt}}=2\). These limitations of the architecture can be traced back to the ambiguity of the expected output given a particular input. That is, given the material parameters for V1, the overstress \(\sigma_{1}^{\text{ov}}\) can be determined solely from the values of \({}^{n}\varepsilon\) and \({}^{n}\sigma\) given in the input. The material state is therefore determined unambiguously. For material V2, on the other hand, only the sum \(\sigma_{1}^{\text{ov}}+\sigma_{2}^{\text{ov}}\) of the overstresses can be determined from these quantities, but not their distinct values. Thus, the inner state is not clearly described and the stress prediction is not possible precisely, as Fig. 11(b) shows. For this reason, the network is provided with further information of an additional time step, see also Fig. 11(b).

\begin{table}
\begin{tabular}{c c c c c c c c} \hline \hline Material & \(N^{\text{pt}}\) & \(N_{2}\) & \(\mathcal{A}_{2}^{l}\) & \(N^{\text{ds}}\) & Iterations & \(\mathcal{L}\) & Val. \\ \hline V1 & 1 & 15 & tanh & 1000 & 3094 & \(2\cdot 10^{-4}\) & Fig. 13a \\ V2 & 1 & 25 & tanh & 1000 & 911 & \(2\cdot 10^{-3}\) & Fig. 13b \\ V2 & 2 & 25 & tanh & 1000 & 2903 & \(1\cdot 10^{-4}\) & Fig. 13b \\ P1 & 1 & 15 & ReLU & 1000 & 356 & \(3\cdot 10^{-6}\) & Fig. 13c \\ P2 & 1 & 15 & ReLU & 1000 & 1240 & \(8\cdot 10^{-3}\) & Fig. 13d \\ \hline \hline \end{tabular}
\end{table}
Table 3: FNN\({}^{\sigma}\) architectures: network hyperparameters with number of preceding time steps \(N^{\text{pt}}\), neurons in the hidden layer \(N_{2}\), and activation function \(\mathcal{A}_{2}^{l}\) as well as characteristic values from training with number of training data \(N^{\text{ds}}\) and loss value \(\mathcal{L}\) after the given number of iterations. The figure numbers with the corresponding validation load case are given in column Val.
The state of the elastoplastic material with kinematic hardening P1 is, analogously to V1, sufficiently described by \({}^{n}\epsilon\) and \({}^{n}\sigma\). However, in contrast to V2, elastoplasticity with mixed kinematic-isotropic hardening as present in P2 cannot be described by adding an additional time step to the input. This would require information about the current position and extent of the elastic region, which cannot be provided by a fixed amount of preceding times steps. Interestingly, the model tries to approximate the material behavior with kinematic hardening as good as possible.

Summarizing, FNN\({}^{\sigma}\) is thus only applicable for certain types of material behavior. Moreover, no further information is obtained apart from the stress prediction. On the other hand, the training process is computationally inexpensive and requires only stress-strain pairs and time steps for the data set.

RNN\({}^{\sigma}\)Using the same reference data generated with V1, V2, P1, and P2, the capability of the RNN\({}^{\sigma}\) model according to Wu et al. [29] is now investigated, where the network hyperparameters given in Tab. 4 have been used. With that, the results given in Fig. 14 could be achieved for the validation test. The RNN\({}^{\sigma}\) model is thus able to predict the stress response almost perfectly for all four test materials. Due to the comparatively large number of weights and the expensive gradient calculations for networks with recurrent cells, this broad applicability comes at the cost of the computational time needed for the training. Similar to FNN\({}^{\sigma}\), only stresses and strains are required for the training data set, but no additional physical information is obtained. However, both models, FNN\({}^{\sigma}\) and RNN\({}^{\sigma}\), do not allow statements to be made about whether the processes described by the respective model are embedded in a meaningful thermodynamic framework.

#### 6.1.2 Neural networks enforcing physics in a weak form

Fnn\({}^{\xi+\psi}\)After the investigation of the pure black box models given above, the first approach belonging to the class of NNs enforcing physics in a weak form is analyzed, the FNN\({}^{\xi+\psi}\) model according to Masi et al. [11]. Again, the same reference data generated with V1, V2, P1, and P2, are used here.

Using the hyperparameters given in Tab. 5, four FNN\({}^{\xi+\psi}\) networks were created. The results for the validation loading case given in Fig. 15 show that FNN\({}^{\xi+\psi}\) is able to produce precise predictions for stress and internal variables, i.e.,

\begin{table}
\begin{tabular}{c c c c c c c c} \hline \hline Material & \(N^{c}\) & \(N_{2}^{\text{FF}}\) & \(\mathcal{A}_{2}^{\text{FF}}\) & \(N^{\text{eqs}}\) & \(N^{\text{us}}\) & Iterations & \(\mathcal{L}\) & Val. \\ \hline V1 & 6 & 10 & tanh & 100 & 100 & 3465 & \(9\cdot 10^{-4}\) & Fig. 14a \\ V2 & 10 & 10 & tanh & 100 & 100 & 2219 & \(8\cdot 10^{-5}\) & Fig. 14b \\ P1 & 10 & 10 & ReLU & 100 & 100 & 5352 & \(7\cdot 10^{-3}\) & Fig. 14c \\ P2 & 12 & 20 & ReLU & 100 & 100 & 5281 & \(1\cdot 10^{-3}\) & Fig. 14d \\ \hline \hline \end{tabular}
\end{table}
Table 4: RNN\({}^{\sigma}\) architectures: network hyperparameters with number of values in the cell state \(N^{\text{c}}\), neurons in the hidden layer of the feedforward network \(N_{2}^{\text{FF}}\) with activations \(\mathcal{A}_{2}^{\text{FF}}\) as well as characteristic values from training with number of sequences \(N^{\text{seq}}\) of length \(N^{\text{ts}}\) and loss value \(\mathcal{L}\) after training for the given number of iterations.

Figure 13: Stress prediction of FNN\({}^{\sigma}\) for the four underlying test materials. The chosen hyperparameters are given in Tab. 3.
the available quantities. In addition to the stress-strain plots, diagrams showing free energy \(\psi\), internal variables \(g^{\alpha}\), and dissipation rate \(D\) are added. It can be seen that \(\psi\) and \(D\) vary from the expected results, although the order of magnitude and rough course mostly coincide. The condition \(D\geq 0\) is satisfied reliably. Thus, in contrast to the black box models, \(\mathrm{FNN}^{\xi+\psi}\) puts the predictions in a more physical framework, which has positive effects on the generalization capability of the network. For a detailed study on this architecture, see Masi et al. [11]. However, \(\mathrm{FNN}^{\xi+\psi}\) requires the provision of information about the internal variables for the training data, which is not a trivial task. This problem is discussed in Masi and Stefanou [47] for multiscale simulations.

Rnn\({}^{\xi+\psi}\)The second analyzed NN-based approach that enforces physics in a weak form is the \(\mathrm{RNNN}^{\xi+\psi}\) approach similar to He and Chen [48]. The network hyperparameters and the information on the training are given in Tab. 6. With that, the results given in Fig. 16 could be achieved for the considered validation path.

One can see, that the stress predictions are very precise for all four materials. Furthermore, the Clausius-Duhem inequality is also complied with, with the exception of few individual time steps. The internal variables, however, vary vastly from the expected results. Since no internal variables are given in the training data set and no further restrictions on their course are made, the network is free to choose any set of internal variables that allow accurate stress predictions and compliance with \(D\geq 0\). The network thus finds another representation of the set of internal variables and the free energy function.

Note, that for the kinematic hardening material P1, the network predicts only a single internal variable (\(N^{\xi}=1\)), although \(N=2\). This is possible since plastic strain \(\varepsilon^{\mathrm{pl}}\) and kinematic hardening variable \(\alpha\) coincide and the material response can be expressed in terms of only the plastic strain. The same simplification could as well be made for material P2 with \(N^{\xi}=2\) despite \(N=3\), but \(N^{\xi}=3\) is chosen.

In contrast to \(\mathrm{FNN}^{\xi+\psi}\), the provision of internal variables is not explicitly necessary to obtain precise, physically consistent stress predictions. This advantage comes at cost of a computational expensive recurrent cell. For both \(\mathrm{FNN}^{\xi+\psi}\) and \(\mathrm{RNN}^{\xi+\psi}\), although \(D\geq 0\) is enforced by a penalty term in the loss function and is also satisfied in this validation path without extrapolation, this cannot be fully guaranteed for this class of models, see Sect. 6.2.

Figure 14: Stress prediction of \(\mathrm{RNNN}^{\sigma}\) for the four underlying test materials. The chosen hyperparameters are given in Tab. 4.
\begin{table}
\begin{tabular}{c c c c c c c c c c c c c} Mat. & \(N^{\mathrm{c}}\) & \(N^{\mathrm{c}}\) & \(N^{\mathrm{c}}\) & \(N_{2}^{\mathrm{FP}}\) & \(N_{2}^{\mathrm{FP}}\) & \(\mathcal{A}_{2}^{\mathrm{FP}^{\mathrm{c}}}\) & \(\mathcal{A}_{2}^{\mathrm{FP}^{\mathrm{c}}}\) & \(\mathcal{A}_{2}^{\mathrm{FP}^{\mathrm{c}}}\) & \(w^{\sigma}\) / \(w^{D2\otimes}\) & \(N^{\mathrm{sq}}\) & \(N^{\mathrm{u}}\) & Iter. & \(\mathcal{L}^{\sigma}\) / \(\mathcal{L}^{D2\otimes}\) & Val. \\ \hline V1 & 6 & 1 & 10 & 10 & tanh & tanh & 1/5 & 100 & 100 & 5079 & \(2\cdot 10^{-3}/4\cdot 10^{-7}\) & Fig. 16a \\ V2 & 6 & 2 & 10 & 10 & tanh & tanh & 1/5 & 100 & 100 & 5244 & \(8\cdot 10^{-3}/3\cdot 10^{-5}\) & Fig. 16b \\ P1 & 10 & 1 & 10 & 10 & ReLU & tanh & 1/5 & 100 & 100 & 5329 & \(7\cdot 10^{-3}/1\cdot 10^{-4}\) & Fig. 16c \\ P2 & 12 & 3 & 15 & 20 & ReLU & tanh & 1/5 & 100 & 100 & 5140 & \(1\cdot 10^{-3}/5\cdot 10^{-6}\) & Fig. 16d \\ \hline \end{tabular}
\end{table}
Table 6: RNN\({}^{\mathrm{s+\psi}}\) architectures: network hyperparameters with number of values in the cell state \(N^{\mathrm{c}}\), number of predicted internal variables \(N^{\mathrm{\xi}}\) neurons in the hidden layer of the feedforward networks \(N_{2}^{\mathrm{FF}^{\mathrm{c}}}\), \(N_{2}^{\mathrm{FP}^{\mathrm{v}}}\) with activations \(\mathcal{A}_{2}^{\mathrm{FF}^{\mathrm{c}}}\), \(\mathcal{A}_{2}^{\mathrm{FF}^{\mathrm{v}}}\) as well as weighting factors of the loss function \(w^{\sigma}\), \(w^{D2\otimes}\), number of training sequences \(N^{\mathrm{seq}}\) of length \(N^{\mathrm{ts}}\) and loss values \(\mathcal{L}^{\sigma}\), \(\mathcal{L}^{D2\otimes}\) after training for the given number of iterations. Validation load cases can be found in the figures given in Val.

Figure 15: Prediction of stress, free energy, internal variables and dissipation rate for the four test materials using \(\mathrm{FNN}^{\mathrm{\xi+\psi}}\).
Figure 16: Prediction of stress, free energy, internal variables and dissipation rate for the four test materials using RNN\({}^{\xi+\psi}\).
#### 6.1.3 Neural networks enforcing physics in a strong form

FNN\({}^{\psi+\phi^{*}}\)Finally, the analysis of approaches belonging to the class of _NNs enforcing physics in a strong form_ is done by using the reference data generated with V1, V2, P1, and P2. The first model is the FNN\({}^{\psi+\phi}\) model according to Huang et al. [49]. Here, the network hyperparameters and the information on the training are given in Tab. 7 and the results achieved for the considered validation path are given in Fig. 17.

One can see, that the stress predictions are very precise for both viscoelastic materials, V1 and V2. In addition to the stress-strain plots, diagrams showing internal variables \(\xi^{a}\), free energy \(\psi\), internal forces \(\tau^{a}\), dissipation rate \(\mathcal{D}\), and dissipation potential \(\phi\) are added. From that one can see that all of these quantities are almost identical to the reference model, although these values do not explicitly appear in the loss, cf. Eq. (31). Only for V1 there is a visible discrepancy in the dissipation potential. Furthermore, in contrast to the former two models, the important condition \(\mathcal{D}\geq 0\) is now guaranteed for all admissible load paths by construction of the network architecture. However, FNN\({}^{\psi+\phi}\) requires the provision of information about the internal variables for the training data when no adapted training is used [7].

Regarding rate independent materials, the FNN\({}^{\psi+\phi}\) model is able to rebuild the response of the elastoplastic material P1. Interestingly, however, the internal forces \(\tau^{a}\) do not match with the reference. This is due to the fact that \(\xi^{1}=\epsilon^{\mathrm{pl}}\) and \(\xi^{2}=\mathbf{a}\) coincide, i.e. \(\mathbf{a}=\epsilon^{\mathrm{pl}}\), cf. Tab. 1. Therefore, there are any number of ways to divide the internal forces, all of which lead to the same, good result, even if the relation between the internal variables is not incorporated into the network, c.f. Remark 1. Furthermore, regarding the plots of \(\mathcal{D}\) and \(\phi\), one can see that the the original non-continuous curve shapes at the zero line are approximated by smooth ones. This process is often called a regularization of the rate independent model, i.e., an approximation by a rate dependent one [59]. Due to the choice of the softplus activation function this results all by itself within the training process. Regarding the results for P2, a poor prediction becomes apparent in Fig. 17 as soon as the plastic regime is reached. The stress response is, similar to FNN\({}^{\psi}\), approximated with only kinematic hardening instead of kinematic and isotropic hardening. According to the authors, this is due to the fact that there is no function \(\phi\), for which Eq. (4) can be applied without explicitly considering the relation between the rates of the plastic strain \(\dot{\epsilon}^{\mathrm{pl}}\) and the isotropic hardening variable \(\dot{\hat{\mathbf{a}}}\), c.f. Remark 1. If this relation, i.e., \(\dot{\hat{\mathbf{a}}}=|\dot{\epsilon}^{\mathrm{pl}}|\), is incorporated into the model architecture, the model is able to make accurate predictions. The authors have already verified this.

FNN\({}^{\psi+\phi}\)Instead of using the dissipation potential \(\phi\), the model FNN\({}^{\psi+\phi^{*}}\) according to As'ad and Farhat [6] makes use of the dual dissipation potential \(\phi^{*}\). Here, the network hyperparameters and the information on the training are given in Tab. 8. With that, the results given in Fig. 18 could be achieved for the considered validation path.

Similar to FNN\({}^{\psi+\phi*}\), the predictions for stress, free energy, dissipation rate, etc. are very precise for the viscoelastic materials V1 and V2. Again, only for V1 there is a visible discrepancy in the dual dissipation potential. FNN\({}^{\psi+\phi^{*}}\) requires the provision of information about the internal variables for the training data when no adapted training is used.

Regarding the prediction of FNN\({}^{\psi+\phi^{*}}\) for the elastoplastic materials P1 and P2, a rather poor result shows up in Fig. 18 even for kinematic hardening. However, this is not surprising, since the dual dissipation potentials \(\phi^{*}\) of the reference models P1 and P2 have a shape that cannot be reasonably represented with the selected activation functions, cf. Tab. 1, and restrictions to the networks weights.

FNN\({}^{\psi+\phi+\xi}\)This architecture models the free energy and dissipation potential without requiring knowledge about the internal variables during training. The architecture is restricted to a single internal variable, i.e., a single FNN for the internal variable with scalar output. Using the network parameters in Tab. 9, the results in Fig. 19 could be achieved. These results show, that the architecture is able to learn reasonable representations of the potentials for V1, V2 and P1, but not P2. All of the predictions are less precise compared to FNN\({}^{\psi+\phi}\). Similar to RNN\({}^{\psi+\psi}\), the internal variables are
not learned exactly as in the reference model. However, normalization would show that the courses are similar. Because other potentials are learned if \(\xi(t)\) differs from the reference model, this leads to equivalent results.

Interestingly, the architecture finds a set of potentials for V2, that allows for surprisingly accurate stress predictions, although only one internal variable is used. The prediction for P1 shows a rather strong rate dependency. This can be attributed to a second regularization mechanism: besides the regularization of the dissipation potential, the sudden increase of \(\dot{e}^{\text{pl}}\) when leaving the elastic region cannot be modeled accurately by the FNN \(\xi^{\text{NN}}(t)\). P2 in turn cannot be modeled at all, since at least two internal variables are necessary to describe the path dependency and, which is more important, due to the reasons discussed for FNN\({}^{\psi+\phi}\) above. The architecture FNN\({}^{\psi+\phi+\xi}\) is thus more suitable for viscoelastic materials.
Figure 17: Prediction of stress, internal variables, free energy, internal forces, dissipation rate and dissipation potential for the four test materials using \(\mathrm{FNN}^{\mathsf{w}+\phi}\).
Figure 18: Prediction of stress, internal variables, free energy, internal forces, dissipation rate and dual dissipation potential for the four test materials using FNNv+\(\phi^{*}\).
Figure 19: Prediction of stress, internal variables, free energy, internal forces, dissipation rate and dissipation potential for the four test materials using the advanced training method \(\mathrm{FNN}^{\nu+\phi+\xi}\).
### Validation path with extrapolation

After investigating which NN-based approaches can reproduce which material behavior well, i.e. V1, V2, P1 and P2, the extrapolation capability of these approaches is now compared. This is done using the path given in Fig. 12b. Since all NN-based approaches have shown that they can reproduce the viscoelastic material V1 well, only this will be considered in the following. The exact same models that were validated with the strain path without extrapolation are now used to carry out the extrapolation study. The results of this final study are given in Fig. 20.

#### 6.2.1 Black box NNs

To start with, the _extrapolation_ capabilities of the first _black box model_ FNN\({}^{\sigma}\) are evaluated. As can be seen, the model is easily able to extrapolate to strains of up to \(\epsilon=3\,\%\) which is outside the training range \(\epsilon^{\text{train}}\in[-2,2]\,\%\). After reentering the training range without notable inaccuracies, the model fails to extrapolate into ranges of higher strain rates and exhibits large errors. However, after a short relaxation, the model is able to produce accurate predictions for strain increments outside the training range and is even able to yield reasonable values for the stress up to \(4\,\%\). Thus, the FNN\({}^{\sigma}\) model is surprisingly good at extrapolating, except for increased strain rates \(\epsilon\). However, it should be noted that there is no possibility to make statements about the thermodynamics of the model.

The second _black box model_, RNN\({}^{\sigma}\), initially shows a similarly good prediction quality, but also fails to predict stresses for larger strain rates and deviates for time increments of \(\Delta t=0.2\,\text{s}\). It should be noted that the extrapolation behavior of RNNs often differs significantly for different training configurations. Even strong oscillations could be observed in some cases. Thus, compared to the FNN\({}^{\sigma}\) model, the RNN\({}^{\sigma}\) is worse at extrapolating. However, it can be used for a broader class of material behavior compared to FNN\({}^{\sigma}\), see Figs. 13 and 14. The missing possibility to make statements about the thermodynamics remains.

#### 6.2.2 NNs enforcing physics in a weak form

Now, the _extrapolation_ behavior of _NNs enforcing physics in a weak form_ is considered. As shown in Fig. 20, the model FNN\({}^{5+\psi}\) agrees very well with the reference V1 up to maximum strains of \(2.5\,\%\), but slightly deviates above. After reentering the training range, precise predictions are made again. Increased strain rates do lead to errors as well, but these are substantially reduced compared to the previous models. The model is able to cope with increased time increments and again precisely predicts the stress up to \(2.5\,\%\). This also applies to the free energy \(\psi\) and the dissipation rate, for which \(D\geq 0\) applies up to that point. Thus, the model does not violate the second law of thermodynamics up to this point of the loading path. However, when the strain is further increased to \(\epsilon>3\,\%\), significant errors in \(\sigma\) occur. Furthermore, one can see that the model predicts negative values for \(D\). Thus, FNN\({}^{5+\psi}\) is very good at extrapolation for the most part of the loading path. However, from a certain level of strains, unphysical predictions can be seen. This is due to the fact that the fulfillment of the second law of thermodynamics is enforced only by a penalty term in the loss and is not fulfilled a priori.

The RNN\({}^{5+\psi}\) model shows similarly good results, but is more precise in strains up to \(3\,\%\). As with the FNN\({}^{5+\psi}\), completely unphysical predictions with \(D<0\) occur for strains greater than \(3\,\%\). Thus, RNN\({}^{5+\psi}\) provides acceptable results when extrapolating up to \(\epsilon=3\,\%\). For this, however, no inner variables are necessary for the training here. Finally, as it is also the case for FNN\({}^{5+\psi}\), unphysical predictions may occur.

#### 6.2.3 NNs enforcing physics in a strong form

Lastly, the _extrapolation_ behavior of the three _NNs enforcing physics in a strong form_, FNN\({}^{\psi+\phi}\), FNN\({}^{\psi+\phi^{*}}\), and FNN\({}^{\psi+\phi}\), is analyzed. As shown in Fig. 20, highly accurate predictions can be achieved with FNN\({}^{\psi+\phi}\) and FNN\({}^{\psi+\phi^{*}}\) for \(\sigma\), \(\psi\), as well as \(D\). This applies to the entire load path, except for the last piece of FNN\({}^{\psi+\phi^{*}}\) with \(\epsilon>4\,\%\). FNN\({}^{\psi+\phi}\) is even capable to produce very precise results for strains of \(\epsilon=6\,\%\). The extrapolation using FNN\({}^{\psi+\phi+\phi}\) is not as accurate as the two former architectures, but still good considering the mediore interpolation results. Summarizing, neither highly increased strains or strain rates nor time increments outside of the training range lead to considerable deviations from the expected material response. Particularly noteworthy here is that in any case \(D\geq 0\) is ensured. Thus, the predictions are always in accordance with the second law given by the CDI (1). All in all, this model class is best suited for extrapolation which is due to the strong physical background inserted here.
Figure 20: Predictions of the six considered NNs tested for an extrapolation path with \(\varepsilon\), \(\dot{\varepsilon}\), and \(\Delta t\) not included in the training data set: **(a)** stress \(\sigma\), **(b)** free energy \(\psi\), and **(c)** dissipation rate \(D\). The viscoelastic model V1 serves as a reference and has been used for generation of training data.
## 7 Conclusions

In this work, a classification of the variety of NN-based approaches to modeling inelastic constitutive behavior with particular attention to the thermodynamic framework as well as a unified formulation of these approaches is provided. To this end, a division of NN-based approaches into _black box NNs, NNs enforcing physics in a weak form_ and _NNs enforcing physics in a strong form_ is made and an application to both 1D elastoplastic and viscoelastic data is made.

After a compact literature review, a short overview on continuum based constitutive modeling including standard viscoelastic and elastoplastic models and a condensed repetition to the basis of FNNs and RNNs is given. Based on this, a total of seven NN-based approaches are presented with a detailed description on training and application. It is shown in which way the second law of thermodynamics is taken into account in the respective model and what data are necessary for training. Thereafter, the generation of training data and the application of the seven NN-based approaches to these data are shown. It can be seen that all considered models are able to represent viscoelasticity, whereas elastoplasticity cannot be represented by all approaches. Furthermore, the models' extrapolation capabilities are analyzed.

In summary, the results of this work show that NN-based models are promising for the description of complex inelastic behavior and prove to be very flexible. They have the potential to replace the time-consuming task of classical constitutive model formulation and calibration piecemeal and to enable automated workflows. _Black box models_, however, do not allow any conclusions to be drawn as to whether the processes described by the respective model are embedded in a meaningful thermodynamic framework and are therefore not recommended for use. Furthermore, when applied to unknown load paths outside the training domain, the poor extrapolation capability of NNs can lead to large errors within stress predictions. _NNs enforcing physics in a weak form_ on the other hand, have a higher content of physics included into the model. Only the class of _NNs enforcing physics in a strong form_ can really ensure that no violation of thermodynamics occurs. Moreover, they also have shown to be characterized by the best extrapolation behavior of the considered models.

Regarding the follow-up of the study presented here, several extensions are planned in the future. For instance, an extension to the general 3D case at finite strains has to be made. Thus, a variety of further physical principles and conditions have to be included in the comparison for this, e.g., principles as objectivity or material symmetry.

## CRediT authorship contribution statement

**Max Rosenkranz:** Conceptualization, Formal analysis, Investigation, Methodology, Visualization, Software, Validation, Writing - original draft, Writing - review and editing. **Karl A. Kalina:** Conceptualization, Formal analysis, Methodology, Writing - original draft, Writing - review and editing. **Jorg Brummund:** Formal analysis, Methodology, Writing - review and editing. **Markus Kastner:** Funding acquisition, Resources, Writing - review and editing.
