Neural networks meet hyperelasticity: A guide to enforcing physics
In the present work, a hyperelastic constitutive model based on neural
networks is proposed which fulfills all common constitutive conditions by
construction, and in particular, is applicable to compressible material
behavior. Using different sets of invariants as inputs, a hyperelastic
potential is formulated as a convex neural network, thus fulfilling symmetry of
the stress tensor, objectivity, material symmetry, polyconvexity, and
thermodynamic consistency. In addition, a physically sensible stress behavior
of the model is ensured by using analytical growth terms, as well as
normalization terms which ensure the undeformed state to be stress free and
with zero energy. In particular, polyconvex, invariant-based stress
normalization terms are formulated for both isotropic and transversely
isotropic material behavior. By fulfilling all of these conditions in an exact
way, the proposed physics-augmented model combines a sound mechanical basis
with the extraordinary flexibility that neural networks offer. Thus, it
harmonizes the theory of hyperelasticity developed in the last decades with the
up-to-date techniques of machine learning. Furthermore, the non-negativity of
the hyperelastic neural network-based potentials is numerically examined by
sampling the space of admissible deformations states, which, to the best of the
authors' knowledge, is the only possibility for the considered nonlinear
compressible models. For the isotropic neural network model, the sampling space
required for that is reduced by analytical considerations. In addition, a proof
for the non-negativity of the compressible Neo-Hooke potential is presented.
The applicability of the model is demonstrated by calibrating it on data
generated with analytical potentials, which is followed by an application of
the model to finite element simulations. In addition, an adaption of the model
to noisy data is shown and its [...]
Lennart Linden, Dominik K. Klein, Karl A. Kalina, Jörg Brummund, Oliver Weeger, Markus Kästner
# Neural networks meet hyperelasticity: A guide to enforcing physics

###### Abstract

In the present work, a hyperelastic constitutive model based on neural networks is proposed which fulfills all common constitutive conditions by construction, and in particular, is applicable to compressible material behavior. Using different sets of invariants as inputs, a hyperelastic potential is formulated as a convex neural network, thus fulfilling symmetry of the stress tensor, objectivity, material symmetry, polyconvexity, and thermodynamic consistency. In addition, a physically sensible stress behavior of the model is ensured by using analytical growth terms, as well as normalization terms which ensure the undeformed state to be stress free and with zero energy. The normalization terms are formulated for both isotropic and transversely isotropic material behavior and do not violate polyconvexity. By fulfilling all of these conditions in an exact way, the proposed physics-augmented model combines a sound mechanical basis with the extraordinary flexibility that neural networks offer. Thus, it harmonizes the theory of hyperelasticity developed in the last decades with the up-to-date techniques of machine learning. Furthermore, the non-negativity of the hyperelastic potential is numerically verified by sampling the space of admissible deformations states, which, to the best of the authors' knowledge, is the only possibility for the considered nonlinear compressible models. The applicability of the model is demonstrated by calibrating it on data generated with analytical potentials, which is followed by an application of the model to finite element simulations. In addition, an adaption of the model to noisy data is shown and its extrapolation capability is compared to models with reduced physical background. Within all numerical examples, excellent and physically meaningful predictions have been achieved with the proposed physics-augmented neural network.

**Key words:** hyperelasticity, physics-augmented neural networks, normalization, anisotropy, constitutive modeling, finite element simulation

## 1 Introduction

The mechanical principles underlying hyperelasticity were extensively discussed in the last decades, but for a long time, fulfilling them all at once could be seen as "the main open problem of the theory of material
behavior" (Truesdell and Noll [51]). For instance, while both the polyconvexity [5, 6] and objectivity condition have a sound mechanical motivation, objective strain measures easily violate the polyconvexity condition [24]. This does not mean that different constitutive conditions contradict each other - it rather shows the big challenge of fulfilling them all at the same time. With an increasing amount of restrictions a model should fulfill, this effort increases considerably, and it took sophisticated approaches to construct analytical models which fulfill all relevant conditions at the same time [8, 44]. In addition, the calibration of such models is also not a trivial task and requires a lot of knowledge [42].

To overcome the time consuming task of formulating classical constitutive models and to improve the restricted functional relationships that most of these analytical models bring along, concepts like the data-driven mechanics approach [23] or modern machine learning methods such as Gaussian process regression [11, 13] or _neural networks (NNs)_[1, 26] represent promising alternatives. For the first time, the idea of applying NNs in constitutive modeling was proposed in the early 1990s by Ghaboussi et al. [16]. However, in this early phase, mostly pure black-box approaches were used, i.e., networks that do not take into account any physical principles and therefore can only reproduce the training dataset, here consisting of stress-strain couples, well but extrapolate poorly. To remedy this weakness, a fairly new trend in NN-based constitutive modeling, and in scientific machine learning in general [41], is to include essential underlying physics in a strong form, e.g., by using adapted network architectures, or in a weak form, e.g., by modifying the loss term for the training [34]. These types of approaches, coined as _physics-informed_[22], _mechanics-informed_[4], _physics-augmented_[25], _physics-constrained_[19], or _thermodynamics-based_[37], enable an improvement of the extrapolation capability and the usage of sparse training data [22, 27], which is particularly important when constitutive models are to be fitted to experimental data.

In the following, a brief overview on the mentioned NN-based approaches applied to _finite strain hyper-elasticity_ modeling is given. Regarding isotropic materials, transferred from analytical models, the works [29, 46] propose to approximate the elastic potential by a feed-forward neural network (FFNN) with three deformation-type invariants as input and thus fulfill several constitutive conditions, e.g., thermodynamic consistency, objectivity, or material symmetry. However, similar to the approaches [28, 31, 43] applied to anisotropic problems, the elastic potential is needed directly for training within [29, 46]. In the meantime, NNs using _invariants_ as inputs and the hyperelastic potential as output, thus also being a priori thermodynamically consistent, have become a fairly established approach [12, 19, 24, 25, 30, 32, 33, 48]. Thereby, a more sophisticated training is applied, which allows the direct calibration of the network by tuples of stress and strain, i.e., the derivative of the energy with respect to the deformation is used in the loss term. This technique is also named Sobolev training [52, 54]. Alternatively, in order to ensure thermodynamic consistency a posteriori, a previously trained network predicting stress coefficients can be used to construct a pseudo-potential [20]. An NN-based approach which is coupled to a specific model, the so-called micro-sphere approach, is presented in [57].

Besides the mentioned requirements, namely thermodynamic consistency, objectivity, and material symmetry, there exist further physical conditions, e.g., ellipticity, which ensures material stability [56]. However, since ellipticity is difficult to verify and ensure, the concept of polyconvexity of the strain energy potential [5, 6], which implies ellipticity and is mathematically linked to the existence and stability of solutions of the elasticity problem, is preferable for the formulation of constitutive models [39]. There are several approaches for building polyconvex NNs [4, 7, 12, 24, 33, 47, 48, 50], with the most notable technique for incorporating this condition being the use of _input convex neural networks (ICNNs)_ originally introduced by Amos et al. [3]. For the fulfillment of the growth condition, a special network architecture may be applied [19], whereas using analytical growth terms is more widely spread [24, 25]. Finally, while several works introduce correction terms which ensure normalization conditions, they are either not polyconvex [9, 50] or restricted to the case of nearly incompressible material behavior [33]. E.g., the model proposed in [33] includes terms of the form \((I_{1}-3)^{2}\), where \(I_{1}\) denotes the first invariant of the right Cauchy-Green deformation tensor. However, in order to preserve polyconvexity of \(I_{1}\), the functions acting on it must in general be convex and _non-decreasing_,
cf. [24, Remark A.10]. For the quadratic function used in [33], this holds only if \(I_{1}\) is bounded from below by \(3\), which is only the case for \(\det\mathbf{F}=1\), cf. [44, Corollary A.11], where \(\mathbf{F}\) denotes the deformation gradient. Thus, for \(\det\mathbf{F}\neq 1\), the polyconvexity of the potential in [33] is not ensured by construction. Overall, to the best of our knowledge, the models found in literature so far only fulfill subsets of constitutive conditions belonging to compressible hyperelasticity in an exact way, while the remaining ones are only fulfilled in an approximate fashion, i.e., they are taken into account by penalty terms in the loss [49, 55]. Becoming more specific, while [24] fulfills the polyconvexity condition in an exact way, the normalization condition is only approximated by learning it through the calibration data. On the other side, while [50] uses a stress correction for the exact fulfillment of the normalization condition, this stress correction term includes non-diagonal components of the right Cauchy-Green deformation tensor, and thus is in general not polyconvex and violates material symmetry.

Concluding on NN-based constitutive models, fulfilling all common constitutive conditions of compressible hyperelasticity in an exact way at the same time has so far remained an open challenge.1 In the present work, such an approach which consequently accounts for _thermodynamic consistency, symmetry of the stress tensor, objectivity, material symmetry, polyconvexity, growth condition_ as well as _normalization of energy and stress_ by construction of the network architecture is systematically derived. Regarding the introduced new family of NN-based hyperelastic models, which fulfill all of the aforementioned conditions in an _exact_ way, we advocate for naming them as _physics-augmented neural networks (PANNs)_. The proposed framework will be very valuable in fields, where highly flexible and at the same time physically sensible constitutive models are required, such as the simulation of microstructured materials [14, 19]. For this purpose, the PANN approach is build up by extending the aforementioned model [24] by polyconvex normalization terms, followed by a detailed analytical and numerical study analyzing the overall model. Thereby, it is clearly explained step by step how all the considered physical conditions are incorporated into the PANN approach, which is particularly applicable to compressible anisotropic material behavior. Our framework is applied to the isotropic and transversely isotropic case, where several descriptive examples including multiaxial stress-strain states and noisy data are considered. The data basis for training is thereby generated by using analytical potentials. After the calibration of the models, they are applied within finite element (FE) computations to demonstrate their usability and general accuracy.

Footnote 1: This is only true for the case of compressible elastic materials. A corresponding proposal is made by Linka and Kuhl [33] for the incompressible case. However, apart from the restriction to incompressible materials, it is less general in some other points compared to the approach presented here.

The outline of the manuscript is as follows: In Sec. 2, the fundamentals of hyperelasticity are discussed, which are then applied to the proposed NN model in Sec. 3. In Sec. 4, numerical examples are presented. Finally, in Sec. 5 we conclude with a short discussion on the importance of augmenting neural networks with physics.

NotationThroughout this work the space of tensors

\[\mathcal{L}_{n}:=\underbrace{\mathbb{R}^{3}\otimes\cdots\otimes\mathbb{R}^{3} }_{\text{$n$-times}}\ \forall n\in\mathbb{N} \tag{1}\]

is used, except for a tensor of rank zero. In Eq. (1), \(\mathbb{R}^{3}\), \(\mathbb{N}\) and \(\otimes\) denote the Euclidean vector space, the set of natural numbers without zero and the dyadic product, respectively. Tensors of rank one and two are given by boldface italic symbols in the following, i. e., \(\mathbf{a}\in\mathcal{L}_{1}\) or \(\mathbf{B},\mathbf{C}\in\mathcal{L}_{2}\). Transpose and inverse of a second order tensor \(\mathbf{B}\) are marked by \(\mathbf{B}^{T}\) and \(\mathbf{B}^{-1}\), respectively. Furthermore, trace, determinant and cofactor are denoted by \(\operatorname{tr}\mathbf{B}\), \(\det\mathbf{B}\) and \(\operatorname{cof}\mathbf{B}:=\det(\mathbf{B})\mathbf{B}^{-T}\). The set of invertible second order tensors with positive determinant is denoted by \(\mathcal{E}\mathcal{L}^{+}(3):=\{\mathbf{A}\in\mathcal{L}_{2}\,|\det\mathbf{A}>0\}\), while the orthogonal group and special orthogonal group in \(\mathbb{R}^{3}\) are denoted by \(\mathcal{O}(3):=\big{\{}\mathbf{A}\in\mathcal{L}_{2}\,|\,\mathbf{A}^{T}\cdot\mathbf{A}= \mathbf{1}\big{\}}\) and \(\mathcal{SO}(3):=\big{\{}\mathbf{A}\in\mathcal{L}_{2}\,|\,\mathbf{A}^{T}\cdot\mathbf{A}= \mathbf{1},\,\det\mathbf{A}=1\big{\}}\), respectively.
Here, \(\mathbf{1}\in\mathcal{L}_{2}\) denotes the second order identity tensor. The space of symmetric second order tensors is denoted as \(\mathcal{S}\mathcal{S}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U }\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U }\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U} \mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}\mathcal{U}
a potential

\[\psi:\mathcal{SD}^{+}(3)\to\mathbb{R},\quad\mathbf{F}\mapsto\psi(\mathbf{F}) \tag{5}\]

is introduced, which corresponds to the strain energy density stored in the body and is equal to the Helmholtz free energy density \(W\). Thus, to satisfy the inequality (4) for arbitrary \(\dot{\mathbf{F}}\), the corresponding first Piola-Kirchhoff stress follows to

\[\mathbf{P}=\frac{\partial\psi}{\partial\mathbf{F}}\,. \tag{6}\]

By this definition, the stress tensor is a gradient field which implies energy conservation and path-independency, thus being **thermodynamically consistent** by construction [19, 30, 33].2

Footnote 2: Note that the constitutive equations must always be physically consistent, i.e., without contradiction to the introduced balance equations. Special importance is attached to the compatibility with the second law of thermodynamics, given here in the form of the Clausius-Duhem inequality (4). Therefore, the thermodynamic consistency is especially highlighted here.

The hyperelastic potential is subjected to further mathematical and physical considerations, which are shortly discussed in the following and also illustrated in Fig. 1. For a more detailed introduction to hyperelasticity, the reader is referred to [17, 18, 42].

First of all, \(\psi(\mathbf{F})\) has to be constructed such that the compatibility with the balance of angular momentum is ensured. By using Eqs. (3) and (6) this requirement is expressed as

\[\frac{\partial\psi}{\partial\mathbf{F}}\cdot\mathbf{F}^{T}=\mathbf{F}\cdot \frac{\partial\psi}{\partial\mathbf{F}^{T}}\,. \tag{7}\]

Using the stress transformations introduced in Sec. 2.1, this results in the requirement for the **symmetry of the stress tensors \(\mathbf{\sigma}\)** and \(\mathbf{T}\), e.g., \(\mathbf{T}=\mathbf{F}^{-1}\cdot\frac{\partial\psi}{\partial\mathbf{F}^{T}}\overset{!}{=} \frac{\partial\psi}{\partial\mathbf{F}^{T}}\cdot\mathbf{F}^{-T}=\mathbf{T}^{T}\,.\)

Secondly, a constitutive model should be independent on the choice of observer, which is referred to as **objectivity**. In hyperelasticity, this is formalized as

\[\psi(\mathbf{Q}\cdot\mathbf{F})=\psi(\mathbf{F})\quad\forall\mathbf{F}\in \mathcal{SD}^{+}(3),\mathbf{Q}\in\mathcal{SO}(3)\,. \tag{8}\]

The constitutive equations should also reflect the material's underlying (an-)isotropy which is expressed as **material symmetry** and mathematically written as

\[\psi(\mathbf{F}\cdot\mathbf{Q}^{T})=\psi(\mathbf{F})\quad\forall\mathbf{F}\in \mathcal{SD}^{+}(3),\mathbf{Q}\in\mathcal{G}\subseteq\mathcal{O}(3)\,, \tag{9}\]

where \(\mathcal{G}\) denotes the symmetry group of the material under consideration.

Furthermore, we consider hyperelastic potentials which are **polyconvex**[5, 6, 8, 44], allowing for a representation

\[\psi(\mathbf{F})=\mathcal{P}\left(\mathbf{F},\,\text{cof}\,\mathbf{F},\,\text {det}\,\mathbf{F}\right)\,, \tag{10}\]

where \(\mathcal{P}\left(\mathbf{F},\,\text{cof}\,\mathbf{F},\,\text{det}\,\mathbf{F}\right)\) is convex in its arguments, with \(\text{cof}\,\mathbf{F}=\text{det}(\mathbf{F})\,\mathbf{F}^{-T}\) denoting the cofactor of the deformation gradient. Polyconvexity stems from a quite theoretical context - however, it implies ellipticity [56], which is of practical importance and far more challenging to include in the model formulation than polyconvexity. The ellipticity (or rank-one convexity) condition [39, 56]

\[(\mathbf{a}\otimes\mathbf{b}):\,\frac{\partial^{2}\psi}{\partial\mathbf{F} \partial\mathbf{F}}\,:\,\,(\mathbf{a}\otimes\mathbf{b})\geq 0\qquad\forall\mathbf{a},\mathbf{b} \in\mathcal{L}_{1}\,, \tag{11}\]
ensures material stability of the model, which leads to a favorable behavior in numerical applications. Note that, in order to ensure polyconvexity, all steps that are made in the construction of the hyperelastic potential must preserve polyconvexity, in particular, polyconvex invariants have to be used, cf. Sec. 2.3.

In addition, a variety of coercivity conditions can be considered, the most common one being the **volumetric growth condition**

\[\psi\left(\mathbf{F}\right)\rightarrow\infty\quad\text{as}\quad\left(J\to 0^{+} \quad\vee\quad J\rightarrow\infty\right)\,, \tag{12}\]

in order to take into account the observation that a material body can not be compressed to a volume of zero or expanded to an infinite volume [18]. In particular, the case of large volumetric compression is important, because it may be in a relevant range of practical engineering applications.

Finally, the functional relationship given in Eq. (5) is subjected to further considerations on a physically

Figure 1: Schematic depiction of the common conditions on elastic potential and stresses: (a) thermodynamic consistency, (b) symmetry of the Cauchy stress, (c) objectivity, (d) material symmetry, (e) polyconvexity, (f) volumetric growth condition, (g) and (h) normalization conditions for energy and stress, as well as (i) non-negativity of energy. Within (e)–(i), solid blue lines denote models accounting for the respective condition, while the dashed red lines mark models which violate it. In (c) and (d), \(\mathrm{d}\mathbf{X}\), \(\mathrm{d}\mathbf{X}^{*}\), \(\mathrm{d}\mathbf{x}\), and \(\mathrm{d}\mathbf{x}^{*}\) denote material line elements for different deformation states, where \(\mathrm{d}\mathbf{X}^{*}\) and \(\mathrm{d}\mathbf{x}^{*}\) are transformed by a rotation of \(\mathrm{d}\mathbf{X}\) and \(\mathrm{d}\mathbf{x}\), respectively.
sensible behavior. In the undeformed configuration, i.e., \(\mathbf{F=1}\), the **normalization** conditions

\[\psi\left(\mathbf{F=1}\right)\overset{!}{=}0\quad\text{and}\quad\mathbf{P}(\mathbf{F=1}) \overset{!}{=}\mathbf{0} \tag{13}\]

for both energy and stress should hold. Besides the normalization, the free energy should increase in any case if deformation appears. Thus, in addition to Eq. (13)\({}_{1}\), the **non-negativity of the strain energy**, i.e., \(\psi\left(\mathbf{F}\right)\geq 0\), is required.

By formulating the potential in terms of invariants following from the right Cauchy-Green deformation tensor \(\mathbf{C}\) and a set of structural tensors, summarized in a set \(\mathcal{S}^{\Box}\), which reflects the material symmetry of the material body under consideration [17, 18], i.e.,

\[\psi\!:\mathbb{R}^{m}\to\mathbb{R},\quad\mathbf{I}\mapsto\psi\left(\mathbf{I}\right)\,, \tag{14}\]

both the **objectivity** and the **material symmetry** condition are fulfilled [42].3 Therein, \(\mathbf{I}\!:=\!\left(I_{1},\ldots,I_{m}\right)\in\mathbb{R}^{m}\) denotes an \(m\)-tuple containing a set of complete and irreducible invariants \(I_{\alpha}(\mathbf{C},\mathcal{S}^{\Box})\) of the symmetry group under consideration. Here, complete means that an arbitrary invariant can be expressed as a function of the elements of \(\mathbf{I}\), e.g., \(\psi\left(\mathbf{I}\right)\). Irreducible means that one invariant \(I_{\alpha}\) of the set \(\mathbf{I}\) cannot be expressed by the other elements of \(\mathbf{I}\), cf. [8]. Note that these invariants have to be chosen such that they do not violate the **polyconvexity** condition.4 The corresponding second Piola-Kirchhoff stress tensor is then given as

Footnote 3: Note that \(\psi\left(\mathbf{F}\right)\) and \(\psi\left(\mathbf{I}\right)\) are different functions, as indicated by the arguments. In favor of a reduced set of symbols, this casual notation is used within this work.

Footnote 4: It should be noted that there are symmetry groups for which no complete set can be found that does not violate polyconvexity [25]. In this case it must be decided which property of the model is preferable. Furthermore, it might be possible that there does not exist a set of invariants to describe the anisotropy of materials with arbitrary microstructures, e.g., 3D-printed composites.

\[\mathbf{T}=2\frac{\partial\psi}{\partial\mathbf{C}}=2\sum_{\alpha=1}^{m}\frac{ \partial\psi}{\partial I_{\alpha}}\frac{\partial I_{\alpha}}{\partial\mathbf{C}} =\mathbf{F}^{-1}\cdot\mathbf{P}\,. \tag{15}\]

Finally, formulating the hyperelastic potential in terms of invariants of the right Cauchy-Green tensor implies **symmetry of the Cauchy stress tensor \(\mathbf{\sigma}\)**, which ensures the conservation of angular momentum [18].

### Special material symmetry classes and specific models

Throughout this work, both isotropic (\(\Box:=\oplus\)) and transversely isotropic (\(\Box:=\parallel\)) material behavior are considered as examples of specific material symmetry groups. In the isotropic case, i.e., when the material response is direction-independent, \(\mathcal{G}=\mathcal{O}(3)\), a complete and irreducible set is given by the three invariants

\[I_{1}:=\operatorname{tr}\mathbf{C},\qquad I_{2}:=\operatorname{tr}(\operatorname {cof}\mathbf{C}),\qquad I_{3}:=\operatorname{det}\mathbf{C}\,. \tag{16}\]

In this case, there are no structural tensors needed, hence \(\mathcal{S}^{\otimes}=\varnothing\) holds. For transversely isotropic materials, a complete functional basis is given by the three isotropic invariants \(I_{1}\), \(I_{2}\), \(I_{3}\) together with

\[I_{4}:=\operatorname{tr}(\mathbf{C}\cdot\mathbf{G}),\qquad I_{5}:=\operatorname{tr}( \operatorname{cof}(\mathbf{C})\cdot\mathbf{G})\,, \tag{17}\]

where \(\mathbf{G}\) denotes the second order transversely isotropic structural tensor [8, 45]. In case the preferred direction is parallel to the \(X_{1}\)-direction, the structural tensor's components are given by

\[\left(G_{KL}\right):=\begin{pmatrix}\beta^{2}&0&0\\ 0&\frac{1}{\beta}&0\\ 0&0&\frac{1}{\beta}\end{pmatrix}, \tag{18}\]

where \(\beta\in\mathbb{R}_{>0}\) is a model parameter [45]. Thus, it holds \(\mathcal{S}^{\parallel}=\left\{\mathbf{G}\right\}\).
#### 2.3.1 Isotropic model

At this point, analytical models have to make an explicit choice of the functional relationship for the hyperelastic potential. While some choices have a strong physical motivation [36], most models are of a heuristic nature - and the reduced flexibility that often goes along with this human choice of functional relationship purely stems from the necessity of an explicit form of the model. To become more specific, we consider the isotropic Neo-Hookean model

\[\psi^{\text{nh}}(I_{1},I_{3})=\frac{1}{2}\left(\mu\left(I_{1}-\ln I_{3}-3 \right)+\frac{\lambda}{2}\left(I_{3}-\ln I_{3}-1\right)\right)\,,\quad\mu=\frac {E}{2(1+\nu)}\,,\quad\lambda=\frac{E\nu}{(1+\nu)(1-2\nu)} \tag{19}\]

with the material parameters \((E,\nu)\) corresponding to the Young's modulus and Poisson's ratio. From this potential, the second Piola-Kirchhoff stress can be derived as

\[\mathbf{T}^{\text{nh}}=\mu\mathbf{1}+\left(\frac{\lambda}{2}-\frac{2\mu+\lambda}{2I_{ 3}}\right)\text{cof}\,\mathbf{C}\,. \tag{20}\]

Note that, e.g., the linear dependency on \(I_{1}\) is a quite restrictive choice. Again, the second law of thermodynamics is fulfilled by constructing the stress as a gradient field, while symmetry of the Cauchy stress tensor, objectivity and material symmetry are fulfilled by using invariants. Then, the normalization conditions are fulfilled by cleverly combining both material parameters and invariants, which becomes evident when setting \(\mathbf{C}=\mathbf{1},I_{1}=3,I_{3}=1\) in Eqs. (19) and (20). By using the term \(-\ln I_{3}\) in combination with the remaining invariants, the volumetric growth condition is also fulfilled. Finally, all mathematical operations preserve the polyconvexity of the invariants.

Surprisingly, to the best of the authors' knowledge, it is not possible to prove analytically that the non-negativity of the elastic potential, i.e., \(\psi^{\text{nh}}(I_{1},I_{3})\geq 0\) for arbitrary physically permissible combinations of invariants \(I_{1}\) and \(I_{3}\), is satisfied by this compressible elastic model in every case. This is due to the fact that \(\psi^{\text{nh}}(I_{1},I_{3})\) is polyconvex, i.e., convex with respect to \(\mathbf{F}\), \(\text{cof}\,\mathbf{F}\) and \(\text{det}\,\mathbf{F}\), but not convex with respect to \(\mathbf{F}\) or \(\mathbf{C}\).5 Consequently, the non-negativity condition has to be numerically checked for admissible deformation states \(\mathbf{C}\), see Sec. 4.

Footnote 5: This is in contrast to the incompressible case [33], where the definition range of the invariants is furthermore significantly more limited since \(J\equiv 1\).

#### 2.3.2 Transversely isotropic model

In order to describe transversely isotropic material behavior, we consider the analytical model

\[\psi^{\text{ii}}(I_{1},I_{2},I_{3},I_{4},I_{5})=\alpha_{1}I_{1}+\alpha_{2}I_{ 2}+\delta_{1}I_{3}-\delta_{2}\ln(\sqrt{I_{3}})+\eta^{*}\left(I_{4}^{\alpha_{4 }}+I_{5}^{\alpha_{4}}\right),\quad\eta^{*}:=\frac{\eta_{1}}{\alpha_{4}\,\text{ tr}(\mathbf{G})^{\alpha_{4}}}\,, \tag{21}\]

as proposed by Schroder et al. [45]. Here, the corresponding second Piola-Kirchhoff stress can be derived as

\[\mathbf{T}^{\text{ti}}=2\bigg{(}\alpha_{1}\mathbf{1}+\alpha_{2}\big{(}I_{1}\mathbf{1}-\mathbf{C }\big{)}+\bigg{(}\delta_{1}I_{3}-\frac{\delta_{2}}{2}\bigg{)}\mathbf{C}^{-1}+2 \eta^{*}I_{4}\mathbf{G}+2\eta^{*}I_{5}\bigg{(}I_{5}\mathbf{C}^{-1}-\text{cof}(\mathbf{C}) \cdot\mathbf{G}\cdot\mathbf{C}^{-1}\bigg{)}\bigg{)}\,. \tag{22}\]

The constitutive model thus includes the material parameters \((\beta,\alpha_{1},\alpha_{2},\delta_{1},\delta_{2},\alpha_{4},\eta_{1})\), with restrictions given in [45]. Herein, all conditions introduced in Sec. 2.2 are satisfied by the analytical model. However, as stated in Sec. 2.3.1, the non-negativity of the strain energy density has to be verified numerically.

To close, while the constitutive conditions on hyperelasticity discussed in this section have a sound physical and mathematical basis, the restricted functional relationships that most analytical models choose have not. This limitation can be circumvented by using PANNs.
Physics-augmented neural network constitutive model

In the previous section, we discussed the constitutive conditions of hyperelasticity. Now, instead of choosing an analytical formulation for the potential, like in Secs. 2.3.1 and 2.3.2, we aim to exploit the excellent approximation properties of FFNNs [1, 26]. In the following, the overall PANN model is introduced, which satisfies all further introduced conditions, except for the non-negativity of the energy, in an exact way, see Sec. 2.2.

### Basic conditions

By constructing the stress as the gradient of a potential predicted directly by a FFNN, **thermodynamic consistency** of the model is fulfilled. If the potential is additionally formulated in terms of invariants, a **symmetric Cauchy stress tensor**6 is automatically obtained and the condition of **objectivity** as well as the condition of **material symmetry** are satisfied. At this stage, the neural network already fulfills _basic conditions_ by construction, i.e., thermodynamic consistency, objectivity, material symmetry and symmetry of the stress tensor.

Footnote 6: It should be noted that a FFNN directly mapping from \(\mathbf{F}\) to \(\mathbf{P}\) does not necessarily yield an elastic model which accounts for \(\mathbf{\sigma}=\mathbf{\sigma}^{T}\). The same restriction holds for a model which maps from \(\mathbf{F}\) to \(\mathbf{\psi}(\mathbf{F})\). E.g., in Naumann and Ihlemann [38], a pseudo elastic model is investigated which violates the integrability condition, so that the stress is no longer a gradient field.

We assume that the anisotropic hyperelastic material behavior can be described by the irreducible and independent set of invariants \(\mathbf{I}:=(I_{1},\ldots,I_{m})\) with \(I_{\mathbf{\beta}}(\mathbf{C},\mathbf{\mathcal{S}}^{\Box})\), \(\beta\in\mathbb{N}_{\leq m}\). However, it may be necessary to incorporate further invariants \(I_{\mathbf{\gamma}}^{*}(\mathbf{I}),\gamma\in\mathbb{N}_{\leq A}\) with \(A\in\mathbb{N}_{\geq 0}\) into the argument list of the predicted potential, e.g., in order to fulfill additional physical conditions or to increase the approximation quality of the predictions [19, 24]. By using the extended set of invariants \(\mathbf{I}^{*}:=(I_{1},\ldots,I_{m},I_{1}^{*},\ldots,I_{A}^{*})\) as inputs of a FFNN with scalar-valued output which is then taken as a hyperelastic potential \(\psi(\mathbf{I}^{*})\), the flexibility of the resulting model exceeds the one of analytical formulations by far. Furthermore, in a lot of practical applications it is sufficient to restrict the network architecture to only one hidden layer containing \(N^{\text{NN}}\) neurons. Using the activation function \(\mathcal{F}:\mathbb{R}\to\mathbb{R}\), which needs to be continuously differentiable twice, the most simple case of an invariant-based potential using only one hidden layer is given by

\[\psi^{\text{NN},\Box}(\mathbf{I}^{*}):=\sum_{\alpha=1}^{N^{\text{NN}}}W_{\alpha} \mathcal{F}\left(\sum_{\beta=1}^{m}w_{\alpha\beta}I_{\beta}+\sum_{\gamma=1}^ {A}w_{\alpha\gamma}^{*}I_{\gamma}^{*}+b_{\alpha}\right)\,, \tag{23}\]

where \(W_{\alpha},w_{\alpha\beta},w_{\alpha\gamma}^{*}\) and \(b_{\alpha}\) denote weights and bias values, respectively, which together form the set of parameters \(\mathbf{\mathcal{P}}\in\mathbb{R}^{P}\), with \(P\in\mathbb{N}\) denoting the total number of parameters, to be optimized in the calibration process to fit a given dataset. Note that none of the methods introduced in this paper is restricted to FFNNs with one hidden layer, but can directly be applied for multilayered network architectures, cf. Appendix A. Nevertheless, we use only one layer here to better illustrate the derivations. Then, the stress prediction is done by applying

\[\mathbf{T}^{\text{NN},\Box}=2\sum_{\alpha=1}^{m}\frac{\partial\psi^{\text{NN}, \Box}}{\partial I_{\alpha}}\frac{\partial I_{\alpha}}{\partial\mathbf{C}}+2\sum_{ \gamma=1}^{A}\sum_{\beta=1}^{m}\frac{\partial\psi^{\text{NN},\Box}}{\partial I _{\gamma}^{*}}\frac{\partial I_{\gamma}^{*}}{\partial I_{\beta}}\frac{\partial I _{\beta}}{\partial\mathbf{C}}. \tag{24}\]

This special choice of _input quantities_ - namely, invariants - is the first way of including physics into the NN.

### Polyconvexity

For ensuring **polyconvexity** of the potential, it is necessary to use polyconvex invariants of the symmetry group under consideration. In addition, the _network architecture_ must be adapted in a certain way. By using
a convex and non-decreasing activation function \(\mathcal{F}\) and non-negative weights, polyconvexity of the overall potential \(\psi^{\mathrm{NN,\Box}}(\mathbf{I}^{*})\) is ensured [24]. Here, the _Softplus_ activation function \(\mathcal{SP}(x):=\log(1+\exp(x))\in\mathcal{C}^{\infty}\) is applied, which is convex and non-decreasing, which overall leads to the conditions

\[W_{\alpha},w_{\alpha\beta},w^{*}_{\alpha\gamma}\in\mathbb{R}_{\geq 0},b_{ \alpha}\in\mathbb{R}\quad\forall\alpha\in\mathbb{N}_{\leq N^{\mathrm{NN}}}, \beta\in\mathbb{N}_{\leq m},\gamma\in\mathbb{N}_{\leq A}\,. \tag{25}\]

**Remark 3.1** (**Input convex neural networks (ICNNs))**.: This special kind of FFNN, namely, one with its output being convex in its input arguments, is referred to as ICNN [3]. In order to understand what makes this kind of NNs _convex_, we take a step back and consider the univariate function

\[f\colon\mathbb{R}\to\mathbb{R},\quad x\mapsto f(x):=(g\circ h)(x)\,, \tag{26}\]

where \(f\) is composed of two functions \(g,h\colon\mathbb{R}\to\mathbb{R}\). Given that all of the above functions are twice continuously differentiable, convexity of \(f(x)\) in \(x\) is equivalent to the non-negativity of its second derivative

\[f^{\prime\prime}(x)=(g^{\prime\prime}\circ h)(x)\,h^{\prime}(x)^{2}+(g^{ \prime}\circ h)(x)\,h^{\prime\prime}(x)\geq 0\,. \tag{27}\]

A sufficient, albeit not necessary condition for this is that the innermost function \(h\) is convex (\(h^{\prime\prime}\geq 0\)), while \(g\) is convex and non-decreasing (\(g^{\prime\prime}\geq 0\) and \(g^{\prime}\geq 0\)). This can be generalized to arbitrarily many function compositions, where the innermost function must be convex, while every following function must be convex and non-decreasing. Transferred to FFNNs, which can be seen as compositions of multiple vector-valued functions, this generalizes to the condition that in an ICNN the first hidden layer must be node-wise convex, while every subsequent layer must be node-wise convex and non-decreasing, cf. [24, Appendix A].

When applied to hyperelasticity, the architecture of ICNNs has to be further adapted: as the invariants are nonlinear functions of the arguments defined in the polyconvexity condition, i.e., \((\mathbf{F},\operatorname{cof}\mathbf{F},\operatorname{det}\mathbf{F})\), cf. Eq. (10), they are the innermost function acting on the arguments of the polyconvexity condition. Thus, already the first hidden layer has to be convex and _non-decreasing_. The only exception is the use of \(J=\sqrt{I_{3}}\) as an additional invariant (instead of \(I_{3}\)), since \(J\) is an argument of the polyconvexity condition. For this reason, the activation function acting on \(J\) must only be convex and not necessarily non-decreasing. This is pragmatically taken into account by including the additional invariant \(I_{1}^{*}:=-2J\) in the set of invariants, which is furthermore essential to represent negative stresses at all [24]. Note that there are different ways of constructing polyconvex neural networks, e.g., [48] or [7]. However, the simple structure and excellent flexibility of ICNNs makes them a very natural choice for this task. For a more extensive introduction to polyconvex neural networks and explicit proofs, see [24].

### Growth and normalization conditions

Finally, **growth and normalization conditions** remain to be included in the neural network, ensuring a physically sensible stress behavior of the model. For this, growth and normalization terms are added to the original potential \(\psi^{\mathrm{NN,\Box}}(\mathbf{I}^{*})\) of Eq. (23). Then, with \(\Box\) denoting the symmetry group under consideration, the overall PANN model, given as

\[\psi^{\mathrm{PANN,\Box}}(\mathbf{I}^{*}):=\psi^{\mathrm{NN,\Box}}(\mathbf{I}^{*})+ \psi^{\mathrm{stress,\Box}}(J,I_{4},\dots,I_{m})+\psi^{\mathrm{energy,\Box}}+ \psi^{\mathrm{growth}}(J)\,, \tag{28}\]

fulfills all constitutive conditions introduced in Sec. 2.2, except for the non-negativity of the energy, in an exact way. In Fig. 2, the overall structure of the PANN model is exemplarily illustrated for one hidden layer. Consequently, the corresponding stress is given by the expression

\[\mathbf{T}^{\mathrm{PANN,\Box}}=2\frac{\partial\psi^{\mathrm{PANN,\Box}}}{ \partial\mathbf{C}}=\mathbf{T}^{\mathrm{NN,\Box}}+\mathbf{T}^{\mathrm{stress,\Box}}+\mathbf{T}^ {\mathrm{energy,\Box}}+\mathbf{T}^{\mathrm{growth}}\,\,. \tag{29}\]
One way to fulfill the volumetric **growth condition** are coercive functions. However, since ICNNs are not necessarily coercive, they are not suited to fulfill this condition and therefore an analytical term

\[\psi^{\text{growth}}(J):=\left(J+\frac{1}{J}-2\right)^{2} \tag{30}\]

is introduced, which is chosen in such a way that polyconvexity is not violated. This leads to the corresponding stress contribution

\[\mathbf{T}^{\text{growth}}=2\left(J+\frac{1}{J}-2\right)\left(1-\frac{1}{J^{2}} \right)J\mathbf{C}^{-1}\,. \tag{31}\]

Another way to fulfill the volumetric growth condition is further adapting the network architecture, cf. [19], however, using an analytical term is more straightforward.

The correction term for the **energy normalization** is given by

\[\psi^{\text{energy},\Box}:=-\psi^{\text{NN},\Box}(\mathbf{T}^{*}) \Big{|}_{\mathbf{C}=\mathbf{1}}\in\mathbb{R}\,. \tag{32}\]

Since \(\psi^{\text{energy},\Box}\) is a constant, it holds \(\mathbf{T}^{\text{energy},\Box}=\mathbf{0}\). Together with the stress normalization term denoted by \(\psi^{\text{stress},\Box}(J,I_{4},\dots,I_{m})\), meaning vanishing gradients of the potential for the undeformed state, we ensure

\[\psi^{\text{PANN},\Box}(\mathbf{T}^{*})\Big{|}_{\mathbf{C}=\mathbf{1}}=0 \tag{33}\]

as well as

\[\mathbf{T}^{\text{PANN},\Box}\Big{|}_{\mathbf{C}=\mathbf{1}}=2\left(\sum_{ \alpha=1}^{m}\frac{\partial\psi^{\text{PANN},\Box}}{\partial I_{\alpha}} \frac{\partial I_{\alpha}}{\partial\mathbf{C}}+\sum_{\gamma=1}^{A}\sum_{\beta=1}^ {m}\frac{\partial\psi^{\text{PANN},\Box}}{\partial I_{\gamma}^{*}}\frac{ \partial I_{\gamma}^{*}}{\partial I_{\beta}}\frac{\partial I_{\beta}}{ \partial\mathbf{C}}\right)\Bigg{|}_{\mathbf{C}=\mathbf{1}}=\mathbf{0}. \tag{34}\]

Thus, the potential has a local minimum exactly for the undeformed state with \(\mathbf{C}=\mathbf{1}\), where for this case we define the invariants as \(I_{\alpha}^{0}:=I_{\alpha}(\mathbf{C}=\mathbf{1})\) and \(I_{\beta}^{*0}:=I_{\beta}^{*}(\mathbf{C}=\mathbf{1})\). However, due to the fact that \(\psi^{\text{PANN},\Box}(\mathbf{T}^{*})\) is polyconvex, i.e., convex with respect to \(\mathbf{F}\), \(\text{col}\,\mathbf{F}\) and \(\text{det}\,\mathbf{F}\), but not convex with respect to \(\mathbf{F}\) or \(\mathbf{C}\), the **non-negativity of the energy**, i.e., \(\psi^{\text{PANN},\Box}(\mathbf{T}^{*})\geq 0\) does not automatically follow. As already stated in Sects. 2.3.1 and 2.3.2, a numerical test for admissible deformation states is needed to prove the fulfillment of this condition, see Sec. 4.

The polyconvex stress correction \(\psi^{\text{stress},\Box}(J,I_{4},\dots,I_{m})\) depends on the symmetry group under consideration and is now further discussed.

Figure 2: Illustration of the PANN based constitutive model for the material symmetry group \(\Box\) under consideration. Note that the hidden-layer (yellow) of the NN may be multilayered.
#### 3.3.1 Isotropic normalization term

In the isotropic case with \(\mathbf{\varGamma}^{*}:=(I_{1},I_{2},I_{3},I_{1}^{*})\) the normalization term is given by

\[\psi^{\text{stress,\tiny\textregistered}}(J):=-\mathfrak{n}(J-1)\, \tag{35}\]

where the constant

\[\mathfrak{n}:=2\left(\frac{\partial\psi^{\text{NN,\tiny\textregistered}}}{ \partial I_{1}}+2\frac{\partial\psi^{\text{NN,\tiny\textregistered}}}{\partial I_{ 2}}+\frac{\partial\psi^{\text{NN,\tiny\textregistered}}}{\partial I_{3}}+\frac{ \partial\psi^{\text{NN,\tiny\textregistered}}}{\partial I_{1}^{*}}\frac{\partial I _{1}^{*}}{\partial I_{3}}\right)\Bigg{|}_{\mathbf{C}=\mathbf{1}}\in\mathbb{R} \tag{36}\]

is a weighted sum of derivatives of the ICNN potential with respect to the invariants for the undeformed state \(\mathbf{C}=\mathbf{1}\). The corresponding stress contribution is given by

\[\mathbf{T}^{\text{stress,\tiny\textregistered}}=-\mathfrak{n}\mathbf{\varGamma}^{-1}. \tag{37}\]

and leads to the fact that the undeformed state is stress-free by construction.

This definition of the correction term comprises two ideas: first of all, this approach preserves poly-convexity of the potential, since the additional term is a linear function in \(J\), which is an invariant quantity included in the arguments of the polyconvexity condition, cf. Eq. (10). Furthermore, the partial derivatives

\[\left.\frac{\partial I_{\alpha}}{\partial\mathbf{C}}\right|_{\mathbf{C}=\mathbf{1}}=\xi \mathbf{1}\quad\text{and}\quad\left.\frac{\partial I_{1}^{*}}{\partial\mathbf{C}} \right|_{\mathbf{C}=\mathbf{1}}=\eta\mathbf{1} \tag{38}\]

of the isotropic invariants with respect to \(\mathbf{C}\) for the undeformed state are multiples of the identity tensor \(\mathbf{1}\) for all \(\alpha\in\mathbb{N}_{\leq 3}\) with constants \(\xi,\eta\in\mathbb{R}\). Hence, it can normalize the stress for the undeformed state to zero, which becomes evident when setting \(\mathbf{C}=\mathbf{1}\) in Eq. (37). Further details are provided in Appendix B, where all necessary tensor derivatives are given in analytical form.

#### 3.3.2 Transversely isotropic normalization term

In the case of transverse isotropy with \(\mathbf{\varGamma}^{*}:=(I_{1},I_{2},I_{3},I_{4},I_{5},I_{1}^{*})\) the partial derivatives of the additional invariants \(I_{4},I_{5}\) with respect to \(\mathbf{C}\) for the undeformed state are not multiples of the identity anymore but include the structural tensor \(\mathbf{G}\), e.g.,

\[\left.\frac{\partial I_{4}}{\partial\mathbf{C}}\right|_{\mathbf{C}=\mathbf{1}}=\mathbf{G} \quad\text{and}\quad\left.\frac{\partial I_{5}}{\partial\mathbf{C}}\right|_{\mathbf{C }=\mathbf{1}}=I_{5}\mathbf{1}-\mathbf{G}. \tag{39}\]

Thus, the correction term

\[\psi^{\text{stress,\tiny\textregistered}}(J,I_{4},I_{5}):=-\mathfrak{o}(J-1)+ \mathfrak{p}(I_{4}-I_{4}^{0})+\mathfrak{q}(I_{5}-I_{5}^{0}) \tag{40}\]

is introduced, with \(I_{4}^{0}=I_{5}^{0}=\mathfrak{tr}\,\mathbf{G}\). Again, the constant

\[\mathfrak{o}:=2\left(\frac{\partial\psi^{\text{NN,\tiny\textregistered}}}{ \partial I_{1}}+2\frac{\partial\psi^{\text{NN,\tiny\textregistered}}}{\partial I _{2}}+\frac{\partial\psi^{\text{NN,\tiny\textregistered}}}{\partial I_{3}}+\frac{ \partial\psi^{\text{NN,\tiny\textregistered}}}{\partial I_{1}^{*}}\frac{\partial I _{1}^{*}}{\partial I_{3}}+\frac{\partial\psi^{\text{NN,\tiny\textregistered}}}{ \partial I_{5}}\,\mathfrak{tr}\,\mathbf{G}+\mathfrak{q}\,\mathfrak{tr}\,\mathbf{G} \right)\Bigg{|}_{\mathbf{C}=\mathbf{1}}\in\mathbb{R} \tag{41}\]

is a weighted sum of the derivatives of the ICNN potential with respect to the invariants for the undeformed state \(\mathbf{C}=\mathbf{1}\). Furthermore, with the ReLU-function denoted as \(\mathcal{R}\mathcal{L}\), the non-negative constants

\[\mathfrak{p}:=\mathcal{R}\mathcal{L}(-x)\in\mathbb{R}_{\geq 0},\quad\mathfrak{q}:= \mathcal{R}\mathcal{L}(x)\in\mathbb{R}_{\geq 0} \tag{42}\]
are defined with the argument

\[x:=\left(\frac{\partial\psi^{\text{NN},\parallel}}{\partial I_{4}}-\frac{\partial \psi^{\text{NN},\parallel}}{\partial I_{5}}\right)\bigg{|}_{\mathbf{C}=\mathbf{1}}\,. \tag{43}\]

Due to the non-negativity of \(\mathsf{o}\), \(\mathsf{p}\) as well as the polyconvexity of \(I_{4},I_{5}\), the correction term is again polyconvex. It should be emphasized again that in the correction terms no assumptions on the number of hidden layers in the neural networks are made, and the approach can directly be applied to multilayered network architectures. Overall, in the transversely isotropic case, the stress contribution of the normalization term is given by

\[\mathbf{T}^{\text{stress},\parallel}=-\mathsf{o}\mathbf{JC}^{-1}+2\mathsf{p}\mathbf{G}+2 \mathsf{a}\left(I_{5}\mathbf{C}^{-1}-\operatorname{cof}(\mathbf{C})\cdot\mathbf{G}\cdot \mathbf{C}^{-1}\right)\,, \tag{44}\]

where all necessary tensor derivatives of the Appendix B are integrated.

Note that in the definitions of all correction terms, even though \(\psi^{\text{stress},\Box}(J,I_{4},I_{5})\) and \(\psi^{\text{energy}}\) depend on the evaluation of \(\psi^{\text{NN},\Box}(\mathbf{T}^{*})\) and its partial derivatives at \(\mathbf{C}=\mathbf{1}\), no assumptions on the number of hidden layers in the ICNN \(\psi^{\text{NN},\Box}(\mathbf{T}^{*})\) are made, and the approaches can directly be applied to multilayered network architectures.

### Model calibration

Finally, the model has to be calibrated to data of a specific material including the set of structural tensors, summarized in \(\mathcal{S}^{\Box}\), which corresponds to the material symmetry group \(\Box\) under consideration. Throughout this work, datasets of the form

\[\mathcal{D}=\left\{\left(\,{}^{1}\mathbf{C},\,{}^{1}\mathbf{T}\,\right),\left(\,{}^{2 }\mathbf{C},\,{}^{2}\mathbf{T}\,\right),\dots\right\} \tag{45}\]

consisting of strain-stress tuples in terms of right Cauchy-Green deformation and second Piola-Kirchhoff stress tensors are used. In order to examine the generalization of a model, i.e., its prediction of general load cases, it is essential to evaluate it on data not seen in the calibration process. Thus, the overall dataset \(\mathcal{D}\) is split into a calibration dataset \(\mathcal{D}_{\text{c}}\) and a test dataset \(\mathcal{D}_{\text{t}}\) with \(\mathcal{D}_{\text{c}}\cap\mathcal{D}_{\text{t}}=\varnothing\). Then, after the calibration of the model on the dataset \(\mathcal{D}_{\text{c}}\), its predictions can be evaluated on \(\mathcal{D}_{\text{t}}\). Only when the model is able to predict the load cases of the test dataset, and only if the load cases included in the test dataset are sufficiently general, it can be assumed that the model generalizes well and can predict the stress for arbitrary deformations.

The next step of the model calibration is the choice of its hyperparameters, here, the number of hidden layers and the nodes they contain. For this, a "sufficiently large" number of layers and nodes should be chosen, so that the model is flexible enough for the material behavior under consideration. For general NNs, overfitting [1, 26], i.e., a too accurate interpolation of the calibration data which results in a bad prediction of general data, is an issue. However, as we will demonstrate, for the PANN proposed in this work, the inclusion of physics provides the model with a pronounced mathematical structure, which makes it less prone to overfitting. Then, with a fixed model architecture, the model parameters \(\mathbf{\mathcal{P}}\) must be optimized, i.e., its weights and biases. To calibrate the model parameters \(\mathbf{\mathcal{P}}\) on the dataset \(\mathcal{D}_{\text{c}}\), the loss function defined in terms of the mean squared error

\[\mathcal{MSE}^{\Box}(\mathbf{\mathcal{P}})=\frac{1}{|\mathcal{D}_{\text{c}}|}\sum _{i=1}^{|\mathcal{D}_{\text{c}}|}\left\|i\mathbf{T}-\mathbf{T}^{\text{model},\Box}(i \mathbf{C};\mathbf{\mathcal{P}})\right\|^{2} \tag{46}\]

is minimized, where \(|\mathcal{D}_{\text{c}}|\) denotes the number of tuples in \(\mathcal{D}_{\text{c}}\) and \(\left\|\cdot\right\|\) denotes the Frobenius norm. Note that the hyperelastic potential \(\psi^{\text{PANN}}(\mathbf{T}^{*})\) is calibrated only through its gradients, i.e., the corresponding stress tensor, which is referred to as Sobolev training [53]. For the optimization process, the _SLSQP optimizer_ (Sequential Least Squares Programming) is applied [40]. According the proposed model, loss function and optimizer have no influence on the underlying physics, and they could also be chosen otherwise. An implementation of the described workflow is realized using _Python, TensorFlow_ and _SciPy_.
**Remark 3.2**.: In the calibration process, it is important to already account for the growth term \(\psi^{\rm growth}(J)\) included in the model \(\psi^{\rm PANN,\Box}(\mathbf{T}^{*})\), cf. Eq. (28), as in particular for calibration data including high volumetric compressions, it can have a considerable influence on the model. On the other side, the _stress and energy normalization_ terms may not necessarily have to be included in the calibration process: the energy normalization term \(\psi^{\rm energy}\) can be added to the PANN model after calibration, as it is a constant it does not influence the stress prediction at all. Then, even when the model is calibrated without the stress correction term \(\psi^{\rm stress,\Box}(J,I_{4},\ldots,I_{m})\), it can learn through suitable data how to approximate the stress-free reference configuration [24]. Given that the approximation of the stress-free reference configuration is good enough, the derivatives of the potential w.r.t. the invariants become very small at the identity \(\mathbf{C}=\mathbf{1}\), and consequently, the constants \(\mathfrak{u},\mathfrak{o},\mathfrak{p}\) on which the stress correction depends become very small, cf. Eqs. (36,41,42). Then, both the correction term \(\psi^{\rm stress,\Box}(J,I_{4},\ldots,I_{m})\) and its gradient \(\mathbf{T}^{\rm stress,\Box}\) approximately vanish, and adding the stress correction after the model calibration has a negligibly small influence on the overall model behavior. While it is indeed possible to fulfill the normalization conditions in good approximation by learning them through suitable data, they are only fulfilled exactly when incorporated in the model formulation, e.g., by normalization terms.

## 4 Numerical examples

After the introduction of the PANN hyperelastic constitutive model in the former section, the ability of our approach is now demonstrated by several numerical examples and compared to approaches with reduced physical foundation. We start with simple stress states and analyze interpolation and extrapolation behavior, also for perturbed data. Thereafter, we go on to complex multiaxial deformation states. Finally, the trained PANN model is applied within a finite element simulation.

### Simple stress-strain states

In this subsection, three types of models, namely a network \(\mathbf{P}^{\rm simple}(\mathbf{F})\) directly mapping from \(\mathbf{F}\) to \(\mathbf{P}\), a network \(\psi^{\rm NN,\emptyset}(\mathbf{I}^{*})\) accounting for the basic conditions, and the PANN \(\psi^{\rm PANN,\emptyset}(\mathbf{I}^{*})\), are compared with respect to their interpolation and extrapolation capabilities for simple stress-strain states. Thereby, the first model, given by

\[P^{\rm simple}_{kL}(\mathbf{F})\coloneqq B_{kL}+\sum_{\alpha=1}^{N^{\rm NN}}W_{ \alpha kL}\,\mathcal{S}\mathcal{P}\left(\sum_{i=1}^{3}\sum_{J=1}^{3}w_{\alpha iJ }F_{iJ}+b_{\alpha}\right)\text{ with }B_{kL},W_{\alpha kL},w_{\alpha iJ},b_{ \alpha}\in\mathbb{R}\,, \tag{47}\]

does not take into account any of the introduced conditions, the second, which is defined by Eq. (23), respects for thermodynamic consistency, symmetry of \(\mathbf{\sigma}\), objectivity and material symmetry, and the PANN fulfills the conditions of the previous model and additionally includes polyconvexity, growth condition, as well as energy and stress normalization, see Eq. (28).

In the following examples, the architectures of the three models are set to one hidden layer with \(N^{\rm NN}\coloneqq 4\) neurons. Furthermore, the softplus activation function is used for all models. The calibration of the model accounting for basic conditions and the PANN is performed according to Sec. 3.4. In contrast, naturally, the training of the \(\mathbf{F}\)-\(\mathbf{P}\) model is performed with data for \(\mathbf{F}\) and \(\mathbf{P}\). Furthermore, the Adam optimizer is used for the training of this model. The training data for this study are generated by using the isotropic Neo-Hooke model, cf. Eq. (20), where the constants \((E,\nu)\coloneqq(1\,\mathrm{MPa},0.3)\) are chosen.

#### 4.1.1 Interpolation behavior

In order to analyze the interpolation behavior of the three models, we investigate three different sets of uniaxial stress states: ideal data, offset data, as well as noisy data.
Ideal dataThe three models are first trained on analytical data from a uniaxial tension and compression test in \(X_{1}\)-direction, cf. Fig. 3(a1). The principal stretch \(\lambda_{1}\in\mathbb{R}_{>0}\) in \(X_{1}\)-direction is prescribed and the principal stretch \(\lambda_{2}\in\mathbb{R}_{>0}\) has to be calculated in such a way that the stresses \(T_{22}^{\oplus},T_{33}^{\oplus}\) vanish, so that the components of the deformation tensor and the stress tensor result in

\[(C_{KL})^{\text{uniaxial}}:=\begin{pmatrix}\lambda_{1}^{2}&0&0\\ 0&\lambda_{2}^{2}&0\\ 0&0&\lambda_{2}^{2}\end{pmatrix},\quad(T_{KL})^{\text{uniaxial},\oplus}= \begin{pmatrix}T_{11}^{\oplus}&0&0\\ 0&0&0\\ 0&0&0\end{pmatrix}. \tag{48}\]

Accordingly, this leads to a uniaxial stress state, which is commonly applied in experimental investigations. Only 30 tuples \((\,{}^{i}\mathbf{C},\mathbf{T}^{\oplus})\) with \(0.8\leq\lambda_{1}\leq 2\) are used as training data for the models. Note that the tuple for \(\lambda_{1}=1\) occurs twice in the data set, since this point is assigned to both the compression and tension ranges. This has no influence on the approximation behavior of the PANN model \(\psi^{\text{PANN},\oplus}(\mathbf{I}^{*})\), since this model is normalized by construction. However, for both the \(\mathbf{F}\)-\(\mathbf{P}\) model and the basic conditions model \(\psi^{\text{NN},\oplus}(\mathbf{I}^{*})\), the undeformed state is weighted more heavily in the data set.

After calibration, the stress-stretch curves \(T_{11}^{\text{simple}}(\lambda_{1})\), \(T_{11}^{\text{NN},\oplus}(\lambda_{1})\),and \(T_{11}^{\text{PANN},\oplus}(\lambda_{1})\) of the three models are available in analytical form and are compared with the training data in Fig. 3(a1). As expected, it can be

Figure 3: Predicted stress-stretch curves \(T_{11}^{\oplus}(\lambda_{1})\) of the PANN, the NN fulfilling the basic conditions, as well as the \(\mathbf{F}\)-\(\mathbf{P}\) model with one hidden layer containing \(N^{\text{NN}}:=4\) neurons. The trained models are compared to data from a uniaxial tension/compression test in \(X_{1}\)-direction for: (a1) ideal isotropic Neo-Hooke data, (b1) data with offset, and (c1) noisy data. In (a2)–(c2), the corresponding energies of the models are shown. For reasons of improved comparability, energy normalization is applied for \(\psi^{\text{NN},\oplus}(\mathbf{I}^{*})\) which, however, has no influence on the stress prediction.
seen that all three NN-based models are able to perfectly approximate the ideal dataset which is also shown by the MSEs given in Tab. 1. The basic conditions model as well as the PANN also reproduce the energy \(\psi^{\otimes}(\lambda_{1})\) with high accuracy although it was not trained directly, see Fig. 3(a2). Thus, incorporating physics does not degrade the model prediction quality. In addition, we would like to point out that only the PANN model really fulfills the normalization condition, i.e., \(T_{11}^{\text{PANN},\oplus}(\lambda_{1}=1)=0\), exactly.

Offset dataThe next step is to investigate the flexibility of the models for non-ideal data. For this, the stress components of the ideal dataset are shifted according to

\[(T_{KL})^{\text{offset},\oplus}:=(T_{KL})^{\text{uniaxial},\oplus}+\begin{pmatrix} 100&0&0\\ 0&0&0\\ 0&0&0\end{pmatrix}\text{kPa}\,, \tag{49}\]

so that the calibration dataset is no longer normalized for the undeformed state, i.e., \(T_{11}^{\text{offset},\oplus}(\lambda_{1}=1)\neq 0\). Thereby, the normalization condition included into the PANN model is of particular interest now.

As can be seen in Fig. 3(b1), the \(\mathbf{F}\)-\(\mathbf{P}\) model \(\mathbf{P}^{\text{simple}}(\mathbf{F})\) perfectly approximates the training data again, see also the MSE given in Tab. 1. However, this is only possible since this simple model does not know about the existence of a potential and is thus not really an hyperelastic one. In contrast, the basic conditions model \(\psi^{\text{NN},\oplus}(\mathbf{I}^{*})\) is not able to reproduce the training data well near \(\lambda=1\), which is due to the fact that it cannot violate the material symmetry and that the stress results from the derivative of an elastic potential. However, for the undeformed state a non-zero stress is predicted, which is in contradiction to the expectation of a standard elastic model.

As can be seen in Fig. 3(b1), due to the model approach according to Sec. 3, the normalization condition is only exactly fulfilled by the PANN model \(\psi^{\text{PANN},\oplus}(\mathbf{I}^{*})\) after model calibration, i.e., \(T_{11}^{\text{PANN},\oplus}(\lambda_{1}=1)=0\). Although it becomes evident from the proofs in Sec. 3, it is now demonstrated that this important condition is fulfilled at all times by construction, even for the case when the training data have an offset. Similar to the basic conditions model, due to the stress normalization and the polyconvexity, the PANN model does not approximate the data points in the neighborhood of \(\lambda_{1}=1\) well - but in a physically meaningful way - resulting in an inflection point within the stress-stretch curve. Outside this neighborhood, the approximation of the data points is close to perfect. Thereby, the non-convexity of \(T_{11}\) in \(\lambda_{1}\) should not be mistaken as a violation of the polyconvexity condition. The proposed NN model is polyconvex by construction, cf. Remark 3.1, and furthermore, polyconvexity does not pose restrictions on the dependency of \(\mathbf{T}\) in \(\mathbf{F}\).

Due to the described physical restrictions, the MSEs of the basic conditions model and the PANN are several orders of magnitude larger compared to the \(\mathbf{F}\)-\(\mathbf{P}\) model for the offset data, cf. Tab. 1. The elastic energies given by \(\psi^{\text{NN},\oplus}(\mathbf{I}^{*})\) and \(\psi^{\text{PANN},\oplus}(\mathbf{I}^{*})\) are shown in Fig. 3(b2). Thereby, no significant difference between both models occurs. Note that the training data are not included in this plot because there is no energy available for the offset dataset under consideration.

Noisy dataIn the last step, the approximation behavior of the three models with respect to noisy data is investigated. For this purpose, the stress components of a larger ideal dataset containing 100 tuples \((\mathbf{{}^{i}}\mathbf{C},\mathbf{{}^{i}}\mathbf{T}^{\otimes})\) are shifted according to

\[(T_{KL})^{\text{noisy},\otimes}:=(T_{KL})^{\text{uniaxial},\oplus}+ \begin{pmatrix}\xi&0&0\\ 0&0&0\\ 0&0&0\end{pmatrix}\quad\text{with}\quad\xi\sim\mathcal{N}(\mu,\sigma^{2})\,, \tag{50}\]

where \(\xi\) describes a normally distributed noise with mean value \(\mu=0\,\text{kPa}\) and standard deviation \(\sigma=50\,\text{kPa}\).

Although the training data are overlaid with a strong Gaussian noise, as it can be seen in Fig. 3(c1), all three calibrated models give a result that appears to be physically reasonable. To be more precise, the
predicted stress-stretch curves are monotonously increasing and the potentials \(\psi^{\text{PANN},\oplus}(\lambda_{1})\) and \(\psi^{\text{NN},\oplus}(\lambda_{1})\) shown in Fig. 3(c2) are convex in \(\lambda_{1}\), which is both to be expected for uniaxial stress, they are smooth, and do not show any markedly oscillations. It can be seen that the stress prediction is approximately zero for the undeformed state in this case, see Fig. 3(c1). However, only the PANN model satisfies all relevant conditions in an exact manner, which also applies to stress-strain states not included in the training. The advantage of this property will be revealed in the next study. The _F-P_ model even violates the balance of angular momentum, since \(\mathbf{T}^{\text{simple}}\) and \(\mathbf{\sigma}^{\text{simple}}\) are not symmetrical.

The MSEs for all three models are again given in Tab. 1. Due to the imperfect data, the values are increased by several orders of magnitude compared to the ideal case.

#### 4.1.2 Extrapolation behavior

Now, the extrapolation behavior of the compared models is analyzed by considering three load cases: uniaxial tension/compression, biaxial tension/compression, as well as simple shear. For the training of \(\mathbf{P}^{\text{simple}}(\mathbf{F})\), \(\psi^{\text{NN},\oplus}(\mathbf{I}^{*})\), and \(\psi^{\text{PANN},\oplus}(\mathbf{I}^{*})\) only uniaxial stress states and the corresponding deformations within the narrow range \(0.8\leq\lambda_{1}\leq 1.1\) have been used here, cf. Eq. (48), which are stored in the dataset \(\mathcal{D}^{\text{uniaxial},\oplus}\) with \(|\mathcal{D}^{\text{uniaxial},\oplus}|=15\). Again, note that the tuple for \(\lambda_{1}=1\) occurs twice in the data set, since this point is assigned to both the compression and tension ranges. Given that the PANN model \(\psi^{\text{PANN},\oplus}(\mathbf{I}^{*})\) is normalized by construction, this does not affect the approximation behavior. In contrast, for both the _F-P_ model and the basic conditions model \(\psi^{\text{NN},\oplus}(\mathbf{I}^{*})\), the undeformed state is weighted more heavily in the data set.

Uniaxial stress statesAs one can see in Fig. 4(a1), all models again perfectly approximate the training data which is also evident from the MSEs given in Tab. 2. However, the _F-P_ model fails immediately when it has to extrapolate. The model \(\psi^{\text{NN},\oplus}(\mathbf{I}^{*})\), in contrast, is very good at extrapolating up to a stretch of \(\lambda_{1}=4\), which is really impressing from the author's point of view. Likewise, the related energy shown in Fig. 4(a2) is very well reproduced. Thus, the advantage of invariant-based approaches [12, 19, 24, 25, 30, 32, 33, 48] which approximate the energy and not directly the stress is particularly evident here. Even more impressive is the result for the PANN model, which is able to reproduce the data almost perfectly despite extrapolation. This is also evident from the MSE for the extrapolated data, which is two orders of magnitude lower for the PANN model than for the basic conditions model, cf. Tab. 2. The observed improvement results from the insertion of the additional physical principles, namely polyconvexity, growth condition, as well as energy and stress normalization, into the PANN model.

Biaxial stress states and simple shearNow, as a next step, we want to evaluate how the models perform if they have to extrapolate for completely unknown load cases, i.e., biaxial stress states

\[(C_{KL})^{\text{biaxial}}:=\begin{pmatrix}\lambda_{1}^{2}&0&0\\ 0&\lambda_{1}^{2}&0\\ 0&0&\lambda_{2}^{2}\end{pmatrix},\quad(T_{KL})^{\text{biaxial},\oplus}= \begin{pmatrix}T_{11}^{\oplus}&0&0\\ 0&T_{11}^{\oplus}&0\\ 0&0&0\end{pmatrix} \tag{51}\]

\begin{table}
\begin{tabular}{c c c c} Model & Ideal data & Data with offset & Noisy data \\ \hline _F-P_ model & \(3.14\,\text{kPa}^{2}\) & \(2.87\,\text{kPa}^{2}\) & \(2.02\cdot 10^{3}\,\text{kPa}^{2}\) \\ Basic conditions NN & \(1.74\cdot 10^{-4}\,\text{kPa}^{2}\) & \(1.22\cdot 10^{3}\,\text{kPa}^{2}\) & \(2.03\cdot 10^{3}\,\text{kPa}^{2}\) \\ PANN model & \(5.92\cdot 10^{-5}\,\text{kPa}^{2}\) & \(2.77\cdot 10^{3}\,\text{kPa}^{2}\) & \(2.02\cdot 10^{3}\,\text{kPa}^{2}\) \\ \end{tabular}
\end{table}
Table 1: MSE according to Eq. (46) achieved with the _F-P_ model, the NN fulfilling the basic conditions and the PANN for training with different uniaxial stress data. The respective stress predictions are shown in Fig. 3.
and simple shear

\[(C_{KL})^{\text{shear}}:=\begin{pmatrix}1&\gamma&0\\ \gamma&\gamma^{2}+1&0\\ 0&0&1\end{pmatrix} \tag{52}\]

with \(0\leq\gamma\leq 2\) denoting the shearing. The results for the stress predictions are given in Fig. 4(b1) and (c1). Again, the simple approach \(\mathbf{P}^{\text{simple}}(\mathbf{F})\) completely fails for the unseen data. In the simple shear load case, even a completely implausible shear stress \(T_{12}^{\text{simple}}\approx 0\) is predicted, since the model has learned in training that the diagonal elements of the stress tensor disappear. In contrast, the basic conditions model \(\psi^{\text{NN},\emptyset}(\mathbf{T}^{\ast})\) is still quite good at physically plausible extrapolation, but shows noticeable deviations for \(\lambda_{1}>1.4\) within the biaxial stress loading. A similar result is obtained for simple shear, where noticeable deviations occur starting at \(\gamma\approx 0.8\). The same holds for the approximation of the energy which is shown in Fig. 4(b2) and (c2). Surprisingly, the PANN model is able to predict the data almost perfectly for biaxial and simple shear, although full extrapolation is required here. The significantly improved extrapolation capability of the PANN model \(\psi^{\text{PANN},\emptyset}(\mathbf{I}^{\ast})\) compared to \(\mathbf{P}^{\text{simple}}(\mathbf{F})\) and \(\psi^{\text{NN},\emptyset}(\mathbf{I}^{\ast})\) becomes also evident from the MSEs given

Figure 4: Predicted stress-stretch curves \(T_{KL}^{\otimes}(\lambda_{1})\) and \(T_{KL}^{\otimes}(\gamma)\) of the PANN, the NN fulfilling the basic conditions, as well as the \(\mathbf{F}\)-\(\mathbf{P}\) model with one hidden layer containing \(N^{\text{NN}}:=4\) neurons. Shown is the extrapolation of the models trained with uniaxial stress data stored in \(\mathcal{D}^{\text{uniaxial},\emptyset}\) from an isotropic Neo-Hooke model for: (a1) uniaxial tension/compression test, (b1) biaxial tension/compression test, and (c1) a simple shear test. The corresponding energies are given in (a2)–(c2).
in Tab. 2. To illustrate into which ranges the models must extrapolate for the considered test cases, the deformation states examined are shown within the invariant space in Fig. 5.

Thus, summarizing the findings of the presented study, the adaption of NN-based models in such a way that they fulfill physical conditions for arbitrary loadings does not necessarily improve the approximation of training data, at least for the simple load cases considered, but it does allow for a significant improvement in the extrapolation capability.

### Complex multiaxial stress-strain states

Now, as we have shown the advantage of _invariant-energy-based NN approaches_, we want to analyze only this model class in the following. Thereby, it is of interest whether the exact fulfillment of the physical principles introduced in Sec. 2.2 could be too restrictive, so that the model is not flexible enough to approximate complex multiaxial stress-strain states given by a non-linear material behavior sufficiently well. Thereby, the four different architectures \(\psi^{\text{(i),\tiny{D}}}(\mathbf{I}^{*})\), \(\psi^{\text{(ii),\tiny{D}}}(\mathbf{I}^{*})\), \(\psi^{\text{(iii),\tiny{D}}}(\mathbf{I}^{*})\), as well as \(\psi^{\text{(iv),\tiny{D}}}(\mathbf{I}^{*})\) satisfying

1. the basic conditions, i.e., thermodynamic consistency, symmetric \(\mathbf{\sigma}\), objectivity and material symmetry,
2. the conditions of (i) + polyconvexity,
3. the conditions of (ii) + growth condition, as well as
4. the conditions of (iii) + energy and stress normalization

are compared to each other for isotropic as well as transversely isotropic behavior, respectively. Note that the first and the fourth model are equal to \(\psi^{\text{NN,\tiny{D}}}(\mathbf{I}^{*})\) and \(\psi^{\text{PANN,\tiny{D}}}(\mathbf{I}^{*})\).

\begin{table}
\begin{tabular}{c c c c c} Model & Training uniaxial stress & Extrapolation uniaxial stress & Biaxial stress & Simple shear stress \\ \hline \(\mathbf{F}\)_-\(\mathbf{P}\)_ model & \(1.95\cdot 10^{-1}\) kPa\({}^{2}\) & \(9.43\cdot 10^{4}\) kPa\({}^{2}\) & \(1.35\cdot 10^{5}\) kPa\({}^{2}\) & \(8.06\cdot 10^{5}\) kPa\({}^{2}\) \\ Basic conditions NN & \(1.30\cdot 10^{-4}\) kPa\({}^{2}\) & \(1.26\cdot 10^{5}\) kPa\({}^{2}\) & \(6.24\cdot 10^{5}\) kPa\({}^{2}\) & \(5.65\cdot 10^{4}\) kPa\({}^{2}\) \\ PANN model & \(3.91\cdot 10^{-5}\) kPa\({}^{2}\) & \(6.21\cdot 10^{2}\) kPa\({}^{2}\) & \(4.11\cdot 10^{3}\) kPa\({}^{2}\) & \(1.58\cdot 10^{-5}\) kPa\({}^{2}\) \\ \end{tabular}
\end{table}
Table 2: MSE according to Eq. (46) achieved with the \(\mathbf{F}\)_-\(\mathbf{P}\)_ model, the NN fulfilling the basic conditions, and the PANN model for training with uniaxial tension/compression (\(0.8\leq\lambda_{1}\leq 1.1\)) as well as extrapolation to unknown states.

Figure 5: Location of the deformation states within the isotropic invariant space for uniaxial tension/compression test, biaxial tension/compression test, and simple shear test. Shown are the sectional planes \(I_{1}\)-\(I_{2}\), \(I_{1}\)-\(I_{3}\), and \(I_{2}\)-\(I_{3}\).
#### 4.2.1 Generation of training data

In a first step, the data basis for the training of the NNs has to be acquired. In the absence of real experiments, it is generated numerically here. To this end, a uniaxial tensile test is applied on a virtual sample in an FE simulation, see Fig. 6(a) and (b). Here, the sample's geometric dimensions are specified by \(L_{x_{1}}\times L_{x_{2}}\times L_{x_{3}}=(100\times 100\times 5)\,\mathrm{mm}\) and within the uniaxial tensile test the displacement boundary condition with a maximum value \(\hat{u}\leq 40\,\mathrm{mm}\) is linearly increased in each increment of the FE simulation. The geometry generation and the meshing was realized using the tool _Gmsh_[15].

Within this virtual experiment, the required dataset \(\mathcal{D}^{\oplus}\) consisting of the tuples \(\mathcal{D}^{\oplus}_{i}:=(^{i}\mathbf{C},^{i}\mathbf{T}^{\oplus})\) is collected at the quadrature points of the finite elements within 30 increments. In order to demonstrate the ability of the proposed NN-based method, the nonlinear stress-strain relation (20) is chosen for the isotropic constitutive behavior of the sample's material. Here the material parameters for the analytical model are chosen according to Tab. 3.

Following the work of Kalina et al. [20] and prescribing a relative tolerance \(\eta:=1\,\%\), the dataset was then filtered with respect to the invariants \((^{i}I_{1},^{i}I_{2},^{i}I_{3})\in\mathbb{R}^{3}\) for the corresponding deformation state \({}^{i}\mathbf{C}\) to obtain a reduced dataset \(\mathcal{D}^{\mathrm{red},\oplus}\) with \(|\mathcal{D}^{\mathrm{red},\oplus}|=963\) for the calibration process of the PANN. This procedure took advantage of the fact that the introduced model lives in the space of invariants rather than in the space of deformations.

Finally, different deformation states are available, which are stored for the isotropic case with the corresponding stresses of the Neo-Hooke model in the dataset \(\mathcal{D}^{\mathrm{red},\oplus}\). For the transversely isotropic case with the preferred direction parallel to the \(X_{1}\)-direction, the same deformation states are chosen as data basis and stored with the corresponding stresses of Schroder's model (22) in the dataset \(\mathcal{D}^{\mathrm{red},\parallel}\) with \(|\mathcal{D}^{\mathrm{red},\parallel}|=963\) as well. The chosen parameters of the analytical transversely isotropic model are also given in Tab. 3.

#### 4.2.2 Overall prediction quality

First, for the reduced dataset \(\mathcal{D}^{\mathrm{red},\Box}\), we will compare the overall prediction quality of the trained NNs (i)-(iv) satisfying different physical constraints, introduced in Sec. 2.2, for the isotropic as well as the transversely isotropic case. Thereby, for the transversely isotropic case \(\beta=2\) is chosen for the structural tensor \(\mathbf{G}\) given

\begin{table}
\begin{tabular}{c c||c c c c c c} \(E\) & \(\nu\) & \(\beta\) & \(\alpha_{1}\) & \(\alpha_{2}\) & \(\delta_{1}\) & \(\delta_{2}\) & \(\alpha_{4}\) & \(\eta_{1}\) \\ \hline \(10^{3}\,\mathrm{kPa}\) & 0.3 & 2 & 8 kPa & 0 kPa & 10 kPa & 56 kPa & 2 & 10 kPa \\ \end{tabular}
\end{table}
Table 3: Material parameters of isotropic Neo-Hooke and transversely isotropic analytical model proposed by Schroder et al. [45]. The models are given in Eqs. (20) and (22), respectively.

Figure 6: Uniaxial tensile test for data generation: (a) applied boundary conditions with prescribed displacement \(\hat{u}\,\mathbf{e}_{1}\) and (b) inhomogeneous specimen geometry with \(L_{x_{1}}\times L_{x_{2}}\times L_{x_{3}}=(100\times 100\times 5)\,\mathrm{mm}\).
in Eq. (18). In the following examples, the network architecture is set to one hidden layer with \(N^{\text{NN}}:=8\) neurons. The NN-based models are trained with respect to the dataset \(\mathcal{D}^{\text{red,\Box}}\), where a random division into calibration (70 %) and test (30 %) data is made once. Within the training process, the weights and bias values \(W_{\alpha},b_{\alpha},w_{\alpha B}\) and \(w^{*}_{\alpha\gamma}\) are then determined according to Sec. 3.4. Within one training run, the respective NN is trained 30 times, and the parameters of the best achieved training state with the lowest MSE (46), see Sec. 3.4, are stored at the end [20].7 In order to evaluate the approximation behavior of the trained NNs, we compute the relative error measure

Footnote 7: Due to local minima within the loss function, the optimization procedure which is applied here depends on the starting values of the weights and biases. Thus, the network is trained several times to overcome this, cf. Kalina et al. [19, 20].

\[\varepsilon^{\Box}:=\frac{\max\limits_{i\in\mathbb{N}_{\leq|\Box^{\Box}}} \left\|i\mathbf{T}^{\Box}-i\mathbf{T}^{\text{model,\Box}}\right\|}{\max\limits_{j\in \mathbb{N}_{\leq|\Box^{\Box}}}\left\|i\mathbf{T}^{\Box}\right\|} \tag{53}\]

for the Frobenius norm \(\left\|\mathbf{T}\right\|\) of the second Piola-Kirchhoff stress. To exclude random effects, a statistical study is performed, i.e., a total of 300 training runs has been performed for each model. Since no uniform distribution with respect to the error measure \(\varepsilon^{\Box}\) can be seen from the results, no underlying distribution is assumed here, see the histogram plots given in App. C. Thus, median \(\varepsilon^{\text{med,\Box}}\) and quantiles are used to compare the interpolation quality of the models with each other. The results of the described study are shown in Fig. 7.

Isotropic modelRegarding the NNs' stress approximation for the isotropic case, an extremely good prediction quality with a median of the errors \(\varepsilon^{\text{med,\phi}}<0.003\,\%\) is achieved for all models, cf. Fig. 7(a). As one can see, the approximation quality of \(\psi^{\text{(ii),\phi}}(\mathbf{I}^{*})\) accounting for the polyconvexity worsens in comparison to the architecture \(\psi^{\text{(i),\phi}}(\mathbf{I}^{*})\) which only fulfills the basic conditions. Thus, the limitation to positive weights reduces the approximation quality of the NN in the statistical sense. However, as already mentioned, the errors are still extremely low. If the growth condition is further added, the error of this NN denoted as \(\psi^{\text{(iii),\phi}}(\mathbf{I}^{*})\), which fulfills basic conditions + polyconvexity + growth condition, stays in a similar range. Surprisingly, when the normalization conditions are finally added, the resulting distribution achieved with \(\psi^{\text{(iv),\phi}}(\mathbf{I}^{*})\) is very similar to the first model, i.e., the basic conditions model. Thus, summarizing the study

Figure 7: Boxplots with the median \(\varepsilon^{\text{med,\Box}}\), 25th and 75th as well as the 1st and 99th percentile of relative error measure \(\varepsilon^{\Box}\) according to Eq. (53) for the reduced dataset \(\mathcal{D}^{\text{red,\Box}}\) and the (a) isotropic NNs or (b) transversely isotropic NNs with (i) basic conditions, (ii) polyconvexity, (iii) growth condition and polyconvexity as well as (iv) the PANN satisfying all conditions including normalization.
carried out, adding all common physical principles of hyperelasticity to the NN-based model does not lead to a deterioration in the approximation for the isotropic case. This is although conditions such as positive weights to account for polyconvexity are reducing the variability of the NN.

Finally, regarding the non-negativity condition on the elastic energy, a numerical test has been applied. To scan only physically meaningful deformation states, the principle stretches \(\lambda_{1},\lambda_{2},\lambda_{3}\in\mathbb{R}_{>0}\) are varied and \(\mathbf{C}^{\text{diag}}(\lambda_{1},\lambda_{2},\lambda_{3}):=(\lambda_{1}^{2} \mathbf{e}_{1}\otimes\mathbf{e}_{1}+\lambda_{2}^{2}\mathbf{e}_{2}\otimes\mathbf{e}_{2}+\lambda_ {3}^{2}\mathbf{e}_{3}\otimes\mathbf{e}_{3})\) is computed as a diagonal tensor which is sufficient for isotropy. The invariants are calculated and the respective energy is determined. Within a range of \(1/10\leq\lambda_{\alpha}\leq 10\) with \(\alpha\in\mathbb{N}_{\leq 3}\), only positive energies have been numerically detected for both the Neo-Hooke model (19) and the trained isotropic PANN with relative error \(\varepsilon^{\text{\tiny\textregistered}}\), which is closest to the median \(\varepsilon^{\text{med},\text{\tiny\textregistered}}\), exemplarily.

Transversely isotropic modelRegarding the NNs' stress approximation for the transversely isotropic scenario given in Fig. 7(b), compared to isotropy, the error is now an order of magnitude higher even for the best model \(\psi^{(\text{i}),\parallel}(\mathbf{I}^{*})\). This can be explained by the increased complexity of the reference model (21). However, with a median of the errors \(\varepsilon^{\text{med},\parallel}<0.4\,\%\), a high level of prediction quality is still attained for all models. Now the four models (i)-(iv) are compared. Applying the model \(\psi^{(\text{ii}),\parallel}(\mathbf{I}^{*})\), which takes polyconvexity into account, leads to errors that are increased by an order of magnitude compared to the basic conditions architecture \(\psi^{(\text{i}),\parallel}(\mathbf{I}^{*})\). Thus, the limitation to positive weights significantly reduces the approximation quality of the NN. If the growth condition is further added, the error of this NN, denoted as \(\psi^{(\text{iii}),\parallel}(\mathbf{I}^{*})\), stays in a similar range. Finally, when the normalization conditions are added, the resulting distribution of errors achieved with \(\psi^{(\text{iv}),\parallel}(\mathbf{I}^{*})\) shifts again to slightly larger values.

We would like to emphasize that the models are based on only one hidden layer with \(N^{\text{NN}}=8\) neurons. The approximation quality of the NN-based models could easily be increased by adapting the network architecture, i.e., using several hidden layers or more neurons, in order to be able to represent the complex material behavior sufficiently well. Summarizing, adding all common physical principles of hyperelasticity to the NN-based model leads to a deterioration of the prediction quality of approximately one order of magnitude for the transversely isotropic case. Since the errors are nevertheless very small, the PANN model should be chosen anyway, especially with regard to the very good extrapolation capability.

Regarding the non-negativity condition on the elastic energy, a numerical test has been applied for transverse isotropy, too. Again, to scan only physically meaningful deformation states, the principle stretches \(\lambda_{1},\lambda_{2},\lambda_{3}\in\mathbb{R}_{>0}\) are varried and the diagonal tensor \(\mathbf{C}^{\text{diag}}(\lambda_{1},\lambda_{2},\lambda_{3})\) is computed. Here, in addition, rotations of these states perpendicular to the preferred direction \(X_{1}\) have to be considered. Thus, we end up with a five parameter space for the deformation to be sampled:

\[\mathbf{C}^{\parallel}(\lambda_{1},\lambda_{2},\lambda_{3},\varphi_{2},\varphi_{3 })=\mathbf{R}(\varphi_{2},\varphi_{3})\cdot\mathbf{C}^{\text{diag}}(\lambda_{1}, \lambda_{2},\lambda_{3})\cdot\mathbf{R}^{T}(\varphi_{2},\varphi_{3}). \tag{54}\]

In the equation above, \(\mathbf{R}(\varphi_{2},\varphi_{3})=\mathbf{R}_{x_{2}}(\varphi_{2})\cdot\mathbf{R}_{x_{3}} (\varphi_{3})\in\mathcal{G}(3)\) is a rotation tensor with \(\mathbf{R}_{x_{2}}(\varphi_{2})\) and \(\mathbf{R}_{x_{2}}(\varphi_{3})\) denoting rotations around the \(X_{2}\)- and \(X_{3}\)-axis, respectively. With Eq. (54) and the structural tensor given in Eq. (18) invariants are calculated and the respective energy is determined. Within a range of \(1/10\leq\lambda_{\alpha}\leq 10\), \(0\leq\varphi_{\beta}\leq\pi/2\) with \(\alpha\in\mathbb{N}_{\leq 3},\beta\in\{2,3\}\), only positive energies have been numerically detected for both the model (21) and the trained transversely isotropic PANN with relative error \(\varepsilon^{\parallel}\), which is closest to the median \(\varepsilon^{\text{med},\parallel}\), exemplarily.

#### 4.2.3 Application of the calibrated PANN within an FE simulation

In order to prove the suitability of invariant-based NNs for the numerical simulation of complex shaped samples and components, a comparison to reference results generated with the isotropic Neo-Hooke model (20), which has been used for the calibration, is shown in the following. Thereby, the model fulfilling the basic
conditions and the PANN are each analyzed in two different scenarios: once trained with multiaxial stress-strain states \(\mathcal{D}^{\text{red,\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny \
material tangents

\[\mathbb{C}^{\text{NN},\oplus}:=4\frac{\partial^{2}\psi^{\text{NN},\oplus}}{ \partial\mathbf{C}\partial\mathbf{C}}\in\mathcal{L}_{4}\quad\text{and}\quad\mathbb{C}^ {\text{PANN},\oplus}:=4\frac{\partial^{2}\psi^{\text{PANN},\oplus}}{\partial\bm {C}\partial\mathbf{C}}\in\mathcal{L}_{4}\;, \tag{55}\]

which are required within the solution via a standard Newton-Raphson scheme, are calculated by means of _automatic differentiation_.

To evaluate the NNs' prediction quality, a comparison to the local stress field \(P_{31}^{\text{nh}}\) is considered in the following, cf. Fig. 9(b). We start with the models trained by complex multiaxial stress-strain states stored in the dataset \(\mathcal{D}^{\text{red},\oplus}\), where the covered domain in the invariant space is shown in Fig. 8. For both NNs, relative errors below \(0.004\,\%\) occur, see the surface plots given in Fig. 9(c), (d), i.e., the predictions are almost perfect with respect to the reference stresses \(P_{31}^{\text{nh}}\). This is not surprising, since a mapping of the deformations of the torsional sample into the invariant space shows that it is almost completely covered by \(\mathcal{D}^{\text{red},\oplus}\). Therefore, nearly no extrapolation is necessary here, see also the discussion in [20]. In the next step, we consider the NN-based models only trained by uniaxial stress states stored in the dataset \(\mathcal{D}^{\text{uniaxial},\oplus}\), where the covered curve in the invariant space is again shown in Fig. 8. Looking now at the results in Figs. 10(c) and (d), we see a significant error of \(3.7\,\%\) for the basic conditions model, which is due to the need for extrapolation. However, despite this, very low errors below \(0.14\,\%\) are achieved with the PANN model.

Consequently, in this example, both the basic conditions model and the PANN approach are very well able to describe the learned nonlinear constitutive behavior within the FE simulation of a comparatively complex load case, whereas the PANN is better, especially if extrapolation is required. Thereby, NNs with \(N^{\text{NN}}=8\) and \(N^{\text{NN}}=4\) neurons in only one hidden layer are used, respectively, which is very small compared to typical NNs applied to problems originating from computational mechanics. Moreover, the implemented basic conditions model and the PANN provide the typical quadratic convergence of Newton iteration and are therefore computationally very efficient.

Figure 10: FE-simulation of torsional sample: (a) loading conditions, (b) macroscopic stress field \(P_{31}^{\text{nh}}\) on the deformed configuration \(\mathcal{B}\) by specifying a distortion of \(\hat{\phi}=45^{\circ}\), and (c) relative error of the PANN-stress field \(P_{31}^{\text{PANN},\oplus}\) as well as (d) relative error of the stress field \(P_{31}^{\text{NN},\oplus}\) of basic conditions model with respect to \(P_{31}^{\text{nh}}\). The NNs trained with only uniaxial stress-strain data \(\mathcal{D}^{\text{uniaxial},\oplus}\), \(|\mathcal{D}^{\text{uniaxial},\oplus}|=15\), and \(N^{\text{NN}}=4\) neurons in only one hidden-layer were implemented as constitutive equations each.
Conclusion

In the present work, an NN-based constitutive model for compressible finite strain hyperelasticity is proposed. This approach denoted as PANN fulfills all common constitutive conditions belonging to the class of hyperelasticity, i.e., _thermodynamic consistency, symmetry of the stress tensor, objectivity, material symmetry, growth condition_, as well as _normalization of energy and stress_, in an exact way. Furthermore, the condition on the _non-negativity of the strain energy_ is numerically validated for the trained models. The proposed model allows the description of highly nonlinear hyperelastic relationships while taking into account the underlying physics and is trainable by using standard machine learning libraries such as Tensorflow.

Starting with a short literature review on NN-based elastic models, an introduction on finite strain hyperelasticity including an overview on general requirements as well as two specific anisotropy classes and models is given. Based on this, the PANN approach is built up step by step: using different sets of invariants as inputs for a convex neural network, the model fulfills the balance of angular momentum, objectivity and material symmetry conditions, as well as thermodynamic consistency and polyconvexity. Then, the volumetric growth condition is fulfilled by using an analytical growth term. Finally, energy and stress normalization are fulfilled by polyconvex normalization terms. The stress normalization terms depend on the material symmetry group and are exemplarily derived for isotropic and transversely isotropic material behavior. However, the procedure for fulfilling physical conditions, e.g., using normalization terms, can also be applied to analytical or other machine learning approaches. Afterwards, the applicability of the PANN models is demonstrated, where a calibration to isotropic and transversely isotropic data generated with analytical potentials is performed. For all cases, even for highly multiaxial deformation states and noisy stress-strain data, a highly accurate and robust prediction quality has been shown. In addition, it has been shown that the PANN is characterized by an extremely good extrapolation capability. Finally, the straightforward application into FE simulations is demonstrated.

Summarizing, the introduced PANN model for compressible finite strain hyperelasticity has shown to be an efficient tool, which can be used in numerous applications stemming from solid mechanics. Thereby, including physics into the NN-based model is the crucial step for several reasons: first of all it leads to reliable, i.e., physically sensible, model predictions. But more than that, it is also essential in order to improve the generalization properties of NNs [22, 27] and allows for extrapolation [19, 25]. Furthermore, only with the pronounced mathematical structure that the inclusion of constitutive conditions provides, it is possible to calibrate the models with small amounts of data which are usually available in engineering applications. Finally, this also enables the use of comparatively small network architectures, cf. [19, 25]. Thus, when constructing NN-based constitutive models, as many constitutive conditions should be fulfilled in an exact way as possible, as this is the only way to ensure their fulfillment with absolute certainty. However, for some applications this might not be possible, e.g., when for the symmetry group under consideration no complete functional basis in invariants is available and using invariants would restrict the model flexibility too much, cf. [24]. Only then the structure of the model should be weakened by fulfilling some constitutive conditions in an approximate fashion, in order to gain more model flexibility. Besides that, as already mentioned, the structure and reliability that the exact fulfillment of constitutive conditions provides should always be prioritized.

In order to generalize our proposed NN-based framework for elastic materials, several extensions are planned in the future. For instance, polyconvex normalization terms for further material symmetry groups [8] have to be derived. Furthermore, in order to allow an automated discovery of type and orientation of the underlying anisotropy, the usage of tensor-basis NNs would be a valuable addition [12]. The application of the PANN model in the identification of material models from experimental is also promising [10, 50]. Finally, an extension to multiphysics problems [21, 25] is feasible to expand the possible field of application.
## Acknowledgements

Dominik K. Klein and Oliver Weeger acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG - German Research Foundation) - Grant No. 492770117 and support by the Graduate School of Computational Engineering within the Centre of Computational Engineering at the TU Darmstadt.

All presented computations were performed on a PC-Cluster at the Center for Information Services and High Performance Computing (ZIH) at TU Dresden. The authors thus thank the ZIH for generous allocations of computer time. Finally, the authors would like to thank Franz Hirsch and Philipp Metsch for providing the serverjob scripts to communicate with the HPC-Cluster.

## CRediT authorship contribution statement

**Lennart Linden:** Conceptualization, Formal analysis, Investigation, Methodology, Visualization, Software, Validation, Visualization, Writing - original draft, Writing - review and editing. **Dominik Klein:** Conceptualization, Formal analysis, Methodology, Visualization, Software, Validation, Writing - original draft, Writing - review and editing. **Karl A. Kalina:** Conceptualization, Formal analysis, Methodology, Visualization, Software, Writing - original draft, Writing - review and editing. **Jorg Brummund:** Formal analysis, Methodology, Writing - review and editing. **Oliver Weeger:** Conceptualization, Funding acquisition, Resources, Writing - review and editing.

## Declarations

**Conflict of interest:** The authors declare that they have no conflict of interest.

## Appendix A Multilayered neural networks

In this work, sets of invariants are used as inputs for FFNNs [1, 26] with scalar-valued output, where the output is used to model a hyperelastic potential. In Sec. 3, a FFNN architecture with only one hidden layer is introduced, which prooves to be flexible enough for a lot of practical applications [19, 25]. Nevertheless, the methods introduced in Sec. 3 are not restricted to network architectures containing only one hidden layer, and for the sake of completeness, multilayered network architectures, cf. Fig. 11, are now introduced.

In a nutshell, FFNNs can be seen as a composition of multiple vector-valued functions, where the components are referred to as nodes or neurons and the functions acting in each node are referred to as activation functions. FFNNs can gain flexibility in two ways: either the number of nodes in each hidden layer is increased, as it is done for one hidden layer in Sec. 3, or the amount of hidden layers is increased. Generalizing the single-layered architecture of Eq. (23) for a network architecture with \(H\) hidden layer and \(N^{\text{NN},H}\) nodes in each hidden layer yields

\[\mathbf{A}_{\alpha}^{[1]} =\mathcal{F}\bigg{(}\sum_{\beta=1}^{m}w_{\alpha\beta}^{[1]}I_{ \beta}+\sum_{\gamma=1}^{A}w_{\alpha\gamma}^{*[1]}I_{\gamma}^{*}+b_{\alpha}^{[1 ]}\bigg{)}\in\mathbb{R}^{N^{\text{NN},1}}\,, \tag{56}\] \[\mathbf{A}_{\alpha}^{[h]} =\mathcal{F}\bigg{(}\sum_{\beta=1}^{N^{\text{NN},h-1}}w_{\alpha \beta}^{[h]}\mathbf{A}_{\beta}^{[h-1]}+b_{\alpha}^{[h]}\bigg{)}\in\mathbb{R}^{N^{ \text{NN},h}}\text{ with }h=2,\ldots,H\,,\] (57) \[\psi^{\text{NN},\Box}(\mathbf{I}^{*}) =\sum_{\alpha=1}^{N^{\text{NN},H}}W_{\alpha}\,\mathbf{A}_{\alpha}^{[ H]}\in\mathbb{R}\,, \tag{58}\]
with the polyconvex, irreducible and independent invariants \(I_{\beta}(\mathbf{C},\mathcal{S}^{\mathbb{D}})\) as well as the additional invariants \(I_{\gamma}^{*}(\mathbf{I})\) as defined in Sec. 3. By the special choice of activation function \(\mathcal{F}\) as convex and non-decreasing in every hidden layer and \(W_{\alpha}\geq 0\), the polyconvexity of the invariants is preserved. Again, the _Softplus_ activation function \(\mathcal{SP}(x):=\log(1+\exp(x))\) is applied, which is convex and non-decreasing for non-negative weights and arbitrary bias values, and overall the conditions

\[W_{\alpha},w_{\alpha\beta}^{[h]},w_{\alpha\gamma}^{*[1]}\in\mathbb{R}_{\geq 0 },b_{\alpha}\in\mathbb{R}\quad\forall h\in\mathbb{N}_{\leq H},\alpha\in \mathbb{N}_{\leq N^{\mathrm{NN},h-1}},\beta\in\mathbb{N}_{\leq m},\gamma\in \mathbb{N}_{\leq A} \tag{59}\]

result in a polyconvex neural network, cf. Remark 3.1, see also [24] for a more extensive discussion and explicit proofs.

## Appendix B Derivatives of invariants

For the convenience of the reader, in Tab. 4, both isotropic and transversely isotropic invariants, as well as their derivatives w.r.t. \(\mathbf{C}\), are provided. From this, the derivative of the determinant \(J=\sqrt{I_{3}}\) follows as

\[\frac{\partial J}{\partial\mathbf{C}}=\frac{J}{2}\mathbf{C}^{-1}\, \tag{60}\]

which implies the derivative of the adapted invariant \(I_{1}^{*}:=-2J\) as

\[\frac{\partial I_{1}^{*}}{\partial\mathbf{C}}=-\mathbf{JC}^{-1}\,. \tag{61}\]

\begin{table}
\begin{tabular}{|l|l|c|} \hline \hline  & \(I_{\alpha}\) & \(\frac{\partial I_{\alpha}}{\partial\mathbf{C}}\) \\ \hline  & \(I_{1}:=\operatorname{tr}\mathbf{C}\) & \(\mathbf{1}\) \\ Isotropy & \(I_{2}:=\operatorname{tr}(\operatorname{cof}\mathbf{C})\) & \(I_{1}\mathbf{1}-\mathbf{C}\) \\  & \(I_{3}:=\operatorname{det}\mathbf{C}\) & \(\operatorname{cof}\mathbf{C}\) \\ \hline  & \(I_{4}:=\operatorname{tr}(\mathbf{C}\cdot\mathbf{G})\) & \(\mathbf{G}\) \\ isotropy & \(I_{5}:=\operatorname{tr}(\operatorname{cof}(\mathbf{C})\cdot\mathbf{G})\) & \(I_{5}\mathbf{C}^{-1}-\operatorname{cof}(\mathbf{C})\cdot\mathbf{G}\cdot\mathbf{C}^{-1}\) \\ \hline \hline \end{tabular}
\end{table}
Table 4: Derivatives of isotropic and transversely isotropic invariants.

Figure 11: Illustration of the multilayered PANN based constitutive model for the material symmetry group \(\square\) under consideration.
Stochastic

In this appended section, the raw data obtained in the statistical study which is discussed in Sec. 4.2.2 are presented. The histograms including median and 25th as well as 75th percentile for the isotropic (\(\Box=\otimes\)) and the transversely isotropic (\(\Box=\|\)) case are given in Fig. 12(a)-(b), respectively. For each model in this statistical investigation, a total of 300 training runs have been completed. Based on the reduced data set \(\mathcal{D}^{\text{red,\Box}}\), the respective NN is trained 30 times in one training run, and the parameters of the optimal training state with the lowest MSE (46), as described in Sec. 3.4, are stored.
Figure 12: Histograms with \(N^{\text{bins}}:=30\) of the error measure \(\varepsilon^{\Box}\) given in Eq. (53) for different anisotropies: (a) isotropic and (b) transversely isotropic invariant-based models. (i) basic conditions, (ii) polyconvexity, (iii) growth condition and polyconvexity, as well as (iv) PANN satisfying all conditions including normalization. The results were generated with 300 training runs, selecting the best of 30 trains each. Calibration has been done with the reduced data set \(\mathcal{D}^{\text{red,\Box}}\).