On the stability of the generalized, finite deformation correspondence model of peridynamics

TL;DR

This paper critically examines the stability of a recent finite deformation peridynamic model, revealing potential instabilities and errors, and discusses stabilization methods for practical use.

Contribution

It provides a detailed stability analysis of the new finite deformation correspondence model, challenging previous claims of improved stability and highlighting issues with homogeneous deformation behavior.

Findings

Homogeneous deformations exhibit unstable behavior at the continuum level.

Discretization effects can induce additional instability.

Numerical examples show large errors due to instability.

Abstract

A class of peridynamic material models known as constitutive correspondence models provide a bridge between classical continuum mechanics and peridynamics. These models are useful because they allow well-established local constitutive theories to be used within the nonlocal framework of peridynamics. A recent finite deformation correspondence theory (Foster and Xu, 2018) was developed and reported to improve stability properties of the original correspondence model (Silling et al., 2007). This paper presents a stability analysis that indicates the reported advantages of the new theory were overestimated. Homogeneous deformations are analyzed and shown to exibit unstable material behavior at the continuum level. Additionally, the effects of a particle discretization on the stability of the model are reported. Numerical examples demonstrate the large errors induced by the unstable behavior. Stabilization strategies and practical applications of the new finite deformation model are discussed.

Figures

Equations

$$\mathcal{H}(\mathbf{X})=\left\{\mathbf{X}+\mathbf{\xi }|\mathbf{X}+\mathbf{\xi }\in {\mathcal{B}}_0,0<|\mathbf{\xi }|\le \delta \right\},$$

$${\rho }_0\ddot{\mathbf{u}}(\mathbf{X},t)={\int }_{\mathcal{H}(\mathbf{X})}\left[\underline{\mathbf{T}}(\mathbf{X},t)⟨\mathbf{\xi }⟩-\underline{\mathbf{T}}(\mathbf{X}+\mathbf{\xi },t)⟨-\mathbf{\xi }⟩\right]\mathrm{d}\mathbf{\xi }+{\rho }_0\mathbf{b}(\mathbf{X},t),$$

$$\bar{\mathbf{F}}=\left({\int }_{\mathcal{H}}\underline{\omega }⟨\mathbf{\xi }⟩\underline{\mathbf{Y}}⟨\mathbf{\xi }⟩\otimes \mathbf{\xi }\mathrm{d}\mathbf{\xi }\right){\mathbf{K}}^{-1},$$

$$\underline{\mathbf{Y}}(\mathbf{X},t)⟨\mathbf{\xi }⟩=\mathbf{x}(\mathbf{X}+\mathbf{\xi },t)-\mathbf{x}(\mathbf{X},t);$$

$$\mathbf{K}={\int }_{\mathcal{H}}\underline{\omega }⟨\mathbf{\xi }⟩\mathbf{\xi }\otimes \mathbf{\xi }\mathrm{d}\mathbf{\xi }.$$

$${\underline{T}}_k=\underline{\omega }⟨\mathbf{\xi }⟩{\sigma }_{kp}{K}_{pq}^{-1}{\xi }_q.$$

$$\bar{\mathbf{C}}_{(m)}=\left[\int_{\mathcal{H}}\underline{\mathbf{\omega}}\langle{\boldsymbol{\xi}}\rangle\,\Bigg{(}\frac{\left|\underline{\mathbf{Y}}\langle{\boldsymbol{\xi}}\rangle\right|}{\left|\boldsymbol{\xi}\right|}\Bigg{)}^{2m}\frac{\boldsymbol{\xi}\otimes\boldsymbol{\xi}}{\left|\boldsymbol{\xi}\right|^{2}}\,{\rm d}\boldsymbol{\xi}\right]:\mathbb{L}^{-1},\quad m\neq 0,$$

$$\mathbb{L}={\int }_{\mathcal{H}}\underline{\omega }⟨\mathbf{\xi }⟩\frac{\mathbf{\xi }\otimes \mathbf{\xi }\otimes \mathbf{\xi }\otimes \mathbf{\xi }}{{\left|\mathbf{\xi }\right|}^4}\mathrm{d}\mathbf{\xi }.$$

$${\bar{\mathbf{E}}}_{(m)}=\frac{1}{2m}\left({\bar{\mathbf{C}}}_m-\mathbf{I}\right),m\ne 0,$$

$$\bar{\mathbf{E}}_{(0)}=\left[\int_{\mathcal{H}}\underline{\mathbf{\omega}}\langle{\boldsymbol{\xi}}\rangle\ln\Bigg{(}\frac{\left|\underline{\mathbf{Y}}\langle{\boldsymbol{\xi}}\rangle\right|}{\left|\boldsymbol{\xi}\right|}\Bigg{)}\frac{\boldsymbol{\xi}\otimes\boldsymbol{\xi}}{\left|\boldsymbol{\xi}\right|^{2}}\,{\rm d}\boldsymbol{\xi}\right]:\mathbb{L}^{-1},$$

$${\underline{T}}_k=\underline{\omega }⟨\mathbf{\xi }⟩S_{(m)ij}\frac{{\xi }_p{\xi }_q}{|\mathbf{\xi }{|}^2}{L}_{pqij}^{-1}{\left(\frac{|\underline{\mathbf{Y}}|}{|\mathbf{\xi }|}\right)}^{2m}\frac{Y_k}{|\underline{\mathbf{Y}}{|}^2},$$

$${\int }_{\mathcal{H}(\mathbf{X})}\left[\underline{\mathbf{T}}(\mathbf{X})⟨\mathbf{\xi }⟩-\underline{\mathbf{T}}(\mathbf{X}+\mathbf{\xi })⟨-\mathbf{\xi }⟩\right]\mathrm{d}\mathbf{\xi }=0.$$

$$\underline{\mathbf{Y}}=\mathbf{F}\cdot \mathbf{\xi }.$$

$${\underline{T}}_k=\underline{\omega }⟨\mathbf{\xi }⟩S_{(m)ij}{L}_{pqij}^{-1}{F}_{kl}{\left(\frac{|\mathbf{F}\cdot \mathbf{\xi }|}{|\mathbf{\xi }|}\right)}^{2m-2}\frac{{\xi }_p{\xi }_q{\xi }_l}{|\mathbf{\xi }{|}^4}.$$

$$\underline{\mathbf{T}}(\mathbf{X})⟨\mathbf{\xi }⟩=-\underline{\mathbf{T}}(\mathbf{X}+\mathbf{\xi })⟨-\mathbf{\xi }⟩.$$

$${\int }_{\mathcal{H}(\mathbf{X})}\left[\underline{\mathbf{T}}(\mathbf{X})⟨\mathbf{\xi }⟩-\underline{\mathbf{T}}(\mathbf{X}+\mathbf{\xi })⟨\mathbf{\xi }⟩\right]\mathrm{d}\mathbf{\xi }=2{\int }_{\mathcal{H}(\mathbf{X})}\underline{\mathbf{T}}(\mathbf{X})⟨\mathbf{\xi }⟩\mathrm{d}\mathbf{\xi }.$$

$$\mathbf{f}(\mathbf{X})={\int }_{\mathcal{H}(\mathbf{X})}\left[\underline{\mathbf{T}}(\mathbf{X})⟨\mathbf{\xi }⟩-\underline{\mathbf{T}}(\mathbf{X}+\mathbf{\xi })⟨-\mathbf{\xi }⟩\right]\mathrm{d}\mathbf{\xi }.$$

$${\underline{J}}_{kl}(\mathbf{X})⟨\mathbf{\eta }⟩=\frac{\partial f_k(\mathbf{X})}{\partial x_l(\mathbf{X}+\mathbf{\eta })},$$

$$\frac{\partial f_k(\mathbf{X})}{\partial x_l(\mathbf{X}+\mathbf{\eta })}={\int }_{\mathcal{H}(\mathbf{X})}\left(\frac{\partial {\underline{T}}_k(\mathbf{X})⟨\mathbf{\xi }⟩}{\partial x_l(\mathbf{X}+\mathbf{\eta })}-\frac{\partial {\underline{T}}_k(\mathbf{X}+\mathbf{\xi })⟨-\mathbf{\xi }⟩}{\partial x_l(\mathbf{X}+\mathbf{\eta })}\right)\mathrm{d}\mathbf{\xi }.$$

$${\underline{T}}_k(\mathbf{X})⟨\mathbf{\xi }⟩=\underline{\omega }⟨\mathbf{\xi }⟩\frac{{\xi }_p{\xi }_q}{|\mathbf{\xi }{|}^{2m+2}}{L}_{pqij}^{-1}{S}_{(m)ij}(\mathbf{X})|{\underline{Y}}_n⟨\mathbf{\xi }⟩{\underline{Y}}_n⟨\mathbf{\xi }⟩{|}^{m-1}{\underline{Y}}_k⟨\mathbf{\xi }⟩,$$

$$\displaystyle\frac{\partial\underline{T}_{k}(\mathbf{X})\langle{\boldsymbol{\xi}}\rangle}{\partial x_{l}(\mathbf{X}+\boldsymbol{\eta})}=\underline{\omega}\langle{\boldsymbol{\xi}}\rangle\frac{\xi_{p}\,\xi_{q}}{|\boldsymbol{\xi}|^{2m+2}}L^{-1}_{pqij}\,\Bigg{[}$$

$$\displaystyle+S_{(m)ij}(\mathbf{X})\left|\underline{\mathbf{Y}}\langle{\boldsymbol{\xi}}\rangle\right|^{2m-2}\frac{\partial\underline{Y}_{k}\langle{\boldsymbol{\xi}}\rangle}{\partial x_{l}(\mathbf{X}+\boldsymbol{\eta})}\Bigg{]}.$$

$$\frac{\partial {\underline{Y}}_k⟨\mathbf{\xi }⟩}{\partial x_l(\mathbf{X}+\mathbf{\eta })}$$

$$=\underline{\delta }⟨\mathbf{\xi },\mathbf{\eta }⟩{\delta }_{kl}-\underline{\delta }⟨0,\mathbf{\eta }⟩{\delta }_{kl},$$

$$\underline{\delta }⟨\mathbf{\xi },\mathbf{\eta }⟩=\{\begin{array}{ll}1& \text{mbox}if\mathbf{\xi }=\mathbf{\eta },\\ 0& \text{mbox}otherwise.\end{array}$$

$$\frac{\partial S_{(m)ij}(\mathbf{X})}{\partial x_l(\mathbf{X}+\mathbf{\eta })}=\frac{\partial S_{(m)ij}}{\partial {\bar{E}}_{(m)rs}}\frac{\partial {\bar{E}}_{(m)rs}(\mathbf{X})}{\partial x_l(\mathbf{X}+\mathbf{\eta })}.$$

$$\frac{\partial {\bar{C}}_{(m)rs}(\mathbf{X})}{\partial x_l(\mathbf{X}+\mathbf{\eta })}=L_{turs}^{-1}{\int }_{\mathcal{H}(\mathbf{X})}\underline{\omega }⟨\mathbf{\xi }⟩m\left(2{\underline{Y}}_n⟨\mathbf{\xi }⟩\frac{\partial {\underline{Y}}_n⟨\mathbf{\xi }⟩}{\partial x_l(\mathbf{X}+\mathbf{\eta })}\right)|\underline{\mathbf{Y}}⟨\mathbf{\xi }⟩{|}^{2m-2}\frac{{\xi }_t{\xi }_u}{|\mathbf{\xi }{|}^{2m+2}}\mathrm{d}\mathbf{\xi },$$

$$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ =L^{-1}_{turs}\int_{\mathcal{H}(\mathbf{X})}2m\,\underline{\omega}\langle{\boldsymbol{\xi}}\rangle\,\bigg{(}\underline{Y}_{n}\langle{\boldsymbol{\xi}}\rangle\left(\underline{\delta}\langle{\boldsymbol{\xi},\boldsymbol{\eta}}\rangle-\underline{\delta}\langle{\mathbf{0},\boldsymbol{\eta}}\rangle\right)\delta_{nl}\bigg{)}\,|\underline{\mathbf{Y}}\langle{\boldsymbol{\xi}}\rangle|^{2m-2}\frac{\xi_{t}\,\xi_{u}}{|\boldsymbol{\xi}|^{2m+2}}\,{\rm d}\boldsymbol{\xi},$$

$$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ =2m\,L^{-1}_{turs}\Bigg{(}\underline{\omega}\langle{\boldsymbol{\eta}}\rangle\,\underline{Y}_{l}\langle{\boldsymbol{\eta}}\rangle\,|\underline{\mathbf{Y}}\langle{\boldsymbol{\eta}}\rangle|^{2m-2}\frac{\eta_{t}\,\eta_{u}}{|\boldsymbol{\eta}|^{2m+2}}\,{\rm d}\boldsymbol{\eta}-\underline{\delta}\langle{\mathbf{0},\boldsymbol{\eta}}\rangle\int_{\mathcal{H}(\mathbf{X})}\underline{\omega}\langle{\boldsymbol{\xi}}\rangle\,\underline{Y}_{l}\,\langle{\boldsymbol{\xi}}\rangle\,|\underline{\mathbf{Y}}\langle{\boldsymbol{\xi}}\rangle|^{2m-2}\frac{\xi_{t}\,\xi_{u}}{|\boldsymbol{\xi}|^{2m+2}}\,{\rm d}\boldsymbol{\xi}\Bigg{)}.$$

$$\frac{\partial {\bar{E}}_{(m)rs}(\mathbf{X})}{\partial x_l(\mathbf{X}+\mathbf{\eta })}$$

$$\displaystyle=L^{-1}_{turs}\Bigg{(}\underline{\omega}\langle{\boldsymbol{\eta}}\rangle\,\underline{Y}_{l}\langle{\boldsymbol{\eta}}\rangle\,|\underline{\mathbf{Y}}\langle{\boldsymbol{\eta}}\rangle|^{2m-2}\frac{\eta_{t}\,\eta_{u}}{|\boldsymbol{\eta}|^{2m+2}}\,{\rm d}\boldsymbol{\eta}-\underline{\delta}\langle{\mathbf{0},\boldsymbol{\eta}}\rangle\int_{\mathcal{H}(\mathbf{X})}\underline{\omega}\langle{\boldsymbol{\xi}}\rangle\,\underline{Y}_{l}\langle{\boldsymbol{\xi}}\rangle\,|\underline{\mathbf{Y}}\langle{\boldsymbol{\xi}}\rangle|^{2m-2}\frac{\xi_{t}\,\xi_{u}}{|\boldsymbol{\xi}|^{2m+2}}\,{\rm d}\boldsymbol{\xi}\Bigg{)}.$$

$$\frac{\partial {\bar{E}}_{(m)rs}(\mathbf{X})}{\partial x_l(\mathbf{X}+\mathbf{\eta })}=L_{turs}^{-1}\underline{\omega }⟨\mathbf{\eta }⟩{\underline{Y}}_l⟨\mathbf{\eta }⟩|\underline{\mathbf{Y}}⟨\mathbf{\eta }⟩{|}^{2m-2}\frac{{\eta }_t{\eta }_u}{|\mathbf{\eta }{|}^{2m+2}}\mathrm{d}\mathbf{\eta }.$$

$$\displaystyle\frac{\partial\underline{T}_{k}(\mathbf{X})\langle{\boldsymbol{\xi}}\rangle}{\partial x_{l}(\mathbf{X}+\boldsymbol{\eta})}=\ \underline{\mathbf{\omega}}\langle{\boldsymbol{\xi}}\rangle\frac{\xi_{p}\,\xi_{q}}{\left|\boldsymbol{\xi}\right|^{2m+2}}L^{-1}_{pqij}\left|\underline{\mathbf{Y}}\langle{\boldsymbol{\xi}}\rangle\right|^{2m-2}\Bigg{[}$$

Full text

1 Introduction

Peridynamics was introduced by Silling (2000) in order to represent the mechanics of continuous and discontinuous media using a single consistent set of equations. The theory has gained notable attention, due to its capabilities in naturally representing material failure, removing the need for complicated numerical treatments in dealing with discontinuities, i.e. cracks. Several peridynamic material models have been proposed since introduction of the theory (Silling et al., 2007; Foster et al., 2010; Mitchell, 2011; Foster and Xu, 2018) and applied to numerous problems, from micro to macro scales, mostly involving brittle fracture (Gerstle et al., 2007; Askari et al., 2008; Ha and Bobaru, 2010; Bobaru et al., 2018; Jafarzadeh et al., 2018; Behzadinasab et al., 2018). A special class of material models was proposed by Silling et al. (2007), called constitutive correspondence models, which provide a convenient framework for using classical constitutive relations in peridynamic material modeling.

While the correspondence model approach offers convenience, it has shortcomings associated with material instability. Stability issues have been associated with the continuum model and its discretization (Littlewood, 2011; Tupek and Radovitzky, 2014; Silling, 2017). Zero-energy mode oscillations appear in meshless simulations that use correspondence material models and result in severe degradation in the accuracy of displacement fields, particularly in areas of steep gradients (Silling, 2017). These oscillations where originally reported as numerical issues; however, Tupek and Radovitzky (2014) demonstrated some kinematically unphysical deformation modes that are permitted in the mathematical formulation of the model. To address the issue, they introduced a nonlocal family of Seth-Hill nonlinear bond strain measures to create what they called an “extended constitutive correspondence formulation”. Recently, Foster and Xu (2018) generalized Tupek’s model and presented a finite deformation correspondence model which does not suffer from any surface effects – a well-known issue in peridynamics (Mitchell et al., 2015).

Presented here is an analysis of the stability, in the sense of Lyapunov, of the generalized finite deformation model and its meshless discretization. It is shown that the new constitutive correspondence theory exhibits zero-energy mode oscillations as a manifestation of a new type of material instability which occurs even in the presence of uniform deformations. The root cause of the problem is investigated to better understand the role of nonlocality on the Seth-Hill strain measures.

The rest of the paper is outlined as follows: In Section 2 a brief background of the peridynamic theory and an overview of the correspondence material modeling is presented. In Section 3 the stability of equilibrium for homogeneous deformations is analyzed at the continuum level. Section 4 presents the effects of particle discretization on the stability. A set of unstable equilibrium points are analytically found in Sections 3 and 4, limited to pure hydrostatic deformations. Numerical examples are provided in Section 5 to demonstrate the stability of the theory in a more general loading scenario. In Section 6 the connection between the stability method presented in this work and the method of Silling (2017) is discussed along with possible stabilization strategies and the source of instability. Section 7 offers concluding remarks.

2 Background on peridynamics and correspondence modeling

The peridynamic theory recasts the mechanics of solid deforamtions in a nonlocal setting. The mechanics of continuous and discontinuous bodies are treated using the same set of equations, which fundamentally allows for discontinuities in the deformation field. Suppose $\mathbf{x}(\mathbf{X},t)$ is the current position of a material point $\mathbf{X}$ in a deformed peridynamic body $\mathcal{B}$. The motion of $\mathbf{x}$ is influenced by other material points through long-range interactions with a nonlocal family. The family of $\mathbf{X}$ is defined in its undeformed configuration by

[TABLE]

where the vector $\boldsymbol{\xi}$ is called a bond, and $\delta$ is the radius of a spherical neighborhood, called the horizon. Note that the point $\mathbf{X}$ itself is excluded from the family. Fig. 1 shows an example of a peridynamic body.

Peridynamics describes the balance of linear momentum using the following integro-differential equation:

[TABLE]

where $\rho_{0}$ is the reference mass density, $\mathbf{u}$ is the displacement, and $\mathbf{b}$ is a prescribed body force density field. $\underline{\mathbf{T}}$ is a vector field (with units of force per square volume), called the force vector state, which is used to describe the constitutive interactions of material points.

Silling et al. (2007) introduced the idea of constitutive correspondence where a nonlocal analog of a kinematic variable is utilized as a bridge between the classical (i.e. local) and peridynamic theories. Silling et al. (2007) defined a nonlocal deformation gradient tensor by

[TABLE]

where $\underline{\mathbf{\omega}}\langle{\boldsymbol{\xi}}\rangle$ is called the influence state, describing the relative degree of interactions between the neighboring points; $\underline{\mathbf{Y}}\langle{\boldsymbol{\xi}}\rangle$ is the deformed image of the bond in its current configuration,

[TABLE]

and $\mathbf{K}$ is a symmetric second-order shape tensor, given by

[TABLE]

The authors showed that in a continuous media, the peridynamic deformation gradient recovers exactly the local quantity under homogeneous deformations. Using this kinematic variable as a metric and incorporating classical theories by means of the first Piola-Kirchhoff stress tensor $\boldsymbol{\sigma}=\boldsymbol{\sigma}(\bar{\mathbf{F}})$, the force vector state was derived

[TABLE]

It was observed (Tupek and Radovitzky, 2014; Silling, 2017) that the nonlocal analog $\bar{\mathbf{F}}$ allows for unphysical deformation modes such as matter interpenetration and can result in loss of material stability. The instability can effectively destroy the accuracy of displacement fields in areas of large displacement gradients. The peridynamic finite deformation correspondence framework (Tupek and Radovitzky, 2014) and its generalization (Foster and Xu, 2018) were proposed to improve the stability of the original correspondence theory by using a family of nonlinear strain measures. As an intermediate step in determining the force vector-state, the generalized theory (Foster and Xu, 2018) utilizes a family of nonlocal Lagrangian Seth–Hill strains. A set of nonlocal right Cauchy-Green deformation tensors is introduced

[TABLE]

where $\mathbb{L}$ is a symmetric fourth-order shape tensor defined by

[TABLE]

The peridynamic Lagrangian Seth–Hill strains are then

[TABLE]

and

[TABLE]

which is shown to recover exactly their local analogs under uniform deformations. Note that $m=0$ and $m=1$ are two important members of the class, which correspond to the widely used Hencky and Green-Lagrange strain measures. A local hyperelastic material model can then be integrated into peridynamics, where the force vector-state is given as

[TABLE]

where $\mathbf{S}_{(m)}$ is the generalized Kirchhoff stress and is a function of $\bar{\mathbf{E}}_{(m)}$.

3 Stability under homogeneous deformation

In this section the stability of the generalized, finite deformation correspondence theory when dealing with uniform deformations is analyzed in a continuum setting. It is first shown that homogeneous deformations constitute a set of equilibria for a peridynamic material under static conditions and in the absence of any external body forces. The stability of these equilibrium points in the context of the finite deformation model are discussed.

3.1 Static equilibrium

For a static problem in the absence of body forces, the equation of motion, Eq. 1, reduces to

[TABLE]

For a static homogeneous deformation field with $\mathbf{F}(\mathbf{X})=\mathbf{F}$, the deformation vector state for a peridynamic bond can be described as

[TABLE]

Then Eq. 6 simplifies to

[TABLE]

Limiting the influence state to be spherical, i.e. $\underline{\mathbf{\omega}}\langle{\boldsymbol{\xi}}\rangle=\underline{\mathbf{\omega}}_{s}(|\boldsymbol{\xi}|)$, $\underline{\mathbf{T}}\langle{\boldsymbol{\boldsymbol{\xi}}}\rangle$ is clearly an odd function of $\boldsymbol{\xi}$. Therefore, for two neighboring points in the bulk of material and away from the boundaries (having the same shape tensor)

[TABLE]

Using Eqs. 7 and 9

[TABLE]

Since $\underline{\mathbf{T}}\langle{\boldsymbol{\xi}}\rangle$ is an odd function of $\boldsymbol{\xi}$, this integral is zero for a full neighborhood and therefore satisfies Eq. 7 automatically in the bulk of material. With appropriate boundary conditions, which equilibrate the internal forces in the vicinity of the boundaries, static equilibrium can be achieved for any homogeneous deformation. A static uniform deformation cannot be obtained if any body force is present in the bulk of material.

3.2 Stability

Analyzing the Jacobian matrix of a mechanical system is a way to understand its stability at equilibrium. If all the eigenvalues of the Jacobian matrix are negative, the system is stable. Presence of one or more positive eigenvalues indicates unstable behavior. A zero eigenvalue can correspond to either stable or unstable behavior.

Let the net force exerted on $\mathbf{X}$ in Eq. 1 be called $\mathbf{f}(\mathbf{X})$. In statics, it is given by

[TABLE]

Using partial differentiation, the Jacobian matrix is computed as

[TABLE]

or

[TABLE]

Rewriting Eq. 6

[TABLE]

and taking the partial derivative

[TABLE]

The partial derivative of the deformation state is calculated

[TABLE]

where $\underline{\delta}\langle\boldsymbol{\xi},\boldsymbol{\eta}\rangle$ is the Kronecker-delta state, i.e.

[TABLE]

Using $\mathbf{S}_{(m)}=\mathbf{S}(\mathbf{\bar{E}}_{(m)})$ and the chain rule

[TABLE]

Taking the derivative of Eq. 4,

[TABLE]

Using Eqs. 5 and 14,

[TABLE]

It can be shown that this equation is valid for $m=0$. Similar to the argument made in deriving Eq. 8, for a uniform deformation, the integrand is odd, and the integral is equal to zero. Therefore,

[TABLE]

Substitute Eqs. 12, 13 and 15 into Eq. 11 and simplify to get

[TABLE]

Noting that $\underline{\mathbf{Y}}\langle{\boldsymbol{\xi}}\rangle=-\underline{\mathbf{Y}}\langle{-\boldsymbol{\xi}}\rangle$, $\mathbf{S}_{(m)}(\mathbf{X})=\mathbf{S}_{(m)}(\mathbf{X}+\boldsymbol{\xi})=\mathbf{S}_{(m)}$, and $\underline{\mathbf{\omega}}\langle{\boldsymbol{\xi}}\rangle=\underline{\mathbf{\omega}}\langle{-\boldsymbol{\xi}}\rangle$, we use Eq. 16 to derive

[TABLE]

Substitute Eqs. 16 and 17 in Eq. 10 and rearrange terms to get

[TABLE]

Bringing the constants out of the integrations and grouping terms to further simplify

[TABLE]

The first integral is zero for uniform deformations. Using Eq. 5, the above equation is further reduced

[TABLE]

The first term is nonzero if $\mathbf{X}+\boldsymbol{\eta}$ is among the neighbors of neighbors, the second term is nonzero if $\mathbf{X}+\boldsymbol{\eta}$ is a neighbor of $\mathbf{X}$, and the third term is nonzero only when $\boldsymbol{\eta}=\mathbf{0}$. The special case $\boldsymbol{\eta}=\mathbf{0}$ and $k=l$ corresponds to the diagonal values of the Jacobian matrix, which need to be negative as a necessary condition for stability. Focusing on the stability of $\mathbf{X}$ when itself is perturbed.

[TABLE]

In the continuum limit the second term dominates as the first one is of $O({\rm d}\mathbf{X})$. Hence, in a continuous media

[TABLE]

This property depends on $m$, from the finite deformation class, and the loading condition.

Remark 1

Using $\boldsymbol{\eta}=\mathbf{0}$ in Eq. 15, and noting that $\underline{\mathbf{\omega}}\langle{\mathbf{0}}\rangle=0$,

[TABLE]

therefore,

[TABLE]

Foster and Xu (2018, Section 2.3)** claimed that the finite deformation correspondence theory eliminates the center-point instability of the original correspondence model and hence improves it. However, it is observed that this type of instability is still present under homogeneous deformations.

3.2.1 Hydrostatic deformations

Consider the pure hydrostatic loading

[TABLE]

with $a>-1$. Using Eqs. 20 and 19

[TABLE]

Using the following hydrostatic constitutive relation:

[TABLE]

where $\kappa$ is a positive material constant (a multiple of bulk modulus, depending on the problem dimension). Using Eqs. 21 and 22,

[TABLE]

where

[TABLE]

A necessary condition for stability is $\dfrac{\partial f_{k}(\mathbf{X})}{\partial x_{k}(\mathbf{X})}<0$; thus,

[TABLE]

is required.

Proposition 1

Depending on the problem dimension, a critical $m_{\rm cr}$ exists such that $\Gamma(m_{\rm cr})=0$.

Proof. Eq. 24 is simplified to find $m_{\rm cr}$ for 1, 2, and 3 dimensions.

1D: $L^{-1}=\dfrac{1}{V_{\underline{\mathbf{\omega}}}}$, where $V_{\underline{\mathbf{\omega}}}$ is the weighted volume

[TABLE]

Then,

[TABLE]

[TABLE]

2D: $L_{pqii}^{-1}=\dfrac{2}{V_{\underline{\mathbf{\omega}}}}\delta_{pq}$ (see Appendix A),

[TABLE]

$\Gamma(m)$ is then evaluated to be

[TABLE]

[TABLE]

3D: $L_{pqii}^{-1}=\dfrac{3}{V_{\underline{\mathbf{\omega}}}}\delta_{pq}$ (cf. (Foster and Xu, 2018, Section 3.2)), then performing calculations similar to 2D results in

[TABLE]

[TABLE]

Remark 2

It is evident that

[TABLE]

Therefore, using Eqs. 25 and 2, there are 3 possibilities:

$m>m_{\rm cr}$: $a<0$ is not favored. Thus, the system is unstable under pure compression loadings.

• 2.

$m=m_{\rm cr}$: $a\Gamma(m)=0$ for any loading, and the unstable loadings cannot be determined yet.

• 3.

$m<m_{\rm cr}$: If $a>0$, the system shows unstable behavior. In other words, pure expansions are unstable equilibria for the system.

Note that the violation of Eq. 25 is indicative of instability. However, even if Eq. 25 holds, stability is not guaranteed; the full Jacobian matrix should be considered for proving stability.

Remark 3

It is realized that in the continuum, for $m\neq m_{\rm cr}$, either the slightest compression or tension will cause instability. It will be shown in Section 4 that discretization would help the stability and reduce the unstable regions.

3.3 Stability of Silling’s correspondence model

The same procedure is repeated for the original correspondence theory (Silling et al., 2007). Similar to Section 3.1, it can be shown that homogeneous deformations constitute a set of equilibria for the static case, in the absence of body forces. To analyze the stability of such equilibria, the Jacobian matrix elements are computed. Taking derivative of Eq. 3,

[TABLE]

Using $\boldsymbol{\sigma}=\boldsymbol{\sigma}(\bar{\mathbf{F}})$,

[TABLE]

Using Eqs. 2 and 12,

[TABLE]

The last integral is zero since the integrand is odd; therefore,

[TABLE]

Combining Eqs. 10, 26, 27 and 28 to obtain the Jacobian element

[TABLE]

Again, paying special attention to the case $\boldsymbol{\eta}=\mathbf{0}$

[TABLE]

For a typical influence state

[TABLE]

the integral resultant equals the shape tensor, which is

[TABLE]

in 3D. Hence,

[TABLE]

If the property of positive-definiteness in $\dfrac{\partial\sigma_{ij}}{\partial F_{kl}}$ is inherited in the underlying classical model, all the eigenvalues of $\bar{J}_{kl}(\mathbf{X})\langle{\mathbf{0}}\rangle$ would be negative and the model will be stable for any homogeneous deformation.

The main takeaway is that the unstable behavior identified in the generalized, finite deformation model is not present in Silling’s model. There are, however, other types of instabilities associated with Silling’s model, such as the zero-energy mode oscillations (Breitenfeld et al., 2014; Silling, 2017; Chowdhury et al., 2019), which are associated to non-uniform deformations and particle discretization.

4 Effect of discretization on the stability

The effect of particle dicretization on the stability of the model is studied in this section. In the discrete form the first term in Eq. 18 is no longer negligible. A general uniform loading is considered first, then it is specialized to the hydrostatic loading condition. Since $m_{\rm cr}$ is different in 1D and higher dimensions (cf. Section 3.2.1), 1D and 2D are discussed separately. In the subsequent sections, to facilitate the notation, $\mathbf{X}$, $\mathbf{X}+\boldsymbol{\xi}$, and $\boldsymbol{\xi}$ are referred to as I, J, and IJ indices, respectively.

4.1 1D

Consider a discretized form of Eq. 18 in a one-dimensional setup and simplify:

[TABLE]

Adopting a uniform discrete representation with $\Delta V^{J}=A\Delta X$ and $\delta=N\Delta X$ and the specific choice of influence function in Eq. 30, the shape tensor is computed to be $L=2AN\Delta X$. Using a homogeneous deformation with $x=(1+a)X$,

[TABLE]

[TABLE]

Using Eqs. 32, 33 and 34,

[TABLE]

Next, the peridynamic Lagrangian Seth–Hill strain is computed to evaluate the stress state.

[TABLE]

and

[TABLE]

A linear 1D constitutive model is used

[TABLE]

Substituting Eq. 36 in Eq. 35,

[TABLE]

for $m\neq 0$, and

[TABLE]

Stability enforces positivity of the terms inside the brackets. In other words, the model is unstable if the following conditions are violated:

[TABLE]

and

[TABLE]

Based on $m$ and $a$, stable and unstable regions are distinguished in Fig. 2, for a typical discretization with 3 nodes per horizon ($N=3$). As discussed in Section 3.2.1, compression is a problem for $m>1/2$, while instability is present in tensile mode for $m<1/2$.

Remark 4

Increasing the nonlocality $N$ is unfavorable to the stability. Consider $m=1$, for instance, where $a_{\rm cr}=-0.184$ for $N=1$ and $a_{\rm cr}=-0.074$ for $N=3$ (the model is unstable for $a<a_{\rm cr}$). In other words, a greater compression is required to make the former unstable.

4.2 2D

A discretized form of Eq. 18 in a two-dimensional setup is considered here.

[TABLE]

Considering a uniform discretization with $\Delta V^{J}=b(\Delta X)^{2}$ and $\delta=N\Delta X$ and an influence function as in Eq. 30, the inverse of the shape tensor is computed first (see Appendix A). Then, Eq. 37 is simplified for a general 2D uniform deformation with $x_{i}=X_{i}+a_{ij}X_{j}$.

[TABLE]

where $Y_{i}^{IJ}=\xi_{i}^{IJ}+a_{ij}\xi_{j}^{IJ}$. It will be tedious, if possible, to analyze a general deformation with all the 4 $a_{ij}$ constants. Let consider a pure hydrostatic loading condition with $a_{11}=a_{22}=a$ and $a_{12}=a_{21}=0$. Eq. 38 is then simplified to

[TABLE]

Using $\boldsymbol{\xi}^{{\rm I}{\rm J}}=p\Delta X\hat{x}_{1}+q\Delta X\hat{x}_{2}$ and involving symmetry,

[TABLE]

where the last summation is approximated by an integral to obtain a closed-form solution. Then,

[TABLE]

Similarly, the following summations are approximated using integration:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Next, by using symmetry and noting that only the even combinations of $\xi^{{\rm I}{\rm J}}_{i}$ produce nonzero values, the summations in Eq. 39 are simplified.

[TABLE]

[TABLE]

[TABLE]

where $\gamma_{ijklrs}$ is the following:

[TABLE]

Using the above calculations to simplify Eq. 39 (detailed in Appendix B),

[TABLE]

Because the deformations are limited to hydrostatic, a simple constitutive law is adopted as

[TABLE]

where $\kappa$ is twice the bulk modulus in a two-dimensional problem. Note that there is no shear for hydrostatic loadings.

For a pure homogeneous hydrostatic loading with $x=(1+a)X$

[TABLE]

and

[TABLE]

Substitute these back in Eq. 41 and simplify (see Appendix C):

[TABLE]

for $m\neq 0$, and

[TABLE]

The eigenvalues of this diagonal matrix are identified as:

[TABLE]

for $m\neq 0$, and

[TABLE]

An equilibrium point is unstable if it does not satisfy the following criteria:

[TABLE]

[TABLE]

Similar to the 1D case, the stable and unstable loading conditions are recognized for different $m$’s, depicted in Fig. 3

Remark 5

Increasing the dimension expands the unstable regions. For example, consider $m=1$ and $N=3$ in 2D, where $a_{cr}=-0.037$. This critical value is half the 1D case (cf. Remark 4), meaning that a lower compression would cause instability. It can be shown that the region of instability is also larger in 3D than 2D.

Remark 6

Comparing the results of the discrete level from Section 4 with the continuum scale from Section 3.2, the stable regions are evidently larger in the discrete form. In other words, numerical representation helps the Lyapunov stability in the case of homogeneous hydrostatic deformations.

Remark 7

For $m=0$ and a small window of $m<0$ (in the discrete form), the model may be stable for any uniform hydrostatic deformation. However, other types of loading conditions may impose instability, which actually is the case and will be discussed in Section 5.3.

While obtaining a closed-form solution for more complex loading scenarios is not easily achievable, stability of the model under such conditions can be examined through numerics.

5 Numerical examples

In this section, stability of the finite deformation model is numerically examined for different problems. The two special members of the material class, $m=0$ (Hencky strain) and $m=1$ (Green-Lagrange strain) are the primary targets.

5.1 1D singular bar

The problem was designed by Breitenfeld et al. (2014) to show the unstable behavior of the original correspondence theory. It involves a 1D bar characterized with a spatially varying Young’s modulus

[TABLE]

fixed at one end and loaded with a stress $\sigma$ at the other end. The local analytical solution gives the normalized axial displacement

[TABLE]

and the normalized strain

[TABLE]

along the bar.

This static problem is solved using the $m=1$ model for an applied tensile stress (with $\sigma/E_{0}=10^{-3}$) and a smaller compressive stress (with $\sigma/E_{0}=-10^{-5}$). Normalized displacement and strain values are shown in Fig. 4, when $\alpha=10$. While the model correctly predicts the response of the bar with no apparent instability, the unstable behavior in compression mode leads to wrong results. The original correspondence model, in contrast, has a similar unstable behavior in both tensile and compressive modes (Breitenfeld et al., 2014).

5.2 Influence of step size in different loading regimes

The stable load step size can be significantly different between loading modes, even in the stable regions. For example, consider a 1-dimensional hyperelastic bar subjected to prescribed displacements at its ends. The equilibrium state is calculated using the $m=1$ material model with 3 nodes per horizon ($N=3$). Note that $a_{cr}=-7.4\%$ for this problem (cf. Remark 4). The bar is initially undeformed and boundary displacements are applied instantly (in one step) over a boundary region that contains 6 nodes at each end. As shown in Fig. 5, under tensile mode, the theory correctly finds the solution for an extremely large step size. However, the stable load step is several orders of magnitude smaller when compressive loading is applied, even when the compression state is well inside the stable region.

5.3 Stability of the logarithmic model ($m=0$)

As discussed in Section 4.2 and Remark 7, it appears that the $m=0$ material model is stable in dealing with uniform hydrostatic loadings in 2D and 3D. However, other types of deformation can cause instability. Consider a simple problem involving a 3D hyperelastic cuboid ($4\times 1\times 1$) under uniaxial loading (along $X$-axis). Explicit time integration is employed to simulate the problem with prescribed displacements at the ends of the material. As shown in Fig. 6, the model becomes unstable under relatively small tension or compression states. It is seen that using a smaller time step has no effect on the predicted macroscopic behavior in this case. While the $X$-displacement values appear in a smooth gradient, instabilities are evident in the $Y$- and $Z$-displacements.

5.4 Waves in a hyperelastic bar

As suggested by Belytschko et al. (2000), stability of a 1D hyperelastic bar is analyzed here. The material is modeled by different members of the finite deformation family, upon introduction of a perturbation deformation in the form of a planar wave

[TABLE]

in which $k$ is the wave number and $\omega$ is a complex angular frequency.

In 1D, the shape tensor $L=\int_{\mathcal{H}}{\rm d}\xi=\omega_{0}$. The force state in Eq. 6 then reads

[TABLE]

and

[TABLE]

where

[TABLE]

Using Eqs. 44, 45 and 46, the 1D peridynamic equation of motion is reduced as

[TABLE]

[TABLE]

Due to the nature of the finite deformation model, this equation is clearly nonlinear and dependent on $X$, $t$, and $u_{0}$. Angular frequency of $X=0$ and $t=0$ is calculated using a Hookean classical model. The dispersion curves are shown in Fig. 7 for two different $u_{0}$. The periodic zero angular frequencies evident in all members of the material model are indicative of rank deficiency as discussed by Belytschko et al. (2000). As expected, for a very small displacement field $u_{0}$ the predicted behavior of all different members of the material class are the same.

6 Discussion

In this section, the connection between the stability analysis presented in this work and a stability condition developed by Silling (2017) is discussed. Additionally, possible stabilization methods and the source of the instability are presented.

6.1 Silling’s material stability test

Using energy minimization, Silling (2017) derived a condition for stability

[TABLE]

The equilibrium is stable if Eq. 49 holds for a given displacement field $\underline{\mathbf{Y}}$.

In the context of the finite deformation theory, using Eq. 6,

[TABLE]

[TABLE]

Specializing to the case of uniform deformations and perturbing a single material point, due to symmetry and uniform loading, the stress state does not change for the point itself (cf. Remark 1), i.e. ${\rm d}S_{(m)ij}=0$. Thus, using Eq. 50,

[TABLE]

Similar to Section 3.1, for a pure hydrostatic loading as in Eq. 20, using the constitutive relation in Eq. 22, and limiting the perturbation to a single coordinate, e.g. ${\rm d}Y_{1}=-\varepsilon$ and ${\rm d}Y_{2}={\rm d}Y_{3}=0$, we obtain

[TABLE]

Note the similarity between Eq. 51 and Eq. 23. The similarity is expected as the energy minimization approach and the Jacobian matrix method are rooted to the same idea. The explicit Jacobian matrix approach allows an investigation into the cause of the instability.

6.2 Possibility of stabilization

Several stabilization methods were suggested for the original correspondence model, whose instability is mostly associated with nonuniform deformations. Silling (2017) added a penalty term to the material model, similar to the method of Littlewood (2010), which resists deviation from homogeneous deformations. Breitzman and Dayal (2018) and Chen (2018) proposed two different bond-level deformation gradients to reduce the influence of the average nodal integration and increase the impact of each bond deformation in evolving its own force state. Chowdhury et al. (2019) somehow removed the instability in some problems by using an arbitrary split of a family of bonds into multiple sub-horizon regions. It comes from the mathematical formulation of the original model that the averaging integral operator allows the cancellation of nonuniform parts of deformation, which is the main reason behind the unstable behavior. All the mentioned stabilization techniques aim to eliminate the spurious unphysical zero-energy modes that appear in a particle discretization of the model. The finite deformation correspondence theory, in contrast, has a different issue as homogeneous deformations can cause instability in its continuum level. Therefore, none of the mentioned stabilization methods can address the observed unstable behavior, and any potentially stabilization method would look artificial.

Another idea is to choose the $m$-member of the material class depending on the nature of the problem. For example, using $m>m_{\rm cr}$ when dealing with tension, and $m<m_{\rm cr}$ for compression. This idea, however, is doubted to be practical. First of all, as discussed in Section 5, instabilities are not limited to pure hydrostatic loadings and can happen in a variety of loading scenarios. It will be tedious, if possible, to distinguish all the stable $m$-models for every possible loading condition. In addition, significant inconsistencies would emerge in dealing with history-dependent problems, such as plasticity.

Pursuing a stabilization method is not recommended for a number of reasons. First, the finite deformation correspondence theory was supposed to naturally improve the stability of the original model. However, a different type of material instability is inherent in the model, which appears even for uniform deformations. Second, the finite deformation theory is computationally more expensive than the original model. In addition, the recent theory is limited to isotropic materials. Therefore, seeking its stabilization seems pointless, and using an existing stabilized version of the original model is more favorable.

6.3 Source of the instability

The origin of the problem is in fact rooted to the underlying classical model. It is well understood that the local Seth-Hill strain family suffers from the observed instability (Ball, 1976), i.e. $m<0$ group has tensile issues and $m>0$ shows compression difficulties. However, the problem is considerably milder in the classical model such that extremely larger tension or compression is required to cause an unstable behavior. The nonlocal integral operator in the correspondence theory largely exaggerates the problem, which results in small compressive or tensile deformations causing instability. The unstable behavior imposes significant difficulties in application of the finite deformation peridynamic model to any practical engineering problem.

Remark 8

Recently, Behzadinasab and Foster (2019), within the context of Sandia Fracture Challenge 2017 (Kramer et al., 2019), used the $m=1$ model to predict the deformations and failure of an additively manufacture metal under tensile loading. While the application of the theory resulted in a qualitatively good blind simulation results, the load-carrying capacity of the structure was underestimated and an early failure was predicted. Although tension was the dominant loading mode of the problem, compression and shear states were inevitable in the material, due to a complicated geometry. While the superiority of the peridynamic Green-Lagrange strain model in dealing with tensile loading was showed (see Section 5.1), even instability in a small region can destroy the accuracy of the macroscopic response.

7 Conclusions

The stability of a recently proposed peridynamic finite deformation correspondence model (Foster and Xu, 2018) is studied in a detailed examination. The finite deformation theory was developed by using a family of nonlinear strain measure and reported to improve the stability of the original correspondence model (Silling et al., 2007) whose instability is attributed to nonuniform deformations and particle discretizations. The offered analysis, however, shows that while the newer model has made some improvements over the original theory (cf. Section 5.1), it fails to maintain stability for homogeneous deformations in the mathematical foundation of the theory. The observed spurious unphysical deformation modes do not emerge because of particle discretization, which, actually interestingly, enhances the stability (cf. Remark 6). The unstable behavior is originated from the underlying local theory, in a less severe manner, and exaggerated largely through the nonlocal integral operator. The identified instability can significantly lower the reliability of the model in practical problems where complicated loading scenarios occur (cf. Remark 8). Since the recent model fails to naturally improve the stability of the original correspondence theory, is more computationally costly, and is limited to isotropic materials, its stabilization is not recommended. A stabilized version of the original model would be more desirable (see Section 6.2).

Acknowledgments

The authors are grateful for the financial support provided by the AFOSR MURI Center for Materials Failure Prediction through Peridynamics: Project NO. ONRBAA12-020. M. Behzadinasab also acknowledges the fellowship funding by The University of Texas at Austin.

Appendix Appendix A Simplified form of the shape tensor in the bulk of a Two-dimensional material

Constraining the influence state to spherical, i.e. $\underline{\mathbf{\omega}}_{s}\langle{\boldsymbol{\xi}}\rangle=\underline{\mathbf{\omega}}_{s}(|\boldsymbol{\xi}|)$, and noting symmetry, integration of all the odd combinations of indices in the shape tensor results in zero. The shape tensor is then evaluated to be (Foster and Xu, 2018, Section 3.2)

[TABLE]

These integrals are then computed in 2D:

[TABLE]

where $V_{\underline{\mathbf{\omega}}}$ is the weighted volume

[TABLE]

therefore,

[TABLE]

Due to its isotropic nature, $\mathbb{L}^{-1}$ can be written as:

[TABLE]

where $A$ and $B$ are unknown constants to be find through the following process:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Comparing the two sides of the equation

[TABLE]

hence,

[TABLE]

Substituting these back into Eq. 52, the inverted shape tensor is obtained

[TABLE]

Appendix Appendix B Simplifying Eq. 39

In this section, Eq. 39 is explicitly simplified by breaking it into different parts.

[TABLE]

[TABLE]

Appendix Appendix C Simplifying Eq. 41

Using Eq. 40 (note that indices can be either 1 or 2, for a two-dimensional setup)

[TABLE]

Then, for the hydrostatic loading condition of Eqs. 42 and 43, Eq. 41 is simplified

[TABLE]