Learning solution of nonlinear constitutive material models using
  physics-informed neural networks: COMM-PINN
We applied physics-informed neural networks to solve the constitutive
relations for nonlinear, path-dependent material behavior. As a result, the
trained network not only satisfies all thermodynamic constraints but also
instantly provides information about the current material state (i.e., free
energy, stress, and the evolution of internal variables) under any given
loading scenario without requiring initial data. One advantage of this work is
that it bypasses the repetitive Newton iterations needed to solve nonlinear
equations in complex material models. Additionally, strategies are provided to
reduce the required order of derivative for obtaining the tangent operator. The
trained model can be directly used in any finite element package (or other
numerical methods) as a user-defined material model. However, challenges remain
in the proper definition of collocation points and in integrating several
non-equality constraints that become active or non-active simultaneously. We
tested this methodology on rate-independent processes such as the classical von
Mises plasticity model with a nonlinear hardening law, as well as local damage
models for interface cracking behavior with a nonlinear softening law. In order
to demonstrate the applicability of the methodology in handling complex path
dependency in a three-dimensional (3D) scenario, we tested the approach using
the equations governing a damage model for a three-dimensional interface model.
Such models are frequently employed for intergranular fracture at grain
boundaries. We have observed a perfect agreement between the results obtained
through the proposed methodology and those obtained using the classical
approach. Furthermore, the proposed approach requires significantly less effort
in terms of implementation and computing time compared to the traditional
methods.
Shahed Rezaei, Ahmad Moeineddin, Ali Harandi
Learning solution of nonlinear constitutive material models using physics-informed neural networks: COMM-PINN

###### Abstract

We applied physics-informed neural networks to solve the constitutive relations for nonlinear, path-dependent material behavior. As a result, the trained network not only satisfies all thermodynamic constraints but also instantly provides information about the current material state (i.e., free energy, stress, and the evolution of internal variables) under any given loading scenario without requiring initial data. One advantage of this work is that it bypasses the repetitive Newton iterations needed to solve nonlinear equations in complex material models. Additionally, strategies are provided to reduce the required order of derivation for obtaining the tangent operator. The trained model can be directly used in any finite element package (or other numerical methods) as a user-defined material model. However, challenges remain in the proper definition of collocation points and in integrating several non-equality constraints that become active or non-active simultaneously. We tested this methodology on rate-independent processes such as the classical von Mises plasticity model with a nonlinear hardening law, as well as local damage models for interface cracking behavior with a nonlinear softening law. Finally, we discuss the potential and remaining challenges for future developments of this new approach.

keywords: Physics-informed neural networks, Constitutive relations, Nonlinear material behavior, path-dependent material models, Finite element analysis. +
Footnote †: journal:

## 1 Introduction

Implementing and evaluating nonlinear material models is a challenging yet essential task for achieving reliable predictions in many applications. For example, advanced models may be needed to represent the complex plastic behavior of metals, soils, and polymers, or to account for the softening of materials as damage and microcracks evolve [1; 2; 3; 4; 5]. Typically, modelers use a thermodynamic framework to derive these models, which helps them to set up all the necessary equations and evolution laws for the internal state of the materials. Once this is done, complex codes must be programmed in a specific programming language, which can be a difficult and time-consuming process. To solve material constitutive equations, the iterative Newton method is often used, which can be highly effective. However, this method can also become computationally expensive, as it requires repeating the algorithm for each
integration point and nonlinear iteration within the finite element framework. In summary, developing and implementing nonlinear material models require significant expertise and effort, and access to appropriate computational resources. In addition to the previous point, material behavior is strongly influenced by lower scales. Multiscale material modeling involves solving a boundary value problem at the microscale. Once the microscale solution is obtained, the homogenized results are transferred to the structural problem at the macroscale [6]. However, the transformation of information between scales is a bottleneck in multiscale analysis, as it involves computationally-intensive tasks that can be time-consuming.

A recent trend in computational material mechanics is to use deep learning (DL) methods such as neural networks (NN) to bypass expensive calculations. The concept involves training a NN on a sufficient amount of data and using the trained NN as a surrogate model for a specific material system [7; 8; 9]. The success of this approach relies heavily on the quality and quantity of available data as well as the machine learning optimizer used. For more information, readers can refer to review articles [10; 11; 12], as well as research and review papers in [13; 14; 15; 16; 17].

In the previous approach, the predictions of the neural network (NN) are solely based on the available data, which raises concerns about the reliability of the NN for cases outside the training range. However, recent studies have shown that it is possible to improve the NN's predictive accuracy for unseen situations by incorporating thermodynamic constraints into the training process [18]. Liu et al. [19] developed a data-driven multiscale material modeling method that integrates principles of homogenization theory and essential physics with machine learning techniques. See also [20] for further extension of the method to reproduce the creep behavior of fiber-reinforced materials. Heider et al. [21] discussed the application of the informed-graph-based neural networks for anisotropic elasto-plastic behavior and proposed remedies to generate a frame-invariant machine-learning constitutive model. Flaschel et al. [22] developed an approach for the unsupervised discovery of material laws belonging to an unknown class of constitutive behavior (such as viscoelasticity and elastoplasticity). The main idea is to obtain the Helmholtz free energy and the dissipation potential from the data on displacements and reaction forces. Fuhg et al. [23] discussed a modular elastoplasticity formulation where different components of the model can be chosen to be either a phenomenological or a data-driven model depending on the amount of available information. Weber et al. [24] improved NN training with physical constraints for hyperelastic materials and showed that the introduced enhancements lead to better learning behavior. Kalina et al. [25] presented a data-driven multiscale framework utilizing physics-constrained NN for finite strain hyperelasticity problems. The idea is to set the Helmholtz free energy density as output and train the NN by using a set of invariants as the input of the NN. As a result, several physical principles are fulfilled automatically. See also investigations by [26] on polyconvex neural network constitutive models.

Masi et al. [27] proposed the thermodynamics-based NN for constitutive modeling. The authors proposed to learn the free energy of the material as a function of strain and internal state variables. See also investigations by [28] on learning effective strain energy of elastic cellular metamaterials. Masi and Stefanou [29] also extended their methodology to allow for identifying, from data and first principles, admissible sets of internal variables in complex materials.

In addition, to feed-forward NNs, other DL methods such as recurrent neural networks (RNN) are also suitable choices for path-dependent material behavior [11; 30; 31]. He and Chen [32] proposed a data-driven constitutive modeling approach with the consideration of thermodynamics principles for path-dependent materials. Bonatti and Mohr [33] considered an elastoplastic example and proposed an RNN architecture that respects self-consistency and performs best when applied to long sequences of small increments. Koeppe et al. [34] proposed
a physics-explaining approach, which interprets RNN. The authors investigated case studies on elastoplasticity, and viscoelasticity constitutive models. Danoun et al. [35] presented a hybrid physics-AI-based model to predict non-linear mechanical behaviors via RNN where specific thermodynamical constraints are considered during the training phase. See also developments by [36], [37] and [38]. In the reviewed articles, one requires enough data to train the NN.

Following the idea of physics-informed neural networks (PINN) [39], one can employ the physical laws as loss functions to enhance the network outcome [40; 41]. Once the underlying physics of the problem is completely known, the PINN can be trained without any initial data and solely based on the given set of equations [42; 43; 44]. Wei et al. [45] established a consistent neural network approach to model the constitutive behavior of interfaces by integrating physical conditions such as positive energy dissipation as additional training constraints in the loss function. Haghighat et al. [46] presented a physics-informed neural network formulation for constitutive modeling where inequality constraints of elastoplasticity theory are embedded in the PINN loss functions. In [47], the authors illustrated an extension for PINN in solid mechanics for the case with von Mises elastoplasticity and show that the PINN model can accurately predict the solution for a wide range of parameters. See also investigation by [48] for the extension to finite strain elastoplasticity. Eghbalian et al. [49] presented a NN architecture that is enriched with the physics of elastoplasticity and show that embedding this aspect enhanced the extrapolation capability for loading regimes outside the training data.

The reviewed works make significant contributions to two main categories. In the first group, NNs are trained on sufficient data obtained from numerical simulations or experimental measurements. To improve predictive accuracy, physical constraints are incorporated into these methods. In the second group, authors applied PINN techniques to study nonlinear material behavior for specific boundary value problems or for the characterization of constitutive models. While there has been previous research on using machine learning methods in material modeling, the specific combination of physics-informed neural networks and the ability to satisfy thermodynamic constraints and predict real-time material states without requiring initial data is a novel contribution to the field (see Fig. 1). It is worth mentioning that a simple FFNN will be used in this work, which can be integrated into any FE computation. Moreover, automatic differentiation will be used to bypass the derivation parts related to the tangent operator.

In the following sections, we will summarize the steps for deriving a consistent nonlinear material model. We will also present case studies that involve elastoplastic behavior with a nonlinear isotropic hardening law, as well as a local damage model with a nonlinear softening law applied to model interface cracking. Finally, we will discuss the results, draw conclusions, and suggest possible future developments.

Figure 1: Utilizing neural networks to substitute the return mapping algorithm in numerical material modeling.
## 2 Obtaining consistent constitutive relations and solving them with PINN

We start with the second law of thermodynamics in terms of the Clausius-Duhem inequality. The local form of dissipation inequality for the case of small deformation, pure mechanical case, and an isothermal process reads [4]

\[-\dot{\psi}(\boldsymbol{\varepsilon},\boldsymbol{\xi}_{k})+\boldsymbol{\sigma} :\dot{\boldsymbol{\varepsilon}}\stackrel{{!}}{{\geq}}0. \tag{1}\]

Here, \(\psi=\psi(\boldsymbol{\varepsilon},\boldsymbol{\xi}_{k})\) is the free energy of the material, \(\boldsymbol{\sigma}\) is the stress tensor, \(\boldsymbol{\xi}_{k}\) represent a set of internal (state) variables and \(\boldsymbol{\varepsilon}\) is the strain tensor. The above inequality is simplified as

\[(\underbrace{\boldsymbol{\sigma}-\partial_{\boldsymbol{\varepsilon}}\psi}_{= \boldsymbol{0}})\ \dot{\boldsymbol{\varepsilon}}-\Sigma_{k}\underbrace{\partial_{\boldsymbol{ \xi}_{k}}\psi}_{=\boldsymbol{q}_{k}}\ \dot{\boldsymbol{\xi}}_{k}\stackrel{{!}}{{\geq}}0. \tag{2}\]

Next, the thermodynamic conjugate forces follow automatically as

\[\boldsymbol{\sigma}=\partial_{\boldsymbol{\varepsilon}}\psi,\quad\boldsymbol{q }_{k}=\partial_{\boldsymbol{\xi}_{k}}\psi. \tag{3}\]

The remaining part of the dissipation inequality (i.e. \(\mathcal{D}=-\Sigma_{k}\boldsymbol{q}_{k}\ \dot{\boldsymbol{\xi}}_{k} \stackrel{{!}}{{\geq}}0\)), still has to be satisfied for an arbitrary process. To assure a positive dissipation rate and model the nonlinear material response, usually one introduces proper yield criterion \(YLD(\boldsymbol{\xi}_{k},...)\) as well as evolution equations \(EVL(\boldsymbol{\xi}_{k},...)\) for the state variables. In the following sections, we will investigate this part in more detail for two practical examples of plasticity and damage models.

At this point, we have all the ingredients to set up a neural network for predicting nonlinear material behavior. In Fig. 2, three architectures are proposed for solving equations in an arbitrary constitutive material model with physics informed neural networks (COMM-PINN).

Next, we will examine the pros and cons of each suggested network. In all cases, the user defines the free energy function and the required state variables. Furthermore, the inputs of the network are the current strain and the history of the material.

For case A, only the state variables are the outputs. All the thermodynamic conjugate forces (\(\boldsymbol{\sigma}^{i+1}\) and \(\boldsymbol{q}_{k}^{i+1}\)) and the consistent tangent operator \(\mathbb{C}^{i+1}=\mathrm{d}\boldsymbol{\sigma}^{i+1}/\mathrm{d}\boldsymbol{ \varepsilon}^{i+1}\) are calculated via automatic differentiation (AD). One possible drawback is the requirement for second-order differentiation, for which proper activation functions must be used [27]. Loss functions in this case are based on the yield function and evolution laws. For case B, the thermodynamic conjugate forces are defined as an additional output, and their definitions are used as additional loss functions (see also [50] and references therein for ideas on breaking the order of derivation in PINN). One advantage of the second architecture is that we only require first-order derivations to obtain the tangent operator. However, the network in this case might be denser to handle multiple outputs. In both of the above cases, the only requirement from the user is to define the free energy function, and all the thermodynamic forces, as well as updated state variables, are natural outcomes of the network. Finally, it should be noted that training the NN in the above cases might be very challenging due to the entangled nature of the equations. For case C, we make training easier by taking advantage of the analytical derivation of the free energy function and feeding the network with the analytical expression for the thermodynamic conjugate forces, such as the stress tensor. Obtaining such derivations can be a drawback for complicated energy functions, but it is usually straightforward to get such partial derivations. The advantage is fewer loss terms, less entanglement of equations, and first-order derivation for obtaining the tangent operator.
Figure 2: Three suggestions as potential neural network architectures and loss functions for solving constitutive material equations via physics-informed neural networks (COMM-PINN).
## 3 A case study on von mises plasticity model

### Model derivation

We limit ourselves to a 1D setting for small deformation. Nevertheless, the current investigations can be extended to higher dimensions and also case studies where more nonlinear behavior is expected. The strain variable \(\varepsilon\) is defined as derivation of displacement \(u\) with respect to the coordinate \(x\) (i.e. \(\varepsilon=\mathrm{d}u/\mathrm{d}x\)). In the case of elastoplastic behavior, the strain variable is additively decomposed into an elastic and a plastic part as

\[\varepsilon=\varepsilon_{e}+\varepsilon_{p}. \tag{4}\]

The Helmholtz free energy \(\psi\) is a contribution of the elastic part (\(\psi_{e}\)) and the plastic part (\(\psi_{p}\)). Therefore, we write the following expressions

\[\psi(\varepsilon,\varepsilon_{p},\xi_{p}) =\psi_{e}(\varepsilon,\varepsilon_{p})+\psi_{p}(\xi_{p}), \tag{5}\] \[\psi_{e}\left(\varepsilon,\varepsilon_{p}\right) =\frac{1}{2}E(\varepsilon-\varepsilon_{p})^{2},\] (6) \[\psi_{p}\left(\xi_{p}\right) =h_{1}\left(\xi_{p}+\frac{e^{-h_{2}\xi_{p}}-1}{h_{2}}\right). \tag{7}\]

In the elastic part of the energy (\(\psi_{e}\)), we have the material Young's modulus \(E\). Moreover, in the plastic part of the energy we introduce the parameter \(\xi_{p}\) which represents the accumulative plastic strain. Only isotropic hardening is considered and the constants \(h_{1}\) and \(h_{2}\) are corresponding plasticity hardening parameters. Using the chain rule of differentiation and considering the additive elastoplastic split of the strain rate, the Clausius-Duhem inequality in its isothermal local form reads

\[\left(\sigma-\partial_{\varepsilon}\psi\right)\,\dot{\varepsilon}-\partial_{ \epsilon_{p}}\psi\,\,\dot{\varepsilon}_{p}-\partial_{\xi_{p}}\psi\,\,\dot{ \xi}_{p}\stackrel{{!}}{{\geq}}0. \tag{8}\]

The above expression must hold for arbitrary thermomechanical processes. A common choice is to set the expressions in the bracket to zero. The state relations of the model as well as the thermodynamic conjugate forces follow automatically as

\[\sigma =\partial_{\varepsilon}\psi=\partial_{\varepsilon}\psi_{e}=E( \varepsilon-\varepsilon_{p}), \tag{9}\] \[q_{p} =\partial_{\xi_{p}}\psi=\partial_{\xi_{p}}\psi_{p}=h_{1}\,\,(1-e^{ -h_{2}\xi_{p}}). \tag{10}\]

The remaining part of the dissipation inequality can now be written as

\[\sigma\,\,\dot{\varepsilon}_{p}-q_{p}\,\,\dot{\xi}_{p}\stackrel{{! }}{{\geq}}0. \tag{11}\]

The plasticity yield criterion (based on the von Mieses model) is written as

\[\phi_{p}=|\sigma|-\left(\sigma_{y0}+q_{p}\right). \tag{12}\]

Moreover, the following evolution laws for the plastic internal variables are used:

\[\dot{\varepsilon}_{p} =\dot{\lambda}_{p}\,\,\frac{\partial\phi_{p}}{\partial\sigma}= \dot{\lambda}_{p}\,\,\mathrm{sgn}(\sigma), \tag{13}\] \[\dot{\xi}_{p} =-\dot{\lambda}_{p}\,\,\frac{\partial\phi_{p}}{\partial q_{p}}= \dot{\lambda}_{p}. \tag{14}\]

Finally, the plastic loading/unloading conditions of the model are taken into account:

\[\dot{\lambda}_{p}\geq 0,\quad\phi_{p}\leq 0,\quad\dot{\lambda}_{p}\,\,\phi_{p}=0. \tag{15}\]

By considering the relations 13, 14 as well as constraints in 15, the nonequality 11 is satisfied.
### Implementation and algorithmic aspects for the plasticity model

Algorithm 1 summarizes the return mapping algorithm commonly used to solve the plasticity equations. In order to solve for the unknowns, the equations for plasticity (or any other nonlinear material routines) must be linearized, which is accomplished using a local Newton solver at each integration point of the finite element method. Additionally, the user must derive equations for the tangent operator to be used in finite element programs. These tasks must be repeated for each global iteration step at the structural level and for each integration point in the domain, making them computationally expensive.

```
0: material parameters, current strain: \(\varepsilon^{i+1}\), history variables: \(\varepsilon^{i}_{p}\), \(\xi^{i}_{p}\)
0: current stress: \(\sigma^{i+1}\), current history variables: \(\varepsilon^{i+1}_{p}\), \(\xi^{i+1}_{p}\), tangent: \(C^{i+1}_{p}=\mathrm{d}\sigma^{i+1}/\mathrm{d}\varepsilon^{i+1}\)
1:\(k=1\)
2: Set trial values: \(\varepsilon^{i+1,k}_{p}=\varepsilon^{i}_{p}\), \(\xi^{i+1,k}_{p}=\xi^{i}_{p}\)
3:loop plasticity solverdo
4:if\(\phi^{tr}_{p}=|E(\varepsilon^{i+1}-\varepsilon^{i}_{p})|-\left(\sigma_{y0}+q_{p} (\xi^{i}_{p})\right)\leq 0\)
5:\(\varepsilon^{i+1}_{p}=\varepsilon^{i}_{p}\) and \(\xi^{i+1}_{p}=\xi^{i}_{p}\) and \(C^{i+1}=E\)
6:else solve for \(\varepsilon^{i+1,k+1}_{p}\) and \(\xi^{i+1,k+1}_{p}=\xi^{i}_{p}+\Delta\lambda^{k+1}_{p}\) (see Eqs 12,13)
7:\(r^{(1)}_{p}=\varepsilon^{i+1,k+1}_{p}-\varepsilon^{i}_{p}-(\xi^{i+1,k+1}_{p}- \xi^{i}_{p})\frac{\partial\phi_{p}}{\partial\sigma}\stackrel{{!}} {{=}}0\)
8:\(r^{(2)}_{p}=\phi^{i+1}_{p}=|\sigma^{i+1,k+1}|-\left(\sigma_{y0}+q_{p}\left( \xi^{i+1}_{p}\right)\right)\stackrel{{!}}{{=}}0\)
9:\([\varepsilon^{i+1,k+1}_{p}\quad\xi^{i+1,k+1}_{p}]^{T}=[\varepsilon^{i+1,k}_{p} \quad\xi^{i+1,k}_{p}]^{T}-\mathbf{K}^{-1}_{p}\mathbf{r}_{p}\)
10:endif
11:if\(\left|\varepsilon^{i+1,k+1}_{p}-\varepsilon^{i+1,k}_{p}\right|>\mathrm{tol}\) OR\(\left|\xi^{i+1,k+1}_{p}-\xi^{i+1,k}_{p}\right|>\mathrm{tol}\)
12:\(k++\)
13:CYCLE loop plasticity
14:else
15:\(g^{i+1}_{p}=g^{i+1,k+1}_{p}\), \(\xi^{i+1}_{p}=\xi^{i+1,k+1}_{p}\)
16:EXIT loop plasticity solver
17:endif
18:end loop
19: Compute stress \(\sigma^{i+1}=E(\varepsilon^{i+1}-\varepsilon^{i+1}_{p})\) see Eq. 33
20: Compute tangent \(C^{i+1}=E(1-\frac{E}{E+h_{1}h_{2}\ \exp(-h_{2}\xi^{i+1}_{p})})\)
```

**Algorithm 1** Solving the elastoplastic governing equations at pseudo time \(t^{i+1}\)

The superindex \(k\) in Algorithm 1 shows the \(k-\)th iteration in the internal loop _plasticity-solver_. Furthermore, we used the notation \(\Delta\lambda_{p}=\Delta t\ \lambda_{p}\), where \(\Delta t\) is the pseudo time step. Once the plasticity is active (\(\phi_{p}>0\)) the plastic residuals (\(r^{(1)}_{p}\) and \(r^{(2)}_{p}\)) are linearized according to the Newton-Raphson method to find the solution of the unknowns (\(\varepsilon^{i+1,k+1}_{p}\) and \(\xi^{i+1,k+1}_{p}\)). The matrix \(\mathbf{K}_{p}=\partial\mathbf{r}_{p}/\partial\mathbf{U}_{p}\) includes the derivatives of the plastic residual vector \(\mathbf{r}_{p}\) with respect to the unknowns (\(\mathbf{U}_{p}=\{\varepsilon^{i+1}_{p},\xi^{i+1}_{p}\}\)). If the changes in the solution of the unknown internal variables are smaller than a certain tolerance (\(\mathrm{tol}=10^{-10}\)), the obtained results are considered as the converged values and will be used to compute the stress value and the tangent operator. Finally, we are able to set up the global residuals at the finite element level [4; 51].
### PINN to solve for elasto-plastic behavior with nonlinear isotropic hardening law

Based on the methodologies proposed by Raissi et al. [39], the input for the network is the location of the collocation points. In the context of material modeling, the collocation points include the information on the current loading state (applied strain/gap) as well as admissible values for history variables (i.e. plastic strain, damage, etc.). Here we propose to use architecture C according to Fig. 2. Therefore, the output layer includes the current (updated) history variables (i.e. new plastic strain, damage, etc.). The other quantities such as stress, energy, and tangent operator are then calculated based on the known free energy equation of the material. In other words, we do not need any additional NN for computing those quantities.

The structure of each neural network takes the standard form where it can be split into a single input layer, several possible hidden layers, and the output layer. Each layer is connected to the next layer for transferring information [52]. In every single layer, there is no connection among its neurons. Therefore, we represent the information bypass from the \(l-1\) layer to \(l\) via the vector \(\mathbf{z}^{l}\). Every component of vector \(\mathbf{z}^{l}\) is computed by

\[z_{m}^{l}=a(\sum_{n=1}^{N_{l}}w_{mn}^{l}z_{n}^{l-1}+b_{m}^{l}),\quad l=1,\dots,L. \tag{16}\]

In Eq. (16), \(z_{n}^{l-1}\), is the \(n\)-th neuron within the \(l-1\)-th layer. The component \(w_{mn}\) shows the connection weight between the \(n\)-th neuron of the layer \(l-1\) and the \(m\)-th neuron of the layer \(l\). Every individual neuron in the \(l\)-th hidden layer owns a bias variable \(b_{m}^{l}\). The number \(N_{I}\) corresponds to the number of neurons in the \(l\)-th hidden layer. The total number of hidden layers is \(L\). The letter \(a\) stands for the activation function in each neuron. The activation function \(a(\cdot)\) is usually a non-linear function. The proper choice of the activation function is problem dependent and shall be obtained based on hyperparameter studies [42], [53].

For the case of plasticity, the input layer contains \(\varepsilon^{i+1}\), \(\varepsilon_{p}^{i}\) and \(\xi_{p}^{i}\). We denote the input layer via vector \(\mathbf{X}=\{\varepsilon^{i+1},\varepsilon_{p}^{i},\varepsilon_{p}^{i}\}\). The output layer is written as \(\mathbf{Y}=\{\varepsilon_{p}^{i+1},\xi_{p}^{i+1}\}\). We use separate fully connected FFNN for each output variable (see also Fig. 3). The trainable set of parameters of the network is represented by \(\mathbf{\theta}=\{\mathbf{W},\mathbf{b}\}\), which are the weights and biases of a neural network, respectively. Considering each neural network structure as \(\mathcal{N}\), the outcome of the constitutive material modeling via neural network reads

\[\varepsilon_{p}^{i+1}=\mathcal{N}_{\varepsilon_{p}}(\mathbf{X};\mathbf{\theta}),\quad \quad\xi_{p}^{i+1}=\mathcal{N}_{\xi_{p}}(\mathbf{X};\mathbf{\theta}). \tag{17}\]

The neural network outputs are functions of the trainable parameters and the training is done by minimizing physical loss functions. Next, we build the residuals and conditions for the introduced plasticity model in terms of the defined input and output layers. To do so, one requires to integrate the yield function as well as evolution equations for the internal variables into loss functions for the neural networks. By denoting the summation of total loss terms for plasticity by \(\mathcal{L}_{pt}\), and based on the Alg. 1, one defines \(\mathcal{L}_{pt}\) as

\[\mathcal{L}_{pt}=\underbrace{\mathcal{L}_{uep}+\mathcal{L}_{uxp}}_{\text{ elastic response}}+\underbrace{\mathcal{L}_{ev}+\mathcal{L}_{yl}}_{\text{plastic evolution}}+\underbrace{\mathcal{L}_{kte}+\mathcal{L}_{kty}}_{\text{KKT conditions}}. \tag{18}\]

In Eq. (18), different loss terms cover all the possible loading, unloading, and reloading scenarios. The first two terms (\(\mathcal{L}_{uep}\) and \(\mathcal{L}_{uxp}\)) guarantee that there is no evolution of plastic strain and accumulative plastic strain when the trial yield function (\(\phi_{p}^{tr}\)) is negative. Once \(\phi_{p}^{tr}>0\), the
main plastic residuals for evolution law as well as the current yield function become active which are denoted by \(\mathcal{L}_{\textit{evp}}\) and \(\mathcal{L}_{\textit{J}\textit{p}}\), respectively. Finally, to make sure that the KKT conditions are always satisfied, we have the last two loss terms (\(\mathcal{L}_{\textit{kep}}\) and \(\mathcal{L}_{\textit{kyp}}\)). All these relevant loss terms are summarized in what follows

\[\mathcal{L}_{\textit{uep}} =\text{MSE}\left((\varepsilon_{p}^{i+1}-\varepsilon_{p}^{i}) \text{Relu}(-\phi_{p}^{tr})\right), \tag{19}\] \[\mathcal{L}_{\textit{uxp}} =\text{MSE}\left((\xi_{p}^{i+1}-\xi_{p}^{i})\text{Relu}(-\phi_{p} ^{tr})\right),\] (20) \[\mathcal{L}_{\textit{evp}} =\text{MSE}\left((\varepsilon_{p}^{i+1}-\varepsilon_{p}^{i}-( \xi_{p}^{i+1}-\xi_{p}^{i})\text{sgn}(E(\varepsilon^{i+1}-\varepsilon_{p}^{i+1} )))\text{Relu}(\phi_{p}^{tr})\right),\] (21) \[\mathcal{L}_{\textit{ylp}} =\text{MSE}\left(\left(\text{abs}(E(\varepsilon^{i+1}-\varepsilon _{p}^{i+1}))-(\sigma_{y0}+h_{1}(1-\exp(-h_{2}\xi_{p}^{i+1})))\right)\text{ Relu}(\phi_{p}^{tr})\right),\] (22) \[\mathcal{L}_{\textit{kep}} =\text{MSE}\left(\text{Relu}(\phi_{p}^{i+1})\right),\] (23) \[\mathcal{L}_{\textit{kyp}} =\text{MSE}\left(\text{Relu}(-\xi_{p}^{i+1}+\xi_{p}^{i})\right). \tag{24}\]

In the above relation, \(\phi_{p}^{tr}\) is the so-called trial yield function which is evaluated by means of the (known) quantities at the input layer (see also Alg. 1):

\[\phi_{p}^{tr}=\text{abs}(E(\varepsilon^{i+1}-\varepsilon_{p}^{i}))-\left( \sigma_{y0}+h_{1}(1-\exp(-h_{2}\xi_{p}^{i}))\right). \tag{25}\]

For completeness, the expressions "abs" and "sgn" represent absolute value and sign function, respectively. The mean squared error is defined as

\[\text{MSE}(\bullet)_{\textit{type}}=\frac{1}{n}\sum_{i=1}^{n}(\bullet)^{2}, \tag{26}\]

where \(n\) is the number of observation points. The final loss term is minimized at every single collocation point. The mathematical optimization problem is written as

\[\mathbf{\theta}^{*}=\arg\min_{\mathbf{\theta}}\mathcal{L}_{pt}(\mathbf{X};\mathbf{\theta}), \tag{27}\]

where \(\mathbf{\theta}^{*}\) are the optimal trainable parameters (weights and biases) of the network.

**Remark 1** In the above mentioned loss terms, the "Relu" function acts as a switch (if condition). In other words, we differentiate between elastic loading/unloading and the evolution of plastic variables via the sign of the \(\phi_{p}^{tr}\). See also the same procedure in classical iterative solvers (i.e. Alg. 1). One can also use other options such as "sgn" function for this purpose. Based on our studies, we realized that utilizing "Relu" functions is more beneficial.

In Fig. 3, we summarized all loss terms and the idea behind constitutive material modeling via physics-informed neural networks. For the sake of clarity, all the input variables are denoted in blue color while the output (unknown) variables are presented in red color. Note that the material properties for the time being are kept constant and represented by green color.

After the training is completed, the network predicts the unknown internal variables (i.e. \(\varepsilon_{p}^{i+1}\) and \(\xi_{p}^{i+1}\) for the case of elastoplasticity). Having the updated internal variables as well as the current strain (\(\varepsilon^{i+1}\)), one can construct the free energy function (\(\psi^{i+1}\)) as well as the stress value (\(\sigma^{i+1}\)). Note that the functionality of stress can be obtained by analytical derivation of the energy function with respect to the strain which is according to the thermodynamical framework [27, 25], see also architecture C in Fig. 2. The final important parameter for the FE calculations is the tangent operator which is defined as the full derivation of the current
stress value with respect to the current strain value. This parameter is obtained via automatic differentiation of the trained network. The main difference of the current work compared to other studies (especially those in [27]), is the fact that here we do not provide any data to the FFNN.

**Remark 2** Although the material properties are kept constant in this study, one can also include them as additional input parameters. As a result, the NN learns the solution to the elastoplastic problem for vast range of given material properties. See investigations by [42, 54] in direction of using transfer learning to include the influence of material properties.

**Remark 3** In the above loss functions, one can add the contribution of the data in such form: \(\mathcal{L}_{data}=\text{MSE}\left(\varepsilon_{p}^{i+1}-\varepsilon_{alg}^{i +1}\right)+\text{MSE}\left(\xi_{p}^{i+1}-\xi_{alg}^{i+1}\right)\), where subindex \(alg\) represent the data coming from the return-mapping algorithm. We examine and compare such network designs in the result section.

**Remark 4** It is possible to derive the terms for the tangent operator via analytical derivation. The advantage is that one can insert the current values of strain and other state variables to find these quantities without running into the problem of vanishing gradients. The disadvantage is the possible difficulties in obtaining the analytical derivation. In the case of complicated differentiation, one can also use recent advances in parametric differentiation to avoid tedious derivations.

Figure 3: Network architecture and loss functions for the COMM-PINNs applied to elasto-plastic material models with isotropic Voce type hardening law.
## 4 A case study on damage model for interface cracking

### Model derivation

The gap value \(g\) is utilized to formulate the interface energy and mechanics. The gap variable is related to the displacement jump \(\left\langle\mathbf{u}\right\rangle\) or the distance between two opposite sides of the interface. Under mode I opening, we write the normal gap variable as

\[g=\left\langle\mathbf{u}\right\rangle\cdot\mathbf{n}, \tag{28}\]

where \(\mathbf{n}\) is the normal vector to the interface. Damage at the interface is represented via a scalar parameter \(d\) which is an internal variable and describes the loss of bonding forces. The (normal) traction parameter \(T\) is related to the gap value \(g\) and other internal variables such as damage \(d\). Similar to the stress-strain behavior for bulk we have the traction-separation relation for the mechanics of a sharp interface. The Helmholtz free energy \(\psi\) at the interface is a contribution of the elastic part (\(\psi_{e}\)) and damage hardening (\(\psi_{d}\)) which lead in total to the following expression:

\[\psi(g,d,\xi_{d}) =\psi_{e}(g,d)+\psi_{d}(\xi_{d}), \tag{29}\] \[\psi_{e}\left(\varepsilon,\varepsilon_{p}\right) =\frac{1}{2}f_{\mathrm{d}}(d)\ K\ g^{2},\] (30) \[\psi_{d}\left(\xi_{d}\right) =h_{1}\left(\xi_{d}+\frac{e^{-h_{2}\xi_{d}}-1}{h_{2}}\right), \tag{31}\]

Here, \(\xi_{d}\) is damage hardening variables. Moreover, \(h_{1}\) and \(h_{2}\) are interface damage hardening parameters and control the softening behavior [51]. The parameter \(K\) is the initial stiffness of the interface. Furthermore, \(f_{\mathrm{d}}(d)\) is a positive scalar value that varies between 1 for an undamaged interface and a critical value close to 0 for a completely damaged interface. A common choice for this function is as follows

\[f_{\mathrm{d}}(d)=(1-d)^{2}. \tag{32}\]

Applying the second law of thermodynamics in terms of the Clausius-Duhem inequality, consistent state relations of the model are derived. The conjugate forces are as follows

\[T =\partial_{g}\psi=\partial_{g}\psi_{e}=f_{\mathrm{d}}(d)K\ g, \tag{33}\] \[Y =-\partial_{d}\psi=-\partial_{d}\psi_{e}=-\frac{\mathrm{d}f_{d}} {\mathrm{d}d}\frac{1}{2}Kg^{2},\] (34) \[q_{d} =\partial_{\xi_{d}}\psi=\partial_{\xi_{d}}\psi_{d}=h_{1}(1-e^{-h_ {2}\xi_{d}}). \tag{35}\]

Similar to plasticity, the damage yield criterion is defined as

\[\phi_{d}(Y,q_{d})=Y-\left(Y_{0}+q_{d}\right). \tag{36}\]

In this model, the damage initiation point is controlled by the interface parameter \(Y_{0}\), whereas \(q_{d}\) accounts for nonlinear damage hardening. Evolution laws for the damage internal variables read:

\[\dot{d} =\dot{\lambda}_{d}\ \frac{\partial\phi_{d}}{\partial Y}=\dot{ \lambda}_{d}, \tag{37}\] \[\dot{\xi}_{d} =-\dot{\lambda}_{d}\ \frac{\partial\phi_{d}}{\partial q_{d}}=\dot{ \lambda}_{d}. \tag{38}\]

Finally, the loading/unloading conditions of the damage model are taken into account:

\[\dot{\lambda}_{d}\geq 0,\quad\phi_{d}\leq 0,\quad\dot{\lambda}_{d}\ \phi_{d}=0. \tag{39}\]
### Implementation and algorithmic aspects for the plasticity model

Similar to the case for plasticity, the superindex \(k\) in Algorithm 2 indicates the \(k\)-th iteration in the internal loop of the damage solver. In this algorithm, \(\Delta\lambda_{d}=\Delta t\ \lambda_{d}\), where \(\Delta t\) is the pseudo time step. When the damage is active (\(\phi_{d}>0\)), the damage residuals (\(r_{d}^{(1)}\) and \(r_{d}^{(2)}\)) are linearized using the Newton-Raphson method to solve for the unknowns (\(d^{k+1}\) and \(\xi_{d}^{k+1}\)). The matrix \(\mathbf{K}_{d}=\partial\mathbf{r}_{d}/\partial\mathbf{U}_{d}\) includes the derivatives of the damage residual vector \(\mathbf{r}_{d}\) with respect to the unknowns (\(\mathbf{U}_{d}=d^{i+1},\xi_{d}^{i+1}\)). If the changes in the solution of the unknown internal variables are smaller than a certain tolerance (\(\mathrm{tol}=10^{-10}\)), the results are considered converged and used to compute the traction values at the interface and the tangent operator. These quantities are used to set up the global residual vector at the finite element level [51].

```
0: interface's parameters, current gap: \(g^{i+1}\), history variables: \(d^{i}\), \(\xi_{d}^{i}\)
0: current traction: \(T^{i+1}\), current history variables \(d^{i+1}\), \(\xi_{d}^{i+1}\), tangent: \(C_{d}=\mathrm{d}T^{i+1}/\mathrm{d}g^{i+1}\)
1:\(k=1\)
2: Set trial values: \(g_{p}^{k}=g_{p}^{i}\), \(d^{k}=\xi_{d}^{i}=d^{i}\)
3:loop damage solverdo
4:if\(\phi_{d}^{tr}=K(g^{i+1})^{2}-(Y_{0}+q_{d}(\xi^{i}))\leq 0\)
5:\(d^{i+1}=d^{i}\), \(\xi_{d}^{i+1}=\xi_{d}^{i}\) and \(C_{d}=(1-d^{i+1})^{2}K\)
6:else solve for \(d^{k+1}\) and \(\xi_{d}^{k+1}=\xi_{d}^{i}+\Delta\lambda_{d}^{k+1}\) (see Eqs 36,37)
7:\(r_{d}^{(1)}=d^{k+1}-d^{i}-(\xi_{d}^{k+1}-\xi_{d}^{i})\frac{\partial\phi_{d}}{ \partial Y}\overset{!}{=}0\)
8:\(r_{d}^{(2)}=\phi_{d}=Y-\left(Y_{0}+q_{d}\left(\xi_{d}^{i+1}\right)\right) \overset{!}{=}0\)
9:\([d^{i+1,k+1}\quad\xi_{d}^{i+1,k+1}]^{T}=[d^{i+1,k}\quad\xi_{d}^{i+1,k}]^{T}- \mathbf{K}_{d}^{-1}\mathbf{r}_{d}\)
10:endif
11:if\(\left|d^{i+1,k+1}-d^{i+1,k}\right|>\mathrm{tol}\ \mathbf{\mathrm{OR}}\ \left|\xi_{d}^{i+1,k+1}-\xi_{d}^{i+1,k}\right|> \mathrm{tol}\)
12:\(k++\)
13:CYCLE loop damage solver
14:else
15:\(d^{i+1}=d^{i+1,k+1}\), \(\xi_{d}^{i+1}=\xi_{d}^{i+1,k+1}\)
16:EXIT loop damage solver
17:endif
18:endloop
19: Compute traction \(T=(1-d^{i+1})^{2}Kg^{i+1}\)
20: Compute tangent \(C_{d}\)
```

**Algorithm 2** Solving the damage governing equations at pseudo time \(t^{i+1}\)

### PINN to solve for local damage model or nonlinear softening behavior at the interface

For the case of local damage evolution (interface fracture), the input layer includes the gap value \(g^{i+1}\), the previous damage value \(d^{i}\), and the damage hardening variable \(\xi_{d}^{i}\). We denote the input layer via vector \(\mathbf{X}=\{g^{i+1},d^{i},\xi_{d}^{i}\}\). The output layer includes \(\mathbf{Y}=\{d^{i+1},\xi_{d}^{i+1}\}\). Similar to the previous case, we intend to use the separate fully connected feed-forward neural networks for each output variable (see also Fig. 4). The outcome of the constitutive material modeling via neural network reads

\[d^{i+1}=\mathcal{N}_{d}(\mathbf{X};\mathbf{\theta}),\ \ \ \ \ \xi_{d}^{i+1}=\mathcal{N}_{\xi_{d}}(\mathbf{X};\mathbf{ \theta}). \tag{40}\]
Next, we build the residuals for the introduced damage model. Here, one requires to build the damage yield function as well as evolution equations for the internal variables. These expressions will be used to construct loss functions for the neural networks. Denoting the summation of total loss terms for damage by \(\mathcal{L}_{dt}\), it is defined based on the Alg. 2, as

\[\mathcal{L}_{dt}=\underbrace{w_{ued}\mathcal{L}_{ued}+w_{uxd} \mathcal{L}_{uxd}}_{\text{elastic response}}+\underbrace{w_{evd}\mathcal{L}_{ evd}+w_{yld}\mathcal{L}_{yld}}_{\text{damage evolution}}+\underbrace{w_{ked}\mathcal{L}_{ked}+w_{kyd}\mathcal{L}_{kyd}}_{\text{KKT conditions}} \tag{41}\]

In Eq. (41), different loss terms cover all the possible loading, unloading, and reloading scenarios. The first two terms (\(\mathcal{L}_{ued}\) and \(\mathcal{L}_{uxd}\)) guarantee that there is no evolution of damage when the trial yield function (\(\phi_{d}^{tr}\)) is negative. Once \(\phi_{d}^{tr}>0\), the main damage residuals for evolution law, as well as the current yield function, become active which is denoted by \(\mathcal{L}_{evd}\) and \(\mathcal{L}_{yld}\), respectively. Finally, to make sure that the KKT conditions are always satisfied, we have the last two loss terms (\(\mathcal{L}_{ked}\) and \(\mathcal{L}_{kyd}\)). All these relevant loss terms are summarized in what follows

\[\mathcal{L}_{ued} =\text{MSE}\left((d^{i+1}-d^{i})\text{Relu}(-\phi_{d}^{tr}) \right), \tag{42}\] \[\mathcal{L}_{uxd} =\text{MSE}\left((\xi_{d}^{i+1}-\xi_{d}^{i})\text{Relu}(-\phi_{ d}^{tr})\right),\] (43) \[\mathcal{L}_{evd} =\text{MSE}\left((d^{i+1}-d^{i}-(\xi_{d}^{i+1}-\xi_{d}^{i})) \text{Relu}(\phi_{d}^{tr})\right),\] (44) \[\mathcal{L}_{yld} =\text{MSE}\left(\left((1-d^{i+1})K(g^{i+1})^{2}-(Y_{0}+h_{1}(1- \exp(-h_{2}\xi_{d}^{i+1}))\right)\text{Relu}(\phi_{d}^{tr})\right),\] (45) \[\mathcal{L}_{ked} =\text{MSE}\left(\text{Relu}(\phi_{d}^{i+1})\right),\] (46) \[\mathcal{L}_{kyd} =\text{MSE}\left(\text{Relu}(-\xi_{d}^{i+1}+\xi_{d}^{i})\right). \tag{47}\]

The so-called trial yield function \(\phi_{d}^{tr}\) is evaluated by means of the quantities at the input layer (see also Alg. 2):

\[\phi_{d}^{tr}=(1-d^{i})Kg^{i+1}-\left(Y_{0}+h_{1}(1-\exp(-h_{2} \xi_{d}^{i}))\right). \tag{48}\]

In Fig. 4, we summarized the main loss terms for the interface fracture model with nonlinear softening behavior. For the sake of clarity, all the input variables are denoted in blue color while the output (unknown) variables are presented in red color. Note that the material (interface) properties are kept constant and represented by green color.

After the training is completed, the network predicts the unknown internal variables (i.e. \(d^{i+1}\) and \(\xi_{d}^{i+1}\)). Having the updated internal variables as well as the current gap value (\(g^{i+1}\)), one can construct the free energy function as well as the traction value (\(T^{i+1}\)). Note that the functionality of traction is obtained by analytical derivation of the energy function with respect to the gap. Finally, the tangent operator which is defined as the full derivation of the current traction value with respect to the current gap value is obtained via automatic differentiation.

The final loss term is minimized at every single collocation point. The mathematical optimization problem is written as

\[\mathbf{\theta}^{*}=\arg\min_{\mathbf{\theta}}\mathcal{L}_{dt}(\mathbf{X}; \mathbf{\theta}), \tag{49}\]

where \(\mathbf{\theta}^{*}\) are the optimal trainable parameters (weights and biases) of the network.
## 5 Generation of collocation points

The location of collocation points plays a crucial role in training [17; 42]. _Collocation points_ are initial inputs for which the governing equations are going to be satisfied (collected). The location of collocation points can be different than _data points_ for which we know the solution in advance. The inputs of the NN are the given strain or gap at step \(i+1\) (i.e. \(\varepsilon^{i+1}\) or \(g^{i+1}\)) as well as state variables from the previous time step \(i\) (i.e. \(\varepsilon^{i}_{p}\) and \(\xi^{i}_{p}\) or \(d^{i}\) and \(\xi^{i}_{d}\)).

In this approach, we propose to assign admissible values to the set of inputs. However, two aspects should be considered when generating these points: 1) the desired range of the inputs and 2) avoiding any violation of important physical aspects when assigning random values. For instance, in the case of elastoplasticity, the accumulative plastic strain is always positive and greater than the plastic strain. Similarly, in the case of interface damage and the local damage model, the damage variable and damage hardening variable are always positive, below one, and identical. These restrictions are typically well-known and straightforward to account for. However, for more complex material models, simpler strategies will be required for future developments. It is essential to ensure that the generated values are physically consistent and within the desired range to prevent issues during training and improve the accuracy of the model.

Another idea is to generate random loading paths and gather the data from material models as initial values and substitute these models with the trained NN. For the current work we have tried both of these options and concluded that with both of these methods, one is able

Figure 4: Network architecture and loss functions for the COMM-PINN applied to interface failure with nonlinear softening law.
to solve the constitutive material equations. For the sake of simplicity, we report the more straightforward approach for the generation of collocation points as explained in the following algorithms.

For the case of elastoplastic behavior, readers are referred to Alg. 3. Here, we are using the following intervals for collocation points generation: \(bg_{e}=0\), \(st_{e}=0.01\), \(en_{e}=+1\), \(bg_{\varepsilon_{p}}=0\), \(st_{\varepsilon_{p}}=0.01\), \(en_{\varepsilon_{p}}=+1\), \(bg_{\xi_{p}}=0\), \(st_{\xi_{p}}=0.01\), and \(en_{\xi_{p}}=+1\). The main outcome of the Alg. 3 is the location of collocation points stored in the matrix \(\mathbf{IN}_{p}\) which has three columns for the three main inputs of the NN. The number of rows depends on the user's choice of intervals and steps.

```
0: step size: \(st\), begin \(bg\) and end of interval \(en\)
0: admissible set of inputs for training \(\mathbf{IN}_{p}\)
1:\(n=1\)
2:Initializing: \(\boldsymbol{\varepsilon}=[bg_{\varepsilon}:st_{\varepsilon}:en_{\varepsilon}]\), \(\boldsymbol{\varepsilon}_{p}=[bg_{\varepsilon_{p}}:st_{\varepsilon_{p}}:en_ {\varepsilon_{p}}]\), \(\boldsymbol{\xi}_{p}=[bg_{\xi_{p}}:st_{\xi_{p}}:en_{\xi_{p}}]\)
3:for l in \(\boldsymbol{\varepsilon}\)do
4:for j in \(\boldsymbol{\varepsilon}_{p}\)do
5:for k in \(\boldsymbol{\xi}_{p}\)do
6:\(\mathbf{IN}_{p}[n,1]=\boldsymbol{\varepsilon}[l]\)
7:\(\mathbf{IN}_{p}[n,2]=\boldsymbol{\varepsilon}_{p}[j]\)
8:if\(\boldsymbol{\xi}_{p}[k]\geq\boldsymbol{\varepsilon}_{p}[j]\)
9:\(\mathbf{IN}_{p}[n,3]=\boldsymbol{\xi}_{p}[k]\)
10:\(n++\)
11:end if
```

**Algorithm 3** Collocation point generation for solving elasto-plastic constitutive relations

For the case of the interface damage model, readers are referred to Alg. 4. We are using following intervals for collocation points generation: \(bg_{g}=0\), \(st_{g}=0.02\), \(en_{g}=+1\), \(bg_{d}=0\), \(st_{d}=0.02\), \(en_{d}=+1\), \(bg_{\xi_{d}}=0\), \(st_{\xi_{p}}=0.02\), and \(en_{\xi_{d}}=+1\). The main outcome of the Alg. 4 is the location of collocation points for training the damage model stored in the matrix \(\mathbf{IN}_{d}\).

```
0: step size: \(st\), begin \(bg\) and end of interval \(en\)
0: admissible set of inputs for training \(\mathbf{IN}_{d}\)
1:\(n=1\)
2:Initializing: \(\boldsymbol{g}=[bg_{g}:st_{g}:en_{g}]\), \(\boldsymbol{d}=[bg_{d}:st_{d}:en_{d}]\), \(\boldsymbol{\xi}_{d}=[bg_{\xi_{d}}:st_{\xi_{d}}:en_{\xi_{d}}]\)
3:for l in \(\boldsymbol{g}\)do
4:for j in \(\boldsymbol{d}\)do
5:for k in \(\boldsymbol{\xi}_{d}\)do
6:\(\mathbf{IN}_{d}[n,1]=\boldsymbol{g}[l]\)
7:\(\mathbf{IN}_{d}[n,2]=\boldsymbol{d}[j]\)
8:if\(\boldsymbol{\xi}_{d}[k]=\boldsymbol{d}[j]\)
9:\(\mathbf{IN}_{d}[n,3]=\boldsymbol{\xi}_{d}[k]\)
10:\(n++\)
11:endif
```

**Algorithm 4** Collocation point generation for solving the interface damage model

For the case of the interface damage model, readers are referred to Alg. 4. We are using following intervals for collocation points generation: \(bg_{g}=0\), \(st_{g}=0.02\), \(en_{g}=+1\), \(bg_{d}=0\), \(st_{d}=0.02\), \(en_{d}=+1\), \(bg_{\xi_{d}}=0\), \(st_{\xi_{p}}=0.02\), and \(en_{\xi_{d}}=+1\). The main outcome of the Alg. 4 is the location of collocation points for training the damage model stored in the matrix \(\mathbf{IN}_{d}\).

```
0: step size: \(st\), begin \(bg\) and end of interval \(en\)
0: admissible set of inputs for training \(\mathbf{IN}_{d}\)
1:\(n=1\)
2:Initializing: \(\boldsymbol{g}=[bg_{g}:st_{g}:en_{g}]\), \(\boldsymbol{d}=[bg_{d}:st_{d}:en_{d}]\), \(\boldsymbol{\xi}_{d}=[bg_{\xi_{d}}:st_{\xi_{d}}:en_{\xi_{d}}]\)
3:for l in \(\boldsymbol{g}\)do
4:for j in \(\boldsymbol{d}\)do
5:for k in \(\boldsymbol{\xi}_{d}\)do
6:\(\mathbf{IN}_{d}[n,1]=\boldsymbol{g}[l]\)
7:\(\mathbf{IN}_{d}[n,2]=\boldsymbol{d}[j]\)
8:if\(\boldsymbol{\xi}_{d}[k]=\boldsymbol{d}[j]\)
9:\(\mathbf{IN}_{d}[n,3]=\boldsymbol{\xi}_{d}[k]\)
10:\(n++\)
11:endif
```

**Algorithm 5** Collocation point generation for solving the interface damage model

For the case of elastoplastic behavior, readers are referred to Alg. 3. Here, we are using the following intervals for collocation points generation: \(bg_{e}=0\), \(st_{e}=0.01\), \(en_{e}=+1\), \(bg_{\xi_{p}}=0\), \(st_{\varepsilon_{p}}=0.01\), \(en_{\varepsilon_{p}}=+1\), \(bg_{\xi_{p}}=0\), \(st_{\xi_{p}}=0.01\), and \(en_{\xi_{p}}=+1\). The main outcome of the Alg. 3 is the location of collocation points stored in the matrix \(\mathbf{IN}_{p}\) which has three columns for the three main inputs of the NN. The number of rows depends on the user's choice of intervals and steps.

```
0: step size: \(st\), begin \(bg\) and end of interval \(en\)
0: admissible set of inputs for training \(\mathbf{IN}_{d}\)
1:\(n=1\)
2:Initializing: \(\boldsymbol{g}=[bg_{g}:st_{g}:en_{g}]\), \(\boldsymbol{d}=[bg_{d}:st_{d}:en_{d}]\), \(\boldsymbol{\xi}_{d}=[bg_{\xi_{d}}:st_{\xi_{d}}:en_{\xi_{d}}]\)
3:for l in \(\boldsymbol{g}\)do
4:for j in \(\boldsymbol{d}\)do
5:for k in \(\boldsymbol{\xi}_{d}\)do
6:\(\mathbf{IN}_{d}[n,1]=\boldsymbol{g}[l]\)
7:\(\mathbf{IN}_{d}[n,2]=\boldsymbol{d}[j]\)
8:if\(\boldsymbol{\xi}_{d}[k]=\boldsymbol{d}[j]\)
9:\(\mathbf{IN}_{d}[n,3]=\boldsymbol{\xi}_{d}[k]\)
10:\(n++\)
11:endif
```

**Algorithm 6** Collocation point generation for solving elasto-plastic constitutive relations
## 6 Results

The algorithms developed in the current work are implemented in the SciANN package [55] and the methodology can be easily transferred to other programming platforms. For all the reported results, the Adam optimizer is employed. We start with the results of the damage model due to its simplicity compared to the case of plasticity. Since the two models share some similarities, the majority of the hyperparameter studies are only reported for one of these models to avoid unnecessary repetition.

### Case studies on the interface damage model

The relevant network parameters for learning the damage model at the interface are summarized in Table 1. Some of these parameters are obtained after extensive (hyper) parameter studies. To make their influence clear to the readers, we will change these values and functions to investigate their impact on the obtained results.

The material parameters reported in Table 2 are chosen in a way to keep the stress values below value one and they do not necessarily represent any realistic material at this point. Nevertheless, one can calibrate the model based on some other measurements (see [51]) and normalize the model variables to keep them in the range between \(-1\) and \(1\). As we checked, the performance of the proposed methodology did not change for other chosen material properties.

The evolution of each loss term is shown in Fig. 5 as a function of epochs. The collocation points for the training are based on the last section. We used equal weightings for all the loss terms at the beginning and it turned out that the prediction of the NN model is acceptable by doing so (i.e. we have \(w_{ued}=w_{uxd}=w_{evd}=w_{ydl}=w_{ked}=w_{kyd}=1.0\)). Furthermore, for the results provided in Fig. 5, we utilize Relu as the activation function and 5 hidden layers with 100 neurons in each. The rest of the parameters are according to Table 1. All the loss functions decay simultaneously and we did not notice any significant improvement after 1000 epochs. We also observed the same response for all the loss terms through different hyper-parameter studies.

\begin{table}
\begin{tabular}{l l} \hline Parameter & Value \\ \hline Input, Output & \(\{g^{i+1},d^{i},\xi_{d}^{i}\}\), \(\{d^{i+1},\xi_{d}^{i+1}\}\) \\ Activation function & Relu, Tanh, Sigmoid \\ Number of layers and neurons per layer (\(L\), \(N_{l}\)) & \((5,\,10),\quad(5,\,50),\quad(5,\,100)\) \\ Batch size & \(500\) \\ Learning rate \(\alpha\), number of epochs & \((10^{-4},10^{3})\) \\ \hline \end{tabular}
\end{table}
Table 1: Summary of the COMM-PINN network parameters for the damage model.

\begin{table}
\begin{tabular}{l l l} \hline  & Unit & Value \\ \hline \hline Normal interface stiffness \(K\) & [MPa/mm] & \(5.0\) \\ Damage initiation criterion \(Y_{0}\) & [MPa mm] & \(0.1\) \\ Damage hardening parameter \(h_{1}\) & [MPa mm] & \(2.0\) \\ Damage hardening parameter \(h_{2}\) & [1/mm] & \(1.0\times 10^{2}\) \\ \hline \end{tabular}
\end{table}
Table 2: Material parameters for the interface damage model described in section 4.3.
In what follows, we will report the performance of the trained NN versus the response we obtained based on the classical material model routine (i.e. return mapping algorithm) which is our reference solution. Since the NN has not trained based on any given data (i.e. known and available solution), any arbitrary loading path we define for the strain value can be seen as a test case to evaluate its performance. Nevertheless, as we go through the results section, we examine the performance of the NN beyond the range of the collocation points. Some early prediction for various loading paths are shown in Fig. 6. Here, the first column represents the input gap value in time (i.e. \(g(t)\)), and the second column represents the obtained traction in time (i.e. \(T(t)\)). For the first and second rows, we applied \(g(t)=t\) and \(g(t)=t^{3}\), respectively. For the third and fourth rows, we have \(g(t)=0.5|t\,\sin(5\pi t)|+0.5|\sin(2\pi t)|\) and \(g(t)=1.0|t\,\sin(3\pi t)|\), respectively.

To evaluate the NN performance, we also over-plot the relative error for the predicted traction values. The error \(err\) is calculated via a point-wise comparison of the traction from NN \(T_{N}\) and calculated traction from the standard return mapping algorithm \(T_{M}\). Therefore, we have \(err=\frac{|T_{N}-T_{M}|}{T_{M}}\times 100\). The average error for different loading scenario is about 1% and in some extreme cases the maximum error goes up to 2.5%. We will discussed possible suggestions to improve the NN predictions through this section.

To avoid repetition, in rest of the paper, we chose a rather complicated loading path for all the studies which includes several loading, unloading, and reloading scenarios. Consequently, we can make sure that all the possible scenarios can become active. For the following studies, we force the gap vector to change through time via the following equation \(g(t)=2.0|t\,\sin(3\pi t)|\). In this loading case, the gap will go beyond the value of 1.5mm while we only train the NN with collocation points up to a gap value of 1.0. Therefore, toward the end of the loading path, we are checking the extrapolation capabilities of the train NN.

Note that even in the case of PINN without any initial data, the location of collocation points has a big impact on the final performance. In other words, by going beyond the range of the collocation points, there is no guarantee about the acceptable predictions of the NN.

Figure 5: Prominent loss terms for the interface damage model with nonlinear softening law utilizing the proposed COMM-PINN algorithm.
#### 6.1.1 Influence of neuron's number in NN architecture

In Fig. 7, we study the influence of various NN architectures. Here collocation points are generated with step size \(0.02\) and the Relu activation function is utilized. By comparing different NNs, we conclude that one requires enough neurons to capture all the complicated nonlinear behavior. At some point, adding extra neurons or layers is not beneficial anymore and might even result in the overfitting of the model. Similar performance was also observed for the number of layers where upon choosing a very small number of layers, the performance is really poor. Interestingly enough, for the red curve where we used a proper architecture, the NN

Figure 6: Predictions of the trained COMM-PINN for the interface damage model with nonlinear softening law.
can accurately capture the nonlinear response for all the loading, unloading, and reloading test cases. Towards the end of the loading where we go beyond the range of the provided collocation points for the training, we observe some deviations in the predicted response. In other words, even thermodynamical constraints can be violated when it comes to evaluating the NN response beyond the range of collocation points. On the positive side, one can always increase the range and density of the collocation points easily, unlike pure data-driven methods for which we need to obtain new data for the new range. It should be noted that using more collocation points will increase the training time. A fair comparison between pure physics-driven and pure data-driven NN training is not so trivial as they are different in the nature of loss functions and initial data sets and computational costs. Nevertheless, we will provide later on some initial comparisons to observe the deviation for a purely data-driven approach when it comes to extrapolation.

Figure 7: Influence of the NN layers on the prediction of the network. (a) the loading (gap over time) is based on \(g(t)=2.0|t\)\(\sin(3\pi t)|\). (b) Predicted traction \(T(t)\) via train NN over time. (c) predicted damage variable \(d(t)\) via train NN over time. (d) predicted damage hardening variable \(\xi_{d}(t)\) via train NN over time. (e) predicted free energy function \(\psi_{d}(t)\) over time. (f) predicted traction versus the gap value.
#### 6.1.2 Influence of activation function

In Fig. 8, we look into the influence of various choices for the activation functions on the predictions. The explanations of the initial training set are similar to those before. Here 5 layers with 100 neurons in each are utilized for all the case studies. We concluded that the Relu function shows the best performance compared to others. Although not reported, our conclusion is the same even by using other activation functions such as Swish and Softplus. The better performance for the Relu function is perhaps due to its nature with a sharp transition that can capture the sudden change of behavior in the unloading, reloading, and yield condition.

#### 6.1.3 Influence of collocation point density

In Fig. 9, we look into the influence of various choices for the density of the collocation points on the predictions. The chosen activation function for this study is the Relu function and the

Figure 8: Influence of the choice of activation function on the prediction of the network. (a) the loading (gap over time) is based on \(g(t)=2.0|t\,\sin(3\pi t)|\). (b) Predicted traction \(T(t)\) via train NN over time. (c) predicted damage variable \(d(t)\) via train NN over time. (d) predicted damage hardening variable \(\xi_{d}(t)\) via train NN over time. (e) predicted free energy function \(\psi_{d}(t)\) over time. (f) predicted traction-separation relation.
network has 5 layers with 100 neurons. The density of the collocation points is varied from 0.1 to 0.02. We concluded that even with a very low number of points, the NN can capture the overall behavior of such a nonlinear loading path. We also observe that the performance of the NN model is improved by feeding more collocation points to the NN for the training process. Furthermore, having more points will improve the extrapolation capabilities of the NN (see the obtained response after time 0.8). On the other hand, more collocation points will significantly increase the computational cost of the training. Training with step size 0.02 takes twice the time compared to step size 0.1. One may have to adapt the architecture of the NN adequately based on the number of collocation points.

It should be noted that in all the cases, the test loading path is generated with 50 points in time (i.e. the time step size is 0.02). As we will discuss this further in the next section, one can choose different step sizes and call the same NN and yet have acceptable performance.

Figure 9: Influence of the step size for generation of collocation points on the predictions. (a) the loading (gap over time) is based on \(g(t)=2.0|t\,\sin(3\pi t)|\). (b) Predicted traction \(T(t)\) via train NN over time. (c) predicted damage variable \(d(t)\) via train NN over time. (d) predicted damage hardening variable \(\xi_{d}(t)\) via train NN over time. (e) predicted free energy function \(\psi_{d}(t)\) over time. (f) predicted traction versus the gap value.
### Case studies on the plasticity model

The relevant network parameters for the case of elastoplasticity are summarized in Table 3.

There are 4 input parameters: \(E\), \(\sigma_{y0}\), \(h_{1}\), and \(h_{2}\) for the plasticity model. The chosen values for the elastoplastic model are summarized in Table 4. The current approach may also be used for reverse analysis and calibrating these material parameters based on given data [46].

The location of collocation points is according to the explanations in the previous section. Before we report the results, we introduce two modifications to enhance the NN model performance. The first modification is about balancing the value of different loss functions. For the second modification, we tend to smoothen the hard switch conditions in the plastic model.

Initially, we used equal weightings for all the loss terms which resulted in an acceptable performance within the training range but does not perform very well when we test the NN beyond the range of the collocation points. After some numerical studies, we concluded that this issue might be due to unbalance nature of different loss terms in the training process. We performed a systematic parameter study and changed the weighting of different loss terms to make them closer to each other. The process is as follows. At first, we start with equal weightings, and then we raise the weighting of those loss functions which are very low in magnitude. By doing so, we observed a huge improvement in the obtained results. This idea is also according to adaptive weighting strategies available in some programming packages. Other adaptive weight techniques should be investigated further in future developments to find more easy-to-use strategies. One of the final optimum results is achieved by selecting \(w_{uep}=100\), \(w_{uxp}=100\), \(w_{evp}=1\), \(w_{ylp}=10\), \(w_{kep}=100\), \(w_{kyp}=10\).

In the next modification, we loosen the condition of Relu and sgn function in Eqs. 19 to 24. Therefore, we replace \(\text{Relu}(x)\) with \(\text{Swish}(x)=x\ \text{Sigmoid}(Rx)\), where \(R\) determines how sharp is the transition zone. Finally, we replace the \(\text{sgn}(x)\) with \(\text{Sigmoid}(Rx)\). In the current study, we used \(R=300\). To make the influence of the introduced modifications clear to the readers, we perform studies without these modifications as well. The evolution of each loss term is shown in Fig. 10, where we use 5 hidden layers with 80 neurons in each. Furthermore, for the results provided in Fig. 10, we considered the above modifications. The rest of the parameters are according to Table 3. Most of the loss functions decay simultaneously, and we did not notice any significant improvement after 1000 epochs.

\begin{table}
\begin{tabular}{l l l} \hline  & Unit & Value \\ \hline \hline Elastic stiffness \(E\) & [MPa] & 3.0 \\ Yield stress \(\sigma_{y0}\) & [MPa] & 0.6 \\ Plasticity hardening parameter \(h_{1}\) & [MPa] & 0.4 \\ Plasticity hardening parameter \(h_{2}\) & [\(-\)] & 10.0 \\ \hline \end{tabular}
\end{table}
Table 4: Material parameters for the plasticity model described in section 3.3.

\begin{table}
\begin{tabular}{l l} \hline Parameter & Value \\ \hline Input, Output & \(\{\varepsilon^{t+1},\varepsilon^{t}_{p},\varepsilon^{t}_{p}\}\), \(\{\varepsilon^{t+1}_{p},\varepsilon^{t+1}_{p}\}\) \\ Activation function & Relu \\ Number of layers and neurons per layer (\(L\), \(N_{l}\)) & (5, 80) \\ Batch size & 100 \\ Learning rate \(\alpha\), number of epochs & \((10^{-4},10^{3})\) \\ \hline \end{tabular}
\end{table}
Table 3: Summary of the COMM-PINN network parameters for the plasticity damage model.
#### 6.2.1 Influence of weightings in the loss function

The studies presented in Fig.11 were conducted using a loading path of \(\varepsilon(t)=2.0|t\mathrm{sin}(3\pi t)|\). The time step size used to generate this loading path was 0.01. The material was initially at rest without any prior history. The results show that the yield point, nonlinear hardening path, and stress saturation for large plastic strains were accurately captured by the NN model. One should note that the NN's response for accumulative plastic strain larger than one can be considered as extrapolation. For the case with equal weights, the performance of the NN is not reliable beyond the range of the collocation points. However, we observed that by tuning and balancing the loss functions, the NN was able to accurately capture the material behavior even beyond the range of the collocation points.

#### 6.2.2 Modifying switch option in the loss functions

In Fig.12, we studied the case where the function Relu is used for switching between different elastic and plastic conditions (see also Eqs.19 to 24). We found that using a smoother version of this function (i.e., \(\mathrm{Relu}(x)\) with \(x\)\(\mathrm{Sigmoid}(Rx)\)) is more beneficial. Our findings indicate that this is particularly important for the condition in \(\mathrm{sgn}(x)\), which is replaced by \(\mathrm{Sigmoid}(Rx)\), where we use \(R=300\). See also investigations by Haghighat et al. [46].

#### 6.2.3 Influence of the time step size

Here, we examine the output by resolving the loading path with different time step sizes. These time steps are chosen in a way to ensure that they do not overlap with the collocation points. In the first row of Fig. 13, the NN is evaluated using a time step size of 0.1. In addition to the stress-strain curve, we also plot the tangent operator, which is computed solely via automatic differentiation. For the second row, we evaluate the model using a step size of 1/80 to ensure that there is no overlap with the collocation points (which are generated using a step size of 1/100). Again, the predictions are acceptable. However, as we use finer time steps (see third row with a step size of 1/200), we observe some deviation from the correct solution, and the tangent operator shows slight oscillations. As a remedy, one can always utilize finer steps in the generation of the collocation points to improve the accuracy of the NN predictions.

Figure 10: Prominent loss terms for the plasticity model utilizing the proposed COMM-PINN algorithm.
Figure 11: Influence of the adapted weightings for the different loss functions. (a) stress versus strain curve. (b) predicted stress \(\sigma(t)\). (c) predicted plastic strain \(\epsilon_{p}(t)\). (d) predicted accumulative plastic strain \(\xi_{p}(t)\).

Figure 12: Effect of modifying the switch option in the loss functions. (a) stress versus strain curve. (b) predicted stress \(\sigma(t)\). (c) predicted plastic strain \(\epsilon_{p}(t)\). (d) predicted accumulative plastic strain \(\xi_{p}(t)\).
#### 6.2.4 Comparison with the pure data-driven approach

Remark 3 suggests that a purely data-driven approach can be used for the same analysis. In this approach, the data on the incremental evolution of the state variables is obtained from a solution of the governing equations on some random loading paths using the return mapping algorithm (see Alg.1). Once the data is obtained, we use the same neural network architecture described in Table 3. In this case, the loss function is based on the difference between the predicted and actual values of plastic and accumulative plastic strain, i.e., \(\mathcal{L}_{data}=\text{MSE}\left(\varepsilon_{p}^{i+1}-\varepsilon_{alg}^{ i+1}\right)+\text{MSE}\left(\xi_{p}^{i+1}-\xi_{alg}^{i+1}\right)\). We generated 625 loading paths for data generation, which are shown in gray in part (a) of Fig.14. All loading paths are limited to a strain of 1, which is similar to the range of collocation points. Furthermore, each loading path is discretized by a time step size of 0.01, which gives us a total of 62500 data points. Although the initial data points and collocation points are similar, it is not straightforward to have an identical initial condition for these two neural networks.

Figure 13: Obtained stress-strain curves and the tangent operator via the train NN and classical material model routine. The time step size for creating the loading path is refined to evaluate the performance of the NN for a strain path that does not fully overlap with the collocation points.
In Fig. 14, we compare the two methodologies and their ability to predict material behavior for complex unseen loading paths. However, in parts (e) and (f) of Fig. 14, we examine the response of these methods in more detail during the first unloading cycle. While the data-driven method shows acceptable predictions, it violates some fundamental physical aspects, such as the decay of plastic strain during elastic unloading, which is unacceptable. Conversely, the results of the COMM-PINN method follow the expected straight line. For the final loading cycle, we evaluate the extrapolation capabilities of the NN models from both the physics-based and purely data-driven methods. Both methods may fail when extrapolating to extreme conditions. However, the physics-based method with physical constraints performs slightly better. In future developments, it would be interesting to explore the combination of these methods. By adding more collocation points outside the training region, we can enhance the prediction of the hybrid NN model, particularly in cases where there are limited data points for a limited strain interval.

Figure 14: Comparison of the proposed approach against pure data-driven methods. (a) applied strain path is shown in blue and all the loading paths used for training the data-driven method are shown in gray. (b) stress versus strain curve. (c) predicted plastic strain \(\epsilon_{p}(t)\). (d) predicted accumulative plastic strain \(\xi_{p}(t)\).
#### 6.2.5 Integration of the COMM-PINN model with finite element method

We tested the performance of the COMM-PINN model by integrating it into a finite element package. Since the model is trained based on a 1D analysis, we examined a truss-based structure (or a simplified meta-material) where each rod is represented by a single 1D finite element mesh using standard linear shape functions. For this purpose, we utilized the package _trusspy_[56]. The geometry and applied boundary conditions are shown in the upper part of Fig. 15, where we fixed the whole structure on the left-hand side and loaded it from the right vertical edge. After the loading phase, the whole system was unloaded.

The obtained reaction force from this process is shown in the lower part of Fig.15. We first performed finite element calculations based on the standard return mapping algorithm and introduced a VOC-type plasticity model. The results are shown in blue. The same boundary value problem was calculated using the trained COMM-PINN model instead of the standard user-defined material model. By comparing the red and blue curves, we conclude that the new method has the potential to replace the classical method in a finite element calculation without any changes to the available codes. Finally, on the right-hand side of Fig.15, we compared the stress distribution from the two methodologies at the end of the loading path and observed identical results.

Figure 15: Comparison of the proposed approach against classical return mapping algorithm when they are integrated into the finite element package (trusspy [56]). Top left: the geometry and boundary conditions for a simple truss-based structure. Bottom: obtained reaction force from two methodologies. Right: distribution of stress in the truss-based system, obtained from the two methodologies.
## 7 Conclusion and outlooks

This work demonstrates the potential of physics-informed neural networks in solving constitutive relations for material mechanics. The approach is applicable when the free energy function and internal variables of the model are known beforehand. We proposed strategies to construct the loss functions for training the network, utilizing expressions for thermodynamic forces, evolution laws, and yield functions. This eliminates the need for labeled data during training. We also suggested techniques to simplify the process and reduce differentiation degrees for evaluating the tangent operator.

Compared to classical return mapping algorithms and other purely data-driven methods, the proposed COMM-PINN method offers several advantages. The model can instantly predict the material response for any loading path within the range of the collocation points. By tuning the weighting of loss functions and performing hyperparameter studies, the model's extrapolation capabilities can be enhanced. The proposed method avoids complex linearizations and yields the tangent operator naturally from the neural network. Labeled data is not necessary for model training, and everything is based on the governing material equations derived from standard thermodynamic principles. Interestingly, the proposed approach showed better performance than purely data-driven methods for a complex test loading path. Finally, the trained material routine can be easily integrated into any finite element package.

In future studies, it will be necessary to extend the framework and train the system for other material properties. Investigations into multi-dimensional cases and more complex material routines are also required to better understand the advantages of this approach. Ideally, the NN can be trained once for each material model and used as a substitute for standard material subroutines, eliminating the need to rewrite codes and material models in different programming languages.

**Data Availability**:

The codes are available at COMM-PINN.

**Author Statement**:

S.R.: Conceptualization, Software, Supervision, Writing - Review & Editing. A. M.: Software, Writing - Review & Editing. A. H.: Software, Writing - Review & Editing.
