A deep learning framework for solution and discovery in solid mechanics

The equations expressing momentum balance, the constitutive model and the kinematic relations are, respectively,

Here, $\sigma_{ij}$ denotes the Cauchy stress tensor. For the two-dimensional problems considered here $i,j=1,2$ (or $i,j=x,y$). We use the summation convention, and an subscript comma denotes partial derivative. The function $f_{i}$ denotes a body force, $u_{i}$ represents the displacements, $\varepsilon_{ij}$ is the infinitesimal stress tensor and $\delta_{ij}$ is the Kronecker delta. The Lamé parameters $\lambda$ and $\mu$ are the quantities to be inferred using PINN.

In this section, we provide an overview of the Physics-Informed Neural Networks (PINN) architecture, with emphasis on their application to model inversion. Let $\mathcal{N}(\mathbf{x};\mathbf{W},\mathbf{b}):\mathbb{R}^{d_{\mathbf{x}}}\rightarrow\mathbb{R}^{d_{\mathbf{y}}}$ be an $L$-layer neural network with input vector $\mathbf{x}$, output vector $\mathbf{y}$, and network parameters $\mathbf{W},\mathbf{b}$. This network is a feed-forward network, meaning that each layer creates data for the next layer through the following nested transformations:

where $\mathbf{z}^{0}\equiv\mathbf{x}$ and $\mathbf{z}^{L}\equiv\mathbf{y}$ are inputs and outputs of the model, $\mathbf{W}^{l},\mathbf{b}^{l}$ are parameters of each layer $l$, known as weights and biases, respectively. The functions $\sigma^{l}$ are called activation functions and make the network nonlinear with respect to the inputs. For instance, an ANN functional of some field variable, such as displacement $u(\mathbf{x})$, with three hidden layers and with $\sigma^{l}=\tanh$ as the activation function for all layers except the last can be written as

This model can be considered as an approximate solution for the field variable $u$ of a partial differential equation.

In the PINN architecture, the network inputs (also known as features) are space and time variables, i.e., $(x,y,z,t)$ in Cartesian coordinates, which makes it meaningful to perform the differentiation of the network’s output with respect to any of the input variables. Classical implementations based on finite difference approximations are not accurate when applied to deep networks (see [35] for a review). Thanks to modern graph-based implementation of the feed-forward network (e.g., Theano [29], Tensorflow [30], MXNet [36]), this can be carried out using Automatic Differentiation at machine precision, therefore allowing for many hidden layers to represent nonlinear response. Hence, evaluation of a partial differential operator $\mathcal{P}$ acting on $u$ is achieved naturally with graph-based differentiation and can then be incorporated in the cost function along with initial and boundary conditions as:

where $\partial\Omega$ is the domain boundary, $u_{0}-u_{0}^{*}$ is the initial condition at $t=t_{0}$, and $0^{*}$ indicates the expected (true) value for the differential relation $\mathcal{P}u$ at any given training point. The norm $|\cdot|$ of a generic quantity $g$ defined in $\Omega$ denotes $\frac{1}{N}\sum_{i=1}^{N}g(\mathbf{x}_{i})^{2}$ where the $\mathbf{x}_{i}$’s are the spatial points where the data is known. The dataset is then fed to the neural network and an optimization is performed to evaluate all the parameters of the model, including the parameters of the PDE.

Different algorithms that can be used to train a neural network. Among the choices available in Keras [37] we use the Adam optimization scheme [38], which we have found to outperform other choices such as Adagrad [39], for this task. Several algorithmic parameters affect the rate of convergence of the network training. Here we adopt the terminology in Keras [37], but the terminology in other modern machine learning packages is similar. The algorithmic parameters include batch-size, epochs, shuffle, and patience. Batch-size controls the number of samples from a dataset used to evaluate one gradient update. A batch-size of 1 would be associated with a full stochastic gradient descent optimization. One epoch is one round of training on a dataset. If a dataset is shuffled, then a new round of training (epoch) would result in an updated parameter set because the batched-gradients are evaluated on different batches. It is common to re-shuffle a dataset many times and perform the back-propagation updates. The optimizer may, however, stop earlier if it finds that new rounds of epochs are not improving the cost function. That is where the last keyword, patience, comes in. This is mainly because we are dealing with non-convex optimization and we need to test the training from different starting points and in different directions to build confidence on the parameters evaluated from minimization of the cost-function on a dataset. Patience is the parameter that controls when the optimizer should stop the training.

There are three ways to train the network: (1) generate a sufficiently large number of datasets and perform a one-epoch training on each dataset, (2) work on one dataset over many epochs by reshuffling the data, and (3) a combination of these. When dealing with synthetic data, all approaches are feasible to pursue. However, strategy (1) above is usually impossible to apply in practice, specially in space, where sensors are installed at fixed and limited locations. In the original work on PINN [22], approach (1) was used to train the model, where datasets are generated on random space discretizations at each epoch. Here, we follow approach (2) to use training data that we could realistically have in practice. For all examples, unless otherwise noted, we use a batch-size of 64, a limit of 10,000 epochs with shuffling, and a patience of 500 to perform the training.

To illustrate the application of the proposed approach, we consider an elastic plane-strain problem on the unit square (Fig. 1), subject to the boundary conditions depicted in the figure. The body forces are:

The exact solution of this problem is

which is plotted in Fig. 2, for parameter values of $\lambda=1$, $\mu=0.5$, and $Q=4$.

Due to the symmetry of the stress and strain tensors, the quantities of interest for a two-dimensional problem are $u_{x}$, $u_{y}$, $\varepsilon_{xx}$, $\varepsilon_{yy}$, $\varepsilon_{xy}$, $\sigma_{xx}$, $\sigma_{yy}$, $\sigma_{xy}$. There are a few potential architectures that we can use to design our network. The input features (variables) are the spatial coordinates $(x,y)$, for all the network choices. For the outputs, a potential design is to have a densely connected network with two outputs as $(u_{x},u_{y})$. Another option is to have two densely connected independent networks with only one output each, associated with $u_{x}$ and $u_{y}$, respectively (Fig. 3). Then, the remaining quantities of interest, i.e., $\sigma_{ij},\varepsilon_{ij}$, can be obtained through differentiation. Alternatively, we may have $(u_{x},u_{y},\sigma_{xx},\sigma_{yy},\sigma_{xy})$ or $(u_{x},u_{y},\varepsilon_{xx},\varepsilon_{yy},\varepsilon_{xy})$ as outputs of one network or multiple independent networks. As can be seen from Fig. 3, these choices affect the number of parameters of the network and how different quantities of interest are correlated. Equation (3) shows that the the feed-forward neural network imposes a special functional form to the network that may not necessarily follow any cross-dependence between variables in the governing equations (1). Our data shows that using separate networks for each variable results in a far more effective strategy. Therefore, we propose to have variables $u_{x},u_{y},\sigma_{xx},\sigma_{yy},\sigma_{xy}$ defined as independent ANNs as our architecture of choice (see Fig. 4), i.e.

The cost function is defined as

The quantities with asterisks represent given data. We will train the networks so that their output values are as close as possible to the data, which may be real field data or, in this paper, synthetic data from the exact solution to the problem or the result of a high-fidelity simulation. The values without asterisk represent either direct outputs of the network (e.g., $u_{x}$ or $\sigma_{xx}$; see Eq. (8)) or quantities obtained through automatic graph-based differentiation [35] of the network outputs (e.g., $\varepsilon_{xx}=u_{x,x}$). In Eq. (9), $f_{x}^{*}$ and $f_{y}^{*}$ represent data on the body forces obtained as $f_{i}^{*}=-\sigma_{ij,j}^{*}$.

The different terms in the cost function represent measures of the error in the displacement and stress fields, the momentum balance, and the constitutive law. This cost function can be used for deep-learning-based solution of PDEs as well as for identification of the model parameters. For the solution of PDEs, $\lambda$ and $\mu$ are treated as fixed numbers in the network. For parameter identification, $\lambda$ and $\mu$ are treated as network parameters that change during the training phase (see Fig. 4). In TensorFlow [30] this can be accomplished defining $\lambda$ and $\mu$ as Constant (PDE solution) or Variable (parameter identification) objects, respectively. We set up the problem using the SciANN [40] framework, a high-level Keras [37] wrapper for physics-informed deep learning and scientific computations. Experimenting with all of the previously mentioned network choices can be easily done in SciANN with minimal coding.¹¹1The code for some of the examples solved here is available at: https://github.com/sciann/examples.

Here, we use PINN to identify the model parameters $\lambda$ and $\mu$. Our data corresponds to the exact solution with parameter values $\lambda=1$, $\mu=0.5$ and $Q=4$. Our default dataset consists of 100$\times$100 sample points, uniformly distributed. We study how the accuracy and the efficiency of the identification process depend on the architecture and functional form of the network; the available data; and whether we use one or several independent networks for the different quantities of interest. To study the impact of the architecture and functional form of the ANN, we use 4 different networks with either 5 or 10 hidden layers, and either 20 or 50 neurons per layer; see Table 1. The role of the network functional form is studied comparing the performance of the two most widely used activation functions, i.e., $\tanh$ and ReLU, where $\textrm{ReLU}(x)=\max(0,x)$ [1].

Studying the impact of the available data on the identification process is crucial because we are interested in identifying the model parameters with as little data as possible. We undertake the analysis considering two scenarios:

1. (a)

   Stress-complete data: In this case, we have data at a set of points for the displacements and their first-order derivatives, that is, $u_{x}^{*}$, $u_{y}^{*}$, $\sigma_{xx}^{*}$, $\sigma_{yy}^{*}$, $\sigma_{xy}^{*}$. Because our cost function (9) involves also data that depends on the stress derivatives ($f_{x}^{*}$ and $f_{y}^{*}$), this approach relies on an additional algorithmic procedure for differentiation of stresses. In this section we compute the stress derivatives using second-order central finite-difference approximations.

2. (b)

   Force-complete data: In this scenario, we have data at a set of points for the displacements, their first derivatives and their second derivatives. The availability of the displacement second derivatives allows us to determine data for the body forces $f_{x}^{*}$ and $f_{y}^{*}$ using the momentum balance equation without resorting to any differentiation algorithm.

In Fig. 5 we compare the evolution of the cost function for stress-complete data (Fig. 5a) and force-complete data (Fig. 5b). Both figures show a comparison of the four network architectures that we study; see Table 1. We find that training on the force-complete data performs slightly better (lower loss) at a given epoch.

The result of convergence of model identification is shown in Fig. 6. The training converges to the true values of parameters, i.e., $\lambda=1$ and $\mu=1/2$, for all cases. We find that the optimization is very quick on the parameters while it takes far more epochs to fit the network on the field variables. Additionally, we observe that deeper networks produce less accurate parameters. We attribute the loss of accuracy as we increase the ANN complexity to over-fitting [1, 3]. Convergence of the individual terms in the loss function (9) is shown in Fig. 7 for Net-ii (see Table 1). We find that all terms in the loss, i.e., data-driven and physics-informed, show oscillations during the optimization. Therefore, no individual term is solely responsible for the oscillations in the total loss (Fig. 5).

The impact of the ANN functional form can be examined comparing the data in Figs. 5b and 8a, which show the evolution of the cost function using the activation functions $\tanh$ and ReLU, respectively. The function ReLU has discontinuous derivatives, which explains its poor performance for physics-informed deep learning, whose effectiveness relies heavily on accurate evaluation of derivatives.

A comparison of Figs. 5b and 8b shows that using independent networks for displacements and stresses is more effective than using a single network. We find that the single network leads to less accurate elastic parameters because the cross-dependencies of the network outputs through the kinematic and constitutive relations may not be adequately represented by the $\tanh$ activation function.

Fig. 9 analyzes the effect of availability of data on the training. We computed the exact solution on four different uniform grids of size $10\times 10$, $40\times 40$, $160\times 160$, and $640\times 640$; and carried out the parameter identification process. We performed the comparison using force-complete data and a network with 10 layers and 20 neurons per layer (network iii). The training process found good approximations to the parameters for all cases, including that with only $10\times 10$ points. The results show that fewer data points require many more epoch cycles, but the overall computational cost is far lower.

Here, we generate synthetic data from FEM solutions, and then perform the training. The domain is discretized with a mesh comprised of $40\times 40$ elements. Four datasets are prepared using quadrilateral bilinear, biquadratic, bicubic, and biquartic Lagrange $C^{0}$ elements using the commercial FEM software COMSOL [41]. We evaluate the FEM displacements, strains, stresses and stress derivatives at the center of each element. Then, we map the data to a $100\times 100$ training grid using SciPy’s griddata module with cubic interpolation. This step is performed as a data-augmentation procedure, which is a common practice in machine learning [1].

To analyze the importance of data satisfying the governing equations of the system, we focus our attention on network ii and we study cases with stress- and force-complete data. The results of training are presented in Fig. 10. As can be seen here, the bilinear element performs poorly on the learning and identification. The performance of training on the other elements is good, comparable to that using the analytical solution. Further analysis shows that this is indeed expected as FEM differentiation of bilinear elements provides a poor approximation of the body forces. The error in the body forces is shown in Fig. 11, which indicates a high error for bilinear elements. We conclude that the standard bilinear elements are not suitable for this problem to generate numerical data for deep learning. Fig. 10a2 confirms that pre-processing the data can remove the error that was present in the numerical solution with bilinear elements, and enable the optimization to successfully complete the identification.

Observing the lowest loss $\mathcal{L}$ on the analytical solution, we decided to study the influence of the global continuity of the numerical solution. We generated a $C^{3}$-continuous dataset using Isogeometric analysis [42]. We, therefore, analyze the system using $C^{3}$ IGA elements with again a grid of $40\times 40$ dimension. The data are then mapped on to a grid of $100\times 100$ and used to train the PINN models. The training results are shown in Fig. 12. The outputs are very similar to the high-order FEM datasets.

Here we explore the applicability of our PINN framework to transfer learning: a neural network that is pre-trained is used to perform identification on a new dataset. The expectation is that since the initial state of neural network is not randomly chosen anymore, training should converge faster to the solution and parameters of the data. This is crucial for many practical aspects including adaptation to new data for online search or purchase history [34] or in geosciences, where we can train a representative PINN in highly-instrumented regions and use them at other locations with limited observational datasets. To this end, we use the pre-trained model on Net-iii (Fig. 5), which was trained on a dataset with $\lambda=1.0$ and $\mu=0.5$ and then we explore how the loss evolves and the training converges when data is generated with different values of $\mu\in\left\{2.0,1.5,1.0,0.1\right\}$.

In Fig. 13 we show the convergence of the model with different datasets. Note that the loss is normalized by the initial value $\mathcal{L}_{0}$ from the pre-trained network on $\mu=0.5$ (Fig. 5). As can be seen here, re-training on new datasets costs only a few hundred epochs with a smaller initial value for the loss. This is pointing to the advantage of deep learning and PINN, where retraining on similar data is much less costly than classical methods that rely on forward simulations.

Performing sensitivity analysis is an expensive task when the analytical solution is not available, since it requires performing many forward numerical simulations. Alternatively, if we can construct a surrogate model to be a function of parameters of interest, then performing sensitivity analysis becomes tractable. However, construction of such a surrogate model is itself an expensive task within classical frameworks. Within PINN, however, this seems to be naturally possible. Let us suppose that the parameter of interest is shear modulus $\mu$. Consider an ANN model with inputs as $(x,y,\mu)$ and outputs as $(u_{x},u_{y},\sigma_{xx},\sigma_{yy},\sigma_{xy})$. We can, therefore, use a similar framework to construct a model that is a function of $\mu$ in addition to the space variables. Again, PINN can constrain the model to adapt to the physics of interest and therefore there is less data needed to construct such a model.

Here, we explore if a PINN model trained on multiple datasets generated with various material parameters, i.e., different values of $\mu$, can be used as a surrogate model to perform sensitivity analysis. The network in Fig. 4 is now slightly adapted to carry $\mu$ as an extra input (in addition to $x,y$). The training set is prepared based on $\lambda=1$ and $\mu\in\left\{1/4,2/3,3/2,4\right\}$. Note that there is no identification in this case, and therefore the parameters $\lambda$ and $\mu$ are known at any given training data. The results of the analysis are shown in Fig. 14. For a wide range of values of $\mu\in(0,9)$, the model performs very well in terms of displacements; it is less accurate, but still very useful, in terms of stresses with a maximum error for near-incompressible conditions, $\mu\approx 0$.

We adopt the classic elastoplasticity postulate of additive decomposition of the strain tensor [43],

The stress tensor is now linearly dependent on the elastic strain tensor:

The plastic part of deformation tensor is evaluated through a plasticity model. The von Mises model implies that the plastic deformation occurs in the direction of normal to a yield surface $\mathcal{F}$ defined as $\mathcal{F}(\sigma_{ij}):=q-\sigma_{Y}$, as

where $\sigma_{Y}$ is the yield stress, $q$ is the equivalent stress defined as $q=\sqrt{3/2s_{ij}s_{ij}}$, with $s_{ij}$ the components of the deviatoric stress tensor, $s_{ij}=\sigma_{ij}-\sigma_{kk}/3\delta_{ij}$. The strain remains strictly elastic as long as the state of stress $\sigma_{ij}$ remains inside the yield surface, $\mathcal{F}<0$. Plastic deformation occurs when the state of stress is on the yield surface, $\mathcal{F}=0$. The condition $\mathcal{F}>0$ is associated with an inadmissible state of stress. Parameter $\gamma$ is the plastic multiplier, subject to the condition $\gamma\geq 0$, and evaluated through a predictor–corrector algorithm by imposing the condition $\mathcal{F}\leq 0$ [43]. In the case of von Mises plasticity, the volumetric plastic deformation is zero, $\varepsilon_{kk}^{p}=0$. It can be shown that the plastic multiplier $\gamma$ is equal to the equivalent plastic strain $\bar{\varepsilon}^{p}=\sqrt{2/3e_{ij}^{p}e_{ij}^{p}}$, where $e_{ij}$ are the components of deviatoric strain tensor, $e_{ij}=\varepsilon_{ij}-\varepsilon_{kk}/3\delta_{ij}$.

Therefore, the elastoplastic relations for a plane-strain problem can be summarized as:

subject to the elastoplasticity conditions:

also known as Karush–Kuhn–Tucker (KKT) conditions. For the von Mises model, the plastic multiplier $\bar{\varepsilon}^{p}$ can be expressed as

where $\bar{\varepsilon}$ is the total equivalent strain, i.e., $\bar{\varepsilon}=\sqrt{2/3e_{ij}e_{ij}}$. Therefore, the parameters of this model are the Lamé elastic parameters $\lambda$ and $\mu$, and the yield stress $\sigma_{Y}$.

The solution variables for a two-dimensional problem are $u_{x}$, $u_{y}$, $\varepsilon_{xx}$, $\varepsilon_{yy}$, $\varepsilon_{xy}$, $\varepsilon_{xx}^{p}$, $\varepsilon_{yy}^{p}$, $\varepsilon_{zz}^{p}$, $\varepsilon_{xy}^{p}$, $\sigma_{xx}$, $\sigma_{yy}$, $\sigma_{zz}$, $\sigma_{xy}$. Since the out-of-plane components are no longer zero, they must be reflected in the choice of independent networks. Following the discussions for the linear elasticity case, we approximate the displacement and stress components $u_{x},u_{y},\sigma_{xx},\sigma_{yy},\sigma_{zz},\sigma_{xy}$ with nonlinear neural networks as:

The associated cost function is then defined as

The KKT positivity and negativity conditions are imposed through a penalty constraint in the loss function. For instance, $\bar{\varepsilon}^{p}\geq 0$ is incorporated in the loss as $(1-\textrm{sign}(\bar{\varepsilon}^{p}))|\bar{\varepsilon}^{p}|$. Therefore, for values of $\bar{\varepsilon}^{p}<0$, the resulting ‘cost’ is $2|\bar{\varepsilon}^{p}|$, which should vanish.

We use a classic example to illustrate our framework: a perforated strip subjected to uniaxial extension [44, 43]. Consider a plate of dimensions $200~{}\text{mm}\times 360~{}\text{mm}$, with a circular hole of diameter $100~{}\text{mm}$ located in the center of the plate. The plate is subjected to extension displacements of $\delta=1~{}\text{mm}$ along the short edge, under plane-strain condition, and without body forces, $f_{i}=0$. The parameters are $\lambda=19.44~{}\textrm{GPa}$, $\mu=29.17~{}\textrm{GPa}$ and $\sigma_{Y}=243.0~{}\textrm{MPa}$. Due to symmetry, only a quarter of the domain needs to be considered in the simulation. The synthetic data is generated from a high-fidelity FEM simulation using COMSOL software [41] on a mesh of 13041 quartic triangular elements (Fig. 15). The plate undergoes significant plastic deformation around the circular hole, as can be seen from $\varepsilon_{ij}^{p}$ contours in Fig. 15. This results in localized deformation in the form of a shear band. While the strain exhibits localization, the stress field remains continuous and smooth—a behavior that is due to the choice of a perfect-plasticity model with no hardening.

We use 2,000 data points from this reference solution, randomly distributed in the simulation domain, to provide the training data. The PINN training is performed using networks with 4 layers, each with 100 neurons, and with a hyperbolic-tangent activation function. The optimization parameters are the same as those used for the linear elasticity problem. The results predicted by the PINN approach match the reference results very closely, as evidenced by: (1) the very small errors in each of the components of the solution, except for the out-of-plane plastic strain components (Fig. 16); and (2) the precise identification of yield stress $\sigma_{Y}$ and relatively accurate identification of elastic parameters $\lambda$ and $\mu$, yielding estimated values $\lambda=18.3~{}\textrm{GPa}$, $\mu=27.6~{}\textrm{GPa}$ and $\sigma_{Y}=243.0~{}\textrm{MPa}$.