Finite element convergence for statebased peridynamic fracture models^{†}^{†}thanks: Funding: This material is based upon work supported by the U. S. Army Research Laboratory and the U. S. Army Research Office under contract/grant number W911NF1610456.
Abstract
We establish the apriori 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 multiwell 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 forcestrain relation, or a multiwell type that can be associated with material degradation and cavitation. We first show the wellposedness of the peridynamic formulation and that peridynamic evolutions exist in the Sobolev space . We show that the finite element approximations converge to the solutions uniformly as measured in the mean square norm. For linear continuous finite elements the convergence rate is shown to be , where is the size of horizon, is the mesh size, and is the size of time step. The constants and are independent of and and may depend on through the norm of the exact solution. We demonstrate the stability of the semidiscrete 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.
Key words. Nonlocal fracture models, peridynamic, state based peridynamic, numerical analysis, finite element approximation
AMS subject classifications. 34A34, 34B10, 74H55, 74S05
1 Introduction
In this work, we study the statebased peridynamic theory and obtain an apriori 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 statebased 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 bondstrain. Let , for , be the material domain. For a displacement field , the bondstrain between two material points is given by
(1) 
Let be the size of horizon and be the neighborhood of a material point . For pairwise interaction, we assume the following form of pairwise interaction potential
(2) 
where is the influence function. We assume where for and for . The potential , see LABEL:fig:f_a, is assumed to be convex for small strains and becomes concave for larger strains. In the widely used prototypical microelastic brittle (PMB) peridynamic material, the strain vs force profile is linear up to some critical strain and is zero for any strain above . In contrast, the peridynamic force given by , is linear near zero strain and as the strain gets larger and reaches the critical strain, () for positive (negative) strain, the bond starts to soften, see LABEL:fig:f_b. For a given potential function , the critical strain is given by and where are the inflection points of the potential function as shown in LABEL:fig:f_a.
The spherical or hydrostatic strain at material point is given by
(3) 
where is the volume of unit ball in dimension . The potential for hydrostatic interaction is of the form
(4) 
where is the potential function associated to hydrostatic strain. Here can be of two types: 1) a quadratic function with only one well at zero strain, and 2) a convexconcave function with a wells at the origin and at , see LABEL:fig:g_a. If is assumed to be quadratic then the force due to spherical strain is linear. If is a multiwell potential, the material softens as the hydrostatic strains exceeds critical value. For the convexconcave type , the critical values are and beyond which the force begins to soften is related to the inflection point and of as follows
(5) 
The critical compressive hydrostatic strain where the force begins to soften for negative hydrostatic strain is chosen much larger in magnitude than , i.e. .
The finite element approximation has been applied to peridynamic fracture, however there remains a paucity of literature addressing the rigorous apriori 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 that are twice differentiable in time, see LABEL:thm:existence. 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 evolutions converge in the limit of vanishing nonlocality 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 and for a bounded time interval.
In this paper we obtain an apriori 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 , see section LABEL:consistent and section LABEL:stable. The bound on the error is uniform in time and is given by , where the constants and are independent of and mesh size , see LABEL:thm:convergence. A more elaborate discussion of the apriori bound is presented in section LABEL:s:convergence. For the linearized model we obtain a stability condition on , LABEL:thm:cflcondition, 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 error for the displacement and velocity that are uniform in time and of the form , where constants the and 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 as opposed to . 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 LABEL:s:numerical. 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 LABEL:sec:_softzone. 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 apriori convergence rates for the FEM given by LABEL:thm:convergence 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 apriori error involves in the denominator and in many cases is chosen small. However typical dynamic fracture experiments last only hundreds of microseconds and the apriori error is controlled by the product of simulation time multiplied by . So for material properties characteristic of Plexiglass and of size , the apriori estimates predict a relative error of for simulations lasting around 100 microseconds. We point out that the apriori 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 , 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 , see [Proposition 6, [Jha and Lipton 2018b]]. Future work will use these observations to focus on adaptive implementation and aposteriori estimates.
Aposteriori 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. Aposteriori 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 LABEL:s:nonlocal_dynamics 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 LABEL:s:fem, we consider the finite element discretization. We prove the stability of a semidiscrete approximation in LABEL:ss:semifem. In LABEL:s:fullfem, we analyze the fully discrete approximation and obtain an apriori bound on errors. The stability of the fully discrete approximation linearized peridynamic force is shown in LABEL:ss:fullfemstability. We present our numerical experiments in LABEL:s:numerical. Proofs of the Lipschitz bound on the peridynamic force and higher temporal regularity of solutions is provided in LABEL:s:proofs. In LABEL:s:conclusions we present our conclusions.
We conclude the introduction by listing the notation used throughout the paper. We denote material domain as where for . Points and vectors in are denoted as bold letters. Some of the key notations are as follows
Time domain  
Size of horizon  
Density  
Horizon of , a ball of radius centered at  
Volume of unit ball in dimension  
Boundary function defined on taking value in the interior and smoothly decaying to as approaches  
Displacement field defined over . We may also use notation to denote field defined over just  
Initial condition on displacement  
Body force defined over  
The unit vector pointing from a point to the point  
Bond strain . We may also use if is a filed defined over just  
Spherical or hydrostatic strain. We may also use if is a filed defined over just  
Critical bond strain  
Critical hydrostatic strain  
Influence function where is integrable with for and for  
Moment of function over with weight  
Potential functions for pairwise and statebased interaction  
Pairwise and statebased potential energy density  
Total peridynamic potential energy at time  
Total dynamic energy at time  
total peridynamic force, pairwise peridynamic force, and statebased peridynamic force respectively  
Nonlinear operator where are vector fields over  
Nonlinear pairwise and statebased operator  
norm over , norm over , and Sobolev norm over (for ) respectively  
Size of mesh and size of time step  
Triangulation of given by triangular/tetrahedral elements  
Continuous piecewise linear interpolation operator on  
Space of functions in such that trace of function is zero on boundary , i.e.  
Space of continuous piecewise linear interpolations on  
Interpolation function of mesh node  
Finite element projection of onto  
Total error in mean square norm at time step  
Approximate displacement and velocity field at time step  
Exact displacement and velocity field at time step 
2 Equation of motion, existence, uniqueness, and higher regularity
We assume to be an open set with boundary. To enforce zero displacement boundary conditions at and to insure a well posed evolution we introduce the boundary function . This function is introduced as a factor into the potentials and . Here the boundary function takes value in the interior of domain and is zero on the boundary. We assume and in our analysis. The hydrostatic strain is modified to include the boundary and is given by
(6) 
The peridynamic potentials LABEL:eq:tensilepot and LABEL:eq:hydropot are modified to see the boundary follows
(7)  
(8) 
We assume that potential function is at least 4 times differentiable and satisfies following regularity condition:
(9) 
If potential function is convexconcave type then we assume that satisfies same regularity condition as . We denote constants , for , similar to above.
The total potential energy at time is given by
(10)  
where potential and are described above. The material is assumed to be homogeneous and the density is given by . The applied body force is denoted by . We define the Lagrangian
here is the velocity given by the time derivative of . 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
(11) 
where
(12) 
is the peridynamic force due to bondbased interaction and is given by
(13) 
and is the peridynamic force due to statebased interaction and is given by
(14) 
The dynamics is complemented with the initial data
(15) 
We prescribe zero Dirichlet boundary condition on the boundary
(16) 
We extend the zero boundary condition outside to whole . In our analysis we will assume mass density without loss of generality.
2.1 Existence of solutions and higher regularity in time
We recall that the space is the closure in the norm of the functions that are infinitely differentiable with compact support in . For suitable initial conditions and body force we show that solutions exist in
(17) 
where is the trace of the function on the boundary of . We will assume that is extended by zero outside . We first exhibit the Lipschitz continuity property and boundedness of the peridynamic force for displacements in . 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 when is a domain.

From Theorem 2.72 of [Demengel, 2012], there exists a constant independent of such that
(18) 
Further application of standard embedding theorems
(e.g., Theorem 2.72 of [Demengel, 2012]), shows there exists a constant independent of such that(19) for any such that when and when .
We have following result which show the Lipschitz continuity property of a peridynamic force .
Theorem 1.
Lipschitz continuity of peridynamic force
Let be a convexconcave function satisfying for and let either be a quadratic function, or be a convexconcave function with for . Also, let boundary function be such that and . Then, for any , we have
(20) 
where constant does not depend on and . Also, for , we have
(21) 
where constant does not depend on and .
Now let be any positive number, a straightforward application of [Theorem 3.2, [Jha and Lipton, 2017]] gives:
Theorem 2.
Existence and uniqueness of solutions over finite time intervals
Let , , and satisfy the hypothesis of LABEL:thm:lipschitzproperty. For any initial condition , time interval , and right hand side continuous in time for such that satisfies , there is a unique solution of peridynamic equation LABEL:eq:equationofmotion. Also, and are Lipschitz continuous in time for .
We can also show higher regularity in time of evolutions under suitable assumptions on the body force:
Theorem 3.
Higher regularity
Suppose the initial data and righthand side satisfy the hypothesis of LABEL:thm:existence and suppose further that exists and is continuous in time for and . Then
and
(22) 
where is a positive constant independent of .
The proofs of Theorems LABEL:thm:lipschitzproperty, and LABEL:thm:higherregularity are given in LABEL:s:proofs. For future reference, we note that for any , we have
(23) 
where constant is given by
(24) 
and .
2.2 Weak form
We multiply LABEL:eq:equationofmotion by a test function in and integrate over to get
(25) 
We have the following integration by parts formula:
Lemma 4.
For any we have
(26) 
where
(27) 
and
(28) 
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
(29) 
Total dynamic energy: We define the total dynamic energy as follows
(30) 
where is defined in LABEL:eq:totalenergy. Time derivative of total energy satisfies
(31) 
Remark. It is readily verified that the peridynamic force and energy are bounded for all functions in . Here the bound on the force follows from the Lipschitz property of the force in , see, LABEL:eq:lipschitzproperty_l2. The peridynamic force is also bounded for functions in . This again follows from the Lipschitz property of the force in using arguments established in LABEL:s:proofs. The boundedness of the energy in both and follows from the boundedness of the bond potential energy and used in the definition of , see LABEL:eq:tensilepotmod and LABEL:eq:hydropotmod. More generally this also shows that for .
We next discuss the spatial and the time discretization of peridynamic equation.
3 Finite element approximation
Let be given by linear continuous interpolations over tetrahedral or triangular elements where denotes the size of finite element mesh. Here we assume the elements are conforming and the finite element mesh is shape regular and .
For a continuous function on , is the continuous piecewise linear interpolant on . It is given by
(32) 
where is the local interpolant defined over finite element and is given by
(33) 
Here is the number of vertices in an element , is the position of vertex , and is the linear interpolant associated to vertex .
Application of Theorem 4.4.20 and remark 4.4.27 in [Brenner and Scott, 2007] gives
(34) 
Let denote the projection of on . For the norm it is defined as
(35) 
and satisfies
(36) 
Since , and LABEL:eq:interpolationerror we see that
(37) 
3.1 Semidiscrete approximation
Let be the approximation of satisfying following for all
(38) 
We have following result:
Theorem 5.
Energy stability of semidiscrete approximation
The semidiscrete scheme is stable and the energy , defined in LABEL:eq:totaldynamicenergy, satisfies the following bound
We note that while proving the stability of semidiscrete scheme corresponding to nonlinear peridynamics, we do not require any assumption on the strain . The proof is similar to [Section 6.2, [Lipton, 2016]].
Proof.
Letting in LABEL:eq:feweak and noting the identity LABEL:eq:energyidentity, we get
(39) 
We also have
where we used the fact that is nonnegative. We substitute above inequality in LABEL:eq:ineq_stab_1 to get
We fix and define as . Then, from the above equation we easily have
Noting that , integrating from to and relabeling as , we get
Proof is complete once we let and take the square of both sides.
4 Central difference time discretization
In LABEL:s:convergence 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 for the linearized peridynamic equation in LABEL:ss:fullfemstability.
At time step , the exact solution is given by where , and their projection onto is given by . The solution of fully discrete problem at time step is given by .
We approximate the initial data on displacement and velocity by their projections and . Let and . For , satisfies, for all ,
(40) 
where we have denoted the projection of , i.e. , as . Combining the two equations delivers central difference equation for . We have
(41) 
For , we have
(42) 
4.1 Implementation details
For completeness we describe the implementation of the time stepping method using FEM interpolants. Let be the shape tensor then are given by
(43) 
where and are dimensional vectors, where is the number of nodal points in the mesh and is the dimension.
From LABEL:eq:central, for all with elements of zero on the boundary, then the following holds for
(44) 
Here the mass matrix and force vector are given by
(45) 
where is defined by
(46) 
We remark that a similar equation holds for .
At the time step we must invert to solve for using
(47) 
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 matrixvector multiplication needs to be carried out at each time step. On the other hand the quadrature error in the computation of the force vector 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 . Here horizon is assumed to be fixed. We first compare the exact solution with its projection in 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 is given by
The error is split into two parts as follows