A nonlocal physics-informed deep learning framework using the peridynamic differential operator

In the PINN framework [24], the solution space is constructed by a deep neural network with the independent variables (e.g., coordinates $\mathbf{x}$) as the network inputs. In this feed-forward network, each layer outputs data as inputs for the next layer through nested transformations. Corresponding to the vector of input variables, $\mathbf{x}$, the output values, $f(\mathbf{x})$ can be mathematically expressed as

where $\mathbf{z}^{0}=\mathbf{x}$, ${z}^{L}={f}(\mathbf{x})$ and $\mathbf{z}^{l}$ represent the inputs, final outputs and hidden layer outputs of the network, and $\mathbf{W}^{\ell}$ and $\mathbf{b}^{\ell}$ represent the weights and biases of each layer, respectively. Note that lowercase and capital boldface letters are used to reflect vector and matrix components while scalars are shown with italic fonts. The activation function is denoted by actf; it renders the network nonlinear with respect to the inputs.

The ‘trained’ ${f}(\mathbf{x})$ can be considered as an approximate solution to the governing PDE. It defines a mapping from inputs $\mathbf{x}$ to the field variable ${f}$ in the form of a multi-layer deep neural network, i.e., ${f}:~{}\mathbf{x}\mapsto\mathcal{N}(\mathbf{x};\mathbf{W},\mathbf{b})$, with $\mathbf{W}$ and $\mathbf{b}$ representing the set of all network parameters. The network inputs $\mathbf{x}$ can be temporal and spatial variables in reference to a Cartesian coordinate system, i.e., $\mathbf{x}=(x,y,t)$ in 2D.

In the PINN framework, the physics, described by a partial differential equation $\mathcal{P}({f})$ with $\mathcal{P}$ as the partial differential operator, is incorporated in the loss or cost function $\mathcal{L}$ along with the training data as

where ${f}^{*}$ is the training dataset (which can be inside the domain or on the boundary), and $0^{*}$ represents the expected (true) value for the differential operation $\mathcal{P}({f})$ at any given training or sampling point. In all modern implementations of the deep-learning framework, such as Theano [17], Tensorflow [18] and MXNet [38], the partial derivatives in $\mathcal{P}$ can be performed using automatic differentiation (AD) [27]—a fundamental aspect of the PINN architecture.

Different optimization algorithms exist for training a neural network; these include Adagrad [39] and Adam [40]. Several algorithmic parameters affect the rate of convergence of network training. 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 is 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) results in an updated parameter set because the batch-gradients are evaluated on different batches. It is common to reshuffle 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 loss function. This situation is described with the keyword patience. It primarily occurs because the loss function is nonconvex. Therefore, the training needs to be tested with different starting points and in different directions to build confidence on the parameters evaluated from minimization of the loss function on a given dataset. Patience is the parameter that controls when the optimizer should stop the training. There are three basic strategies 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. In the original work on PINN [24], the first strategy was used to train the model, with datasets being generated at random space locations at each epoch. This strategy, however, is often not applicable in real-world applications, especially when measurements are collected from sensors that are installed at fixed and limited spatial locations.

In this paper, we rely on SciANN [20], a recent implementation of PINN as a high-level Keras [41] wrapper for physics-informed deep learning and scientific computations. Experimenting with all of the previously mentioned network choices can be easily done, with minimal coding, in SciANN [20, 14].

In the absence of body forces and neglecting inertia, the balance of linear momentum takes the form:

for $i,j=x,y,z$, where $\sigma_{ij}$ is the Cauchy stress tensor, the subscript after a comma denotes differentiation, and a repeated subscript implies summation.

The linear elastic response of the material can be described by the stress–strain relations as

where the pressure or volumetric stress is

and the deviatoric stress tensor is

in which

with the strain tensor defined as

where $u_{i}$ are the components of the displacement field.

The nonlinear material response follows the classical description of elastoplastic behavior [42]. In particular, we adopt the von Mises flow theory with a yield surface $\mathcal{F}$ defined as

in which $\sigma_{Y0}$ is the initial yield stress, $H_{p}$ is the work hardening parameter, $\bar{e}^{p}$ is the equivalent plastic strain, and $\sigma_{e}$ is the effective stress. The plastic deformation occurs in the direction normal to the yield surface, i.e., $n_{ij}=\partial\mathcal{F}/\partial\sigma_{ij}$. We decompose the deviatoric strain tensor into its elastic and plastic components,

To account for plastic deformation, the equations describing the linear elastic material response, Eq. (2.6) can be rewritten as

and

The effective stress $\sigma_{e}$ is defined as

with

The equivalent plastic strain, $e_{ij}^{p}$ can be obtained from Eq. (2.9), Eq. (2.11) and Eq. (2.12) as

where

For linear elasticity under plane-strain conditions, the transverse component of strain, $\varepsilon_{zz}$, is identically equal to zero, and the transverse normal component of stress is evaluated as $\sigma_{zz}=\nu(\sigma_{xx}+\sigma_{yy})$. For elastoplasticity, however, $\sigma_{zz}$ is not predefined while $\varepsilon_{zz}$ remains identically equal to zero.

Here we apply the PINN framework to the solution and inference of two-dimensional quasi-static mechanics. The input variables to the feed-forward neural network are the coordinates, $x$ and $y$, and the output variables are the components of the displacement, $u_{x}$, $u_{y}$, strain tensor, $\varepsilon_{xx}$, $\varepsilon_{yy}$, $\varepsilon_{xy}$, and stress tensor, $\sigma_{xx}$, $\sigma_{yy}$, $\sigma_{xy}$. We define the loss function for linear elasticity as:

where the $\circ$ and $\circ^{*}$ components refer to predicted and true values, respectively. The set $I$ contains all sampling nodes. The set $I_{\square}$ contains all sampling nodes for variable $\square$ where actual data exist. The terms in the loss function represent measures of the error in the displacement, strain and stress fields, the equilibrium equations, and the constitutive relations.

Similarly, the loss function for elastoplasticity is:

These loss functions are used for deep-learning-based solution of the governing PDEs as well as for identification of the model parameters. The constitutive relations and governing equations are tested at all sampling (collocation) points, while data can be selectively imposed. The material parameters are treated as constant values in the network for the solution of governing PDEs. However, they are treated as network parameters, which change during the training phase, during model identification (see Fig. 1). TensorFlow [18] permits such variables to be defined as Constant (PDE solution) or Variable (parameter identification) objects, respectively.

We simulate the deformation of a square domain under plane-strain conditions, as a result of the indentation by a rigid punch (Fig. 4). The body is constrained by roller support conditions along the lateral boundaries, and subject to fixed zero displacement along the bottom boundary. The dimensions of the domain are $W=L=1$ m, and thickness $h=1$ m. The width of the rigid punch is $a=0.2$ m, which indents the body at the top boundary a prescribed vertical displacement $\Delta=1$ mm. These boundary conditions can be expressed as

The material exhibits elastic or elastic-plastic deformation with strain hardening. The elastic modulus, Poisson’s ratio, yield stress and hardening parameter of the material are specified as $E=70$ GPa, $\nu=0.3$, $\sigma_{Y0}=0.1$ GPa and $H_{0}=0.5$ GPa, respectively. The Lamé elastic constants, therefore, have values $\lambda=40.385$ GPa and $\mu=26.923$ GPa.

To generate synthetic (simulated) data to be used in the deep learning frameworks, we simulate the problem described above with the finite element method using COMSOL [43]. The domain is discretized with a uniform mesh of size $100\times 100$ elements of quartic Lagrange polynomials.

The simulated displacement ($u_{x}$, $u_{y}$), strain ($\varepsilon_{xx}$, $\varepsilon_{yy}$, $\varepsilon_{zz}$, $\varepsilon_{xy}$) and stress ($\sigma_{xx}$, $\sigma_{yy}$, $\sigma_{zz}$, $\sigma_{xy}$) is computed for a purely linear elastic response (Fig. 5) and for elastic-plastic deformation (Fig. 6). It is apparent that the distribution of strain and stress components for the elastoplastic case are significantly different from those of the elastic case, with more localized deformation underneath the rigid punch. As expected, the plastic-strain components are zero in most of the domain, except in the vicinity of the corners of the punch, where it exhibits sharp gradients—a feature that, as we will see, poses a challenge for the approximation of the solution with a neural network.

We first apply the established (local) PINN framework for solution and parameter identification of the indentation problem described above [24, 14]. Training of the neural networks is performed with 10,000 training points (nodal solutions of the FEM solution). The convergence of the network training is sensitive to the choice of data-normalization and network size. After a trial-and-error approach, we selected the network architectures and parameters that led to the lowest value of the loss function and the highest accuracy of the physical model parameters. The selected feed-forward neural network has 4 hidden layers, each with 100 neuron units, and employs the hyperbolic-tangent activation function between layers. We adopt batch-training with a total number of 20,000 epochs and a batch-size of 64. We use the Adam optimizer with a learning rate initialized to 0.0005 and decreased gradually to 0.000001 for the last epoch.

The local PINN predictions for elastic deformation do capture the high-gradient regions near the corners of punch, but they are significantly diffused. The differences between the local PINN predictions and the true data are shown in Fig. 7. The Lamé coefficients identified by the network are $\lambda=40.2$ GPa and $\mu=26.2$ GPa—an error of less than 3% (Fig. 8).

In the case of elastic-plastic deformation, $\sigma_{zz}$ depends on the plastic-strain components. Thus, the PINN architecture is defined with networks for $u_{x}$, $u_{y}$, $\sigma_{xx}$, $\sigma_{yy}$, $\sigma_{xy}$ and $\sigma_{zz}$. The error between the local PINN predictions and the true data are shown in Fig. 9. In contrast with the local PINN predictions for the elastic case, the predictions for this case show poor quantitative agreement with the exact solution. The material parameters identified by the method are: $\lambda=40.4$ GPa, $\mu=26.4$ GPa, $\sigma_{Y0}=0.0992$ GPa and $H_{p}=0.00$ GPa. While the elastic Lamé coefficients and the yield stress are identified accurately, the method fails to identify the hardening parameter $H_{p}$ (Fig. 10). We speculate that this is due to the localized plastic deformation in narrow regions in the vicinity of the corners of the rigid punch (Fig. 6). Therefore, there are very few sampling points that contribute to the loss function with the local PINN network.

Given the relative success of the local PINN framework, but also the challenges faced in capturing the sharp gradients in the solution of the elastoplastic deformation problem, here we investigate the application of nonlocal PINN to this problem.

The selected feed-forward neural networks are identical to those used in local PINN: they have 4 hidden layers, each with 100 neuron units, and with hyperbolic-tangent activation functions. ]In the construction of the PD functions, $g_{2}^{p_{1}p_{2}}(\boldsymbol{\xi})$, the TSE is truncated after the second-order derivatives ($N=2$, see  A). The number of family members for each point depends on the order of approximation in the TSE; it is $(N+1)$ points in each dimension, resulting in $(2N+3)\times(2N+3)$ for a square horizon in 2D [31]. Therefore, we choose a maximum number of 49 members as the nonlocal input features. Depending on the location of $\mathbf{x}$, the influence (degree of interaction) of some of these points (family members) might be zero. However, they are all incorporated in the construction of the nonlocal neural network to simplify the implementation procedure.

In what follows we present the results for both AD-PDDO-PINN and PDDO-PINN architectures to the indentation problem with elastic-plastic deformation.