Finite element convergence for state-based peridynamic fracture models

TL;DR

This paper proves the convergence rate of finite element methods for state-based peridynamic fracture models, including stability and numerical validation for dynamic crack propagation.

Contribution

It establishes the a-priori convergence rate for finite element approximations of nonlocal nonlinear fracture models, including stability analysis and numerical simulations.

Findings

Finite element solutions converge uniformly in the mean square norm.

Convergence rate for linear finite elements is C_t Δt + C_s h^2/ε^2.

Numerical simulations support the theoretical convergence rate.

Abstract

We establish the a-priori convergence rate for finite element approximations of a class of nonlocal nonlinear fracture models. We consider state based peridynamic models where the force at a material point is due to both the strain between two points and the change in volume inside the domain of nonlocal interaction. The pairwise interactions between points are mediated by a bond potential of multi-well type while multi point interactions are associated with volume change mediated by a hydrostatic strain potential. The hydrostatic potential can either be a quadratic function, delivering a linear force-strain relation, or a multi-well type that can be associated with material degradation and cavitation. We first show the well-posedness of the peridynamic formulation and that peridynamic evolutions exist in the Sobolev space $H^2$. We show that the finite element approximations converge to the $H^2$ solutions uniformly as measured in the mean square norm. For linear continuous finite elements the convergence rate is shown to be $C_t \Delta t + C_s h^2/\epsilon^2$, where $\epsilon$ is the size of horizon, $h$ is the mesh size, and $\Delta t$ is the size of time step. The constants $C_t$ and $C_s$ are independent of $\Delta t$ and $h$ and may depend on $\epsilon$ through the norm of the exact solution. We demonstrate the stability of the semi-discrete approximation. The stability of the fully discrete approximation is shown for the linearized peridynamic force. We present numerical simulations with dynamic crack propagation that support the theoretical convergence rate.

Equations

$$S(\mathbf{y},\mathbf{x},t;\mathbf{u})=\frac{\mathbf{u}(\mathbf{y},t)-\mathbf{u}(\mathbf{x},t)}{|\mathbf{y}-\mathbf{x}|}\cdot \frac{\mathbf{y}-\mathbf{x}}{|\mathbf{y}-\mathbf{x}|}.$$

$${\mathcal{W}}^{ϵ}(S(\mathbf{y},\mathbf{x},t;\mathbf{u}))=\frac{J^{ϵ}(|\mathbf{y}-\mathbf{x}|)}{ϵ|\mathbf{y}-\mathbf{x}|}f(\sqrt{|\mathbf{y}-\mathbf{x}|}S(\mathbf{y},\mathbf{x},t;\mathbf{u})),$$

$$\theta (\mathbf{x},t;\mathbf{u})=\frac{1}{{ϵ}^d{\omega }_d}{\int }_{H_{ϵ}(\mathbf{x})}{J}^{ϵ}(|\mathbf{y}-\mathbf{x}|)S(\mathbf{y},\mathbf{x},t;\mathbf{u})|\mathbf{y}-\mathbf{x}|d\mathbf{y},$$

$${\mathcal{V}}^{ϵ}(\theta (\mathbf{x},t;\mathbf{u}))=\frac{g(\theta (\mathbf{x},t;\mathbf{u}))}{{ϵ}^2},$$

$${\theta }_c^{+}=r_{\ast }^{+},{\theta }_c^{-}=r_{\ast }^{-}.$$

$$\theta (\mathbf{x},t;\mathbf{u})=\frac{1}{{ϵ}^d{\omega }_d}{\int }_{H_{ϵ}(\mathbf{x})}\omega (\mathbf{y})J^{ϵ}(|\mathbf{y}-\mathbf{x}|)S(\mathbf{y},\mathbf{x},t;\mathbf{u})|\mathbf{y}-\mathbf{x}|d\mathbf{y}.$$

$${\mathcal{W}}^{ϵ}(S(\mathbf{y},\mathbf{x},t;\mathbf{u}))$$

$${\mathcal{V}}^{ϵ}(\theta (\mathbf{x},t;\mathbf{u}))$$

$$C_0^f:=\underset{r}{\sup }|f(r)|<\infty ,C_i^f:=\underset{r}{\sup }|f^{(i)}(r)|<\infty ,\forall i=1,2,3,4.$$

$$P{D}^{ϵ}(\mathbf{u}(t))=\frac{1}{{ϵ}^d{\omega }_d}{\int }_D{\int }_{H_{ϵ}(\mathbf{x})}|\mathbf{y}-\mathbf{x}|{\mathcal{W}}^{ϵ}(S(\mathbf{y},\mathbf{x},t;\mathbf{u}))d\mathbf{y}d\mathbf{x}$$

$$+{\int }_D{\mathcal{V}}^{ϵ}(\theta (\mathbf{x},t;\mathbf{u}))d\mathbf{x},$$

$$\mathrm{L}(\mathbf{u},{\partial }_t\mathbf{u},t)=\frac{\rho }{2}||\dot{\mathbf{u}}|{|}^2-P{D}^{ϵ}(\mathbf{u}(t))+{\int }_D\mathbf{b}(t)\cdot \mathbf{u}(t)d\mathbf{x},$$

$$\rho \ddot{\mathbf{u}}(\mathbf{x},t)={\mathcal{L}}^{ϵ}(\mathbf{u})(\mathbf{x},t)+\mathbf{b}(\mathbf{x},t),for\mathbf{x}\in D,$$

$${\mathcal{L}}^{ϵ}(\mathbf{u})(\mathbf{x},t)={\mathcal{L}}_T^{ϵ}(\mathbf{u})(\mathbf{x},t)+{\mathcal{L}}_D^{ϵ}(\mathbf{u})(\mathbf{x},t),$$

$${\mathcal{L}}_T^{ϵ}(\mathbf{u})(\mathbf{x},t)$$

$$=\frac{2}{{ϵ}^d{\omega }_d}{\int }_{H_{ϵ}(\mathbf{x})}\omega (\mathbf{x})\omega (\mathbf{y})\frac{J^{ϵ}(|\mathbf{y}-\mathbf{x}|)}{ϵ|\mathbf{y}-\mathbf{x}|}{\partial }_{S}f(\sqrt{|\mathbf{y}-\mathbf{x}|}S(\mathbf{y},\mathbf{x},t;\mathbf{u})){\mathbf{e}}_{\mathbf{y}-\mathbf{x}}d\mathbf{y},$$

$${\mathcal{L}}_D^{ϵ}(\mathbf{u})(\mathbf{x},t)$$

$$=\frac{1}{{ϵ}^d{\omega }_d}{\int }_{H_{ϵ}(\mathbf{x})}\omega (\mathbf{x})\omega (\mathbf{y})\frac{J^{ϵ}(|\mathbf{y}-\mathbf{x}|)}{{ϵ}^2}\left[{\partial }_{\theta }g(\theta (\mathbf{y},t;\mathbf{u}))+{\partial }_{\theta }g(\theta (\mathbf{x},t;\mathbf{u}))\right]{\mathbf{e}}_{\mathbf{y}-\mathbf{x}}d\mathbf{y}.$$

$$\mathbf{u}(\mathbf{x},0)={\mathbf{u}}_0(\mathbf{x}),{\partial }_t\mathbf{u}(\mathbf{x},0)={\mathbf{v}}_0(\mathbf{x}).$$

$$\mathbf{u}(\mathbf{x})=0\forall \mathbf{x}\in \partial D.$$

$$W=H^2(D;{\mathbb{R}}^d)\cap H_0^1(D;{\mathbb{R}}^d)=\{v\in H^2(D;{\mathbb{R}}^d):\gamma v=0,on\partial D\}$$

$$||\mathbf{u}|{|}_{\infty }\le C_{e_1}||\mathbf{u}|{|}_2.$$

$$||\nabla \mathbf{u}|{|}_{L^q(D;{\mathbb{R}}^{d\times d})}\le C_{e_2}||\nabla \mathbf{u}|{|}_1\le C_{e_2}||\mathbf{u}|{|}_2,$$

$$||{\mathcal{L}}^{ϵ}(\mathbf{u})-{\mathcal{L}}^{ϵ}(\mathbf{v})|{|}_2$$

$$||{\mathcal{L}}^{ϵ}(\mathbf{u})|{|}_2$$

$$||{\partial }_{ttt}^3\mathbf{u}(\mathbf{x},t)|{|}_2$$

$$||{\mathcal{L}}^{ϵ}(\mathbf{u})-{\mathcal{L}}^{ϵ}(\mathbf{v})||$$

$$L$$

$$(\ddot{\mathbf{u}}(t),\tilde{\mathbf{u}})=({\mathcal{L}}^{ϵ}(\mathbf{u}(t)),\tilde{\mathbf{u}})+(\mathbf{b}(t),\tilde{\mathbf{u}}).$$

$$({\mathcal{L}}^{ϵ}(\mathbf{u}),\mathbf{v})=-a^{ϵ}(\mathbf{u},\mathbf{v}),$$

$$a^{ϵ}(\mathbf{u},\mathbf{v})=a_T^{ϵ}(\mathbf{u},\mathbf{v})+a_D^{ϵ}(\mathbf{u},\mathbf{v})$$

$$a_T^{ϵ}(\mathbf{u},\mathbf{v})$$

Full text

{AMS}

34A34, 34B10, 74H55, 74S05

1 Introduction

In this work, we study the state-based peridynamic theory and obtain an a-priori error bound for the finite element approximation. The peridynamic theory is a reformulation of classical continuum mechanics carried out in the work of Silling in [Silling, 2000, Silling et al., 2007]. The strain inside the medium is expressed in terms of displacement differences as opposed to the displacement gradients. Acceleration of a point is now due to the sum of the forces acting on the point from near by points. The new kinematics bypasses the difficulty incurred by juxtaposing displacement gradients and discontinuities as in the case of classical fracture theories. The nonlocal fracture theory has been applied numerically to model the complex fracture phenomenon in materials, see [Weckner and Abeyaratne, 2005], [Silling and Bobaru, 2005], [Silling and Lehoucq, 2008],

[Silling et al., 2010], [Foster et al., 2011], [Ha and Bobaru, 2010], [Agwai et al., 2011],

[Bobaru and Hu, 2012], [Ghajari et al., 2014], [Du et al., 2013a], [Lipton et al., 2016]. Every point interacts with its neighbors inside a ball of fixed radius called the horizon. The size of horizon sets the length scale of nonlocal interaction. When the forces between points are linear and nonlocal length scale tends to zero, it is seen that peridynamic models converge to the classic model of linear elasticity, [Emmrich et al., 2013, Silling and Lehoucq, 2008, Aksoylu and Unlu, 2014, Mengesha and Du, 2015]. The work of [Tian and Du, 2014] provides an analytic framework for analyzing FEM for linear bond and state-based peridynamics. For nonlinear forces associated with double well potentials the peridynamic evolution converges in the small horizon limit to an evolution with a sharp evolving fracture set and the evolution is governed by the classic linear elastic wave equation away from the fracture set, see [Lipton, 2016, Lipton, 2014, Jha and Lipton, 2018a]. A recent review of the state of the art can be found in [Bobaru et al., 2016] and [Du, 2018a].

In this work we assume small deformation and work with linearized bond-strain. Let $D\subset\mathbb{R}^{d}$, for $d=2,3$, be the material domain. For a displacement field $\boldsymbol{u}:D\times[0,T]\to\mathbb{R}^{d}$, the bond-strain between two material points $\boldsymbol{x},\boldsymbol{y}\in D$ is given by

[TABLE]

Let $\epsilon>0$ be the size of horizon and $H_{\epsilon}(\boldsymbol{x})=\{\boldsymbol{y}\in\mathbb{R}^{d}:|\boldsymbol{y}-\boldsymbol{x}|<\epsilon\}$ be the neighborhood of a material point $\boldsymbol{x}$. For pairwise interaction, we assume the following form of pairwise interaction potential

[TABLE]

where $J^{\epsilon}(|\boldsymbol{y}-\boldsymbol{x}|)$ is the influence function. We assume $J^{\epsilon}(|\boldsymbol{y}-\boldsymbol{x}|)=J(|\boldsymbol{y}-\boldsymbol{x}|/\epsilon)$ where $0\leq J(r)\leq M$ for $r<1$ and $J(r)=0$ for $r\geq 1$. The potential $f$, see 1(a), is assumed to be convex for small strains and becomes concave for larger strains. In the widely used prototypical micro-elastic brittle (PMB) peridynamic material, the strain vs force profile is linear up to some critical strain $S_{c}$ and is zero for any strain above $S_{c}$. In contrast, the peridynamic force given by $\partial_{S}\mathcal{W}^{\epsilon}$, is linear near zero strain and as the strain gets larger and reaches the critical strain, $S_{c}^{+}$ ($S_{c}^{-}$) for positive (negative) strain, the bond starts to soften, see 1(b). For a given potential function $f$, the critical strain is given by $S_{c}^{+}=\frac{r^{+}}{\sqrt{|\boldsymbol{y}-\boldsymbol{x}|}}$ and $S_{c}^{-}=\frac{r^{-}}{\sqrt{|\boldsymbol{y}-\boldsymbol{x}|}}$ where $r^{+}>0,r^{-}<0$ are the inflection points of the potential function $f$ as shown in 1(a).

The spherical or hydrostatic strain $\theta(\boldsymbol{x},t;\boldsymbol{u})$ at material point is given by

[TABLE]

where $\omega_{d}$ is the volume of unit ball in dimension $d=2,3$. The potential for hydrostatic interaction is of the form

[TABLE]

where $g$ is the potential function associated to hydrostatic strain. Here $g$ can be of two types: 1) a quadratic function with only one well at zero strain, and 2) a convex-concave function with a wells at the origin and at $\pm\infty$, see 2(a). If $g$ is assumed to be quadratic then the force due to spherical strain is linear. If $g$ is a multi-well potential, the material softens as the hydrostatic strains exceeds critical value. For the convex-concave type $g$, the critical values are $0<\theta_{c}^{+}$ and $\theta^{-}_{c}<0$ beyond which the force begins to soften is related to the inflection point $r^{+}_{\ast}$ and $r_{\ast}^{-}$ of $g$ as follows

[TABLE]

The critical compressive hydrostatic strain where the force begins to soften for negative hydrostatic strain is chosen much larger in magnitude than $\theta_{c}^{+}$, i.e. $\theta^{+}_{c}<<|\theta_{c}^{-}|$.

The finite element approximation has been applied to peridynamic fracture, however there remains a paucity of literature addressing the rigorous a-priori convergence rate of the finite element approximation to peridynamic problems in the presence of material failure. This aspect provides the motivation for the present work. In this paper we first prove existence of peridynamic evolutions taking values in $H^{2}(D;\mathbb{R}^{d})\cap H^{1}_{0}(D;\mathbb{R}^{d})$ that are twice differentiable in time, see Theorem 2.2. We note that as these evolutions will become more fracture like as the region of nonlocal interaction decreases. These evolutions can be thought of as inner approximations to fracture evolutions. On passing to subsequences it is possible to show that the $H^{2}(D;\mathbb{R}^{d})\cap H^{1}_{0}(D;\mathbb{R}^{d})$ evolutions converge in the limit of vanishing non-locality to a limit solution taking values in the space of special functions of bounded deformation SBD. Here the limit evolution has a well defined Griffith fracture energy bounded by the initial data, see [Lipton, 2016] and [Jha and Lipton, 2019]. We show here that higher temporal regularity can be established if the body force changes smoothly in time. Motivated by these considerations we develop finite element error estimates for solutions that take values in $H^{2}(D;\mathbb{R}^{d})\cap H^{1}_{0}(D;\mathbb{R}^{d})$ and for a bounded time interval.

In this paper we obtain an a-priori $L^{2}$ error bound for the finite element approximation of the displacement and velocity using a central in time discretization. Due to the nonlinear nature of the problem we get a convergence rate by using Lax Richtmyer stability together with consistency. Both stability and consistency are shown to follow from the Lipschitz continuity of the peridynamic force in $L^{2}(D;\mathbb{R}^{d})$, see section 4.2.1 and section 4.2.2. The bound on the $L^{2}$ error is uniform in time and is given by $C_{t}\Delta t+C_{s}h^{2}/\epsilon^{2}$, where the constants $C_{t}$ and $C_{s}$ are independent of $\Delta t$ and mesh size $h$, see Theorem 4.1. A more elaborate discussion of the a-priori bound is presented in section subsection 4.2. For the linearized model we obtain a stability condition on $\Delta t$, Theorem 4.4, that is of the same form as those given for linear local and nonlocal wave equations [Karaa, 2012, Guan and Gunzburger, 2015]. We demonstrate stability for the linearized model noting that for small strains the material behaves like a linear elastic material and that the stability of the linearized model is necessary for the stability of nonlinear model. We believe a more constructive CFL stability condition is possible for the linear case and will pursue this in future work.

Previous work [Jha and Lipton, 2018a] treated spatially Lipschitz continuous solutions and addressed the finite difference approximation and obtained bounds on the $L^{2}$ error for the displacement and velocity that are uniform in time and of the form $C_{t}\Delta t+C_{s}h/\epsilon^{2}$, where constants the $C_{t}$ and $C_{s}$ are as before. For finite elements the convergence rate is seen to be slower than for the FEM model introduced here and is of order $h/\epsilon^{2}$ as opposed to $h^{2}/\epsilon^{2}$. On the other hand the FEM method increases the computational work due to the inversion of the mass matrix.

We carry out numerical experiments for dynamic crack propagation and obtain convergence rates for Plexiglass that are in line with the theory, see section 5. We also compare the Griffith’s fracture energy with the peridynamic energy of the material softening zone we show good agreement between the two energies, see subsection 5.2. Finite difference methods are less expensive than finite element approximations for nonlocal problems, however the latter offers more control on the accuracy of solution see, [Macek and Silling, 2007, Littlewood, 2010, Du et al., 2013c, Gerstle et al., 2007, Du, 2018b].

Here the a-priori $L^{2}$ convergence rates for the FEM given by Theorem 4.1 include the effects of material degradation through the softening of material properties. The FEM simulations presented in this paper show that the material develops localized softening zones (region where bonds exceed critical tensile strain) as it deforms. This is in contrast to linear peridynamic models which are incapable of developing softening zones. For nonlinear peridynamic models with material failure the localization of zones of softening and damage is the hallmark of peridynamic modeling [Silling, 2000, Silling et al., 2007], [Ha and Bobaru, 2010], [Foster et al., 2011]. One notes that the a-priori error involves $\epsilon$ in the denominator and in many cases $\epsilon$ is chosen small. However typical dynamic fracture experiments last only hundreds of microseconds and the a-priori error is controlled by the product of simulation time multiplied by $h^{2}/\epsilon^{2}$. So for material properties characteristic of Plexiglass and $\epsilon$ of size $4mm$, the a-priori estimates predict a relative error of $1/10$ for simulations lasting around 100 microseconds. We point out that the a-priori error estimates assume the appearance of nonlinearity anywhere in the computational domain. On the other hand numerical simulation and independent theoretical estimates show that the nonlinearity concentrates along “fat” cracks of finite length and width equal to $\epsilon$, see [Lipton, 2016, Lipton, 2014]. Moreover the remainder of the computational domain is seen to behave linearly and to leading order can be modeled as a linear elastic material up to an error proportional to $\epsilon$, see [Proposition 6, [Jha and Lipton 2018b]]. Future work will use these observations to focus on adaptive implementation and a-posteriori estimates. A-posteriori convergence for FEM models of peridynamics with material degradation can be seen in the work [Macek and Silling, 2007], [Chen and Gunzburger, 2011],

[Ren et al., 2017]. For other nonlinear and nonlocal models adaptive mesh refinement within FE framework for nonlocal models has been explored in [Du et al., 2013c] and convergence of the adaptive FE approximation is rigorously shown. A-posteriori error analysis of linear nonlocal models is carried out in [Du et al., 2013b].

The paper is organized as follows. We introduce the equation of motion in section 2 and present the Lipschitz continuity of the force, existence of peridynamic solution, and the higher temporal regularity necessary for the finite element error analysis. In section 3, we consider the finite element discretization. We prove the stability of a semi-discrete approximation in subsection 3.1. In section 4, we analyze the fully discrete approximation and obtain an a-priori bound on errors. The stability of the fully discrete approximation linearized peridynamic force is shown in subsection 4.3. We present our numerical experiments in section 5. Proofs of the Lipschitz bound on the peridynamic force and higher temporal regularity of solutions is provided in section 6. In section 7 we present our conclusions.

We conclude the introduction by listing the notation used throughout the paper. We denote material domain as $D$ where $D\subset\mathbb{R}^{d}$ for $d=2,3$. Points and vectors in $\mathbb{R}^{d}$ are denoted as bold letters. Some of the key notations are as follows

[TABLE]

2 Equation of motion, existence, uniqueness, and higher regularity

We assume $D$ to be an open set with $C^{1}$ boundary. To enforce zero displacement boundary conditions at $\partial D$ and to insure a well posed evolution we introduce the boundary function $\omega(\boldsymbol{x})$. This function is introduced as a factor into the potentials $\mathcal{W}^{\epsilon}$ and $\mathcal{V}^{\epsilon}$. Here the boundary function takes value $1$ in the interior of domain and is zero on the boundary. We assume $\sup_{\boldsymbol{x}}|\nabla\omega(\boldsymbol{x})|<\infty$ and $\sup_{\boldsymbol{x}}|\nabla^{2}\omega(\boldsymbol{x})|<\infty$ in our analysis. The hydrostatic strain is modified to include the boundary and is given by

[TABLE]

The peridynamic potentials Equation 2 and Equation 4 are modified to see the boundary follows

[TABLE]

We assume that potential function $f$ is at least 4 times differentiable and satisfies following regularity condition:

[TABLE]

If potential function $g$ is convex-concave type then we assume that $g$ satisfies same regularity condition as $f$. We denote constants $C^{g}_{i}$, for $i=0,1,..,4$, similar to $C^{f}_{i}$ above.

The total potential energy at time $t$ is given by

[TABLE]

where potential $\mathcal{W}^{\epsilon}$ and $\mathcal{V}^{\epsilon}$ are described above. The material is assumed to be homogeneous and the density is given by $\rho$. The applied body force is denoted by $\boldsymbol{b}(\boldsymbol{x},t)$. We define the Lagrangian

[TABLE]

here $\dot{\boldsymbol{u}}$ is the velocity given by the time derivative of $\boldsymbol{u}$. Applying the principal of least action together with a straight forward calculation (see for example [Lipton et al., 2018a] for detailed derivation) gives the nonlocal dynamics

[TABLE]

where

[TABLE]

$\mathcal{L}^{\epsilon}_{T}(\boldsymbol{u})$ is the peridynamic force due to bond-based interaction and is given by

[TABLE]

and $\mathcal{L}^{\epsilon}_{D}(\boldsymbol{u})$ is the peridynamic force due to state-based interaction and is given by

[TABLE]

The dynamics is complemented with the initial data

[TABLE]

We prescribe zero Dirichlet boundary condition on the boundary $\partial D$

[TABLE]

We extend the zero boundary condition outside $D$ to whole $\mathbb{R}^{d}$. In our analysis we will assume mass density $\rho=1$ without loss of generality.

2.1 Existence of solutions and higher regularity in time

We recall that the space $H^{n}_{0}(D;\mathbb{R}^{d})$ is the closure in the $H^{n}$ norm of the functions that are infinitely differentiable with compact support in $D$. For suitable initial conditions and body force we show that solutions exist in

[TABLE]

where $\gamma$ is the trace of the function $v$ on the boundary of $D$. We will assume that $\boldsymbol{u}\in W$ is extended by zero outside $D$. We first exhibit the Lipschitz continuity property and boundedness of the peridynamic force for displacements in $W$. We will then apply [Theorem 3.2, [Jha and Lipton, 2017]] to conclude the existence of unique solutions.

We note the following Sobolev embedding properties of $H^{2}(D;\mathbb{R}^{d})$ when $D$ is a $C^{1}$ domain.

• •

From Theorem 2.72 of [Demengel, 2012], there exists a constant $C_{e_{1}}$ independent of $\boldsymbol{u}\in H^{2}(D;\mathbb{R}^{d})$ such that

[TABLE]

• •

Further application of standard embedding theorems

(e.g., Theorem 2.72 of [Demengel, 2012]), shows there exists a constant $C_{e_{2}}$ independent of $\boldsymbol{u}$ such that

[TABLE]

for any $q$ such that $2\leq q<\infty$ when $d=2$ and $2\leq q\leq 6$ when $d=3$.

We have following result which show the Lipschitz continuity property of a peridynamic force $\mathcal{L}^{\epsilon}$.

Theorem 2.1**.**

**Lipschitz continuity of peridynamic force

Let $f$ be a convex-concave function satisfying $C^{f}_{i}<\infty$ for $i=0,..,4$ and let $g$ either be a quadratic function, or $g$ be a convex-concave function with $C^{g}_{i}<\infty$ for $i=0,..,4$. Also, let boundary function $\omega:D\to[0,1]$ be such that $\sup_{\boldsymbol{x}\in D}|\nabla\omega(\boldsymbol{x})|<\infty$ and $\sup_{\boldsymbol{x}\in D}|\nabla^{2}\omega(\boldsymbol{x})|<\infty$. Then, for any $\boldsymbol{u},\boldsymbol{v}\in W$, we have**

[TABLE]

where constant $\bar{L}_{1}$ does not depend on $\epsilon$ and $\boldsymbol{u},\boldsymbol{v}$. Also, for $\boldsymbol{u}\in W$, we have

[TABLE]

*where constant $\bar{L}_{2}$ does not depend on $\epsilon$ and $\boldsymbol{u}$. *

Now let $T>0$ be any positive number, a straight-forward application of [Theorem 3.2, [Jha and Lipton, 2017]] gives:

Theorem 2.2**.**

**Existence and uniqueness of solutions over finite time intervals

Let $f$, $g$, and $\omega$ satisfy the hypothesis of Theorem 2.1. For any initial condition $\boldsymbol{u}_{0},\boldsymbol{v}_{0}\in W$, time interval $I_{0}=(-T,T)$, and right hand side $\boldsymbol{b}(t)$ continuous in time for $t\in I_{0}$ such that $\boldsymbol{b}(t)$ satisfies $\sup_{t\in I_{0}}||\boldsymbol{b}(t)||_{2}<\infty$, there is a unique solution $\boldsymbol{u}(t)\in C^{2}(I_{0};W)$ of peridynamic equation 11. Also, $\boldsymbol{u}(t)$ and $\dot{\boldsymbol{u}}(t)$ are Lipschitz continuous in time for $t\in I_{0}$. **

We can also show higher regularity in time of evolutions under suitable assumptions on the body force:

Theorem 2.3**.**

**Higher regularity

Suppose the initial data and righthand side $\boldsymbol{b}(t)$ satisfy the hypothesis of Theorem 2.2 and suppose further that $\dot{\boldsymbol{b}}(t)$ exists and is continuous in time for $t\in I_{0}$ and  $\sup_{t\in I_{0}}||\dot{\boldsymbol{b}}(t)||_{2}<\infty$. Then $\boldsymbol{u}\in C^{3}(I_{0};W)$ and**

[TABLE]

where $C$ is a positive constant independent of $\boldsymbol{u}$.

The proofs of Theorems 2.1, and 2.3 are given in section 6. For future reference, we note that for any $\boldsymbol{u},\boldsymbol{v}\in L^{2}_{0}(D;\mathbb{R}^{d})$, we have

[TABLE]

where constant $L$ is given by

[TABLE]

and $\bar{J}^{\alpha}=(1/\omega_{d})\int_{H_{1}(\mathbf{0})}J(|{\boldsymbol{\xi}}|)/|{\boldsymbol{\xi}}|^{\alpha}d{\boldsymbol{\xi}}$ .

2.2 Weak form

We multiply Equation 11 by a test function $\tilde{\boldsymbol{u}}$ in $H^{1}_{0}(D;\mathbb{R}^{d})$ and integrate over $D$ to get

[TABLE]

We have the following integration by parts formula:

Lemma 2.4**.**

For any $\boldsymbol{u},\boldsymbol{v}\in L^{2}_{0}(D;\mathbb{R}^{d})$ we have

[TABLE]

where

[TABLE]

and

[TABLE]

The proof of above lemma is identical to the proof of Lemma 4.2 in [Lipton et al., 2018a].

Using the above Lemma, the weak form of the peridynamic evolution is given by

[TABLE]

*Total dynamic energy: *We define the total dynamic energy as follows

[TABLE]

where $PD^{\epsilon}$ is defined in Equation 10. Time derivative of total energy satisfies

[TABLE]

Remark. It is readily verified that the peridynamic force and energy are bounded for all functions in $L^{2}(D;\mathbb{R}^{d})$. Here the bound on the force follows from the Lipschitz property of the force in $L^{2}(D;\mathbb{R}^{d})$, see, Equation 23. The peridynamic force is also bounded for functions $\boldsymbol{u}$ in $H^{1}(D;\mathbb{R}^{d})$. This again follows from the Lipschitz property of the force in $H^{1}(D;\mathbb{R}^{d})$ using arguments established in section 6. The boundedness of the energy $PD^{\epsilon}(\boldsymbol{u})$ in both $L^{2}(D;\mathbb{R}^{d})$ and $H^{1}(D;\mathbb{R}^{d})$ follows from the boundedness of the bond potential energy $\mathcal{W}^{\epsilon}(S(\boldsymbol{y},\boldsymbol{x},t;\boldsymbol{u}))$ and $\mathcal{V}^{\epsilon}(\theta(\boldsymbol{x},t;\boldsymbol{u}))$ used in the definition of $PD^{\epsilon}(\boldsymbol{u})$, see Equation 7 and Equation 8. More generally this also shows that $PD^{\epsilon}(\boldsymbol{u})<\infty$ for $\boldsymbol{u}\in L^{1}(D;\mathbb{R}^{d})$.

We next discuss the spatial and the time discretization of peridynamic equation.

3 Finite element approximation

Let $V_{h}$ be given by linear continuous interpolations over tetrahedral or triangular elements $\mathcal{T}_{h}$ where $h$ denotes the size of finite element mesh. Here we assume the elements are conforming and the finite element mesh is shape regular and $V_{h}\subset H^{1}_{0}(D;\mathbb{R}^{d})$.

For a continuous function $\boldsymbol{u}$ on $\bar{D}$, $\mathcal{I}_{h}(\boldsymbol{u})$ is the continuous piecewise linear interpolant on $\mathcal{T}_{h}$. It is given by

[TABLE]

where $\mathcal{I}_{T}(\boldsymbol{u})$ is the local interpolant defined over finite element $T$ and is given by

[TABLE]

Here $n$ is the number of vertices in an element $T$, $\boldsymbol{x}_{i}$ is the position of vertex $i$, and $\phi_{i}$ is the linear interpolant associated to vertex $i$.

Application of Theorem 4.4.20 and remark 4.4.27 in [Brenner and Scott, 2007] gives

[TABLE]

Let $\boldsymbol{r}_{h}(\boldsymbol{u})$ denote the projection of $\boldsymbol{u}\in W$ on $V_{h}$. For the $L^{2}$ norm it is defined as

[TABLE]

and satisfies

[TABLE]

Since $\mathcal{I}_{h}(\boldsymbol{u})\in V_{h}$, and Equation 34 we see that

[TABLE]

3.1 Semi-discrete approximation

Let $\boldsymbol{u}_{h}(t)\in V_{h}$ be the approximation of $\boldsymbol{u}(t)$ satisfying following for all $t\in[0,T]$

[TABLE]

We have following result:

Theorem 3.1**.**

**Energy stability of semi-discrete approximation

The semi-discrete scheme is stable and the energy $\mathcal{E}^{\epsilon}(\boldsymbol{u}_{h})(t)$, defined in Equation 30, satisfies the following bound**

[TABLE]

We note that while proving the stability of semi-discrete scheme corresponding to nonlinear peridynamics, we do not require any assumption on the strain $S(\boldsymbol{y},\boldsymbol{x},t;\boldsymbol{u}_{h})$. The proof is similar to [Section 6.2, [Lipton, 2016]].

Proof 3.2**.**

Letting $\tilde{\boldsymbol{u}}=\dot{\boldsymbol{u}}_{h}(t)$ in Equation 38 and noting the identity Equation 31, we get

[TABLE]

We also have

[TABLE]

where we used the fact that $PD^{\epsilon}(\boldsymbol{u})(t)$ is nonnegative. We substitute above inequality in Equation 39 to get

[TABLE]

We fix $\delta>0$ and define $A(t)$ as $A(t)=\mathcal{E}^{\epsilon}(\boldsymbol{u}_{h})(t)+\delta$. Then, from the above equation we easily have

[TABLE]

Noting that $\frac{1}{\sqrt{a(t)}}\frac{da(t)}{dt}=2\frac{d}{dt}\sqrt{a(t)}$, integrating from $t=0$ to $\tau$ and relabeling $\tau$ as $t$, we get

[TABLE]

*Proof is complete once we let $\delta\to 0$ and take the square of both sides. *

4 Central difference time discretization

In subsection 4.2 we calculate the convergence rate for the central difference time discretization of the fully nonlinear problem. We then present a CFL like condition on the time step $\Delta t$ for the linearized peridynamic equation in subsection 4.3.

At time step $k$, the exact solution is given by $(\boldsymbol{u}^{k},\boldsymbol{v}^{k})$ where $\boldsymbol{v}^{k}=\partial\boldsymbol{u}^{k}/\partial t$, and their projection onto $V_{h}$ is given by $(\boldsymbol{r}_{h}(\boldsymbol{u}^{k}),\boldsymbol{r}_{h}(\boldsymbol{v}^{k}))$. The solution of fully discrete problem at time step $k$ is given by $(\boldsymbol{u}^{k}_{h},\boldsymbol{v}^{k}_{h})$.

We approximate the initial data on displacement $\boldsymbol{u}_{0}$ and velocity $\boldsymbol{v}_{0}$ by their projections $\boldsymbol{r}_{h}(\boldsymbol{u}_{0})$ and $\boldsymbol{r}_{h}(\boldsymbol{v}_{0})$. Let $\boldsymbol{u}^{0}_{h}=\boldsymbol{r}_{h}(\boldsymbol{u}_{0})$ and $\boldsymbol{v}^{0}_{h}=\boldsymbol{r}_{h}(\boldsymbol{v}_{0})$. For $k\geq 1$, $(\boldsymbol{u}^{k}_{h},\boldsymbol{v}^{k}_{h})$ satisfies, for all $\tilde{\boldsymbol{u}}\in V_{h}$,

[TABLE]

where we have denoted the projection of $\boldsymbol{b}(t^{k})$, i.e. $\boldsymbol{r}_{h}(\boldsymbol{b}(t^{k}))$, as $\boldsymbol{b}^{k}_{h}$. Combining the two equations delivers central difference equation for $\boldsymbol{u}^{k}_{h}$. We have

[TABLE]

For $k=0$, we have $\forall\tilde{\boldsymbol{u}}\in V_{h}$

[TABLE]

4.1 Implementation details

For completeness we describe the implementation of the time stepping method using FEM interpolants. Let $\boldsymbol{N}$ be the shape tensor then $\boldsymbol{u}^{k}_{h},\tilde{\boldsymbol{u}}\in V_{h}$ are given by

[TABLE]

where $\boldsymbol{U}^{k}$ and $\tilde{\boldsymbol{U}}$ are $Nd$ dimensional vectors, where $N$ is the number of nodal points in the mesh and $d$ is the dimension.

From Equation 41, for all $\tilde{\boldsymbol{U}}\in\mathbb{R}^{Nd}$ with elements of $\tilde{\boldsymbol{U}}$ zero on the boundary, then the following holds for $k\geq 1$

[TABLE]

Here the mass matrix $\boldsymbol{M}$ and force vector $\boldsymbol{F}^{k}$ are given by

[TABLE]

where $\boldsymbol{F}^{k}_{pd}$ is defined by

[TABLE]

We remark that a similar equation holds for $k=0$.

At the time step $k$ we must invert ${\boldsymbol{M}}$ to solve for ${\boldsymbol{U}}^{k+1}$ using

[TABLE]

As is well known this inversion amounts to an increase of computational complexity associated with discrete approximation of the weak formulation of the evolution. Further, the matrix-vector multiplication $\boldsymbol{M}^{-1}\boldsymbol{F}^{k}$ needs to be carried out at each time step. On the other hand the quadrature error in the computation of the force vector $\boldsymbol{F}^{k}_{pd}$ is reduced when using the weak form.

We next show the convergence of approximation.

4.2 Convergence of approximation

In this section, we prove the uniform bound on the error and show that the approximate solution converges to the exact solution with rate given by $C_{t}\Delta t+C_{s}h^{2}/\epsilon^{2}$. Here horizon $\epsilon>0$ is assumed to be fixed. We first compare the exact solution with its projection in $V_{h}$ and then compare the projection with the approximate solution. We further divide the calculation of error between the projection and the approximate solution in two parts, namely consistency analysis and error analysis.

Error $E^{k}$ is given by

[TABLE]

The error is split into two parts as follows

[TABLE]

where the first term is the error between the exact solution and projection, and the second term is the error between the projection and approximate solution. Let

[TABLE]

and

[TABLE]

Using Equation 37, we have

[TABLE]

where

[TABLE]

We have the following a-priori convergence rate given by

Theorem 4.1**.**

**Convergence of Central difference approximation

Let $(\boldsymbol{u},\boldsymbol{v})$ be the exact solution of the peridynamic equation 11. Let $(\boldsymbol{u}^{k}_{h},\boldsymbol{v}^{k}_{h})$ be the FE solution of Equation 40. If $\boldsymbol{u},\boldsymbol{v}\in C^{2}([0,T];W)$, then the scheme is consistent and the error $E^{k}$ satisfies following bound**

[TABLE]

*where the constants $C_{p}$, $C_{t}$, and $C_{s}$ are given by Equation 51 and Equation 58. The constant $L/\epsilon^{2}$ is the Lipschitz constant of the peridynamic force $\mathcal{L}^{\epsilon}(\boldsymbol{u})$ in $L^{2}$, see Equation 23. If the error in initial data is zero then $E^{k}$ is of the order of $C_{t}\Delta t+C_{s}h^{2}/\epsilon^{2}$. *

In Theorem 2.3 we have shown that $\boldsymbol{u},\boldsymbol{v}\in C^{2}([0,T];W)$ for righthand side $\boldsymbol{b}\in C^{1}([0,T];W)$. In section 7 we discuss the behavior of the exponential constant appearing in Theorem 4.1 for evolution times seen in fracture experiments. Since we are approximating the solution of an ODE on a Banach space the proof of Theorem 4.1 will follow from the Lipschitz continuity of the force $\mathcal{L}^{\epsilon}(\boldsymbol{u})$ with respect to the $L^{2}$ norm. The proof is given in the following two sections.

4.2.1 Truncation error analysis and consistency

The results in this section follows the same steps as in [Jha and Lipton, 2017] and therefore we will just highlight the major steps. We can write the discrete evolution equation for $(\boldsymbol{e}^{k}_{h}(\boldsymbol{u})=\boldsymbol{r}_{h}(\boldsymbol{u}^{k})-\boldsymbol{u}^{k}_{h},\boldsymbol{e}^{k}_{h}(\boldsymbol{v})=\boldsymbol{r}_{h}(\boldsymbol{v}^{k})-\boldsymbol{v}^{k}_{h})$ as follows

[TABLE]

where consistency error terms ${\boldsymbol{\tau}}^{k}_{h}(\boldsymbol{u}),{\boldsymbol{\tau}}^{k}_{h}(\boldsymbol{v}),{\boldsymbol{\sigma}}^{k}_{h}(\boldsymbol{u})$ are given by

[TABLE]

When $\boldsymbol{u},\boldsymbol{v}$ are $C^{2}$ in time, we easily see that

[TABLE]

To estimate ${\boldsymbol{\sigma}}^{k}_{h}(\boldsymbol{u})$, we recall the Lipschitz continuity property of the peridynamic force in the $L^{2}$ norm, see Equation 23. This leads us to

[TABLE]

where the constant $L$ is defined in Equation 24.

We now state the consistency of this approach.

Lemma 4.2**.**

**Consistency

Let $\tau$ be given by**

[TABLE]

then the approach is consistent in that

[TABLE]

where

[TABLE]

4.2.2 Stability analysis

In equation for $\boldsymbol{e}^{k}_{h}(\boldsymbol{u})$, we take $\tilde{\boldsymbol{u}}=\boldsymbol{e}^{k+1}_{h}(\boldsymbol{u})$. We have

[TABLE]

which implies

[TABLE]

Similarly, we can show

[TABLE]

We have from Equation 23

[TABLE]

After adding Equation 59 and subsubsection 4.2.2, and substituting Equation 61, we get

[TABLE]

where $\tau$ is defined in Equation 56. Since $e^{k}=||\boldsymbol{e}^{k}_{h}(\boldsymbol{u})||+||\boldsymbol{e}^{k}_{h}(\boldsymbol{v})||$, we can show, assuming $L/\epsilon^{2}\geq 1$,

[TABLE]

Substituting for $e^{k}$ recursively in the equation above, we get

[TABLE]

Noting that

[TABLE]

and $(1+a\Delta t)^{k}\leq\exp[ka\Delta t]\leq\exp[Ta]$ for $a>0$, we have

[TABLE]

This implies

[TABLE]

By substituting above equation in Equation 50, we get stability of the scheme:

Lemma 4.3**.**

**Stability

**

[TABLE]

After taking sup over $k\leq T/\Delta t$ and substituting the bound on $\tau$ from Lemma 4.2, we get the desired result and proof of Theorem 4.1 is complete.

We now consider a stronger notion of stability for the linearized peridynamics model.

4.3 Linearized peridynamics and energy stability

In this section, we linearize the peridynamics model and obtain a CFL like stability condition. For problems where strains are small, the stability condition for the linearized model is expected to apply to the nonlinear model. The slope of peridynamics potential $f$ and $g$ are constant for sufficiently small strain and therefore for small strain the nonlinear model behaves like a linear model.

In section 2, linearization gives

[TABLE]

The corresponding bilinear form is denoted as $a^{\epsilon}_{T,l}$ and is given by

[TABLE]

Similarly, linearization of $\mathcal{L}^{\epsilon}_{D}$ in section 2 gives

[TABLE]

The associated bilinear form is given by

[TABLE]

The total force after linearization is

[TABLE]

and the bilinear operator associated to $\mathcal{L}^{\epsilon}_{l}$ is given by

[TABLE]

We have

[TABLE]

We now discuss the stability of the FEM approximation to the linearized problem. Let $\boldsymbol{u}^{k}_{l,h}$ denote the approximate solution satisfying, for $k\geq 1$,

[TABLE]

and, for $k=0$,

[TABLE]

The following notation will be used to define the discrete energy at each time step $k$

[TABLE]

We also define

[TABLE]

We introduce the discrete energy associated with $\boldsymbol{u}^{k}_{l,h}$ at time step $k$ as follows

[TABLE]

Following [Theorem 4.1, [Karaa, 2012]], the stability of central difference scheme is given by

Theorem 4.4**.**

**Energy Stability of the Central difference approximation of linearized peridynamics

Let $\boldsymbol{u}^{k}_{l,h}$ be the approximate solution of Equation 70 and Equation 71. In the absence of body force $\boldsymbol{b}(t)=0$ for all $t$, if $\Delta t$ satisfies the CFL like condition**

[TABLE]

then the discrete energy is positive and we have the stability

[TABLE]

We skip the proof of above theorem as it is straight-forward extension of Theorem 5.2 in [Jha and Lipton, 2017].

5 Numerical experimants

In this section, we present numerical simulations that are consistent with the theoretical a-priori bound on the convergence rate. We also compare the peridynamic energy of the material softening zone and the classic Griffith’s fracture energy of linear elastic fracture mechanics.

We consider Plexiglass at room temperature and specify the density $\rho=1200\,kg/m^{3}$, bulk modulus $K=25\,GPa$, Poisson’s ratio $\nu=0.245$, and critical energy release rate $G_{c}=500\,Jm^{-2}$. The pairwise interaction and the hydrostatic interaction are characterized by potentials $f(r)=c(1-\exp[-\beta r^{2}])$ and $g(r)=\bar{C}r^{2}/2$ respectively. Here we have used a quadratic hydrostatic interaction potential. The influence function is $J(r)=1-r$. Since pairwise potential $f$ is symmetric for positive and negative strain the critical strain is given by $S_{c}(\boldsymbol{y},bx)=\frac{\pm\bar{r}}{\sqrt{\boldsymbol{y}-\boldsymbol{x}}}$, where $\pm\bar{r}$ is the inflection point of $f(r)$ given by $\bar{r}=\frac{1}{\sqrt{\beta}}$. Following equations 94, 95, and 97 of [Lipton et al., 2018b], the relation between peridynamic material parameters and Lamé constants $(\lambda,\mu)$ and critical energy release rate $G_{c}$ can be written as (for 2-d)

[TABLE]

where $M_{J}$ is given by

[TABLE]

By solving Equation 75, we get $c=4712.4$, $\bar{C}=-1.7349\times 10^{11}$, $\beta=1.5647\times 10^{8}$.

We consider a 2-d domain $D=[0,0.1\,m]^{2}$ (with unit thickness in third direction) with vertical crack of length $0.02\,m$. The boundary conditions are described in Figure 3. The simulation time is $T=40\,\mu s$ and the time step is $\Delta t=0.004\,\mu s$. We consider two horizons $8\,mm$ and $4\,mm$. We run simulations for mesh sizes $h=2,1,0.5\,mm$. We consider the central difference time discretization described by Equation 41 on a uniform mesh consisting of linear triangle elements. Second order quadrature approximation is used in the simulation for each triangle element. To reduce the load on memory and to avoid matrix-vector multiplication at each time step, we approximate the mass matrix by diagonal mass matrix using lumping (row-sum) technique. Suppose exact mass matrix is $\boldsymbol{M}=[m_{ij}]$ where $m_{ij}$ is the element of $\boldsymbol{M}$ corresponding to $i^{\text{th}}$ row and $j^{\text{th}}$ column, then we approximate $\boldsymbol{M}$ by diagonal matrix $\hat{\boldsymbol{M}}=[\hat{m}_{ij}]$ where $\hat{m}_{ii}=\sum_{j}m_{ij}$ and $\hat{m}_{ij}=0$ if $j\neq i$.

5.1 Convergence rate

To compute convergence rate numerically we proceed as follows: consider a fixed horizon $\epsilon$ and three different mesh sizes $h_{1},h_{2},h_{3}$ such that $r=h_{1}/h_{2}=h_{2}/h_{3}$. Let $\boldsymbol{u}_{1},\boldsymbol{u}_{2},\boldsymbol{u}_{3}$ be approximate solutions corresponding to meshes of size $h_{1},h_{2},h_{3}$, and let $\boldsymbol{u}$ be the exact solution. We write the error as $||\boldsymbol{u}_{h}-\boldsymbol{u}||=Ch^{\alpha}$ for some constant $C$ and $\alpha>0$, to get

[TABLE]

From above two equations, it is easy to see that the rate of convergence $\alpha$ is

[TABLE]

The convergence result for horizons $\epsilon=8\,mm$ and $\epsilon=4\,mm$ is shown in Figure 4. In the simulation we have considered second order approximation of integration using quadrature points. The simulations show a rate of convergence that agrees with the a-priori estimates given in Theorem 4.1.

5.2 Fracture energy of crack zone

The extent of damage at material point $\boldsymbol{x}$ is given by the function $Z(\boldsymbol{x})$

[TABLE]

The crack zone is defined as set of material points which have $Z>1$. We compute the peridynamic energy of crack zone and compare it with the Griffith’s fracture energy. For a crack of length $l$, the Griffith’s fracture energy (G.E.) will be $G.E.=G_{c}\times l$. The peridynamic fracture energy (P.E.) associated with the material softening zone is given by

[TABLE]

where $\mathcal{W}^{\epsilon}(S(\boldsymbol{y},\boldsymbol{x};\boldsymbol{u}))$ is the bond-based potential, see Equation 2 and $\mathcal{V}^{\epsilon}(\theta(\boldsymbol{x},t;\boldsymbol{u}))$ is the hydrostatic interaction potential, see Equation 4.

In Figure 5 classical fracture energy and peridynamic fracture energy is shown at different crack length. The error in both energies at different times is shown in Figure 6. The agreement between two energies is good. The damage profile at time $30\,\mu s$ and $40\,\mu s$ is shown in Figure 7. At each node, the damage function $Z$ is computed by treating edges between mesh nodes as bonds. In addition to the damage plots, we show velocity profile at $30\,\mu s$ and $40\,\mu s$ in Figure 8. In Figure 9 we show the plot of the $xx$ component of symmetric gradient of the displacement. Here the region for which the magnitude of the strain is greater than a multiple of the critical strain is the yellow region. It is seen that the high strain region surrounds the crack.

As the crack is propagating vertically it is seen that the high strain region is next to the crack i

6 Lipschitz continuity of peridynamic force and higher temporal regularity of solutions

In this section, we prove Theorem 2.1 and Theorem 2.3. Here $\boldsymbol{u}\in W\subset H^{2}(D;\mathbb{R}^{d})$ and the $||\boldsymbol{u}||_{2}$ norm is given by

[TABLE]

6.1 Proof of Lipschitz continuity with respect to the $\|\cdot\|_{2}$ norm

We assume that the potential function $f$ satisfies $C^{f}_{i}<\infty$ for $i=0,1,2,3,4$. Recall that $C^{f}_{0}=\sup_{r}|f(r)|$ and $C^{f}_{i}=\sup_{r}|f^{(r)}(r)|$ for $i=1,...,4$. $C^{g}_{i}$ is defined similarly for $i=0,1,...,4$. If the potential function $g$ is a convex-concave function then we can assume $C^{g}_{i}<\infty$ for $i=0,1,2,3,4$. In what follows, we will prove the Theorem 2.1 for convex-concave type $g$. If $g$ is a purely a quadratic function the proof follows easily using only a subset of the estimates proved in this section.

Let $\boldsymbol{u},\boldsymbol{v}\in W$. Using the triangle inequality, we get

[TABLE]

where $\mathcal{L}^{\epsilon}_{T}$ and $\mathcal{L}^{\epsilon}_{D}$ is given by section 2 and section 2.

We first write the peridynamic force $\mathcal{L}^{\epsilon}_{T}(\boldsymbol{u})(\boldsymbol{x})$ as follows

[TABLE]

where we substituted $\partial_{S}f(\sqrt{|\boldsymbol{y}-\boldsymbol{x}|}S(\boldsymbol{y},\boldsymbol{x};\boldsymbol{u}))=\sqrt{|\boldsymbol{y}-\boldsymbol{x}|}f^{\prime}(\sqrt{|\boldsymbol{y}-\boldsymbol{x}|}S(\boldsymbol{y},\boldsymbol{x};\boldsymbol{u}))$. The form of the peridynamic force described above is the same as the one given in [Section 6, [Jha and Lipton, 2017]]. We apply Theorem 3.1 in [Jha and Lipton, 2017]] to show

[TABLE]

and

[TABLE]

Next we analyze $||\mathcal{L}^{\epsilon}_{D}(\boldsymbol{u})-\mathcal{L}^{\epsilon}_{D}(\boldsymbol{v})||_{2}$. We define new terms to simplify the calculations. For ${\boldsymbol{\xi}}\in H_{1}(\mathbf{0})$, we set

[TABLE]

Similar notations hold if we exchange $\boldsymbol{x}$, ${\boldsymbol{\xi}}\in H_{1}(\mathbf{0})$, and $\boldsymbol{u}\in W$ by $\boldsymbol{y}$, $\boldsymbol{\eta}\in H_{1}(\mathbf{0})$, and $\boldsymbol{v}\in W$ respectively We will also encounter various moments of the influence function $J$ therefore we define following moments

[TABLE]

Recall that $J(|{\boldsymbol{\xi}}|)=0$ for ${\boldsymbol{\xi}}\notin H_{1}(\mathbf{0})$ and $0\leq J(|{\boldsymbol{\xi}}|)\leq M$ for ${\boldsymbol{\xi}}\in H_{1}(\mathbf{0})$. The boundary function $\omega$ is assumed to satisfy

[TABLE]

We choose finite constants $C_{\omega_{1}}$ and $C_{\omega_{2}}$ such that

[TABLE]

We now collect the following estimates which will be used to estimate $||\mathcal{L}^{\epsilon}_{D}(\boldsymbol{u})-\mathcal{L}^{\epsilon}_{D}(\boldsymbol{v})||_{2}$.

Lemma 6.1**.**

Let $\boldsymbol{u},\boldsymbol{v}\in W$, for any $\boldsymbol{\eta}\in H_{1}(\mathbf{0})$ and $\delta\leq 2\epsilon$, we have

[TABLE]

*Here $\nabla$ in all the equations above is with respect to $\boldsymbol{x}$. The constants $C_{e_{1}},C_{e_{2}}$ are the constants associated to the Sobolev embedding property of space $H^{2}(D;\mathbb{R}^{d})$, see Equation 18 and Equation 19. *

Proof 6.2**.**

Using the notation given in subsection 6.1, we write $\theta(\boldsymbol{x};\boldsymbol{u})$ as

[TABLE]

On noting that $|\bar{\boldsymbol{u}}_{\boldsymbol{\xi}}(\boldsymbol{x})|\leq 2||\boldsymbol{u}||_{\infty}$ and $||\boldsymbol{u}||_{\infty}\leq C_{e_{1}}||\boldsymbol{u}||_{2}$, we easily see that

[TABLE]

In the rest of the proof, we will let $\boldsymbol{y}=\boldsymbol{x}+\delta\boldsymbol{\eta}$, where $0\leq\delta\leq 2\epsilon$ and $\boldsymbol{\eta}\in H_{1}(\mathbf{0})$.

To show Equation 88, we first introduce an important identity which will be used frequently. Let $p({\boldsymbol{\xi}})$ be some function of ${\boldsymbol{\xi}}$, and $\alpha,C\in\mathbb{R}$ then

[TABLE]

where we used the inequality $ab\leq a^{2}/2+b^{2}/2$ in the first step, and definition of $\bar{J}_{\alpha}$ and symmetry of terms in second step.

From the expression of $\theta(\boldsymbol{y};\boldsymbol{u})$ we can show

[TABLE]

We now prove the bound Equation 89. Taking the gradient of $\theta(\boldsymbol{y};\boldsymbol{u})$, with respect to $\boldsymbol{x}$, noting that $\boldsymbol{y}=\boldsymbol{x}+\delta\boldsymbol{\eta}$, we get

[TABLE]

We can show using inequality 6.2 and the estimates $\int_{D}|\nabla\bar{\boldsymbol{u}}_{\boldsymbol{\xi}}(\boldsymbol{y})|^{2}d\boldsymbol{x}\leq 4||\nabla\boldsymbol{u}||^{2}\leq 4||\boldsymbol{u}||_{2}^{2}$, $|\nabla\omega(\boldsymbol{y}+\epsilon{\boldsymbol{\xi}})|\leq C_{\omega_{1}}$, $\int_{D}|\bar{\boldsymbol{u}}_{\boldsymbol{\xi}}(\boldsymbol{y})|^{2}d\boldsymbol{x}\leq 4||\boldsymbol{u}||_{2}^{2}$, to conclude

[TABLE]

We now show Equation 90. We will use Equation 87 and Equation 89, and proceed as follows

[TABLE]

To prove Equation 91 we note expression of $\nabla\theta(\boldsymbol{y};\boldsymbol{u})$ in 6.2 and inequality $(a+b)^{4}\leq 8a^{4}+8b^{4}$ and 6.2 to get

[TABLE]

Application of Fubini’s theorem gives

[TABLE]

Using the Sobolev embedding property, $||\boldsymbol{u}||_{\infty}\leq C_{e_{1}}||\boldsymbol{u}||_{2}$ and $||\nabla\boldsymbol{u}||_{L^{4}}\leq C_{e_{2}}||\boldsymbol{u}||_{2}$, we obtain

[TABLE]

The estimate Equation 92 follows by combining estimates Equation 87 and Equation 91. It now remains to show Equation 93. From expression of $\nabla\theta(\boldsymbol{y};\boldsymbol{u})$ in 6.2, we have

[TABLE]

Using the equation above we can show

[TABLE]

The terms $|\bar{\boldsymbol{u}}_{\boldsymbol{\xi}}(\boldsymbol{y})|^{2}$, $|\nabla\bar{\boldsymbol{u}}_{\boldsymbol{\xi}}(\boldsymbol{y})|^{2}$, and $|\nabla^{2}\bar{\boldsymbol{u}}_{\boldsymbol{\xi}}(\boldsymbol{y})|^{2}$ are bounded by $2(|\boldsymbol{u}(\boldsymbol{y}+\epsilon{\boldsymbol{\xi}})|^{2}+|\boldsymbol{u}(\boldsymbol{y})|^{2})$, $2(|\nabla\boldsymbol{u}(\boldsymbol{y}+\epsilon{\boldsymbol{\xi}})|^{2}+|\nabla\boldsymbol{u}(\boldsymbol{y})|^{2})$, and $2(|\nabla^{2}\boldsymbol{u}(\boldsymbol{y}+\epsilon{\boldsymbol{\xi}})|^{2}+|\nabla^{2}\boldsymbol{u}(\boldsymbol{y})|^{2})$ respectively. Therefore, we have

[TABLE]

*and this completes the proof of lemma. *

**Estimating $||\mathcal{L}^{\epsilon}_{D}(\boldsymbol{u})-\mathcal{L}^{\epsilon}_{D}(\boldsymbol{v})||$: **We apply the notation described in subsection 6.1, and write $\mathcal{L}^{\epsilon}_{D}(\boldsymbol{u})(\boldsymbol{x})$ as follows

[TABLE]

Using the formula above and from the expression for $\theta$ we can easily show

[TABLE]

where $L_{1}=4C^{g}_{2}\bar{J}_{0}^{2}$.

**Estimating $||\nabla\mathcal{L}^{\epsilon}_{D}(\boldsymbol{u})-\nabla\mathcal{L}^{\epsilon}_{D}(\boldsymbol{v})||$: **Taking the gradient of Equation 107 gives

[TABLE]

where we have denoted the first and second term as $G_{1}(\boldsymbol{u})(\boldsymbol{x})$ and $G_{2}(\boldsymbol{u})(\boldsymbol{x})$ for convenience. On using the triangle inequality we get

[TABLE]

From the expression of $G_{1}(\boldsymbol{u})$, we have

[TABLE]

Let

[TABLE]

and we get

[TABLE]

Note that

[TABLE]

Therefore, from Equation 110,

[TABLE]

where we have added and subtracted $g^{\prime\prime}(\theta(\boldsymbol{y};\boldsymbol{u}))\nabla\theta(\boldsymbol{y};\boldsymbol{v})$ and used the fact that $g^{\prime\prime}(r)\leq C^{g}_{2}$ and $|g^{\prime\prime}(r_{1})-g^{\prime\prime}(r_{2})|\leq C^{g}_{3}|r_{1}-r_{2}|$. We use the estimate on $p_{1}$ and proceed as follows

[TABLE]

where we denote $\boldsymbol{x}+\epsilon{\boldsymbol{\xi}}$ as $\boldsymbol{y}$. We apply inequality Equation 89 and Equation 90 of Lemma 6.1 to obtain

[TABLE]

where we have grouped all the constant factors together and denote their product by $L_{2}$. Substituting these estimates in to Equation 111 gives

[TABLE]

where we have introduced the new constant $L_{3}$.

The formula for $G_{2}(\boldsymbol{u})$ is similar to $\mathcal{L}^{\epsilon}_{D}(\boldsymbol{u})$ and therefore we have

[TABLE]

Collecting results, we have shown

[TABLE]

where we have introduced new constant $L_{4}$.

**Estimating $||\nabla^{2}\mathcal{L}^{\epsilon}_{D}(\boldsymbol{u})-\nabla^{2}\mathcal{L}^{\epsilon}_{D}(\boldsymbol{v})||$: **Taking the gradient of subsection 6.1, gives

[TABLE]

It is easy to see that estimate on $||H_{2}(\boldsymbol{u})-H_{2}(\boldsymbol{v})||$ and $||H_{3}(\boldsymbol{u})-H_{3}(\boldsymbol{v})||$ is similar to the estimate for $||G_{1}(\boldsymbol{u})-G_{1}(\boldsymbol{v})||$. Thus, from subsection 6.1, we have

[TABLE]

Also the estimate for $||H_{4}(\boldsymbol{u})-H_{4}(\boldsymbol{v})||$ is similar to the estimate for $||G_{2}(\boldsymbol{u})-G_{2}(\boldsymbol{v})||$ and we conclude

[TABLE]

We now work on $||H_{1}(\boldsymbol{u})-H_{1}(\boldsymbol{v})||$. From expression of $H_{1}(\boldsymbol{u})(\boldsymbol{x})$ in subsection 6.1, we can easily get following

[TABLE]

Let $p_{2}(\boldsymbol{y})$, where $\boldsymbol{y}=\boldsymbol{x}+\epsilon{\boldsymbol{\xi}}$ and $\nabla$ is with respect to $\boldsymbol{x}$, is given by

[TABLE]

We then have

[TABLE]

Note that

[TABLE]

We add and subtract terms to the equation above to get

[TABLE]

Using inequalities $|g^{\prime\prime}(r)|\leq C^{g}_{2}$, $|g^{\prime\prime\prime}(r)|\leq C^{g}_{3}$, $|g^{\prime\prime}(r_{1})-g^{\prime\prime}(r_{2})|\leq C^{g}_{3}|r_{1}-r_{2}|$, $|g^{\prime\prime\prime}(r_{1})-g^{\prime\prime\prime}(r_{2})|\leq C^{g}_{4}|r_{1}-r_{2}|$, and $|\boldsymbol{a}\mathbf{\otimes}\boldsymbol{a}-\boldsymbol{c}\mathbf{\otimes}\boldsymbol{c}|\leq(|\boldsymbol{a}|+|\boldsymbol{c}|)|\boldsymbol{a}-\boldsymbol{c}|$, and the fact that $\theta(\boldsymbol{y};\boldsymbol{u})-\theta(\boldsymbol{y};\boldsymbol{v})=\theta(\boldsymbol{y};\boldsymbol{u}-\boldsymbol{v})$, we have

[TABLE]

Taking the square of the above equation and using $(\sum_{i=1}^{5}a_{i})^{2}\leq 5\sum_{i=1}^{5}a^{2}_{i}$ gives

[TABLE]

We now estimate each term using Lemma 6.1 as follows. Applying Hölder inequality and inequality Equation 91 of Lemma 6.1 to get

[TABLE]

Similarly,

[TABLE]

Using Equation 92 of Lemma 6.1, we get

[TABLE]

For $I_{4}$, we use inequality Equation 93 to get

[TABLE]

In $I_{5}$, we use Equation 87 and Equation 93 to get

[TABLE]

After collecting results, we can find a constant $L_{5}$ such that we have

[TABLE]

We substitute Equation 120 into Equation 119 to show

[TABLE]

where we have introduced the new constant $L_{6}$.

We combine the estimates on $H_{1},H_{2},H_{3},H_{4}$, introducing a new constant $L_{7}$, and get

[TABLE]

On adding the estimates Equation 108, Equation 113, Equation 122 it is evident that the proof of Theorem 2.1 is complete.

6.2 Proof of higher temporal regularity

In this section we prove that the peridynamic evolutions have higher regularity in time for body forces that that are differentiable in time. To see this we take the time derivative of Equation 11 to get a second order differential equation in time for $\boldsymbol{v}=\dot{\boldsymbol{u}}$ given by

[TABLE]

where $Q(\boldsymbol{v};\boldsymbol{u})$ is an operator that depends on the solution $\boldsymbol{u}$ of Equation 11 and acts on $\boldsymbol{v}$. It is given by,

[TABLE]

where

[TABLE]

and

[TABLE]

Clearly, for $\boldsymbol{u}$ fixed, the form $Q(\boldsymbol{v};\boldsymbol{u})$ acts linearly on $\boldsymbol{v}$ which implies that the equation for $\boldsymbol{v}$ is a linear nonlocal equation. The linearity of $Q(\boldsymbol{v};\boldsymbol{u})$ implies Lipschitz continuity for $\boldsymbol{v}\in W$ as stated below.

Theorem 6.3**.**

**Lipschitz continuity of $Q$

Let $\boldsymbol{u}\in W$ be any given field. Then for all $\boldsymbol{v},\boldsymbol{w}\in W$, we have**

[TABLE]

where the constant $L_{8}$ does not depend on $\boldsymbol{u},\boldsymbol{v},\boldsymbol{w}$. This gives for all $\boldsymbol{v}\in W$ the upper bound,

[TABLE]

The proof follows the same steps used in the proof of Theorem 2.1.

If $\boldsymbol{u}$ is a peridynamic solution such that $\boldsymbol{u}\in C^{2}(I_{0};W)$ then we have for all $t\in I_{0}$ the inequality

[TABLE]

Note that Lipschitz continuity of $\dot{\boldsymbol{u}}(t)$ stated in Theorem 2.2 implies $\lim_{t\rightarrow 0^{\pm}}\partial^{2}_{tt}\boldsymbol{u}(\boldsymbol{x},t)=\partial^{2}_{tt}\boldsymbol{u}(\boldsymbol{x},0)$. We now demonstrate that $\boldsymbol{v}(\boldsymbol{x},t)=\partial_{t}\boldsymbol{u}(\boldsymbol{x},t)$ is the unique solution of following initial boundary value problem.

Theorem 6.4**.**

**Initial value problem for $\boldsymbol{v}(\boldsymbol{x},t)$

Suppose the initial data and righthand side $\boldsymbol{b}(t)$ satisfy the hypothesis of Theorem 2.2 and we suppose further that $\dot{\boldsymbol{b}}(t)$ exists and is continuous in time for $t\in I_{0}$ and $\sup_{t\in I_{0}}||\dot{\boldsymbol{b}}(t)||_{2}<\infty$. Then $\boldsymbol{v}(\boldsymbol{x},t)$ is the unique solution to the initial value problem $\boldsymbol{v}(\boldsymbol{x},0)=\boldsymbol{v}_{0}(\boldsymbol{x})$, $\partial_{t}\boldsymbol{v}(\boldsymbol{x},0)=\partial^{2}_{tt}\boldsymbol{u}(\boldsymbol{x},0)$,**

[TABLE]

$\boldsymbol{v}\in C^{2}(I_{0};W)$* and*

[TABLE]

Theorem 2.3 now follows immediately from Theorem 6.4 noting that $\partial_{t}\boldsymbol{u}(\boldsymbol{x},t)=\boldsymbol{v}(\boldsymbol{x},t)$ together with Equation 129 and Equation 131. The proof of Theorem 6.4 follows from the Lipschitz continuity Equation 127 and the Banach fixed point theorem as in [Brezis, 1983].

7 Conclusions

In this article, we have provided a-priori error estimates for finite element approximations to nonlocal state based peridynamic fracture models. We have shown that the convergence rate applies even over time intervals for which the material is softening over parts of the computational domain. The results are established for two different classes of state-based peridynamic forces. The convergence rate of the approximation is of the form $C(\Delta t+h^{2}/\epsilon^{2})$ where the constant $C$ depends on $\epsilon$ and the $H^{2}$ norm of the solution and its time derivatives. For fixed $\Delta t$ numerical simulations for Plexiglass show that the error decreases at the rate of $h^{2}$ at $40\mu$-sec into the simulation. The simulations were carried out in parallel using 20 threads on an workstation with single Intel Xeon processor and with 32 GB of RAM. We anticipate similar convergence rates for longer times on bigger parallel machines.

We reiterate that the a-priori error estimates account for the possible appearance of nonlinearity anywhere in the computational domain. On the other hand numerical simulation and independent theoretical estimates show that the nonlinearity concentrates along “fat” cracks of finite length and width equal to $\epsilon$, see [Lipton, 2016, Lipton, 2014]. Moreover the remainder of the computational domain is seen to behave linearly and to leading order can be modeled as a linear elastic material up to an error proportional to $\epsilon$, see [Proposition 6, [Jha and Lipton 2018b]]. Future work will use these observations to focus on adaptive implementation and a-posteriori estimates.