Advanced discretization techniques for hyperelastic physics-augmented
  neural networks
In the present work, advanced spatial and temporal discretization techniques
are tailored to hyperelastic physics-augmented neural networks, i.e., neural
network based constitutive models which fulfill all relevant mechanical
conditions of hyperelasticity by construction. The framework takes into account
the structure of neural network-based constitutive models, in particular, that
their derivatives are more complex compared to analytical models. The proposed
framework allows for convenient mixed Hu-Washizu like finite element
formulations applicable to nearly incompressible material behavior. The key
feature of this work is a tailored energy-momentum scheme for time
discretization, which allows for energy and momentum preserving dynamical
simulations. Both the mixed formulation and the energy-momentum discretization
are applied in finite element analysis. For this, a hyperelastic
physics-augmented neural network model is calibrated to data generated with an
analytical potential. In all finite element simulations, the proposed
discretization techniques show excellent performance. All of this demonstrates
that, from a formal point of view, neural networks are essentially mathematical
functions. As such, they can be applied in numerical methods as
straightforwardly as analytical constitutive models. Nevertheless, their
special structure suggests to tailor advanced discretization methods, to arrive
at compact mathematical formulations and convenient implementations.
Marlon Franke, Dominik K. Klein, Oliver Weeger, Peter Betsch
# Advanced discretization techniques for hyperelastic physics-augmented neural networks

###### Abstract

In the present work, advanced spatial and temporal discretization techniques are tailored to hyperelastic physics-augmented neural networks, i.e., neural network based constitutive models which fulfill all relevant mechanical conditions of hyperelasticity by construction. The framework takes into account the structure of neural network-based constitutive models, in particular, that their derivatives are more complex compared to analytical models. The proposed framework allows for convenient mixed Hu-Washizu like finite element formulations applicable to nearly incompressible material behavior. The key feature of this work is a tailored energy-momentum scheme for time discretization, which allows for energy and momentum preserving dynamical simulations. Both the mixed formulation and the energy-momentum discretization are applied in finite element analysis. For this, a hyperelastic physics-augmented neural network model is calibrated to data generated with an analytical potential. In all finite element simulations, the proposed discretization techniques show excellent performance. All of this demonstrates that, from a formal point of view, neural networks are essentially mathematical functions. As such, they can be applied in numerical methods as straightforwardly as analytical constitutive models. Nevertheless, their special structure suggests to tailor advanced discretization methods, to arrive at compact mathematical formulations and convenient implementations.

**Keywords:** finite element analysis, dynamic simulations, energy momentum scheme, mixed methods, hyperelasticity, physics-augmented neural networks

## 1 Introduction

In continuum mechanics, the motion of solid bodies is described by field equations, comprised of kinematic relations and balance laws which both are valid in general, and material specific constitutive models [103]. These field equations reflect our understanding of how material bodies behave and carefully fulfill important physical conditions. However, when applying standard numerical methods such as the widely used Newmark scheme or midpoint rule, which are well-established for linear problems, the discretization usually influences the physics in the nonlinear regime and is inclined
to numerical instabilities. This is particularly the case in dynamic simulations, where the temporal discretization does not preserve necessary physical properties of their continuous counterparts [32, 35].

To address this problem, an abundant number of structure-preserving methods, also known as geometric integrators, were proposed. These methods are designed to preserve the physics of the continuous system during the time discretization process and thereby enhance the numerical stability. There is a great variety of structure-preserving methods in the literature, which roughly can be divided into symplectic methods and energy-momentum schemes (EMS) [38, 96]. These methods can be distinguished by preserving different invariants of the system, where for instance the simultaneous preservation of the invariants, i.e. total energy, momentum maps and simplecticity, is in general not possible [68, 111]. Here we restrict our consideration to the established EMS, in particular to the discrete gradient versions thereof, which preserve both linear and angular momentum maps, and the total energy of the system. The mechanism behind the EMS can be understood as an explicit [99] or, in case of the discrete gradients, implicit [32, 33] momentum-preserving projection of the midpoint rule onto the surface of constant energy [96, 98]. EMS are promising since they fulfill most design criteria for a practical time-stepping scheme as suggested by [8]. In that regard, the EMS is applicable for elastic and inelastic problems, it does not use additional variables (e.g. Lagrange multipliers) or parameters to be provided by the analyst, it is second-order accurate, and exhibits similar stability properties in large deformations and long term simulations for coarse time step sizes as energy-conserving integrators in linear analysis. The only drawbacks associated with the EMS are that the method leads to unsymmetric tangent stiffness matrices and that the implementation is more involved by different time point evaluations of stresses and remaining parts of the equations.

The development of EMS can be traced back to the early works of [35, 63, 64, 65], which provide energy and momentum conserving methods for particle dynamics. Especially the work of [35] provided a formula which conserves both momentum maps and total energy for energy densities with only scalar dependencies. In the pioneering work [99] the so-called discrete energy-momentum method has been proposed, which has been improved and properly reformulated in [67]. A milestone in the development of EMS, which circumvents the use of the mean value algorithm and its associated solution of a nonlinear equation on quadrature point level used in [99], is related to the work [32, 33] where the discrete gradient EMS has been proposed. It is applicable to arbitrary nonlinear hyperelastic material models, leading to a projection-based algorithmic stress formula, and as a special case contains both the midpoint averaging of the strains for quadratic potentials [99] and the Greenspan formula [35] for scalar dependencies of the energy density. Other areas of successful developments of EMS are located in the fields of structural mechanics for nonlinear beams and shells [10, 11, 90, 98], inelasticity [37, 44, 76, 79], elastodynamical contact problems [18, 45, 66], and systems with (non-)holonomic constraints [17, 34], to name but a few. An extension of EMS to dissipative systems is based on the general equation for non-equilibrium reversible irreversible coupling (GENERIC) formulation, which provides a systematic framework facilitating the design of energy-momentum and even entropy consistent schemes [12, 36, 60, 75, 84, 92].

Recently, in [9] a new family of partioned discrete gradients has been proposed which is inspired by polyconvex strain energies [6, 7], leading to a new algorithmic stress formula based on three partioned discrete gradients which represent the algorithmic counterparts of the work conjugates of the right Cauchy-Green strain tensor, its cofactor, and its determinant. The formulation in [9] benefits from the so-called tensor cross-product, first published in [21] and kind of rediscovered in the field of computational nonlinear elasticity in [13], which is based on a Hu-Washizu type mixed variational formulation and leads to a remarkably simple structure for the algorithmic stress. This work has given rise to several polyconvexity-inspired works in the field of coupled problems, in particular for finite thermo-elastodynamics [28, 83], nonlinear electro-elastodynamics [26, 82], and
nonlinear thermo-electro-elastodynamics [27, 29].

In addition to the above mentioned sophisticated temporal discretization methods, significant progress has been made in the development of advanced spatial discretization techniques over the past few decades. In partiuclar, an abundant number of mixed finite element (FE) formulations in elasticity based on Hu-Washizu variational functionals [107] have been developed [4, 9, 13, 95, 97]. These formulations aim to address the issue of locking which standard displacement-based formulations are prone to.

Shifting the focus to constitutive modeling, the mechanical conditions underlying hyperelasticity were extensively discussed in the last decades [20, 103], and sophisticated analytical models were formulated to fulfill all constitutive conditions by construction [22, 94]. While the constitutive conditions that these models fulfill have a sound mechanical motivation, their explicit choice of functional relationship usually does not. Rather, it is heuristically motivated, such as the polynomial form of the Ogden model [81]. This human choice of functional relationship easily restricts the function space that a model can represent, thus introducing an undesired _model uncertainty_[50]. This is a major drawback of analytical constitutive models. On the other side, artificial neural networks (NN) are highly flexible functions, in fact, they are universal approximators [47]. However, when applied to mechanical problems, standard NNs easily violate physical conditions and require large amounts of calibration data. Hyperelastic physics-augmented neural networks (PANNs)1 aim to combine the best of both worlds. Here, NN-based models are formulated which fulfill all mechanical conditions of hyperelasticity by construction. This yields highly flexible yet reliable constitutive models. Furthermore, the constitutive conditions included in the model formulation serve as an _inductive bias_[41], which highly improves the models generalization properties and allows for a calibration with moderately sized datasets that could stem from experimental material characterization.

Footnote 1: Following [69], we denote physics-augmented neural networks (PANNs) as NN-based constitutive models which fulfill all relevant mechanical conditions by construction.

Basically, NNs are applied to represent hyperelastic potentials, which benefit from the excellent flexibility of NNs. In the same manner as analytical constitutive models [22, 94], NN-based potentials can be formulated in terms of invariants [52, 55, 56, 71, 102], thus ensuring several mechanical conditions by construction. Other approaches formulate potentials directly in terms of the components of strain tensors [5, 23, 55, 100, 105]. For the formulation of polyconvex potentials, several approaches exist [19, 55, 72, 101, 102], where the most noteworthy approaches are based on input-convex neural networks (ICNNs). Proposed in a pioneering work by [3], this special network architecture is attractive for, e.g., physical applications which require convexity [48] and convex optimization [14, 15]. While the aforementioned NN-based models fulfill some constitutive conditions in an exact way, other conditions are only approximated. E.g., the model proposed by [55] fulfills several conditions including polyconvexity in an exact way, but only approximates the stress-free reference configuration. Indeed, in the framework of NN-based constitutive models, it is possible to fulfill constitutive conditions only in an approximate fashion [78, 108, 109], e.g., the NN can learn to approximate a condition through data augmentation [70], or the loss function can include physical conditions [53]. For some applications this is even necessary, e.g., when a formulation in invariants leads to a loss of information, which is the case for some anisotropy classes [55]. However, since fulfilling a mechanical condition by construction is the only way to be absolutely sure about its fulfillment, if possible, this should be preferred over approximations.

Finally, in the works of [72] and [69], NN-based models were proposed which fulfill all common mechanical conditions of hyperelasticity by construction. While the work of [72] is restricted to incompressible material behavior, the work of [69] is formulated for general, compressible material behavior. In the present work, the model proposed by [69] is applied, which is an extension of the
polyconvex model proposed by [55]. The model fulfills mechanical conditions in several ways, first of all, by a cautious _choice of inputs_ for the NN (i.e., polyconvex invariants). Furthermore, the _special choice of network architecture_ (i.e., ICNNs) preserves the polyconvexity of the invariants. In addition, _analytical growth terms_ and _invariant-based normalization terms_ ensure a physically sensible energy and stress behavior. In that light, one could argue that the PANN model proposed by [69] is actually quite close to analytical constitutive models, with the only difference being the choice of ansatz function for the constitutive model. By using ICNNs as ansatz functions, the flexibility of the model can be increased to an arbitrary amount, while even for large number of model parameters, their calibration is still straightforward.

The most obvious reason to apply PANN models is their extraordinary flexibility in conjunction with their physically well motivated basis. In particular, PANN constitutive models can represent the highly challenging homogenized behavior of microstructured materials, for which analytical models may not be flexible enough [55]. This allows for very efficient simulation of microstructured materials [31, 52, 78]. The flexibility has another implication, namely, PANNs can represent a variety of analytical constitutive models. In fact, for a wide range of applications, analytical constitutive models are flexible enough. However, there is a vast number of analytical constitutive models, and finding a suitable one can be very demanding and requires a lot of expert knowledge in both the choice and calibration of the model [89]. This can be bypassed by NN-based constitutive models, which automatically find the best function out of a very rich function space. Another noteworthy approach for this is the data-driven framework proposed by [24, 25], where an approach for the automated discovery of analytical constitutive models is introduced. On the downside, the flexibility of PANNs goes along with a lack of traceability. Indeed, some analytical constitutive models base their choice of functional relationship on physical reasoning, such as the hyperelastic Hencky model [42, 43, 74]. For such models, the material parameters have a clear physical interpretation. This is in contrast to NN-based models, which have a large amount of parameters which are generally inaccessible to a clear interpretation [55, Sect. 6]. However, it should be noted that this is also the case for moderately flexible analytical models such as the Ogden model [81]. Finally, with data-driven approaches it is also possible to circumvent the formulation of an explicit, analytical constitutive model [16, 54, 88]. To close, in the era of artificial intelligence and big data, it is indeed possible to largely abandon the formulation of models and let ML methods purely work on data [87]. However, it is broadly accepted that including scientific knowledge in ML methods has a multitude of advantages in many scientific fields [53, 61, 73, 85, 106]. Overall, the scientific knowledge developed in the last centuries is not to be omitted, instead, the use of ML methods allows to apply them in a new light, and thereby, to exploit the advantages of up-to-date techniques.

From a formal point of view, NNs are nothing more but mathematical functions. As such, they can be applied in standard numerical methods as straightforwardly as analytical constitutive models. However, compared to analytical constitutive models, NN-based constitutive models comprise far more complicated functional relationships. Thus, while applying standard methods is indeed possible, it potentially leads to demanding mathematical formulations and implementations. This can be circumvented by tailoring numerical methods to the special structure of NN-based constitutive models.

In the present work, advanced spatial and temporal discretization techniques are tailored to the special structure of hyperelastic PANN constitutive models. In particular, inspired by [9, 13, 59], the polyconvex structure of PANN models is taken into account. Furthermore, the numerical methods are formulated to conveniently handle the - possibly very complex - functional relationships of NN-based constitutive models. This is done by setting the focus on invariant-based formulations. Based on these ideas, a new Hu-Washizu like mixed formulation is proposed which uses scalar-valued strain
invariants as additional unknowns. Furthermore, a novel discrete gradient EMS is designed, which mostly relies on scalar-valued formulations in terms of strain invariants. To the best knowledge of the authors, both mixed spatial methods and discrete gradient integration schemes have not yet been tailored to PANN constitutive models so far.

The outline of the manuscript is as follows. In Sect. 2, the basics of hyperelasticity are introduced, which are relevant for the PANN constitutive model introduced in Sect. 3. In Sect. 4, continuum formulation and variation tailored to PANN models are introduced. This is followed by the energy- and momentum consistent time discretization introduced in Sect. 5 and the spatial FE discretization discussed in Sect. 6. Stability and robustness of the proposed framework are examined in representative numerical examples in Sect. 7. After the conclusion in Sect. 8, the appendices provide some additional information on the proposed EMS.

NotationThroughout this work, tensor composition and contraction are denoted by \(\left(\mathbf{A}\,\mathbf{B}\right)_{ij}=A_{ik}B_{kj}\) and \(\mathbf{A}:\mathbf{B}=A_{ij}B_{ij}\), respectively, with vectors \(\mathbf{a}\), \(\mathbf{b}\) and second order tensors \(\mathbf{A}\), \(\mathbf{B}\). Note that the Einstein summation convention is applied if not otherwise stated. The tensor product is denoted by \(\otimes\), the second order identity tensor by \(\mathbf{I}\). Transpose and inverse of a second order tensor \(\mathbf{A}\) are denoted as \(\mathbf{A}^{\mathrm{T}}\) and \(\mathbf{A}^{-1}\), respectively. The tensor cross product as introduced by [21] is defined as \(\mathbf{A}\bigstar\mathbf{B}=\varepsilon_{ijk}\,\varepsilon_{\alpha\beta \gamma}\,A_{j\beta}\,B_{k\gamma}\,\mathbf{e}_{i}\otimes\mathbf{e}_{\alpha}\), where \(\varepsilon_{ijk}\) denotes the third-order permutation tensor. For a discussion of the properties of the tensor cross product, see [13]. Furthermore, trace, determinant and cofactor are denoted by \(\mathrm{tr}\,\mathbf{A}\), \(\det\mathbf{A}\) and \(\mathrm{cof}\,\mathbf{A}:=\det\left(\mathbf{A}\right)\,\mathbf{A}^{-\mathrm{T}}\), respectively. The space of second order tensors is denoted as \(\mathbb{T}^{3}:=\left\{\mathbf{A}\in\mathcal{L}(\mathbb{R}^{3},\mathbb{R}^{3})\right\}\), the set of invertible second order tensors with positive determinant is denoted by \(\mathbb{T}^{3}_{+}:=\left\{\mathbf{A}\in\mathbb{T}^{3}\,|\,\det(\mathbf{A})>0\right\}\), and the set of symmetric invertible second order tensors with positive determinant is denoted by \(\mathbb{S}^{3}_{+}:=\left\{\mathbf{A}\in\mathbb{T}^{3}_{+}\,|\,\mathbf{A}= \mathbf{A}^{\mathrm{T}}\right\}\). The orthogonal and special orthogonal group in \(\mathbb{R}^{3}\) are denoted by \(\mathrm{O}(3):=\left\{\mathbf{A}\in\mathbb{T}^{3}\,|\,\,\mathbf{A}^{\mathrm{T }}\mathbf{A}=\mathbf{I}\right\}\) and \(\mathrm{SO}(3):=\left\{\mathbf{A}\in\mathbb{T}^{3}\,|\,\,\mathbf{A}^{\mathrm{ T}}\mathbf{A}=\mathbf{I},\,\,\det\mathbf{A}=1\right\}\), respectively. The space of positive and negative real numbers are denoted by \(\mathbb{R}_{+}:=\left\{x\in\mathbb{R}\,\,|\,\,x>0\right\}\) and \(\mathbb{R}_{-}:=\left\{x\in\mathbb{R}\,\,|\,\,x<0\right\}\), respectively. The first Frechet derivative of a function \(f\) w.r.t. \(\mathbf{X}\) is denoted by \(\partial_{\mathbf{X}}f\), while discrete gradients are denoted by \(D_{\mathbf{X}}f\).

## 2 Basics of hyperelasticity

In this section, the kinematic relationships of nonlinear solid mechanics are introduced in Sect. 2.1, followed by the constitutive conditions of hyperelasticity in Sect. 2.2. In Sect. 2.3, the formulation of hyperelastic potentials in terms of strain invariants is discussed.

### Kinematics

Let us consider the motion of a solid body. The reference configuration \(\mathcal{B}_{0}\subset\mathbb{R}^{3}\) and the current, deformed configuration \(\mathcal{B}\subset\mathbb{R}^{3}\) are connected via the bijective deformation mapping \(\boldsymbol{\varphi}:\mathcal{B}_{0}\to\mathbb{R}^{3}\), linking material particles \(\mathbf{X}\in\mathcal{B}_{0}\) to \(\mathbf{x}=\boldsymbol{\varphi}(\mathbf{X})\in\mathcal{B}\). The material gradient of the deformation field \(\boldsymbol{\varphi}\), which is denoted as the deformation gradient \(\mathbf{F}:\mathcal{B}_{0}\to\mathbb{T}^{3}_{+}\), is defined as

\[\mathbf{F}=\partial_{\mathbf{X}}\boldsymbol{\varphi}(\mathbf{X})\,. \tag{1}\]

Then, its determinant \(J:\mathcal{B}_{0}\to\mathbb{R}_{+}\) and its cofactor \(\mathbf{H}:\mathcal{B}_{0}\to\mathbb{T}^{3}_{+}\) follow as

\[\begin{split} J&=\det(\mathbf{F})=\frac{1}{3}\, \mathbf{H}:\mathbf{F}=\frac{1}{6}\left(\mathbf{F}\boldsymbol{\Phi}\mathbf{F} \right):\mathbf{F}\,,\\ \mathbf{H}&=\mathrm{cof}(\mathbf{F})=\frac{1}{2}\, \mathbf{F}\boldsymbol{\Phi}\mathbf{F}\,.\end{split} \tag{2}\]
Let us furthermore introduce the symmetric Cauchy-Green strain tensor \(\mathbf{C}:\mathcal{B}_{0}\rightarrow\mathbb{S}_{+}^{3}\) as

\[\mathbf{C}=\mathbf{F}^{\mathrm{T}}\,\mathbf{F}\,, \tag{3}\]

with the corresponding determinant \(C:\mathcal{B}_{0}\rightarrow\mathbb{R}_{+}\) and cofactor \(\mathbf{G}:\mathcal{B}_{0}\rightarrow\mathbb{S}_{+}^{3}\) defined as

\[\begin{split} C&=\det(\mathbf{C})=\frac{1}{3}\, \mathbf{G}:\mathbf{C}=\frac{1}{6}\left(\mathbf{C}\boldsymbol{\Psi}\mathbf{C} \right):\mathbf{C}\,,\\ \mathbf{G}&=\operatorname{cof}(\mathbf{C})=\frac{1}{ 2}\,\mathbf{C}\boldsymbol{\times}\mathbf{C}\,.\end{split} \tag{4}\]

Based on the above symmetric kinematic set, we introduce the isotropic invariants \(I_{\mathbf{C}}:\mathcal{B}_{0}\rightarrow\mathbb{R}_{+}\), \(II_{\mathbf{C}}:\mathcal{B}_{0}\rightarrow\mathbb{R}_{+}\) and \(III_{\mathbf{C}}:\mathcal{B}_{0}\rightarrow\mathbb{R}_{+}\) as

\[I_{\mathbf{C}}=\operatorname{tr}(\mathbf{C})\,,\qquad II_{\mathbf{C}}= \operatorname{tr}(\mathbf{G})\,,\qquad III_{\mathbf{C}}=C\,. \tag{5}\]

### Constitutive conditions

The specific mechanical behavior of a hyperelastic material can be expressed by a potential

\[W^{\mathrm{F}}\colon\mathbb{T}_{+}^{3}\rightarrow\mathbb{R}\,,\qquad\mathbf{F }\mapsto W^{\mathrm{F}}\left(\mathbf{F}\right)\,, \tag{6}\]

which corresponds to the strain energy density stored in the body due to deformation [20, 40, 46]. With the first Piola-Kirchhoff (PK1) stress given as the gradient field

\[\mathbf{P}=\partial_{\mathbf{F}}W^{\mathrm{F}}\left(\mathbf{F}\right)\,, \tag{7}\]

the **(i) second law of thermodynamics** is fulfilled by construction. The second Piola-Kirchhoff (PK2) stress is given by

\[\mathbf{S}=\mathbf{F}^{-1}\mathbf{P}\,. \tag{8}\]

The **(ii) balance of angular momentum** implies that \(\mathbf{S}\) is symmetric. Furthermore, the principle of **(iii) objectivity** states that a model should be independent on the choice of observer, which is formalized as

\[W^{\mathrm{F}}\left(\mathbf{Q}\,\mathbf{F}\right)=W^{\mathrm{F}}\left(\mathbf{ F}\right)\quad\forall\,\mathbf{Q}\in\mathrm{SO}(3)\,. \tag{9}\]

Throughout this work, isotropic material behavior is considered, which is taken into account by the **(iv) material symmetry condition**

\[W^{\mathrm{F}}(\mathbf{F}\,\mathbf{Q}^{\mathrm{T}})=W^{\mathrm{F}}\left( \mathbf{F}\right)\quad\forall\,\mathbf{Q}\in\mathrm{O}(3)\,. \tag{10}\]

Furthermore, we consider **(v) polyconvex** potentials which allow for a representation

\[W^{\mathrm{F}}\left(\mathbf{F}\right)=\hat{W}^{\mathrm{F}}\left(\mathbf{F}, \,\mathbf{H},\,J\right)\,, \tag{11}\]

with the kinematical quantities as introduced in Eq. (2) and the function \(\hat{W}^{\mathrm{F}}\) which is convex in the components of \(\mathbf{F}\), \(\mathbf{H}\), and \(J\). The notion of polyconvexity as introduced by Ball [6, 7] stems from a quite theoretical context and is linked to the existence of solutions in finite elasticity. However, it is also the most straightforward way of fulfilling the ellipticity condition [80, 110]

\[\left(\mathbf{a}\otimes\mathbf{b}\right):\frac{\partial^{2}W^{\mathrm{F}}\left( \mathbf{F}\right)}{\partial\mathbf{F}\partial\mathbf{F}}\colon\left(\mathbf{a }\otimes\mathbf{b}\right)\geq 0\quad\forall\,\mathbf{a},\mathbf{b}\in \mathbb{R}^{3}\,. \tag{12}\]
Also known as material stability, this condition leads to a favorable behavior in numerical applications. Furthermore, a physically sensible energy and stress behavior requires fulfillment of the **(vi) growth condition**

\[W^{\rm F}\to\infty\quad\mbox{as}\quad J\to 0^{+}\,, \tag{13}\]

as well as that the reference configuration \({\bf F}={\bf I}\) is energy- and stress-free, i.e.,

\[W^{\rm F}({\bf I})=0\,,\qquad{\bf P}({\bf I})={\bf 0}\,, \tag{14}\]

also referred to as **(vii) energy normalisation** and **(viii) stress normalisation**.

### Invariant-based modeling

To conclude on the previous section, while in hyperelasticity the **(i) second law of thermodynamics** is fulfilled by construction, the remaining conditions require further considerations in the model formulation. In a first step, by formulating the potential in terms of invariants of the right Cauchy-Green tensor \({\bf C}\), cf. Eq. (5), conditions **(ii-iv)** can be fulfilled. In particular, as isotropic material behavior is assumed, the hyperelastic potential is formulated in terms of isotropic invariants, allowing for the representation

\[W^{\rm I}\colon\mathbb{R}_{+}\times\mathbb{R}_{+}\times\mathbb{R}_{+}\to \mathbb{R}\,,\qquad{\boldsymbol{\mathcal{I}}}^{0}\mapsto W^{\rm I}\left({ \boldsymbol{\mathcal{I}}}^{0}\right) \tag{15}\]

with

\[{\boldsymbol{\mathcal{I}}}^{0}:=(I_{\bf C},\,II_{\bf C},\,J). \tag{16}\]

The isotropic invariants introduced in Eq. (5) are polyconvex, in particular, \(I_{\bf C}\) is convex in \({\bf F}\) and \(II_{\bf C}\) is convex in \({\bf H}\). To preserve the polyconvexity of the invariants, the potential \(W^{\rm I}\) is chosen as a convex function in \({\boldsymbol{\mathcal{I}}}^{0}\) and furthermore as _non-decreasing_ in \(I_{\bf C},\,II_{\bf C}\), which overall ensures fulfillment of the **(v) polyconvexity condition**. The non-decreasing condition is owed to the fact that \(I_{\bf C}\), \(II_{\bf C}\) are already non-linear functions of the arguments of the polyconvexity condition, cf. Eq. (11) and Eq. (5), see also [55, Rem. A.10] for an explicit proof and [57, Sect. 2] for a simple 1D example. In contrast to that, \(J\) is the only invariant quantity in the arguments of the polyconvexity condition, meaning that the potential must be convex in \(J\), but may indeed be non-decreasing in \(J\). This additional flexibility of the model in \(J\) is also the reason to not use \(III_{\bf C}\), which would again be a non-linear function in \(J\), and is furthermore essential for the model to represent negative stress values at all, cf. [55]. The subsequent NN-based model formulation is simplified with the slightly adapted version

\[\hat{W}^{\rm I}\colon\mathbb{R}_{+}\times\mathbb{R}_{+}\times\mathbb{R}_{+} \times\mathbb{R}_{-}\to\mathbb{R}\,,\qquad{\boldsymbol{\mathcal{I}}}\mapsto \hat{W}^{\rm I}\left({\boldsymbol{\mathcal{I}}}\right)\,, \tag{17}\]

with

\[{\boldsymbol{\mathcal{I}}}:=(I_{\bf C},\,II_{\bf C},\,J,\,J^{*})\,\qquad J^{*}:=-J\,. \tag{18}\]

By this adaption, the function \(\hat{W}^{\rm I}\) can be chosen as a convex function in \({\boldsymbol{\mathcal{I}}}\) and a non-decreasing function in all of its arguments. Note that this general form of the potential does not yet fulfill conditions **(vi-viii)**, which ensure a physically sensible energy and stress behavior of the model. For a formulation of the potential in isotropic invariants, the second Piola-Kirchhoff stress follows as

\[{\bf S}=2\partial_{\bf C}\,W^{\rm I}=2\partial_{I_{\bf C}}W^{\rm I}\,{\bf I}+2 \partial_{II_{\bf C}}W^{\rm I}\,{\bf I}\bigstar{\bf C}+\partial_{J}W^{\rm I}\, J^{-1}{\bf G}\,, \tag{19}\]

where we made use of the tensor generators [51]

\[\partial_{\bf C}I_{\bf C}={\bf I}\,,\qquad\partial_{\bf C}II_{\bf C}={\bf I} \bigstar{\bf C}\,,\qquad\partial_{\bf C}J=\frac{1}{2}J^{-1}{\bf G}\,. \tag{20}\]
n analytical example for hyperelastic potentials is the Mooney-Rivlin model, given by

\[W^{\rm MR}=a\left(I_{\bf C}-3\right)+b\left(II_{\bf C}-3\right)+\frac{c}{2}\left( J-1\right)^{2}-d\,\log(J)\,,\qquad a,b,c\geq 0\,, \tag{21}\]

which is formulated such that it fulfills the remaining conditions **(vi-viii)**. In particular, the stress normalization condition is ensured by the additional constraint \(d=2(a+2b)\).

## 3 Physics-augmented neural network constitutive model

In this section, the constitutive equations of the neural network (NN) constitutive model applied throughout this work are introduced in Sect. 3.1, followed by a description of the model calibration process in Sect. 3.2. Afterwards, in Sect. 3.3, the derivatives of the NN-based constitutive model required throughout this work are introduced.

### Constitutive equations

As discussed in Sect. 2.3, hyperelastic potentials are usually formulated in terms of invariants of the right Cauchy-Green tensor. When polyconvex invariants are considered and the hyperelastic potential is a convex and non-decreasing function of the invariants, the overall potential is polyconvex. The simple structure and recursive definition of feed-forward neural networks (FFNNs) [2, 58] make them a very natural choice for the construction of convex and non-decreasing functions [3], which suggests to use them for the formulation of polyconvex hyperelastic potentials [55]. Still, the NN is only a part of the overall physics-augmented neural network (PANN) constitutive model proposed by [69], which is given by

\[W^{\rm PANN}\left(\boldsymbol{\mathcal{I}}^{0}\right)=W^{\rm NN}\left( \boldsymbol{\mathcal{I}}\right)+W^{\rm stress}\left(J\right)+W^{\rm energy }+W^{\rm growth}\left(J\right)\,, \tag{22}\]

with \(\boldsymbol{\mathcal{I}}^{0}\), \(\boldsymbol{\mathcal{I}}\) as defined in Sect. 2.3. Here, \(W^{\rm NN}\) denotes the NN part of the model, which provides the model with its excellent flexibility, while the remaining terms ensure a physically sensible stress and energy behavior of the model, cf. Sect. 2.2. The overall flow and structure of the PANN model is visualized in Fig. 1. In the present work, without loss of generality, neural networks with only one hidden layer are applied, which has shown to be sufficient for the representation of most isotropic hyperelastic potentials [56, 69]. For an introduction to multilayered PANN architectures, the reader is referred to [69]. In the present work, the NN takes the form

\[W^{\rm NN}(\boldsymbol{\mathcal{I}})=\mathbf{w}^{2}\cdot\mathcal{SP}\left( \mathbf{w}^{1}\cdot\boldsymbol{\mathcal{I}}+\mathbf{b}\right)\,, \tag{23}\]

Figure 1: Illustration of the PANN based constitutive model. Note that the hidden layer (yellow) of the NN may be multilayered.
where \(\mathcal{SP}(x)=\log(1+e^{x})\) denotes the _Softplus_ function, which is applied component-wise. This non-linear function is referred to as activation function, while \(\mathbf{w}^{1}\in\mathbb{R}_{\geq 0}^{n\times 4},\,\mathbf{w}^{2}\in\mathbb{R}_{ \geq 0}^{n}\) are referred to as weights, and \(\mathbf{b}\in\mathbb{R}^{n}\) is referred to as bias. When all weights are greater or equal to zero, and the activation function is a convex and non-decreasing function, the overall function Eq. (23) is convex and non-decreasing in \(\boldsymbol{\mathcal{I}}\), thus leading to a polyconvex potential [55]. This is the motivation for using the _Softplus_ activation function, which is convex and non-decreasing, where it should be noted that also other suitable activation functions exist [3, 5]. In general, NNs which are convex in the input are referred to as input-convex neural networks (ICNNs) [3].

Together, weights and bias form the set of model parameters, here denoted as \(\boldsymbol{\mathcal{P}}\), which are optimized in the calibration process to fit the model to a given dataset. For a given input dimension, the size of the weights and bias matrices is determined by \(n\), which specifies the number of nodes / neurons in the NN. This is a hyperparameter of the NN, which is fixed before the model calibration. By increasing \(n\), the model flexibility is increased.

The polyconvex stress normalization term as proposed by [69] is given by

\[W^{\text{stress}}(J):=-\mathfrak{n}\left(J-1\right)\,, \tag{24}\]

where the constant

\[\mathfrak{n}:=2\left(\partial_{I\mathbf{C}}\,W^{\text{NN}}+2\partial_{II \mathbf{C}}\,W^{\text{NN}}+\frac{1}{2}\partial_{J}\,W^{\text{NN}}-\frac{1}{2} \partial_{J^{*}}\,W^{\text{NN}}\right)\bigg{|}_{\mathbf{F}=\mathbf{I}}\in \mathbb{R} \tag{25}\]

is a weighted sum of derivatives of the NN potential with respect to the invariants for the undeformed state \(\mathbf{F}=\mathbf{I}\). The stress correction term preserves polyconvexity as it is linear in \(J\). The energy normalization term is given by

\[W^{\text{energy}}:=-W^{\text{NN}}\Big{|}_{\mathbf{F}=\mathbf{I}}\in\mathbb{R}\,. \tag{26}\]

Finally, the volumetric growth term is chosen as the analytical, polyconvex function

\[W^{\text{growth}}(J)=(J+J^{-1}-2)^{2}\,. \tag{27}\]

Setting aside the nomenclature of machine learning, one could consider Eq. (22) as a classical constitutive model, which uses linear transformations in combination with the _Softplus_ function as an ansatz for hyperelastic potentials, cf. Eq. (23). Then, two benefits compared to other hyperelastic constitutive models become obvious. First of all, the function \(W^{\text{NN}}\) allows for a strong interrelation between the different invariants, which becomes even more clear when considering the derivatives of the potential, cf. Sect. 3.3. This is in contrast to most analytical polyconvex constitutive models. For such models, the strain energy function is usually additively decomposed in the invariants [22], see e.g. Eq. (21). This is due to the fact that multiplication of polyconvex invariants does not preserve polyconvexity [39]. The NN-potential in Eq. (23), however, allows for a strong interrelation between the invariants, while also preserving their polyconvexity. The second benefit is the immediate possibility to increase the model flexibility, basically to an arbitrary amount [47]. While analytical constitutive models generally can become also very flexible, e.g., by considering enough terms in the polynomial ansatz of the Ogden model [81], their calibration quickly becomes infeasible. In contrast, the NN ansatz in Eq. (23) has proven to be very stable in the calibration, even for large number of parameters [55]. In this regard, it should be mentioned that the PANN model does not suffer from phenomena such as overfitting, since the inclusion of physics provides the model with a pronounced mathematical structure [69].
### Model calibration

After fixing the model architecture, here, the number of nodes \(n\), the model parameters \(\mathbf{\mathcal{P}}\) can be optimized to fit the model to a given dataset. Throughout this work, datasets of the form

\[\mathcal{D}=\left\{\left({}^{1}\mathbf{C},\,{}^{1}\mathbf{S}\right),\ldots, \left({}^{m}\mathbf{C},\,{}^{m}\mathbf{S}\right)\right\}\,, \tag{28}\]

are considered, consisting of \(m\) strain-stress tuples in the right Cauchy-Green deformation tensor \(\mathbf{C}\) and the second Piola-Kirchhoff stress \(\mathbf{S}\). Thereby, two distinct datasets are created. A calibration dataset \(\mathcal{D}_{c}\), which is used to calibrate the model, and a test dataset \(\mathcal{D}_{t}\), which is used to examine the models generalization, i.e., its capability to predict load cases not included in the calibration process. To this end, the test dataset \(\mathcal{D}_{t}\) should consist of fairly general load cases. Then, the loss function given as the mean squared error

\[\mathcal{MSE}\left(\mathbf{\mathcal{P}}\right)=\frac{1}{9\text{Pa}^{2}}\frac{1}{m_ {c}}\sum_{i=1}^{m_{c}}\left\|{}^{i}\mathbf{S}-\mathbf{S}^{\text{PANN}}\left({} ^{i}\mathbf{C};\,\mathbf{\mathcal{P}}\right)\right\|^{2}\,, \tag{29}\]

is minimized, where \(m_{c}=|\mathcal{D}_{c}|\) denotes the number of tuples in \(\mathcal{D}_{c}\) and \(\left\|\cdot\right\|\) denotes the Frobenius norm. Note that the model proposed by [69] is independent of the choice of loss function and also other loss functions could be applied for calibration. As can be seen from Eq. (29), the PANN model is calibrated only through its gradients \(\mathbf{S}^{\text{PANN}}=\partial_{\mathbf{C}}W^{\text{PANN}}\), which is referred to as Sobolev training [104, 105]. In general, also values of the potential itself could be included in the loss function. However, while synthetic datasets can yield information about energies, this is generally not the case for real-world experiments. Thus, including energy values in the loss function would lead to a less general formulation and not even lead to significant benefits in the calibrated models [23].

### Derivatives of the neural network potential

For the FE implementation of the PANN model, both first and second derivatives of the potential w.r.t. the invariants are required. For FFNN architectures with a large number of hidden layers, calculating derivatives in an explicit way quickly leads to huge expressions. However, for the single-layered FFNN used in this work, the explicit derivatives are very manageable. For the convenience of the reader, the derivatives of the NN part of the PANN model are provided here. Note that for the calibration process, cf. Sect. 3.2, the derivatives of the NN are not implemented in an explicit way, but the automatic differentiation provided by TensorFlow is applied.

The first derivative is given by

\[\partial_{(\mathbf{\mathcal{I}})_{i}}\,W^{\text{NN}}=\sum_{a=1}^{n}w_{a}^{2}\, \frac{e^{h_{a}}}{1+e^{h_{a}}}\,w_{ai}^{1}\,,\qquad\text{with}\ \ \mathbf{h}:=\mathbf{w}^{1}\cdot\mathbf{\mathcal{I}}+\mathbf{b}\,, \tag{30}\]

where we made use of the derivative of the _Softplus_ function

\[\partial_{x}\,\mathcal{SP}(x)=\frac{e^{x}}{1+e^{x}}\,. \tag{31}\]

Actually, \(\partial_{x}\,\mathcal{SP}(x)\) is the widely used _Sigmoid_ activation function. This shows another benefit of applying the _Softplus_ activation function for the representation of hyperelastic potentials. Since also its derivative is a commonly used activation function, the derivatives of the NN potential benefit from the properties of activation functions [48], which are specifically chosen for the application in NNs. Note that in Eq. (30), the derivatives of the NN potential w.r.t. each single invariant also depends on the remaining invariants. The second derivative follows as

\[\partial_{(\mathbf{\mathcal{I}})_{i}(\mathbf{\mathcal{I}})_{j}}^{2}\,W^{\text{NN}}= \sum_{a=1}^{n}w_{a}^{2}\,\frac{e^{h_{a}}}{(1+e^{h_{a}})^{2}}\,w_{ai}^{1}\,w_{aj }^{1}\,. \tag{32}\]
Continuum formulation of finite deformation elasticity

While in the previous section, material specific constitutive models were discussed, in this section, general balance laws of finite deformation elasticity are introduced. For this, we consider the potential energy \(\Pi\), which is comprised of an internal part \(\Pi^{\rm int}\) and an external part \(\Pi^{\rm ext}\), respectively, such that

\[\Pi(\boldsymbol{\varphi})=\Pi^{\rm int}+\Pi^{\rm ext}\,, \tag{33}\]

where the external part is given by

\[\Pi^{\rm ext}=-\int_{\mathcal{B}_{0}}\bar{\mathbf{B}}\cdot\boldsymbol{\varphi }\ {\rm d}V-\int_{\partial\mathcal{B}_{0}^{T}}\bar{\mathbf{T}}\cdot\boldsymbol{ \varphi}\ {\rm d}A\,. \tag{34}\]

Therein, \(\bar{\mathbf{B}}:\mathcal{B}_{0}\to\mathbb{R}^{3}\) is the prescribed body force and \(\bar{\mathbf{T}}:\partial\mathcal{B}_{0}^{T}\to\mathbb{R}^{3}\) denotes the prescribed Piola-Kirchoff traction vector on \(\partial\mathcal{B}_{0}^{T}\subset\partial\mathcal{B}_{0}\), where the boundary is comprised of a Dirichlet boundary \(\partial\mathcal{B}_{0}^{\boldsymbol{\varphi}}\) with prescribed deformations \(\bar{\boldsymbol{\varphi}}:\partial\mathcal{B}_{0}^{\boldsymbol{\varphi}}\to \mathbb{R}^{3}\) and a Neumann boundary \(\partial\mathcal{B}_{0}^{T}\) with prescribed tractions \(\bar{\mathbf{T}}:\partial\mathcal{B}_{0}^{T}\) such that \(\partial\mathcal{B}_{0}=\partial\mathcal{B}_{0}^{T}\cup\partial\mathcal{B}_{0 }^{\boldsymbol{\varphi}}\) and \(\partial\mathcal{B}_{0}^{T}\cap\partial\mathcal{B}_{0}^{\boldsymbol{\varphi}}=\emptyset\). Furthermore, in Eq. (33), \(\Pi^{\rm int}\) basically denotes the elastic energy stored internally by the deformed body. This can be a purely displacement-based version

\[\Pi^{\rm int}=\int_{\mathcal{B}_{0}}W^{\rm F}(\mathbf{F})\ {\rm d}V\,, \tag{35}\]

which itself depends on the hyperelastic potential introduced in Eq. (6). In the above, the solution function \(\boldsymbol{\varphi}\in\mathbb{V}_{\boldsymbol{\varphi}}\) is subject to the space

\[\mathbb{V}_{\boldsymbol{\varphi}}=\{\boldsymbol{\varphi}:\mathcal{B}_{0}\to \mathbb{R}^{3}\,|\,\varphi_{i}\in H^{1}(\mathcal{B}_{0})\wedge\boldsymbol{ \varphi}=\bar{\boldsymbol{\varphi}}\ {\rm on}\ \partial\mathcal{B}_{0}^{ \boldsymbol{\varphi}}\wedge\det(\mathbf{F}_{\boldsymbol{\varphi}})>0\ {\rm in}\ \mathcal{B}_{0}\}\,, \tag{36}\]

where \(H^{1}\) denotes the Sobolev space. The principle of stationary potential energy requires the satisfaction of the stationary condition, such that

\[\delta\Pi=\delta\Pi^{\rm int}+\delta\Pi^{\rm ext}\stackrel{{!}}{{= }}0\,, \tag{37}\]

which has to hold for arbitrary \(\delta\boldsymbol{\varphi}\in\mathbb{V}_{\boldsymbol{\varphi}}^{0}\) with the space

\[\mathbb{V}_{\boldsymbol{\varphi}}^{0}=\{\delta\boldsymbol{\varphi}:\mathcal{ B}_{0}\to\mathbb{R}^{3}\,|\,\delta\varphi_{i}\in H^{1}(\mathcal{B}_{0})\wedge \delta\boldsymbol{\varphi}=\mathbf{0}\ {\rm on}\ \partial\mathcal{B}_{0}^{ \boldsymbol{\varphi}}\}\,, \tag{38}\]

where the variation of \(\Pi^{\rm ext}\) leads to

\[\delta\Pi^{\rm ext}=-\int_{\mathcal{B}_{0}}\bar{\mathbf{B}}\cdot\delta \boldsymbol{\varphi}\ {\rm d}V-\int_{\partial\mathcal{B}_{0}^{T}}\bar{\mathbf{T}}\cdot\delta \boldsymbol{\varphi}\ {\rm d}A\,. \tag{39}\]

The variation of the internal potential \(\Pi^{\rm int}\) will be dealt with subsequently after discussing different displacement-based and mixed formulations.

### Displacement-based formulations

In solid mechanics, the primary strain measure is the deformation gradient \(\mathbf{F}\). As already pointed out in Sect. 2.2, the notion of polyconvexity introduced by Ball [6, 7] suggests to not only consider the deformation gradient \(\mathbf{F}\) in advanced numerical methods, but also its cofactor \(\mathbf{H}\) and determinant \(J\), cf. Eq. (11). In particular, hyperelastic potentials are usually formulated in this extended set of arguments of \(\mathbf{F}\), i.e., they are formulated in (\(\mathbf{F}\), \(\mathbf{H}\), \(J\)). Considering this suggests to use the tensor cross product firstly introduced by [21], which facilitates both the analytical expressions and the
implementation of advanced numerical methods, as it leads to neat derivatives of (**F**, **H**, \(J\)) w.r.t. **F**. This has already been realized in purely displacement-based and mixed methods [13], as well as energy-momentum schemes [9]. Starting with the former, [13] proposed a new computational framework based on polyconvex hyperelastic potentials of the form

\[\hat{W}^{\mathrm{F}}\colon\mathbb{T}_{3}^{+}\times\mathbb{T}_{3}^{+}\times \mathbb{R}_{+}\to\mathbb{R}\,,\qquad(\mathbf{F},\mathbf{H},J)\mapsto\hat{W}^{ \mathrm{F}}(\mathbf{F},\mathbf{H},J)\,. \tag{40}\]

The above stored energy function inspired the objective version thereof, given by

\[\hat{W}^{\mathrm{C}}\colon\mathbb{S}_{+}^{3}\times\mathbb{S}_{+}^{3}\times \mathbb{R}_{+}\to\mathbb{R}\,,\qquad(\mathbf{C},\mathbf{G},C)\mapsto\hat{W}^{ \mathrm{C}}(\mathbf{C},\mathbf{G},C)\,. \tag{41}\]

Introduced by [9], this formulation facilitates the design of an energy-momentum scheme (EMS).

Furthermore, as has been outlined in Sect. 2.3, hyperelastic potentials are usually formulated in terms of strain invariants, which basically are nonlinear functions of (**F**, **H**, \(J\)). Thus, the invariant-based formulation

\[W^{\mathrm{I}}\colon\mathbb{R}_{+}\times\mathbb{R}_{+}\times\mathbb{R}_{+}\to \mathbb{R}\,,\qquad(I_{\mathbf{C}},II_{\mathbf{C}},J)\mapsto W^{\mathrm{I}}(I_ {\mathbf{C}},II_{\mathbf{C}},J)\,, \tag{42}\]

can be seen as a consecutive development of established polyconvexity-inspired numerical methods [9, 13], now setting the focus directly on invariants. Variation of the displacement-based potential energy \({}^{W^{\mathrm{I}}}\Pi^{\mathrm{int}}\), which depends on the strain energy \(W^{\mathrm{I}}\) as defined in Eq. (42), yields

\[\begin{split}\delta(^{W^{\mathrm{I}}}\Pi^{\mathrm{int}})& =\int_{\mathcal{B}_{0}}\delta W^{\mathrm{I}}(I_{\mathbf{C}},\,II_{ \mathbf{C}},\,J)\,\,\mathrm{d}V\\ &=\int_{\mathcal{B}_{0}}2\,\left(\partial_{I_{\mathbf{C}}}W^{ \mathrm{I}}\,\mathbf{I}+\partial_{II_{\mathbf{C}}}W^{\mathrm{I}}\,\mathbf{I} \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!
To circumvent this issue, a variety of mixed formulations have been proposed. The purpose of the present study is not only to demonstrate the ability of the PANN constitutive model to be employed within a mixed formulation, but also to develop a mixed method which is tailored for the PANN model. As in the previous section, the mixed formulations discussed here are inspired by the notion of polyconvexity, in the sense that they consider the arguments of the polyconvexity condition as additional unknowns, cf. Eq. (11).

In this regard, the Hu-Washizu like 7-field mixed formulation from [13] corresponding to the displacement-based version Eq. (40) is worth mentioning. There, the consistency of \((\mathbf{F},\mathbf{H},J)\) with \(\boldsymbol{\varphi}\) is enforced by means of the Lagrange multipliers \((\boldsymbol{\Lambda}^{\mathbf{F}},\boldsymbol{\Lambda}^{\mathbf{H}},\Lambda^ {J})\), respectively, and the total set of unknowns is \((\boldsymbol{\varphi},\mathbf{F},\mathbf{H},J,\boldsymbol{\Lambda}^{\mathbf{F} },\boldsymbol{\Lambda}^{\mathbf{H}},\Lambda^{J})\). A further notable version is proposed by [9], which corresponds to the displacement-based version Eq. (41). There, the consistency of \((\mathbf{C},\mathbf{G},C)\) is enforced by means of the Lagrange multipliers \((\boldsymbol{\Lambda}^{\mathbf{C}},\boldsymbol{\Lambda}^{\mathbf{G}},\Lambda^ {C})\), with the total set of unknowns \((\boldsymbol{\varphi},\mathbf{C},\mathbf{G},C,\boldsymbol{\Lambda}^{\mathbf{ C}},\boldsymbol{\Lambda}^{\mathbf{G}},\Lambda^{C})\) employed. The extension to enhanced assumed strain (EAS) variants of the above mentioned 7-field mixed formulation should be mentioned and can also be employed in a straight forward manner, cf. [86]. Here, we restrict our consideration without loss of generality to the 5-field mixed formulation \({}^{W}\Gamma_{5}^{\text{int}}\) proposed by [59], which the authors in [59] claim to be robust, efficient and locking-free. It is given by2

Footnote 2: In this subsection, the subscript \(\boldsymbol{\varphi}\) indicates that a quantity is calculated through the displacement field rather than being an additional unknown.

\[{}^{W}{}^{\text{G}}\Pi_{5}^{\text{int}}=\int_{\mathcal{B}_{0}}W^{\text{G}}( \mathbf{C}_{\boldsymbol{\varphi}},\mathbf{G},J)+\boldsymbol{\Lambda}^{\mathbf{ G}}:(\operatorname{cof}(\mathbf{C}_{\boldsymbol{\varphi}})-\mathbf{G})+ \Lambda^{J}\left(\operatorname{det}(\mathbf{F}_{\boldsymbol{\varphi}})-J \right)\,\mathrm{d}V\,, \tag{44}\]

where \(\boldsymbol{\varphi}\in\mathds{V}_{\boldsymbol{\varphi}}\). Furthermore, \(\mathbf{G}\in\mathds{V}_{\boldsymbol{\Lambda}}\) and \(J\in\mathds{V}_{A}\) are enforced by means of Lagrange multipliers \(\boldsymbol{\Lambda}^{\mathbf{G}}\in\mathds{V}_{\boldsymbol{\Lambda}}\) and \(\Lambda^{J}\in\mathds{V}_{A}\) with generalized spaces

\[\mathds{V}_{\boldsymbol{\Lambda}}=\mathds{V}_{\boldsymbol{\Lambda}}^{0}=\{ \mathbf{A}:\mathcal{B}_{0}\to\mathbb{S}_{+}^{3}\,|\,A_{ij}\in\mathds{L}_{2}( \mathcal{B}_{0})\},\quad\mathds{V}_{A}=\mathds{V}_{A}^{0}=\{A:\mathcal{B}_{0} \to\mathbb{R}_{+}\,|\,A\in\mathds{L}_{2}(\mathcal{B}_{0})\}\,. \tag{45}\]

The corresponding variation yields

\[\begin{split}\delta({}^{W}{}^{\text{G}}\Pi_{5}^{\text{int}})=\int _{\mathcal{B}_{0}}& 2\,(\partial\mathbf{C}_{\boldsymbol{\varphi}}W^{ \text{G}}+\boldsymbol{\Lambda}^{\mathbf{G}}\,\boldsymbol{\Re}\,\mathbf{C}_{ \boldsymbol{\varphi}}+\tfrac{1}{2}\,\Lambda^{J}\,J_{\boldsymbol{\varphi}}^{-1 }\,\mathbf{G}_{\boldsymbol{\varphi}}):\tfrac{1}{2}\,\delta\mathbf{C}_{ \boldsymbol{\varphi}}\\ &+(\partial_{\mathbf{G}}W^{\text{G}}-\boldsymbol{\Lambda}^{ \mathbf{G}}):\delta\mathbf{G}\\ &+(\partial_{J}W^{\text{G}}-\Lambda^{J})\,\delta J\\ &+(\operatorname{cof}(\mathbf{C}_{\boldsymbol{\varphi}})-\mathbf{ G}):\delta\boldsymbol{\Lambda}^{\mathbf{G}}\\ &+(\operatorname{det}(\mathbf{F}_{\boldsymbol{\varphi}})-J)\, \delta\Lambda^{C}\,\,\mathrm{d}V\,,\end{split} \tag{46}\]

which has to hold for arbitrary \(\delta\boldsymbol{\varphi}\in\mathds{V}_{\boldsymbol{\varphi}}^{0}\), \(\delta\mathbf{G},\delta\boldsymbol{\Lambda}^{\mathbf{G}}\in\mathds{V}_{ \boldsymbol{\Lambda}}^{0}\), and \(\delta J,\delta\Lambda^{J}\in\mathds{V}_{A}^{0}\). Based on the above newly proposed displacement formulation in Eq. (42), we also suggest a new Hu-Washizu like 7-field mixed formulation \({}^{W}{}^{\text{I}}\Pi_{7}^{\text{int}}\) given by3

Footnote 3: The initial idea to use a strain energy depending solely on the invariants of \(\mathbf{C}\), i.e. \(W=W(I_{\mathbf{C}},II_{\mathbf{C}},III_{\mathbf{C}})\) resulted in the occurrence of hourglassing phenomena for nearly incompressible materials under shear loading conditions. Consequently, this approach was deemed unsuitable and subsequently discarded.

\[\begin{split}{}^{W}{}^{\text{I}}\Pi_{7}^{\text{int}}=\int_{ \mathcal{B}_{0}}W^{\text{I}}(I_{\mathbf{C}},II_{\mathbf{C}},J)&+ \lambda^{I_{\mathbf{C}}}\left(\operatorname{tr}\mathbf{C}_{\boldsymbol{\varphi}} -I_{\mathbf{C}}\right)\\ &+\lambda^{II_{\mathbf{C}}}\left(\operatorname{tr}(\operatorname{ cof}\mathbf{C}_{\boldsymbol{\varphi}})-II_{\mathbf{C}}\right)\\ &+\lambda^{J}\left(\operatorname{det}\mathbf{F}_{\boldsymbol{ \varphi}}-J\right)\,\mathrm{d}V\,,\end{split} \tag{47}\]
with \(\boldsymbol{\varphi}\in\mathds{V}_{\boldsymbol{\varphi}}\) and \(I_{\mathbf{C}},II_{\mathbf{C}},J,\lambda^{I_{\mathbf{C}}},\lambda^{I_{\mathbf{I }}}\mathrm{C},\lambda^{J}\in\mathds{V}_{A}\). The above formulation can also be regarded as tailor made for the employed PANN constitutive model. Variation of this Hu-Washizu like 7-field mixed formulation yields

\[\begin{split}\delta(^{W^{\mathrm{I}}}\Pi_{7}^{\mathrm{int}})= \int_{\mathcal{B}_{0}}& 2\,(\lambda^{I_{\mathbf{C}}}\,\mathbf{I}+\lambda^{I \mathrm{I}\mathrm{C}}\,(\mathbf{I}\boldsymbol{\$}\,\mathbf{\$}\,\mathbf{C}_{ \boldsymbol{\varphi}})+\tfrac{1}{2}\,\lambda^{J}\,J_{\boldsymbol{\varphi}}^{- 1}\,\mathbf{G}_{\boldsymbol{\varphi}}):\tfrac{1}{2}\,\delta\mathbf{C}_{ \boldsymbol{\varphi}}\\ &+(\partial_{I_{\mathbf{C}}}W^{\mathrm{I}}-\lambda^{I_{\mathbf{C} }})\,\delta I_{\mathbf{C}}\\ &+(\partial_{II_{\mathbf{C}}}W^{\mathrm{I}}-\lambda^{I_{\mathbf{ C}}})\,\delta II_{\mathbf{C}}\\ &+(\partial_{J}W^{\mathrm{I}}-\lambda^{J})\,\delta J\\ &+(\mathrm{tr}\,\mathbf{C}_{\boldsymbol{\varphi}}-I_{\mathbf{C} })\,\delta\lambda^{I_{\mathbf{C}}}\\ &+(\mathrm{tr}\,\mathbf{G}_{\boldsymbol{\varphi}}-II_{\mathbf{C} })\,\delta\lambda^{I_{\mathbf{C}}}\\ &+(\mathrm{det}\,\mathbf{F}_{\boldsymbol{\varphi}}-J)\,\delta \lambda^{J}\,\,\mathrm{d}V\,,\end{split} \tag{48}\]

which has to hold for arbitrary \(\delta\boldsymbol{\varphi}\in\mathds{V}_{\boldsymbol{\varphi}}^{0}\) and \(\delta I_{\mathbf{C}},\delta II_{\mathbf{C}},\delta J,\delta\lambda^{I_{ \mathbf{C}}},\delta\lambda^{II_{\mathbf{C}}},\delta\lambda^{J}\in\mathds{V}_{ A}^{0}\).

### Dynamic formulation

For the dynamic formulation, we restrict our consideration to the displacement-based formulation given in Eq. (43). In so doing, we introduce the time interval of interest with range \(\mathcal{R}=[0,T]\), where \(T\in\mathbb{R}_{+}\). The motion of the deformation field \(\boldsymbol{\varphi}_{t}\in\mathds{V}_{\boldsymbol{\varphi}}\) denotes the placement of a particle at \(\mathbf{X}\in\mathcal{B}_{0}\) at time \(t\in\mathcal{R}\). We introduce the velocity field \(\mathds{V}_{\mathbf{v}}\ni\mathbf{v}_{t}:\mathcal{B}_{0}\times\mathcal{R} \rightarrow\mathbb{R}^{3}\) with space

\[\mathds{V}_{\mathbf{v}}=\left\{\mathbf{v}:\mathcal{B}_{0}\times\mathcal{R} \rightarrow\mathbb{R}^{3}\,|\,v_{i}\in\mathds{L}_{2}(\mathcal{B}_{0})\right\}, \tag{49}\]

and the mass density in the reference configuration \(\rho_{0}:\mathcal{B}_{0}\rightarrow\mathbb{R}^{+}\). With that we are able to introduce the action functional

\[S(\boldsymbol{\varphi}_{t},\mathbf{v}_{t})=\int_{t_{0}}^{t_{e}}\int_{\mathcal{ B}_{0}}(\boldsymbol{\dot{\varphi}}_{t}-\tfrac{1}{2}\,\mathbf{v}_{t})\cdot \mathbf{v}_{t}\,\rho_{0}\,\,\mathrm{d}V-\left({}^{W^{\mathrm{I}}}\Pi^{\mathrm{ int}}+\Pi^{\mathrm{ext}}\right)\,\mathrm{d}V\,\,\mathrm{d}t\,. \tag{50}\]

By means of Hamilton's principle, the stationary conditions of the action \(S\) w.r.t. the variations yield

\[\begin{split}\int_{\mathcal{B}_{0}}\delta\mathbf{v}_{t}\cdot( \boldsymbol{\dot{\varphi}}_{t}-\mathbf{v}_{t})\,\rho_{0}\,\,\mathrm{d}V& =0\\ \int_{\mathcal{B}_{0}}\delta\boldsymbol{\varphi}\cdot\rho_{0}\, \dot{\mathbf{v}}_{t}+\mathbf{S}_{t}:\tfrac{1}{2}\,\delta\mathbf{C}\,\, \mathrm{d}V+\Pi^{\mathrm{ext}}(\delta\boldsymbol{\varphi})&=0\,, \end{split} \tag{51}\]

which have to hold for arbitrary \(\delta\mathbf{v}\in\mathds{V}_{\mathbf{v}}^{0}\) and \(\delta\boldsymbol{\varphi}\in\mathds{V}_{\boldsymbol{\varphi}}\) with space

\[\mathds{V}_{\mathbf{v}}^{0}=\left\{\mathbf{v}:\mathcal{B}_{0}\to \mathbb{R}^{3}\,|\,v_{i}\in\mathds{L}_{2}(\mathcal{B}_{0})\right\}. \tag{52}\]

In the above equation, due to Eq. (43), the PK2 stress tensor is defined as

\[\mathbf{S}_{t}=2\,\left(\partial_{I_{\mathbf{C}}}W^{\mathrm{I}}\,\mathbf{I}+ \partial_{II_{\mathbf{C}}}W^{\mathrm{I}}\,\mathbf{I}\boldsymbol{\$}\,\mathbf{ \$}\,\mathbf{\$}\,\mathbf{C}_{t}+\tfrac{1}{2}\,\partial_{J}W^{\mathrm{I}}\,J_{ t}^{-1}\,\mathbf{G}_{t}\right)\,. \tag{53}\]

The verification of the balance laws for the time-continuous system is detailed in Appx. A, which serves as a comprehensive reference for ensuring the validity and consistency of the sytem's fundamental principles.
Time discretization with consistent energy-momentum scheme

Regarding the time discretization, the use of the midpoint rule in nonlinear continuum mechanics typically introduces a violation of the energy consistency solely due to the discretization of the variation of the (strain) energy within the weak form, while all other terms involved do not give rise to any issue. In order to circumvent this issue and to obtain an energy and momentum consistent time integration scheme, a midpoint-type discretization with application of the concept of the discrete gradient, as defined in [33], is employed. Specifically, for the weak form in Eq. (51) the following steps are pursued to ensure energy and momentum consistency:

* substitution of time rates (\(\dot{\bullet}\)) with \(\frac{(\bullet)_{n+1}-(\bullet)_{n}}{\Delta t}=\frac{\Delta(\bullet)}{\Delta t}\),
* algorithmic evaluation of the variation of the (strain) energy \(W^{\rm I}\), which leads to an algorithmic PK2 stress tensor \(\mathbf{S}_{\rm algo}\) based on the concept of the discrete gradients in the sense of [33],
* midpoint evaluation of remaining terms, i.e. \(\big{(}\bullet\big{)}_{n+\frac{1}{2}}=\frac{1}{2}\left((\bullet)_{n}+(\bullet )_{n+1}\right)\).

Employing the above steps to the weak form Eq. (51) yields the semi-discrete weak form as

\[\begin{split}\int_{\mathcal{B}_{0}}\delta\mathbf{v}\cdot\frac{ \Delta\boldsymbol{\varphi}}{\Delta t}\,\rho_{0}\ \mathrm{d}V&=\int_{\mathcal{B}_{0}}\delta\mathbf{v}\cdot \mathbf{v}_{n+\frac{1}{2}}\,\rho_{0}\ \mathrm{d}V\\ \int_{\mathcal{B}_{0}}\delta\boldsymbol{\varphi}\cdot\rho_{0} \,\frac{\Delta\mathbf{v}}{\Delta t}\ \mathrm{d}V&=-\int_{\mathcal{B}_{0}}\mathbf{S}_{\rm algo}: \frac{1}{2}\left(\delta\mathbf{F}^{\rm T}\,\mathbf{F}_{n+\frac{1}{2}}+ \mathbf{F}_{n+\frac{1}{2}}^{\rm T}\,\delta\mathbf{F}\right)\ \mathrm{d}V-\Pi_{n+\frac{1}{2}}^{\rm ext}(\delta \boldsymbol{\varphi})\,.\end{split} \tag{54}\]

Therein, the particular form of the algorithmic version of the PK2 stress tensor \(\mathbf{S}_{\rm algo}\) follows from the so-called directionality property, which is a sufficient condition for energy consistency. The directionality property itself serves as a time-discrete representation of the continuum stress power \(\dot{W^{\rm I}}=\mathbf{S}:\frac{1}{2}\,\dot{\mathbf{C}}\), acting as a sufficient criterion for energy conservation, cf. [91]. In that regard, the directionality property for the present system yields

\[\begin{split}\mathbf{S}_{\rm algo}:\frac{1}{2}\,\Delta\mathbf{C}& =\mathrm{D}_{I_{\mathbf{C}}}W^{\rm I}\,\Delta I_{\mathbf{C}}+ \mathrm{D}_{II_{\mathbf{C}}}W^{\rm I}\,\Delta II_{\mathbf{C}}+\mathrm{D}_{J}W ^{\rm I}\,\Delta J_{\boldsymbol{\varphi}}\\ &=W^{\rm I}(I_{\mathbf{C}_{n+1}},II_{\mathbf{C}_{n+1}},J_{n+1})-W ^{\rm I}(I_{\mathbf{C}_{n}},II_{\mathbf{C}_{n}},J_{n})\,.\end{split} \tag{55}\]

The satisfaction of this condition motivates for proposing the algorithmic PK2 stress tensor as

\[\mathbf{S}_{\rm algo}=2\,\left(\mathrm{D}_{I_{\mathbf{C}}}W^{\rm I}\,\mathbf{ I}+\mathrm{D}_{II_{\mathbf{C}}}W^{\rm I}\,\mathbf{I}\boldsymbol{\Phi}\, \mathbf{C}_{\rm algo}+\frac{1}{2}\,\mathrm{D}_{J}W^{\rm I}\,(J_{\rm algo})^{-1 }\mathbf{G}_{\rm algo}\right)\,, \tag{56}\]

which serves as consistent time-discrete variant of Eq. (53) and ensures energy consistency while maintaining the second-order accuracy of the midpoint rule. In case of midpoint time integration, the only difference compared to the EM scheme is the evaluation of the PK2 stress tensor. In particular, instead of Eq. (56) we would evaluate the PK2 stress tensor as follows

\[\begin{split}\mathbf{S}_{n+\frac{1}{2}}&=2\,(\partial _{I_{\mathbf{C}}}W^{\rm I}(\mathbf{F}_{n+\frac{1}{2}})\,\mathbf{I}+\partial _{II_{\mathbf{C}}}W^{\rm I}(\mathbf{F}_{n+\frac{1}{2}})\,\mathbf{I} \boldsymbol{\Phi}\mathbf{C}((\mathbf{F}_{n+\frac{1}{2}})\\ &\quad+\frac{1}{2}\,\partial_{J}W^{\rm I}(\mathbf{F}_{n+\frac{1}{ 2}})\,(J((\mathbf{F}_{n+\frac{1}{2}}))^{-1}\mathbf{G}(\mathbf{F}_{n+\frac{1}{ 2}}))\,.\end{split} \tag{57}\]

The algorithmic form of the PK2 stress tensor in Eq. (56) is the basic feature of the energy-momentum method and contains the subsequently introduced partitioned discrete derivatives \(\mathrm{D}_{I_{\mathbf{C}}}W^{\rm I},\mathrm{D}_{II_{\mathbf{C}}}W^{\rm I}\), \(\mathrm{D}_{J}W^{\rm I}\) in replacement of the partial derivatives of the strain energy from Eq. (53). The construction
of these particular discrete gradients is provided by Appx. C, which eventually yields the final Greenspan-like versions

\[\begin{split} D_{I\mathbf{C}}W^{1}=&\left(W^{\mathrm{I} }(I_{\mathbf{C}_{n+1}},II_{\mathbf{C}_{n}},J_{n})-W^{\mathrm{I}}(I_{\mathbf{C} _{n}},II_{\mathbf{C}_{n}},J_{n})\right.\\ &\left.+W^{\mathrm{I}}(I_{\mathbf{C}_{n+1}},II_{\mathbf{C}_{n+1}}, J_{n+1})-W^{\mathrm{I}}(I_{\mathbf{C}_{n}},II_{\mathbf{C}_{n+1}},J_{n+1})\right)/(2\, \Delta I_{\mathbf{C}})\,,\\ D_{II\mathbf{C}}W^{1}=&\left(W^{\mathrm{I}}(I_{ \mathbf{C}_{n+1}},II_{\mathbf{C}_{n+1}},J_{n})-W^{\mathrm{I}}(I_{\mathbf{C}_{n +1}},II_{\mathbf{C}_{n}},J_{n})\right.\\ &\left.+W^{\mathrm{I}}(I_{\mathbf{C}_{n}},II_{\mathbf{C}_{n+1}}, J_{n+1})-W^{\mathrm{I}}(I_{\mathbf{C}_{n}},II_{\mathbf{C}_{n}},J_{n+1})\right)/(2\, \Delta II_{\mathbf{C}})\,,\\ D_{J}W^{1}=&\left(W^{\mathrm{I}}(I_{\mathbf{C}_{n+1 }},II_{\mathbf{C}_{n+1}},J_{n+1})-W^{\mathrm{I}}(I_{\mathbf{C}_{n+1}},II_{ \mathbf{C}_{n+1}},J_{n})\right.\\ &\left.+W^{\mathrm{I}}(I_{\mathbf{C}_{n}},II_{\mathbf{C}_{n}},J_{n +1})-W^{\mathrm{I}}(I_{\mathbf{C}_{n}},II_{\mathbf{C}_{n}},J_{n})\right)/(2\, \Delta J)\,.\end{split} \tag{58}\]

As shown in Eq. (56), algorithmic versions of the kinematic relations \(\mathbf{C}_{\mathrm{algo}}\), \(\mathbf{G}_{\mathrm{algo}}\), \(J_{\mathrm{algo}}\) are employed, which directly follows from the directionality condition Eq. (C.17). Alternatively, the algorithmic PK2 tensor from Eq. (56) with its algorithmic kinematic relations can be kind of naturally achieved when deriving it from the corresponding dynamic Hu-Washizu like mixed formulation given based on Eq. (41), see [9] for more details. In any case, the algorithmic kinematic relations \(\mathbf{C}_{\mathrm{algo}}\), \(\mathbf{G}_{\mathrm{algo}}\), \(J_{\mathrm{algo}}\) are chosen as

\[\begin{split}\mathbf{C}_{\mathrm{algo}}&=\frac{1}{2} \left(\mathbf{C}_{n+1}+\mathbf{C}_{n}\right)\\ \mathbf{G}_{\mathrm{algo}}&=\frac{1}{3}\left(\mathbf{ C}_{\mathrm{algo}}\boldsymbol{\Psi}\mathbf{C}_{\mathrm{algo}}+\widetilde{ \mathbf{G}}_{\mathrm{algo}}\right)\text{,\quad with \ }\widetilde{\mathbf{G}}_{\mathrm{algo}}=\frac{1}{2}\left(\mathbf{G}_{n+1}+ \mathbf{G}_{n}\right),\\ J_{\mathrm{algo}}&=\frac{1}{2}\left(J_{n+1}+J_{n} \right),\end{split} \tag{59}\]

and fulfill the necessary directionality property of Eq. (55). It is important to remark that the discrete gradients possess the following properties, see [33] for more information. In particular, the discrete gradients

* satisfy the directionality property Eq. (C.17), which serves as sufficient condition for energy consistency,
* are well defined in the limit \(\|\Delta\Pi^{i}\|\to 0\) where \(\|\Delta\Pi^{i}\|=\sqrt{\langle\Delta\Pi^{i},\Delta\Pi^{i}\rangle}\), \(\Delta\Pi^{i}=\Pi^{i}_{n+1}-\Pi^{i}_{n}\), here \(\boldsymbol{\Pi}=(\Pi^{1},\Pi^{2},\Pi^{3})=(I_{\mathbf{C}},II_{\mathbf{C}},J)\), cf. Appx. C, and \(\langle\bullet,\bullet\rangle\) denotes the inner product,
* are a second-order approximation of the corresponding midpoint evaluations of the partial derivatives,
* and are undefined for handling zero magnitude \(\|\Delta\Pi^{i}\|=0\), which regarding the implementation occurs e.g. within a very first iteration of Newton's method, and can be avoided by midpoint-type evaluation of the partial derivatives \(\mathrm{D}_{\Pi^{i}}W^{\mathrm{I}}=\partial_{\Pi^{i}}W^{\mathrm{I}}(\Pi^{i}_{n +\frac{1}{2}})\),
* results with the aforementioned algorithmic stress evaluation in the generation of a tangent stiffness matrix that exhibits non-symmetry.

**Remark 1**.: _As becomes obvious in Appx. D, the choice of the strain energy \(\hat{W}^{C}(\mathbf{C},\mathbf{G},C)\) in Eq. (41) leads to far more complex versions of the discrete gradients \(\mathrm{D}_{\mathbf{C}}\hat{W}^{C}\), \(\mathrm{D}_{\mathbf{G}}\hat{W}^{C}\) and \(\mathrm{D}_{C}\hat{W}^{C}\) when compared to the choice of strain energy \(W^{I}(I_{\mathbf{C}},II_{\mathbf{C}},J)\) in Eq. (42) with just Greenspan versions of the partitioned discrete gradients \(\mathrm{D}_{I\mathbf{C}}W^{I}\), \(\mathrm{D}_{II\mathbf{C}}W^{I}\) and \(\mathrm{D}_{J}W^{I}\) in Eq. (58). Also the traditional, projection-based discrete gradient \(\mathrm{D}_{\mathbf{C}}W^{C}(\mathbf{C})\) associated with strain energy \(W^{C}(\mathbf{C})\) from Eq. (41) provided by Appx. E can be employed. However, the consistent linearization of the tensor-valued
discrete gradient \(\mathrm{D}_{\mathbf{C}}W^{C}(\mathbf{C})\) and the application of the chain rule become more complex when employing the PANN constitutive model with its inputs \(\boldsymbol{\mathcal{I}}\) in comparison to the proposed version outlined in Eq. (58)._

Standard integrators like the 2nd order accurate midpoint rule usually violate the conserving properties, particularly the balance of angular momentum and total energy. However, the proposed EM scheme ensures that these balance laws are preserved in the time-discrete system. The verification and detailed analysis of the balance laws for the time-discrete system can be found in Appx. B.

## 6 Spatial discretization with the finite element method

In the following, we focus on the semi-discrete displacement-based system from Sect. 5, which makes a variety of spatial discretizations possible. Here, we employ the finite element (FE) method for the spatial discretization of Eq. (54). To this end, the body \(\mathcal{B}_{0}\) is subdivided into \(n_{e}\in\mathbb{N}\) finite elements

\[\mathcal{B}_{0}\approx\mathcal{B}_{0}^{\mathrm{h}}=\bigcup_{e=1}^{n_{e}} \mathcal{B}_{0}^{e}\,, \tag{60}\]

where \(e=1,...,n_{e}\) denotes the respective element. In particular, we use the isoparametric concept, cf. [49], by using the same shape functions for the approximation of the geometry with

\[\mathbf{X}^{\mathrm{h}}\big{|}_{\mathcal{B}_{0}^{e}}=\sum_{a=1}^{n_{\mathrm{ node}}}N^{a}\,\mathbf{X}^{a}\,, \tag{61}\]

and the deformation field \(\boldsymbol{\varphi}_{t}^{h}\in\mathds{V}_{\boldsymbol{\varphi}}^{\mathrm{h} }\subset\mathds{V}_{\boldsymbol{\varphi}}\) and the velocity field \(\mathbf{v}_{t}^{h}\in\mathds{V}_{\mathbf{v}}^{\mathrm{h}}\subset\mathds{V}_ {\mathbf{v}}\) with spaces

\[\mathds{V}_{\boldsymbol{\varphi}}^{\mathrm{h}}=\left\{\boldsymbol{\varphi}_{t }\in\mathds{V}_{\boldsymbol{\varphi}}\,:\,\boldsymbol{\varphi}_{t}^{\mathrm{ h}}\big{|}_{\mathcal{B}_{0}^{e}}=\sum_{a=1}^{n_{\mathrm{node}}}N^{a} \boldsymbol{\varphi}_{t}^{a}\right\},\quad\mathds{V}_{\mathbf{v}}^{\mathrm{h} }=\left\{\mathbf{v}_{t}\in\mathds{V}_{\mathbf{v}}\,:\,\mathbf{v}_{t}^{\mathrm{ h}}\big{|}_{\mathcal{B}_{0}^{e}}=\sum_{a=1}^{n_{\mathrm{node}}}N^{a}\mathbf{v}_{t}^{a} \right\}\,, \tag{62}\]

where \(N^{a}:\mathcal{B}_{0}\to\mathbb{R}\) denote nodal shape functions with corresponding nodal values \(\mathbf{X}^{a}\in\mathbb{R}^{3}\), \(\boldsymbol{\varphi}_{a}:\mathcal{R}\to\mathbb{R}^{3}\), and \(\mathbf{v}_{a}:\mathcal{R}\to\mathbb{R}^{3}\). According to a Bubnov-Galerkin FE ansatz, we employ the same Ansatz functions for solution and variations \(\delta\boldsymbol{\varphi}\in\mathds{V}_{\boldsymbol{\varphi}}^{0,\mathrm{h} }\subset\mathds{V}_{\boldsymbol{\varphi}}^{0}\), and \(\delta\mathbf{v}\in\mathds{V}_{\mathbf{v}}^{0,\mathrm{h}}\subset\mathds{V}_{ \mathbf{v}}^{0}\) with spaces

\[\mathds{V}_{\boldsymbol{\varphi}}^{0,\mathrm{h}}=\left\{\delta\boldsymbol{ \varphi}\in\mathds{V}_{\boldsymbol{\varphi}}^{0}\,:\,\delta\boldsymbol{ \varphi}^{\mathrm{h}}\big{|}_{\mathcal{B}_{0}^{e}}=\sum_{a=1}^{n_{\mathrm{ node}}}N^{a}\,\delta\boldsymbol{\varphi}^{a}\right\},\quad\mathds{V}_{ \mathbf{v}}^{0,\mathrm{h}}=\left\{\delta\mathbf{v}\in\mathds{V}_{\mathbf{v}}^{ 0}\,:\,\delta\mathbf{v}^{\mathrm{h}}\big{|}_{\mathcal{B}_{0}^{e}}=\sum_{a=1}^{ n_{\mathrm{node}}}N^{a}\delta\mathbf{v}^{a}\right\}\,. \tag{63}\]

We insert the above discretized solution and variations into the semi-discrete weak form Eq. (54) and obtain

\[\begin{split}\boldsymbol{\bigwedge}_{e=1}^{n_{e}}\sum_{a=1}^{n_{ \mathrm{node}}}\delta\mathbf{v}^{a}\cdot\mathbf{R}_{\mathbf{v}}^{a,e}& =0\,,\\ \boldsymbol{\bigwedge}_{e=1}^{n_{e}}\sum_{a=1}^{n_{\mathrm{ node}}}\delta\boldsymbol{\varphi}^{a}\cdot\mathbf{R}_{\boldsymbol{\varphi}}^{a,e}& =0\,,\end{split} \tag{64}\]

where \(\mathbf{A}_{e=1}^{n_{e}}\) denotes the assembly operator. In the above, the nodal residuals \(\mathbf{R}_{\mathbf{v}}^{a,e}\in\mathbb{R}^{3}\) and \(\mathbf{R}_{\boldsymbol{\varphi}}^{a,e}\in\mathbb{R}^{3}\) are defined as

\[\begin{split}\mathbf{R}_{\mathbf{v}}^{a,e}&=\mathbf{ M}^{ab,e}\,\left(\frac{1}{\Delta t}\Delta\boldsymbol{\varphi}^{b}-\mathbf{v}_{n+ \frac{1}{2}}^{b}\right)\,,\\ \mathbf{R}_{\boldsymbol{\varphi}}^{a,e}&=\mathbf{M}^{ ab,e}\,\frac{1}{\Delta t}\Delta\mathbf{v}^{b}+\int_{\mathcal{B}_{0}^{e}} \mathbf{B}a^{\mathrm{T}}(\boldsymbol{\varphi}_{n+\frac{1}{2}}^{b})\mathbf{S}_{ \mathrm{algo}}^{\mathrm{h,V}}\,\mathrm{d}V-\int_{\mathcal{B}_{0}^{e}}N^{a} \bar{\mathbf{B}}^{\mathrm{h}}(\boldsymbol{\varphi}_{n+\frac{1}{2}}^{b})\, \mathrm{d}V\,,\end{split} \tag{65}\]
where \(\mathbf{B}^{a}\in\mathbb{R}^{6\times 3}\) denotes the standard nodal operator matrix, \((\bullet)^{\mathrm{V}}\) symbolizes Voigt's vector notation and \(\mathbf{M}^{ab,e}\) denotes the elemental mass matrix defined as

\[\mathbf{M}^{ab,e}=\int_{\mathcal{B}^{e}_{0}}\rho_{0}\,N^{a}\,N^{b}\;\mathrm{d}V \,\mathbf{I}\,. \tag{66}\]

In order to pursue an efficient implementation, cf. [9], Eq. (65)\({}_{1}\) can be rewritten as

\[\mathbf{v}^{b}_{n+1}=\frac{2}{\Delta t}\Delta\mathbf{\varphi}^{b}_{n+1}-\mathbf{v }^{b}_{n}\,, \tag{67}\]

and inserting it into Eq. (65)\({}_{2}\) yields the nodal residual vector \(\widehat{\mathbf{R}}^{a,e}_{\mathbf{\varphi}}\in\mathbb{R}^{3}\) given by

\[\widehat{\mathbf{R}}^{a,e}_{\mathbf{\varphi}}=\frac{2}{\Delta t}\,\left(\frac{1}{ \Delta t}\,\mathbf{M}^{ab,e}\,\Delta\mathbf{\varphi}^{b}-\mathbf{M}^{ab,e}\, \mathbf{v}^{b}_{n}\right)+\mathbf{B}^{a^{\mathrm{T}}}_{n+\frac{1}{2}}\mathbf{ S}^{\mathrm{h,V}}_{\mathrm{algo}}\,\mathrm{d}V-\int_{\mathcal{B}^{e}_{0}}N^{a} \bar{\mathbf{B}}^{\mathrm{h}}(\mathbf{\varphi}^{h}_{n+\frac{1}{2}})\,\mathrm{d}V\,. \tag{68}\]

Accordingly, the velocity field is condensed out and we have a pure displacement-based formulation. Note that, for the midpoint rule, the algorithmic PK2 stress tensor \(\mathbf{S}^{\mathrm{h}}_{\mathrm{algo}}\) is simply replaced with the midpoint evaluated PK2 stress tensor from Eq. (57).

It is important to remark that the discretization of the mixed approaches in Sect. 4.2 is straight forward, where of course the additional fields have to be discretized with suitable spaces considering the inf-sup conditions for the particular choice of approximation formulas, see [59] for more information.

## 7 Numerical examples

In this section, stability and robustness of the proposed framework are examined in representative finite element simulations. In Sect. 7.1, the PANN constitutive model is calibrated to data generated with an analytical potential. Afterwards in Sect. 7.2 and in Sect. 7.3, the PANN constitutive model is applied in both static and dynamic numerical examples.

### Preparation of the neural network-based constitutive model

The PANN constitutive model is now calibrated to data generated with an analytical potential. Details on the data generation are presented in Sect. 7.1.1, while details on the model calibration are presented in Sect. 7.1.2.

#### 7.1.1 Data generation

For the following investigations, the PANN model introduced in Sect. 3 is calibrated to data generated with the analytical Mooney-Rivlin model, cf. Eq. (21). For this, two sets of material parameters are considered, which lead to a compressible and a nearly incompressible model, respectively, cf. Table 1. If not stated otherwise, the parameter set for a compressible material behavior is applied. In the following, the Mooney-Rivlin model is referred to as the ground truth (GT) model, to which the results obtained with the calibrated PANN model are compared.

Following [23], the calibration dataset consists of (i) a uniaxial tension stress state, (ii) an equibiaxial tension stress state, and (iii) a simple shear deformation state, which are applied on the analytical Mooney-Rivlin model. The uniaxial tension is applied in \(x\)-direction with \(F_{11}\in[0.75,\,1.75]\), the equibiaxial tension is applied in \(x-y\)-direction with \(F_{11}=F_{22}\in[0.75,\,1.755]\), and
simple shear is applied with \(F_{12}\in[-0.25,\,0.75]\). In addition, for the test dataset, a mixed shear-tension deformation is applied ("test 3" in [23]), which represents a fairly general deformation mode, cf. Fig. 2. Each load case consists of 100 datapoints, which makes 300 strain-stress tuples in the calibration dataset. It should be mentioned that besides these experimentally motivated load cases, it is also possible to sample the deformation states either by analytical considerations on physically admissible deformation states [62] or by investigating the deformation states which occur in the engineering component under investigation [52].

#### 7.1.2 Model calibration

For the PANN model calibration, six different network architectures with \(n\in\{4,\,8,\,16,\,32,\,64,\,128\}\) nodes in the hidden layer are considered for the compressible material behavior, which corresponds to \(|\mathbf{\mathcal{P}}|=24,\,48,\,96,\,192,\,384,\,768\) NN model parameters. In this way, the computational performance of the algorithms can be examined, and furthermore, it is demonstrated that the choice of nodes \(n\) in the hidden layer is generic for the approach introduced in this work. In general, this also holds for the number of hidden layers, of which the proposed framework is also independendent, with the caveat that the explicit derivatives of the NN, cf. Sect. 3.3, must be adapted when more than one hidden layer is considered. For the nearly incompressible material behavior, \(n=8\) nodes are considered. At this point, it should be noted that even the smallest NN architecture with 4 nodes is flexible enough to perfectly represent the Mooney-Rivlin model, and, furthermore, calibrating a PANN model to such a simple analytical potential should be considered as an academic example. In the following examples, PANN\(n\) denotes a PANN model with \(n\) nodes in the hidden layer.

For the calibration, the models are implemented in TensorFlow 2.5.0 using Python 3.9.9. The optimization is carried out using the ADAM algorithm with default settings. Therefore, the full batch of calibration data is used for 5000 calibration epochs. Each model architecture is calibrated one time. For the reproducibility of the results, the calibrated model parameters are provided in the GitHub repository [https://github.com/CPShub/sim-data](https://github.com/CPShub/sim-data).

All model architectures perform excellent, both on the calibration and the test dataset, and both for the compressible and the nearly incompressible material behavior, cf. Table 2, where the losses of the calibrated models are provided. While there is a general trend for a decreasing loss for larger model architectures, it should be noted that even the small architectures are able to nearly perfectly represent the GT models, which becomes evident in Fig. 2, where the test case is evaluated for the PANN8 models.

\begin{table}
\begin{tabular}{l|l l l} compressible case & \(a\) & 831.25 & Pa \\  & \(b\) & 166.25 & Pa \\  & \(c\) & 10000 & Pa \\  & \(d\) & 2327.5 & Pa \\ mass density (static/transient examples) & \(\rho_{0}\) & 0/100 & kg m\({}^{-3}\) \\ \hline \hline nearly incompressible case & \(a\) & 126 & Pa \\  & \(b\) & 252 & Pa \\  & \(c\) & 81512 & Pa \\  & \(d\) & 1260 & Pa \\ mass density & \(\rho_{0}\) & 0 & kg m\({}^{-3}\) \\ \end{tabular}
\end{table}
Table 1: Mooney-Rivlin material parameters employed for the numerical examples.
### Static deformation of Cook's membrane

Here, a static example is investigated. In Sect. 7.2.1, compressible material behavior is assumed, and the computation time of the PANN constitutive model compared to the GT model is examined. In Sect. 7.2.2, nearly incompressible material behavior is assumed, and the performance of the mixed formulations is demonstrated.

#### 7.2.1 Compressible case

As a first example, a three-dimensional version of the well-known Cook's membrane with compressible material behavior is considered, see e.g. [9]. With this example, the spatial performance of the PANN constitutive model compared to the GT model in a purely static, large deformation solid mechanical scenario is examined. Furthermore, a detailed investigation of the computation times for both PANN and GT models is carried out.

The geometry of the membrane is provided by Fig. 3 (left). The membrane is subject to a clamping on its left side while a pressure load is applied to shear the membrane on the right side, as illustrated in Fig. 3 (right). The resultant force of the pressure load has the magnitude of \(p=200\) Pa and shows into the positive \(\mathbf{e}_{3}\)-direction. For the discretization, displacement-based trilinear hexahedron elements based on to Eq. (42) are employed.

In Fig. 4, typical deformed configurations for both the GT and the PANN8 model are provided. The mesh convergence study in Fig. 5 (left) shows a practically identical behavior for the GT model and the PANN model for all choices of nodes. Unsurprisingly, the PANN model is able to perfectly

\begin{table}
\begin{tabular}{c|c c|c c}  & \multicolumn{2}{c|}{compressible model} & \multicolumn{2}{c}{nearly incompressible model} \\  & calibration & test & calibration & test \\ \hline PANN & 4 & 1.70 & 0.89 & \\  & 8 & 0.84 & 0.61 & 1.05 & 0.23 \\  & 16 & 0.93 & 0.61 & \\  & 32 & 0.25 & 0.15 & \\  & 64 & -0.23 & 0.17 & \\  & 128 & -0.11 & 0.10 & \\ \end{tabular}
\end{table}
Table 2: \(\log_{10}(\mathcal{MSE})\) according to Eq. (29) for the calibrated PANN models.

Figure 2: Mixed shear-tension test evaluated for the GT models (marks) and the calibrated PANN8 models (continous lines) for both the compressible (left) and the nearly incompressible (right) case. Stress in MPa.
represent the behavior of the GT model, even for small numbers of nodes. This is also reflected in the Newton iterations required to solve the system of nonlinear equations, which are visualized in Fig. 5 (right) for a typical mesh and load step. For all choices of nodes in the PANN model, the required Newton iterations and the corresponding norm of the residual are equal to the GT model. This shows that the PANN model is not only able to excellently predict the stress behavior of the GT model, which is closely related to the first derivative of the hyperelastic potential, but also excellently predicts the second derivative, which is essential for the solution of the nonlinear FEA, although the second derivative was not included in the calibration process, cf. Sect. 3.2.

In Fig. 6 and Fig. 7, the percentage computation times for the simulation of Cook's membrane with 14406 degree of freedom (DOFs) on one single core are presented for the GT model and the PANN128 model, respectively. On the element level of the GT model, the evaluation of the derivatives of the analytical Mooney-Rivlin potential ("GT derivatives") only makes up 5% of the computation time and takes even less time than operations such as the evaluation of the tensor cross product. For the PANN128 model, however, the evaluation of the NN potential is significantly more expensive and takes over 50% of the computation time on element level. Thus, the overall simulation time is increased. In Fig. 8, the computation time for the GT model and the PANN model with different numbers of nodes on one single core are visualized. Depending on the number

Figure 4: Deformed configuration with von Mises stress distribution \(\sigma_{vM}\) [Pa] with a typical mesh and load step for compressible GT model (left) and the corresponding PANN8 model (mid).

Figure 3: 2D projection of the geometry (left) and 3D illustration of the boundary conditions (right) of the Cook’s membrane.
of nodes in the NN, the computation time of the PANN model is more than twice as long compared to the GT model, when the simulation is carried out on a single core. However, when the FEA is parallelised on 24 cores, the computation time of the PANN model is equal to the computation time of the analytical GT model, and, within the examined bounds, even independent of the number of nodes in the hidden layer. The parallelization is employed within the research code moofeKIT [30] written in Matlab using the Parallel Computing Toolbox. In particular, the element routine is parallelized within a parallel for loop environment parfor. Obviously, the computation time is highly influenced by many factors, such as computer hardware and programming language. Still, our findings suggest that with a suitable implementation of the NN-based constitutive model, computation times comparable to analytical constitutive models can be achieved.

#### 7.2.2 Nearly incompressible case

In the next step, Cook's membrane is simulated with a nearly incompressible material behavior of both the GT and the PANN model. Thereby, the performance of the PANN constitutive model

Figure 5: Convergence plot with vertical displacement of point A (see Fig. 3) vs. mesh refinement (left) and Newton iterations for a typical mesh and load step (right).

Figure 6: Percentage computation time of the GT model on global level (left) with 1.6% assembly of \(\mathbf{K}\) included in the slice “other”, and percentage computation time on element level (right), where the operator matrix \(\mathcal{D}(\mathbf{A})\) is employed (see [9, Appx. B]). Evaluated for the simulation of Cook’s membrane with 14406 DOFs on one core.
when applied in mixed formulations, which avoid the locking effect, is investigated.

The resultant force of the pressure load has the magnitude of \(p=100\) Pa and shows into the positive \(\mathbf{e}_{3}\)-direction, see Fig. 3. Typical deformed configurations with von Mises stress are provided in Fig. 9. A convergence study is provided in Fig. 10 (left), where different material models and element types are examined. The solutions for both the GT and the PANN8 model are almost indistinguishable for all elements employed. As to be expected in this nearly incompressible, shear-loaded problem, the locking effect becomes obvious for both displacement-based models with trilinear hexahedron elements (GT - H1 and PANN8 - H1). In contrast to that, the 5-field mixed formulation presented in Eq. (44) with trilinear / constant hexahedron elements (G-J) and the new Hu-Washizu like 7-field mixed formulation in Eq. (47) with trilinear / constant hexahedron elements (I-II-J) exhibit fast, locking-free convergence. In particular, although the I-II-J element only considers the _scalar_ quantities \(I_{\mathbf{C}}\), \(II_{\mathbf{C}}\), and \(J\) as additional kinematic unknowns, it shows the same excellent spatial performance as the G-J element, which, in contrast, includes an additional _tensor-valued_ unknown.

Figure 8: Computation times with 1 and 24 core setup for GT and different PANN models, evaluated for the simulation of Cook’s membranan with 36015 DOFs.

Figure 7: Percentage computation time of the PANN128 model on global level (left) with 0.5% assembly of \(\mathbf{K}\) included in the slice “other”, and percentage computation time on element level (right). Evaluated for the simulation of Cook’s membrane with 14406 DOFs on one core.
### Dynamics of L-shaped body

The objective of this investigation is to demonstrate the discrete conservation properties of the proposed EM integrator, particularly applied in the context of the PANN model. In addition to that, the robustness and long-term stability of the midpoint (MP) and the EM integrator are compared for the PANN model.

The initial boundary value problem of an L-shaped body is considered, which is inspired by [9, 99] (in case of pure mechanics) and [26, 27] (in case of thermo-mechanics and electro-thermo-mechanics, respectively). The geometry of the L-shaped body and the initial setting with boundary conditions are depicted in Fig. 11. Therein the applied dead load pressures \(\mathbf{P}_{1}:\mathcal{T}\rightarrow\mathds{R}^{3}\) and \(\mathbf{P}_{2}:\mathcal{T}\rightarrow\mathds{R}^{3}\)

Figure 10: Convergence plot with displacement vs. mesh refinement for the nearly incompressible Cook’s membrane. Different material models (GT and PANN8) and elements (H1: displacement-based trilinear hexahedron, G-J – H1H0: 5-field mixed formulation according to Eq. (44) with trilinear / constant hexahedron elements, and I-II-J – H1H0: new Hu-Washizu like 7-field mixed formulation according to Eq. (47) with trilinear / constant hexahedron elements) are employed. Note that the two alternative constitutive model formulations under investigation essentially lead to the same results for each finite element formulation.

Figure 9: Deformed configuration with von Mises stress distribution \(\sigma_{vM}\) [Pa] with a typical mesh and load step for the nearly incompressible GT model (left) and the corresponding PANN8 model (mid).
follow the functions

\[\mathbf{P}_{1}(t)=f(t)\begin{pmatrix}256/9\\ 512/9\\ 768/9\end{pmatrix}\frac{\mathrm{N}}{\mathrm{m}^{2}}=-\mathbf{P}_{2}(t),\quad \text{with}\quad f(t)=\begin{cases}t&\text{for}\quad t\leq 2.5\ \mathrm{s}\\ 5-t&\text{for}\quad 2.5\ \mathrm{s}\leq t\leq 5\ \mathrm{s}\\ 0&\text{for}\quad t>5\ \mathrm{s}\end{cases}. \tag{69}\]

In Fig. 12 and Fig. 13 snapshots of the motion of the L-shaped body with von Mises stress distribution for the compressible GT and PANN8 models, respectively, are shown. As can be observed, both stress plots show a remarkable similarity. The total energy and the energy difference of both the GT model and corresponding PANN8 model are shown in Fig. 14 and Fig. 15, respectively. Apparently, the midpoint rule aborts due to energy blow ups for both GT and PANN8 model, even when the time step size is reduced. In contrast, the EMS is stable for both the GT and PANN8 model during this long-term simulation and preserves the total energy after the loading phase (\(t>5\) s). Furthermore, it is important to note that the energy results obtained for the GT and the PANN8 model are practically identical. Furthermore, in Fig. 16 the components of the angular momentum of the L-shaped body for both the GT and the PANN8 model in case of the EMS are shown. As evident from the observation, the components of the angular momentum are preserved after the loading phase (\(t>5\) s). Similar to the energy results, the components of the angular momentum exhibit practically identical results for both models.

## 8 Conclusion

In this contribution, advanced spatial and temporal discretization techniques are tailored to hyperelastic PANN constitutive models, i.e., NN-based constitutive models which fulfill all relevant mechanical conditions of hyperelasticity by construction [69]. One of the benefits of PANN-based constitutive models is their extraordinary flexibility, which is achieved by using NNs as highly flexible functions. Actually, the feed-forward neural networks applied for this have a quite simple structure, which is based on only a few ingredients. Furthermore, from a formal point of view, NNs are nothing more but mathematical functions, and thus, constitutive models based on NNs can be applied in numerical methods just like their analytical counterparts. Still, the explicit mathematical formulation of the NN, and even more its derivatives often required for numerical applications, are quite extensive

Figure 11: 2D projection of the geometry (left) and 3D illustration of the boundary conditions (right) of the L-shaped body with \(t=3m\).
Figure 12: Snapshots with von Mises stress distribution \(\sigma_{vM}\) [Pa] of the GT model at times \(t\in\{4.0,8.0,12.0,16.0,20.0,24.0\}s\).

Figure 13: Snapshots with von Mises stress distribution \(\sigma_{vM}\) [Pa] of the PANN8 model at times \(t\in\{4.0,8.0,12.0,16.0,20.0,24.0\}s\).
Figure 16: Components of angular momentum of the L-shaped body for the GT model (marks) and for the PANN8 model (continuous lines), both in case of the energy momentum scheme.

Figure 14: Total energy of the L-shaped body (left) for the energy-momentum scheme (EM) and the midpoint rule (MP), and energy difference (right) for the EM scheme. Evaluated for the GT model.

Figure 15: Total energy of the L-shaped body (left) for the energy-momentum scheme (EM) and the midpoint rule (MP), and energy difference (right) for the EM scheme. Evaluated for the PANN8 model.
compared to analytical constitutive models. This suggests to tailor advanced discretization methods for NN-based constitutive models, to arrive at compact mathematical formulations and convenient implementations with superior stability and robustness in static and dynamic analysis.

In the present work, this is achieved by setting the focus of the numerical methods on the scalar-valued strain invariants that hyperelastic stored energy functions depend on, rather than considering the full tensorial strain measures. Thereby, the evaluation of the hyperelastic potential and its derivatives w.r.t. the strain invariants is separated from the derivatives of the strain invariants w.r.t. the primal tensorial strain measures. This is a very natural choice, since for NN-based constitutive models, the derivatives of the hyperelastic potential w.r.t. the invariants become more complex, while the derivatives of the invariants remain the same as in the analytical theory. In this regard, an elegant 7-field Hu-Washizu like mixed formulation is proposed, where the additional unknowns are only strain invariants. This novel mixed formulation is applicable for nearly incompressible material behavior, in particular, to avoid locking phenomena. Furthermore, a novel energy-momentum scheme (EMS) is designed, which is applicable in dynamical simulations to preserve energy and momentum in the discretized version of the weak form. At the same time, the EMS highly improves the numerical stability in transient simulations. For this, a new algorithmic stress formula is proposed, which, again, only depends on strain invariants. Here, the required discrete gradients are formulated in scalar-valued invariants, therefore only consisting of Greenspan formulas. The employed Greenspan formulas are by far less complex for both consistent tangent computation and implementation when compared to the usually employed projection-based discrete gradients.

Both the mixed method and the EMS are applied in finite element analysis. For this, the hyperelastic PANN model proposed by [69] is calibrated to data generated with an analytical Mooney-Rivlin potential. Unsurprisingly, due to the simple nature of the analytical potential, the results obtained with the PANN model are practically identical to the results obtained with the Mooney-Rivlin model in all simulations. The newly proposed mixed 7-field Hu-Washizu like formulation showed excellent performance in a scenario where purely displacement-based formulations extensively suffer from locking phenomena. Also, it showed similar results when compared to a state-of-the-art mixed, polyconvex 5-field approach [59] which considers tensor-valued quantities as additional unknowns. Furthermore, in a benchmark example, it could be shown that with a suitable parallelization strategy, the computation time of the NN-based model is practically the same as for the analytical constitutive model. In a dynamical simulation, the tailored EMS showed a remarkable stability and robustness even for coarse time steps when compared to the symplectic midpoint-rule. Overall, the numerical examples demonstrate the superior stability and robustness of the proposed numerical methods.

Finally, future work should be concerned with the extension of the proposed methods towards multiphysics [56, 112], inelasticity [1, 77, 93], and structural mechanics. Furthermore, while the present work only considers isotropic material behavior, the extension to further symmetry groups is straightforward [22, 69].

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

**Software.** The source code used for the finite element computations is implemented in Matlab under MIT license and is available at [https://github.com/kit-ifm/moofeKIT](https://github.com/kit-ifm/moofeKIT). Version 1.XX of the code, used in this paper, is archived at [30].

**Acknowledgement.** M. Franke acknowledges the financial support provided by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation, project number 443238377). D.K. Klein and O. Weeger acknowledge the financial support provided by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation, project number 492770117), the Graduate School of Computational Engineering at TU Darmstadt, and Hessian.AI.

## Appendix A Continous balance laws

The verification of the time-continuous balance laws, i.e., balance of angular momentum and balance of total energy for the weak form given in Eq. (51), will be dealt with in the following.

### Balance of angular momentum

In order to verify the balance of angular momentum, we introduce the angular momentum

\[\mathbf{J}_{t}=\int_{\mathcal{B}_{0}}\boldsymbol{\varphi}_{t}\times\rho_{0} \,\mathbf{v}_{t}\,\,\mathrm{d}V\,,\] (A.1)

and the external momentum

\[\mathbf{M}_{t}^{\mathrm{ext}}=\int_{\mathcal{B}_{0}}\boldsymbol{\varphi}_{t} \times\mathbf{B}_{t}\,\,\mathrm{d}V+\int_{\partial\mathcal{B}_{0}^{T}} \boldsymbol{\varphi}_{t}\times\mathbf{T}_{t}\,\,\mathrm{d}A\,.\] (A.2)

Furthermore, we choose the admissible variations

\[\delta\boldsymbol{\varphi}=\boldsymbol{\xi}\times\boldsymbol{\varphi}_{t}\,, \quad\delta\mathbf{v}=\boldsymbol{\xi}\times\dot{\boldsymbol{\varphi}}_{t}\,,\] (A.3)

where \(\mathds{R}^{3}\ni\boldsymbol{\xi}=\mathrm{const}\). For further considerations, we anticipate

\[\mathbf{F}(\delta\boldsymbol{\varphi}=\boldsymbol{\xi}\times\boldsymbol{ \varphi}_{t})=\boldsymbol{\hat{\xi}}\,\mathbf{F}_{t}\,,\] (A.4)

where \(\boldsymbol{\hat{\xi}}\) is skew-symmetric, such that

\[\boldsymbol{\hat{\xi}}\,\mathbf{a}=\boldsymbol{\xi}\times\boldsymbol{\varphi }_{t}\quad\forall\mathbf{a}\in\mathds{R}^{3}\,.\] (A.5)

Recall Eq. (51)\({}_{1}\) and substitute Eq. (A.3)\({}_{2}\), which finally yields

\[\boldsymbol{\xi}\cdot\int_{\mathcal{B}_{0}}(\dot{\boldsymbol{\varphi}}_{t} \times\mathbf{v}_{t})\,\rho_{0}\,\,\mathrm{d}V=0\,,\] (A.6)

where therein the identity \(\dot{\boldsymbol{\varphi}}\times\dot{\boldsymbol{\varphi}}=\mathbf{0}\) has been employed. Consider the integrand of the 2nd term in Eq. (51)\({}_{2}\) and substitute Eq. (A.3)\({}_{1}\)

\[\mathbf{S}_{t}:\tfrac{1}{2}\,\delta\mathbf{C}_{t}=\mathbf{S}_{t}:\frac{1}{2} \left(\mathbf{F}_{t}^{\mathrm{T}}\left(\boldsymbol{\hat{\xi}}^{\mathrm{T}}+ \boldsymbol{\hat{\xi}}\right)\mathbf{F}_{t}\right)=0\,,\] (A.7)

which is due to skew-symmetry of \(\boldsymbol{\hat{\xi}}\), i.e. \(\boldsymbol{\hat{\xi}}^{\mathrm{T}}=-\boldsymbol{\hat{\xi}}\). Recall Eq. (51)\({}_{2}\) and substitute Eq. (A.3)\({}_{1}\)

\[\boldsymbol{\xi}\cdot\int_{\mathcal{B}_{0}}\boldsymbol{\varphi}_{t}\times \dot{\mathbf{v}}_{t}\,\rho_{0}\,\,\mathrm{d}V+\Pi^{\mathrm{ext}}(\boldsymbol {\xi}\times\boldsymbol{\varphi}_{t})=0\,,\] (A.8)

for the above we finally obtain with Eq. (A.6) the desired balance of angular momentum

\[\boldsymbol{\xi}\cdot\left(\frac{\,\mathrm{d}}{\,\mathrm{d}t}\, \int_{\mathcal{B}_{0}}\boldsymbol{\varphi}_{t}\times\mathbf{v}_{t}\,\rho_{0} \,\,\mathrm{d}V-\int_{\mathcal{B}_{0}}\boldsymbol{\varphi}_{t}\times\mathbf{ B}_{t}\,\,\mathrm{d}V-\int_{\partial\mathcal{B}_{0}^{T}}\boldsymbol{\varphi}_{t} \times\mathbf{T}_{t}\,\,\mathrm{d}V\right) =0\] \[\Leftrightarrow\quad\boldsymbol{\xi}\cdot\left(\frac{\,\mathrm{d} }{\,\mathrm{d}t}\,\mathbf{J}_{t}-\mathbf{M}_{t}^{\mathrm{ext}}\right) =0\,.\]
### Balance of total energy

With the definition of the kinetic energy

\[T_{t}=\frac{1}{2}\,\int_{\mathcal{B}_{0}}\rho_{0}\,\mathbf{v}_{t}\cdot\mathbf{v}_ {t}\ \mathrm{d}V\,,\] (A.9)

at hand, we choose the admissible variations

\[\delta\boldsymbol{\varphi}=\dot{\boldsymbol{\varphi}}_{t},\qquad\delta\mathbf{v }=\dot{\mathbf{v}}_{t}\,,\] (A.10)

and substitute Eq. (A.10)\({}_{2}\) into Eq. (51)\({}_{1}\), such that

\[\int_{\mathcal{B}_{0}}\dot{\boldsymbol{\varphi}}_{t}\cdot\dot{\mathbf{v}}_{t} \,\rho_{0}\ \mathrm{d}V=\int_{\mathcal{B}_{0}}\mathbf{v}_{t}\cdot\dot{\mathbf{v}}_{t}\, \rho_{0}\ \mathrm{d}V=\frac{\mathrm{d}}{\mathrm{d}t}\,\left(\frac{1}{2}\,\int_{ \mathcal{B}_{0}}\mathbf{v}\cdot\mathbf{v}_{t}\,\rho_{0}\ \mathrm{d}V\right)=\dot{T}_{t}\,.\] (A.11)

Consider the derivative w.r.t. time of the internal potential

\[\frac{\mathrm{d}}{\mathrm{d}t}\,w^{1}\Pi_{t}^{\mathrm{int}} =\frac{\mathrm{d}}{\mathrm{d}t}\,\int_{\mathcal{B}_{0}}W^{1}(I_{ \mathbf{C}_{t}},II_{\mathbf{C}_{t}},III_{\mathbf{C}_{t}})\ \mathrm{d}V\] \[=\int_{\mathcal{B}_{0}}\partial_{I\mathbf{C}}W^{1}\mathbf{I}: \dot{\mathbf{C}}_{t}+\partial_{II_{\mathbf{C}}}W^{1}\mathbf{I}:\dot{\mathbf{G} }_{t}+\partial_{J}W^{1}\dot{J}_{t})\ \mathrm{d}V\] \[=\int_{\mathcal{B}_{0}}2\,\left(\partial_{I_{\mathbf{C}}}W^{1} \mathbf{I}+\partial_{II_{\mathbf{C}}}W^{1}\mathbf{I}\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!
Considering the integrand of the first term of the right-hand side of Eq. (54)\({}_{2}\) and plugging in Eq. (B.1)\({}_{1}\) yields

\[\mathbf{S}_{\text{algo}}:\frac{1}{2}\,\mathbf{F}_{n+\frac{1}{2}}{}^{\text{T}}\,( \hat{\boldsymbol{\xi}}+\hat{\boldsymbol{\xi}}{}^{\text{T}})\,\mathbf{F}_{n+ \frac{1}{2}}{}^{\text{T}}=0\,.\] (B.4)

With the above in mind, plugging Eq. (B.1)\({}_{1}\) into the whole Eq. (54)\({}_{2}\) and adding Eq. (B.3) yields after some algebra the desired semi-discrete balance of angular momentum

\[\boldsymbol{\xi}\cdot\left(\frac{1}{\Delta t}\,\int_{\mathcal{B}_{ 0}}(\Delta\boldsymbol{\varphi}\times\mathbf{v}_{n+\frac{1}{2}}+\boldsymbol{ \varphi}_{n+\frac{1}{2}}\times\Delta\mathbf{v})\,\rho_{0}\,\,\mathrm{d}V \right)=\boldsymbol{\xi}\cdot\left(\int_{\mathcal{B}_{0}}\boldsymbol{\varphi} _{n+\frac{1}{2}}\times\mathbf{B}_{n+\frac{1}{2}}\,\,\mathrm{d}V\] \[\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\int_{\partial \mathcal{B}_{0}^{\text{T}}}\boldsymbol{\varphi}_{n+\frac{1}{2}}\times \mathbf{T}_{n+\frac{1}{2}}\,\,\mathrm{d}A\right)\] \[\Leftrightarrow\qquad\boldsymbol{\xi}\cdot\frac{1}{\Delta t}\, \int_{\mathcal{B}_{0}}(\rho_{0}\,\boldsymbol{\varphi}_{n+1}\times\mathbf{v}_ {n+1}-\rho_{0}\,\boldsymbol{\varphi}_{n}\times\mathbf{v}_{n})\,\,\mathrm{d}V =\boldsymbol{\xi}\cdot\mathbf{M}_{n,n+1}^{\text{ext}}\] \[\Leftrightarrow\qquad\boldsymbol{\xi}\cdot\left(\frac{1}{\Delta t }\,\Delta\mathbf{J}-\mathbf{M}_{n,n+1}^{\text{ext}}\right)=0\,.\] (B.5)

Accordingly, the proposed EM scheme consistently reproduces the balance of angular momentum for the time-discrete system which is e.g. in contrast to the 2nd order accurate trapezoidal rule.

### Balance of total energy

To verify the balance of energy in the semi-discrete weak form Eq. (54), we choose admissible variations

\[\delta\boldsymbol{\varphi}=\boldsymbol{\varphi}_{n+1}-\boldsymbol{\varphi}_{n }=\Delta\boldsymbol{\varphi},\quad\delta\mathbf{v}=\mathbf{v}_{n+1}-\mathbf{v} _{n}=\Delta\mathbf{v}\,,\] (B.6)

plug Eq. (B.6)\({}_{2}\) into Eq. (54)\({}_{1}\):

\[\int_{\mathcal{B}_{0}}\Delta\boldsymbol{\varphi}\,\frac{1}{\Delta t }\,(\boldsymbol{\varphi}_{n+1}-\boldsymbol{\varphi}_{n})\,\rho_{0}\,\,\mathrm{d}V =\int_{\mathcal{B}_{0}}\Delta\mathbf{v}\cdot\frac{1}{2}\,(\mathbf{ v}_{n+1}-\mathbf{v}_{n})\,\rho_{0}\,\,\mathrm{d}V\] \[=\frac{1}{2}\,\int_{\mathcal{B}_{0}}\rho_{0}\,(\mathbf{v}_{n+1} \cdot\mathbf{v}_{n+1}-\mathbf{v}_{n}\cdot\mathbf{v}_{n})\,\,\mathrm{d}V\] \[=T_{n+1}-T_{n}=\Delta T\,.\] (B.7)

and plug Eq. (B.6)\({}_{1}\) into the integrand of the first term of the right-hand side of Eq. (54)\({}_{2}\):

\[\mathbf{S}_{\text{algo}}:\text{sym}(\Delta\mathbf{F}^{\text{T}}\,\mathbf{F}_ {n+\frac{1}{2}})=\mathbf{S}_{\text{algo}}:(\Delta\mathbf{F}^{\text{T}}\, \mathbf{F}_{n+\frac{1}{2}})=\mathbf{S}_{\text{algo}}:\frac{1}{2}(\mathbf{C}_{n +1}-\mathbf{C}_{n})=\mathbf{S}_{\text{algo}}:\frac{1}{2}\Delta\mathbf{C}\,.\] (B.8)

Here, the symmetry property \(\mathbb{S}^{3}\cap\mathds{K}^{3}=\{0\}\) has been taken into account, where \(\mathbb{S}^{3}=\{\mathbf{A}\in\mathbb{T}^{3}\,|\,\mathbf{A}=\mathbf{A}^{\text{ T}}\}\) and \(\mathds{K}^{3}=\{\mathbf{A}\in\mathbb{T}^{3}\,|\,\mathbf{A}=-\mathbf{A}^{\text{ T}}\}\) denote spaces of symmetric and skew-symmetric, quadratic second order tensors, respectively. Furthermore, for Eq. (B.8) we insert Eq. (56), which yields after some algebra, cf. Eq. (C.17):

\[\mathbf{S}_{\text{algo}}:\frac{1}{2}\Delta\mathbf{C} =\Big{(}\mathrm{D}_{I\mathbf{C}}W^{1}\mathbf{I}+\mathrm{D}_{II \mathbf{C}}W^{1}\mathbf{I}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!
Plugging Eq. (B.6)\({}_{1}\) into the whole Eq. (54)\({}_{2}\), considering the above results Eq. (B.7) and Eq. (B.9), and employing the directionality property Eq. (C.17) yields the desired semi-discrete balance of total energy

\[\Delta T=-\int_{\mathcal{B}_{0}}W^{\mathrm{I}}(I_{\mathbf{C}_{n+1} },II_{\mathbf{C}_{n+1}},J_{\mathbf{C}_{n+1}})-W^{\mathrm{I}}(I_{\mathbf{C}_{n} },II_{\mathbf{C}_{n}},J_{\mathbf{C}_{n}})\ \mathrm{d}V+W^{\mathrm{ext}}_{n,n+1}\] (B.10) \[\Leftrightarrow \Delta T+\Delta^{W^{\mathrm{I}}}\Pi^{\mathrm{int}}=W^{\mathrm{ ext}}_{n,n+1}\,,\]

where

\[W^{\mathrm{ext}}_{n,n+1} =-\left(\Pi^{\mathrm{ext}}_{n+\frac{1}{2}}(\mathbf{\varphi}_{n+1})- \Pi^{\mathrm{ext}}_{n+\frac{1}{2}}(\mathbf{\varphi}_{n})\right)\] (B.11) \[=-\Pi^{\mathrm{ext}}_{n+\frac{1}{2}}(\Delta\mathbf{\varphi})=\int_{ \mathcal{B}_{0}}\Delta\mathbf{\varphi}\cdot\mathbf{B}_{n+\frac{1}{2}}\ \mathrm{d}V+\int_{\partial\mathcal{B}_{0}^{T}}\Delta\mathbf{\varphi}\cdot \mathbf{T}_{n+\frac{1}{2}}\ \mathrm{d}A\,.\]

As shown above the proposed EM scheme consistently reproduces the balance of total energy for the time-discrete system. This achievement can be attributed to the algorithmic stress evaluation, as demonstrated in Eq. (B.9). In contrast, standard integrators such as the 2nd order accurate midpoint rule fail to maintain this discrete energy balance.

## Appendix C Partioned discrete gradients of \(W^{\mathrm{I}}\)

Introducing the tuple \(\mathbf{\Pi}=(\Pi^{1},\Pi^{2},\Pi^{3})=(I_{\mathbf{C}},II_{\mathbf{C}},J)\), the partitioned discrete gradients \(\mathrm{D}_{I_{\mathbf{C}}}W^{\mathrm{I}}\), \(\mathrm{D}_{II_{\mathbf{C}}}W^{\mathrm{I}}\), and \(\mathrm{D}_{J}W^{\mathrm{I}}\) with final presentation in Eq. (58) are computed with

\[D_{\Pi^{i}}W^{\mathrm{I}}=\frac{1}{2}\left(\mathrm{D}_{\Pi^{i}}W ^{\mathrm{I}}|_{\Pi^{J}_{n+1},\Pi^{K}_{n}}+\mathrm{D}_{\Pi^{i}}W^{\mathrm{I}} |_{\Pi^{J}_{n},\Pi^{K}_{n+1}}\right),\] (C.1) \[\forall i\in Y,\quad j\in J<Y>K\ni k\,:\,J\ni j<i<k\in K\,,\]

where in general for tensor-valued quantities \(\Pi^{i}\) we have

\[\mathrm{D}_{\Pi^{i}}W^{\mathrm{I}}|_{\Pi^{J}_{n+1},\Pi^{K}_{n}}= \partial_{\Pi^{i}}W^{\mathrm{I}}(\Pi^{i}_{n+\frac{1}{2}})|_{\Pi^{ J}_{n+1},\Pi^{K}_{n}}\] (C.2) \[+\frac{W^{\mathrm{I}}(\Pi^{i}_{n+1})|_{\Pi^{J}_{n+1},\Pi^{K}_{n}} -W^{\mathrm{I}}(\Pi^{i}_{n})|_{\Pi^{J}_{n+1},\Pi^{K}_{n}}-<\partial_{\Pi^{i}}W ^{\mathrm{I}}(\Pi^{i}_{n+\frac{1}{2}})|_{\Pi^{J}_{n+1},\Pi^{K}_{n}},\Delta\Pi^ {i}>}{<\Delta\Pi^{i},\Delta\Pi^{i}>}\Delta\Pi^{i}\,,\]

and

\[\mathrm{D}_{\Pi^{i}}W^{\mathrm{I}}|_{\Pi^{J}_{n},\Pi^{K}_{n+1}}= \partial_{\Pi^{i}}W^{\mathrm{I}}(\Pi^{i}_{n+\frac{1}{2}})|_{\Pi^ {J}_{n},\Pi^{K}_{n+1}}\] (C.3) \[+\frac{W^{\mathrm{I}}(\Pi^{i}_{n+1})|_{\Pi^{J}_{n},\Pi^{K}_{n+1}} -W^{\mathrm{I}}(\Pi^{i}_{n})|_{\Pi^{J}_{n},\Pi^{K}_{n+1}}-<\partial_{\Pi^{i}}W ^{\mathrm{I}}(\Pi^{i}_{n+\frac{1}{2}})|_{\Pi^{J}_{n},\Pi^{K}_{n+1}},\Delta\Pi^ {i}>}{<\Delta\Pi^{i},\Delta\Pi^{i}>}\Delta\Pi^{i}\,,\]

whereas for scalar-valued invariants \(\Pi^{i}\) the above simplifies to the so-called Greenspan formulas (cf. [35])

\[\mathrm{D}_{\Pi^{i}}W^{\mathrm{I}}|_{\Pi^{J}_{n+1},\Pi^{K}_{n}}=\frac{W^{ \mathrm{I}}(\Pi^{i}_{n+1})|_{\Pi^{J}_{n+1},\Pi^{K}_{n}}-W^{\mathrm{I}}(\Pi^{i}_ {n})|_{\Pi^{J}_{n+1},\Pi^{K}_{n}}}{\Delta\Pi^{i}}\,,\] (C.4)

and

\[\mathrm{D}_{\Pi^{i}}W^{\mathrm{I}}|_{\Pi^{J}_{n},\Pi^{K}_{n+1}}=\frac{W^{ \mathrm{I}}(\Pi^{i}_{n+1})|_{\Pi^{J}_{n},\Pi^{K}_{n+1}}-W^{\mathrm{I}}(\Pi^{i} _{n})|_{\Pi^{J}_{n},\Pi^{K}_{n+1}}}{\Delta\Pi^{i}}\,.\] (C.5)

This finally leads to
1. the partioned discrete gradient \(\mathrm{D}_{I_{\mathbf{C}}}W^{\mathrm{I}}\) \[D_{I_{\mathbf{C}}}W^{\mathrm{I}}=\frac{1}{2}\left(\mathrm{D}_{I_{\mathbf{C}}}W^{ \mathrm{I}}|_{\Pi_{n+1}^{J},\Pi_{n}^{K}}+\mathrm{D}_{I_{\mathbf{C}}}W^{\mathrm{I }}|_{\Pi_{n}^{J},\Pi_{n+1}^{K}}\right)\] (C.6) with being \(i=1\), \(J=\emptyset\), \(K=\{2,3\}\) we have \(\Pi^{1}=I_{\mathbf{C}}\), \(\Pi^{J}=\emptyset\), \(\Pi^{K}=(II_{\mathbf{C}},III_{\mathbf{C}})\) and thus \[\mathrm{D}_{I_{\mathbf{C}}}W^{\mathrm{I}}|_{\Pi_{n+1}^{J},\Pi_{n}^{K}} =\frac{W^{\mathrm{I}}(I_{\mathbf{C}_{n+1}},II_{\mathbf{C}_{n}},J_ {n})-W^{\mathrm{I}}(I_{\mathbf{C}_{n}},II_{\mathbf{C}_{n}},J_{n})}{\Delta I_{ \mathbf{C}}}\] (C.7) \[\mathrm{D}_{I_{\mathbf{C}}}W^{\mathrm{I}}|_{\Pi_{n}^{J},\Pi_{n+1}^{K}} =\frac{W^{\mathrm{I}}(I_{\mathbf{C}_{n+1}},II_{\mathbf{C}_{n+1}}, J_{n+1})-W^{\mathrm{I}}(I_{\mathbf{C}_{n}},II_{\mathbf{C}_{n+1}},J_{n+1})}{ \Delta I_{\mathbf{C}}}\]
2. the partioned discrete gradient \(\mathrm{D}_{II_{\mathbf{C}}}W^{\mathrm{I}}\) \[D_{II_{\mathbf{C}}}W^{\mathrm{I}}=\frac{1}{2}\left(\mathrm{D}_{II_{\mathbf{C}}}W^{ \mathrm{I}}|_{\Pi_{n+1}^{J},\Pi_{n}^{K}}+\mathrm{D}_{II_{\mathbf{C}}}W^{ \mathrm{I}}|_{\Pi_{n}^{J},\Pi_{n+1}^{K}}\right)\] (C.8) with being \(i=2\), \(J=1\), \(K=3\) we have \(\Pi^{2}=II_{\mathbf{C}}\), \(\Pi^{J}=I_{\mathbf{C}}\), \(\Pi^{K}=III_{\mathbf{C}}\) and thus \[\mathrm{D}_{II_{\mathbf{C}}}W^{\mathrm{I}}|_{\Pi_{n+1}^{J},\Pi_{n}^{K}} =\frac{W^{\mathrm{I}}(I_{\mathbf{C}_{n+1}},II_{\mathbf{C}_{n+1}}, J_{n})-W^{\mathrm{I}}(I_{\mathbf{C}_{n+1}},II_{\mathbf{C}_{n}},J_{n})}{ \Delta II_{\mathbf{C}}}\] (C.9) \[\mathrm{D}_{II_{\mathbf{C}}}W^{\mathrm{I}}|_{\Pi_{n}^{J},\Pi_{n+1}^{K}} =\frac{W^{\mathrm{I}}(I_{\mathbf{C}_{n}},II_{\mathbf{C}_{n+1}}, J_{n+1})-W^{\mathrm{I}}(I_{\mathbf{C}_{n}},II_{\mathbf{C}_{n}},J_{n+1})}{ \Delta II_{\mathbf{C}}}\]
3. the partioned discrete gradient \(\mathrm{D}_{J}W^{\mathrm{I}}\) \[D_{J}W^{\mathrm{I}}=\frac{1}{2}\left(\mathrm{D}_{J}W^{\mathrm{I}}|_{\Pi_{n+1}^{J },\Pi_{n}^{K}}+\mathrm{D}_{J}W^{\mathrm{I}}|_{\Pi_{n}^{J},\Pi_{n+1}^{K}}\right)\] (C.10) with being \(i=3\), \(J=\{1,2\}\), \(K=\emptyset\) we have \(\Pi^{3}=J\), \(\Pi^{J}=(I_{\mathbf{C}},II_{\mathbf{C}})\), \(\Pi^{K}=\emptyset\) and thus \[\mathrm{D}_{J}W^{\mathrm{I}}|_{\Pi_{n+1}^{J},\Pi_{n}^{K}} =\frac{W^{\mathrm{I}}(I_{\mathbf{C}_{n+1}},II_{\mathbf{C}_{n+1}}, J_{n+1})-W^{\mathrm{I}}(I_{\mathbf{C}_{n+1}},II_{\mathbf{C}_{n+1}},J_{n})}{ \Delta J}\] (C.11) \[\mathrm{D}_{J}W^{\mathrm{I}}|_{\Pi_{n}^{J},\Pi_{n+1}^{K}} =\frac{W^{\mathrm{I}}(I_{\mathbf{C}_{n}},II_{\mathbf{C}_{n}},J_{n +1})-W^{\mathrm{I}}(I_{\mathbf{C}_{n}},II_{\mathbf{C}_{n}},J_{n})}{\Delta J}\]

For the above partioned discrete gradients \((\mathrm{D}_{I_{\mathbf{C}}}W^{\mathrm{I}},\mathrm{D}_{II_{\mathbf{C}}}W^{ \mathrm{I}},\mathrm{D}_{J}W^{\mathrm{I}})\), conveniently shown in Eq. (58), the directionality property

\[\mathrm{D}_{I_{\mathbf{C}}}W^{\mathrm{I}}\,\Delta I_{\mathbf{C}}+ \mathrm{D}_{II_{\mathbf{C}}}W^{\mathrm{I}}\,\Delta II_{\mathbf{C}}+\mathrm{D} _{J}W^{\mathrm{I}}\,\Delta III_{\mathbf{C}}\] (C.12) \[= \frac{1}{2}\left(W^{\mathrm{I}}(I_{\mathbf{C}_{n+1}},II_{\mathbf{C} _{n}},J_{n})-W^{\mathrm{I}}(I_{\mathbf{C}_{n}},II_{\mathbf{C}_{n}},J_{n})\right.\] \[+W^{\mathrm{I}}(I_{\mathbf{C}_{n+1}},II_{\mathbf{C}_{n+1}},J_{n+1} )-W^{\mathrm{I}}(I_{\mathbf{C}_{n}},II_{\mathbf{C}_{n+1}},J_{n+1})\right)\] (C.13) \[+W^{\mathrm{I}}(I_{\mathbf{C}_{n}},II_{\mathbf{C}_{n+1}},J_{n+1} )-W^{\mathrm{I}}(I_{\mathbf{C}_{n}},II_{\mathbf{C}_{n}},J_{n+1})\] (C.14) \[+\frac{1}{2}\left(W^{\mathrm{I}}(I_{\mathbf{C}_{n+1}},II_{\mathbf{C} _{n+1}},J_{n+1})-W^{\mathrm{I}}(I_{\mathbf{C}_{n+1}},II_{\mathbf{C}_{n+1}},J_{n})\right.\] (C.15) \[+W^{\mathrm{I}}(I_{\mathbf{C}_{n}},II_{\mathbf{C}_{n}},J_{n+1})-W^{ \mathrm{I}}(I_{\mathbf{C}_{n}},II_{\mathbf{C}_{n}},J_{n})\right)\] (C.16) \[= \quad W^{\mathrm{I}}(I_{\mathbf{C}_{n+1}},II_{\mathbf{C}_{n+1}},J_{n +1})-W^{\mathrm{I}}(I_{\mathbf{C}_{n}},II_{\mathbf{C}_{n}},J_{n})\] (C.17)

is obviously fulfilled.
Partioned discrete gradients of \(\hat{W}^{\text{C}}\)

In analogy to Eq. (C.1)-Eq. (C.5) the partitioned discrete gradients \(\text{D}_{\text{C}}\hat{W}^{\text{C}}\), \(\text{D}_{\text{G}}\hat{W}^{\text{C}}\), and \(\text{D}_{C}\hat{W}^{\text{C}}\) of the strain energy \(\hat{W}^{\text{C}}(\text{C},\,\text{G},\,C)\) eventually are

1. the partioned discrete gradient \(\text{D}_{\text{C}}W^{\text{C}}\) \[\begin{split}&(\text{D}_{\text{C}}\hat{W}^{\text{C}})_{ij}=\frac{1} {2}(\partial_{\text{C}}\hat{W}^{\text{C}}(\text{C}_{n+\frac{1}{2}},\text{G}_{n },C_{n})\\ &+\partial_{\text{C}}\hat{W}^{\text{C}}(\text{C}_{n+\frac{1}{2}}, \text{G}_{n},C_{n}))_{kl}\,(\delta_{ki}\delta_{lj}-\frac{(\Delta\text{C})_{kl} (\Delta\text{C})_{ij}}{(\Delta\text{C})_{mn}\,(\Delta\text{C})_{mn}})\\ &+\frac{\hat{W}^{\text{C}}(\text{C}_{n+1},\text{G}_{n},C_{n})+ \hat{W}^{\text{C}}(\text{C}_{n+1},\text{G}_{n+1},C_{n+1})-\hat{W}^{\text{C}}( \text{C}_{n},\text{G}_{n},C_{n})-\hat{W}^{\text{C}}(\text{C}_{n},\text{G}_{n+ 1},C_{n+1})}{2\,(\Delta\text{C})_{mn}\,(\Delta\text{C})_{mn}}(\Delta\text{C})_{ ij}\end{split}\] (D.1)
2. the partioned discrete gradient \(\text{D}_{\text{G}}\hat{W}^{\text{C}}\) \[\begin{split}&(\text{D}_{\text{G}}\hat{W}^{\text{C}})_{ij}= \frac{1}{2}(\partial_{\text{G}}\hat{W}^{\text{C}}(\text{C}_{n+1},\text{G}_{n+ \frac{1}{2}},C_{n})\\ &+\partial_{\text{C}}\hat{W}^{\text{C}}(\text{C}_{n},\text{G}_{n +\frac{1}{2}},C_{n+1}))_{kl}\,(\delta_{ki}\delta_{lj}-\frac{(\Delta\text{G})_{ ij}(\Delta\text{G})_{kl}}{(\Delta\text{G})_{mn}\,(\Delta\text{G})_{mn}})\\ &+\frac{\hat{W}^{\text{C}}(\text{C}_{n+1},\text{G}_{n+1},C_{n})+ \hat{W}^{\text{C}}(\text{C}_{n},\text{G}_{n+1},C_{n+1})-\hat{W}^{\text{C}}( \text{C}_{n+1},\text{G}_{n},C_{n})-\hat{W}^{\text{C}}(\text{C}_{n},\text{G}_{ n},C_{n+1})}{2\,(\Delta\text{G})_{mn}\,(\Delta\text{G})_{mn}}(\Delta\text{G})_{ ij}\end{split}\] (D.2)
3. the partioned discrete gradient \(\text{D}_{C}\hat{W}^{\text{C}}\) \[\begin{split}\text{D}_{C}\hat{W}^{\text{C}}=& \frac{1}{2}\,\left(\frac{\hat{W}^{\text{C}}(\text{C}_{n+1},\text{G}_{n+1},C_{n+ 1})-\hat{W}^{\text{C}}(\text{C}_{n+1},\text{G}_{n+1},C_{n})}{\Delta C}\right. \\ &+\left.\frac{\hat{W}^{\text{C}}(\text{C}_{n},\text{G}_{n},C_{n+ 1})-\hat{W}^{\text{C}}(\text{C}_{n},\text{G}_{n},C_{n})}{\Delta C}\right)\,. \end{split}\] (D.3)

In comparison to an analytical Mooney-Rivlin model as proposed in [9] the above discrete gradients are pretty challenging to handle for the consistent linearization, since for the PANN constitutive model employed the discrete gradients \(\text{D}_{\text{C}}\hat{W}^{\text{C}}\), \(\text{D}_{\text{G}}\hat{W}^{\text{C}}\) are not constant and furthermore are dependend from all its dependencies \((\text{C},\text{G},C)\).

## Appendix E Discrete gradient of \(W^{\text{C}}\)

The projection-based discrete Gradient for \(W^{\text{C}}(\text{C})\) is given by

\[\text{D}_{\text{C}}W^{\text{C}}=\partial_{\text{C}}W^{\text{C}}(\text{C}_{n+ \frac{1}{2}})+\frac{W^{\text{C}}(\text{C}_{n+1})-W^{\text{C}}(\text{C}_{n})- \partial_{\text{C}}W^{\text{C}}(\text{C}_{n+\frac{1}{2}}):\Delta\text{C}}{ \Delta\text{C}:\Delta\text{C}}\,\Delta\text{C}\] (E.1)

where despite the chain rule which is necessary for the PANN constitutive model with its inputs \(\boldsymbol{\mathcal{I}}\) the discrete gradient is tensor-valued and therefore is by far complex when compared to the proposed version in Eq. (58).
