Numerical analysis of nonlocal fracture models in Hölder spaceSubmitted to the editors at SIAM Journal on Numerical Analysis on January 18, 2017. \fundingThis material is based upon work supported by the Army Research Office under award W911NF-16-1-0456.

# Numerical analysis of nonlocal fracture models in Hölder space††thanks: Submitted to the editors at SIAM Journal on Numerical Analysis on January 18, 2017. \fundingThis material is based upon work supported by the Army Research Office under award W911NF-16-1-0456.

Prashant K. Jha Department of Mathematics, Louisiana State University, Baton Rouge, LA ().    Robert Lipton Department of Mathematics, Louisiana State University, Baton Rouge, LA ().
###### Abstract

In this work, we calculate the convergence rate of the finite difference approximation for a class of nonlocal fracture models. We consider two point force interactions characterized by a double well potential. We show the existence of a evolving displacement field in Hölder space with Hölder exponent . The rate of convergence of the finite difference approximation depends on where gives the length scale of nonlocal interaction and is the discretization length. It is shown that this convergence rate holds for both the forward Euler scheme as well as general single step implicit schemes. The Hölder continuous evolutions converge to a brittle fracture evolution in the limit of vanishing nonlocality.

N
\headers

Numerical analysis of nonlocal fracture models in Hölder spacePrashant K. Jha and Robert Lipton

onlocal fracture models, peridynamics, cohesive dynamics, numerical analysis, finite difference approximation

{AMS}

34A34, 34B10, 74H55, 74S20

## 1 Introduction

Nonlocal formulations have been proposed to describe the evolution of deformations which exhibit loss of differentiability and continuity, see [26] and [29]. These models are commonly referred to as peridynamic models. The main idea is to define the strain in terms of displacement differences and allow nonlocal interactions between material points. This generalization of strain allows for the participation of a larger class of deformations in the dynamics. Numerical simulations based on peridynamic modeling exhibit formation and evolution of sharp interfaces associated with phase transformation and fracture [9], [30], [25], [16], [1], [12], [22], [6], [18], [28], [32], [17]. A recent summary of the state of the art can be found in [15].

In this work, we provide a numerical analysis for the class of nonlocal models introduced in [20] and [21]. These models are defined by a double well two point potential. Here one potential well is centered at zero and associated with elastic response while the other well is at infinity and associated with surface energy. The rational for studying these models is that they are shown to be well posed over the class of square integrable non-smooth displacements and, in the limit of vanishing non-locality, their dynamics recover features associated with sharp fracture propagation see, [20] and [21]. The numerical simulation of prototypical fracture problems using this model is carried out in [22]. More recently calibration of the model to experimentally determined elastic properties of material samples is undertaken in [10]. In order to develop an approximation theory, we show the nonlocal evolution is well posed over a more regular space of functions. To include displacement fields which have no well-defined derivatives, we consider displacement fields in Hölder space with Hölder exponent taking any value in . We first show that a unique evolution exists in for initial data and body force. With existence and uniqueness proved we then develop an approximation theory for the forward Euler scheme. It is then shown that these ideas can be easily extended to the backward Euler scheme as well other implicit one step time discretization schemes. We show that the discrete approximation converges to the exact solution in the norm uniformly over finite time intervals with the rate of convergence given by , where is the size of time step, is the size of spatial mesh discretization, and is the length scale of nonlocal interaction. We then show, using the methods developed in [20] and [21], that in the limit , the Hölder continuous evolutions converge to a limiting sharp fracture evolution with bounded Griffiths fracture energy. Here the limit evolution is differentiable off the crack set and satisfies the linear elastic wave equation.

In the language of nonlocal operators, the integral kernel associated with the nonlocal model studied here is Lipschitz continuous guarantying global stability of the finite difference approximation. This is in contrast to PDE based evolutions where stability can be conditional. In addition, we examine local stability criteria for choosing the actual size of time steps in an implementation. Our results show that care must be taken when implementing the finite difference approach in the presence of dynamic instabilities.

There is now a large body of contemporary work addressing the numerical approximation of singular kernels with application to nonlocal diffusion, advection, and mechanics. Numerical formulations and convergence theory for nonlocal -Laplacian formulations are developed in [14], [24]. Numerical analysis of nonlocal steady state diffusion is presented in [31] and [23], and [8]. The use of fractional Sobolev spaces for nonlocal problems is investigated and developed in [13]. Quadrature approximations and stability conditions for linear peridynamics are analyzed in [33] and [27]. The interplay between nonlocal interaction length and grid refinement for linear peridynamic models is presented in [7]. Analysis of adaptive refinement and domain decomposition for linearized peridynamics are provided in [3], [19], and [2]. This list is by no means complete and the literature on numerical methods and analysis continues to grow.

The paper is organized as follows. In section 2, we describe the nonlocal model. In subsection 2.2, we state theorems which show Lipschitz continuity of the nonlocal force (Proposition 1) and the existence and uniqueness of an evolution over any finite time interval (Theorem 2). In section 3, we compute the convergence rate of the forward Euler scheme as well as implicit one step methods. In subsection 3.5, we discuss local stability criteria for choosing the size of time steps. The results of this section show that caution must be used when implementing these methods. In section 4, we give the proof of Proposition 1, Theorem 5, and Theorem 2. The convergence of Hölder continuous evolutions to sharp fracture evolutions as is shown in section 5. We provide conclusions in section 6.

## 2 Double well potential and existence of a solution

In this section, we present the nonlinear nonlocal model. Let , be the material domain with characteristic length-scale of unity. Let be the size of horizon across which nonlocal interaction between points takes place. The material point interacts nonlocally with all material points inside a horizon of length . Let be the ball of radius centered at containing all points that interact with . After deformation the material point assumes position . In this treatment we assume infinitesimal displacements and the strain is written in terms of the displacement as

 S=S(y,x;u) :=u(y)−u(x)|y−x|⋅y−x|y−x|. (1)

Let be the nonlocal potential between material point and . The energy density at is given by

 Wϵ(S,x)=1ϵdωd∫Hϵ(x)Wϵ(S,y−x)dy, (2)

where is the volume of a unit ball in -dimension and is the volume of the ball of radius . The potential energy is written as

 PDϵ(u) =∫DWϵ(S(u),x)dx, (3)

and the displacement field satisfies following equation of motion

 ρ∂2ttu(t,x) =−∇PDϵ(u)+b(t,x) (4)

for all . Here we have

 −∇PDϵ(u)(x)=2ϵdωd∫Hϵ(x)∂SWϵ(S,y−x)y−x|y−x|dy (5)

where is the body force, is the density and is the derivative of potential with respect to the strain.

Let be the layer of thickness surrounding . We prescribe the nonlocal boundary condition on as follows

 u(x)=0∀x∈Dc. (6)

The peridynamic equation, boundary values, together with the initial conditions

 u(0,x)=u0(x)∂tu(0,x)=v0(x) (7)

determine the peridynamic evolution .

### 2.1 Nonlocal potential

We consider the nonlocal two point interaction potential of the form

 Wϵ(S,y−x) =Jϵ(|y−x|)ϵf(|y−x|S2) (8)

where is assumed to be positive, smooth and concave with following properties

 limr→0+f(r)r=f′(0),limr→∞f(r)=f∞<∞ (9)

The potential is of double well type and convex near the origin where it has one well and concave and bounded at infinity where it has the second well. models the influence of separation between points and . We define by rescaling , i.e. . Here is zero outside the ball and satisfies for all .

The potential described in Equation 8 gives the convex-concave dependence, see Figure 1, of on the strain for fixed . Initially the deformation is elastic for small strains and then softens as the strain becomes larger. The critical strain where the force between and begins to soften is given by and the force decreases monotonically for

 (10)

Here is the inflection point of and is the root of following equation

 f′(r2)+2r2f′′(r2)=0. (11)

### 2.2 Existence of solution

Let be the Hölder space with exponent . The norm of is given by

 ∥u∥C0,γ(D;Rd) :=supx∈D|u(x)|+[u]C0,γ(D;Rd), (12)

where is the Hölder semi norm and given by

 [u]C0,γ(D;Rd) :=supx≠y,x,y∈D|u(x)−u(y)||x−y|γ. (13)

Let be the Banach space of functions in such that they satisfy the nonlocal boundary condition described in Equation 6. We write the evolution Equation 4 as an equivalent first order system with and with . Let where and let such that

 F1(y,t) :=y2 (14) F2(y,t) :=−∇PDϵ(y1)+b(t). (15)

The initial boundary value associated with the evolution Equation 4 is equivalent to the initial boundary value problem for the first order system given by

 ddty=Fϵ(y,t), (16)

with initial condition given by .

The function satisfies the Lipschitz continuity given by the following theorem.

###### Proposition 1

Lipschitz continuity and bound
Let . The function , as defined in Equation 14 and Equation 15, is Lipschitz continuous in any bounded subset of . We have, for any ,

 ∥Fϵ(y,t)−Fϵ(z,t)∥X ≤(L1+L2(∥y∥X+∥z∥X))ϵ2+α(γ)∥y−z∥X (17)

where are independent of and depend on peridynamic potential function and influence function and the exponent is given by

 α(γ)={0if γ≥1/21/2−γif γ<1/2. (18)

Furthermore for any and any , we have the bound

 ∥Fϵ(y,t)∥X ≤L3ϵ2(1+∥y∥X)+b (19)

where and is independent of .

In Theorem 6.1 of [21], the Lipschitz property of a peridynamic force is shown in . It is given by

 ∥Fϵ(y,t)−Fϵ(z,t)∥X ≤Lϵ2∥y−y∥X∀y,z∈X,∀t∈[0,T] (20)

for all . For this case does not depend on . We now state the existence theorem.

The following theorem gives the existence and uniqueness of solution in any given time domain .

###### Theorem 2

Existence and uniqueness of Hölder solutions of cohesive dynamics over finite time intervals
For any initial condition , interval , and right hand side continuous in time for and , there is a unique solution of

 y(t)=x0+∫t0Fϵ(y(τ),τ)dτ, (21)

or equivalently

 y′(t)=Fϵ(y(t),t),with y(0)=x0, (22)

where and are Lipschitz continuous in time for .

The proof of this theorem is given in section 4. We now describe the finite difference scheme and analyze its convergence to Hölder continuous solutions of cohesive dynamics.

## 3 Finite difference approximation

In this section, we present the finite difference scheme and compute the rate of convergence. Let be the size of a mesh and be the size of time step. We will keep fixed and assume that . Let be the discretization of material domain and be the discretization of time domain. Let be the index such that . Recall that the exact solution is written , and we will denote the finite difference approximation at the grid point, , as . The exact solution evaluated at grid points is denoted by . We enforce boundary conditions by assuming that discrete set is extended to all and satisfies if .

We begin with the forward Euler time discretization and the finite difference scheme for is written

 ^uk+1i−^ukiΔt =^vki (23) ^vk+1i−^vkiΔt =−∇PDϵ(^uk)(xi)+bki (24)

The scheme is complemented with the discretized initial conditions and . Here we have assumed, without loss of generality, .

Let be the unit cell of a volume corresponding to grid point . The piecewise constant extensions of the discrete sets and are given by

 ^uk(x) :=∑i,xi∈D^ukiχUi(x) (25) ^vk(x) :=∑i,xi∈D^vkiχUi(x) (26)

In this way we represent the finite difference solution as a piecewise constant function. We will show this function provides an approximation of the exact solution.

### 3.1 Convergence results

In this section we provide upper bounds on the rate of convergence of the discreet approximation to the solution of the peridynamic evolution. The approximation error at time , for is defined as

 Ek (27)

The upper bound on the convergence rate of the approximation error is given by the following theorem.

###### Theorem 3

Convergence of finite difference approximation (forward Euler time discretization)
Let be fixed. Let be the solution of peridynamic equation Equation 16. We assume . Then the finite difference scheme given by Equation 23 and Equation 24 is consistent in both time and spatial discretization and converges to the exact solution uniformly in time with respect to the norm. We assume the error at the initial step is zero then the error at time is bounded and satisfies

 sup0≤k≤T/ΔtEk=O(Δt+hγϵ2). (28)

Here we have assumed the initial error to be zero for ease of exposition only.

An identical convergence rate can be established for the general one step scheme and we state it below.

###### Theorem 4

Convergence of finite difference approximation (General single step time discretization)
Let us assume that the hypothesis of Theorem 3 holds. Fix , and let be the solution of following finite difference equation

 ^uk+1i−^ukiΔt =(1−θ)^vki+θ^vk+1i (29) ^vk+1i−^vkiΔt =(1−θ)(−∇PDϵ(^uk)(xi)+bki)+θ(−∇PDϵ(^uk+1)(xi)+bk+1i). (30)

Then, for any fixed , there exists a constant independent of and , such that for the finite difference scheme given by Equation 29 and Equation 30 is consistent in time and spatial discretization. We assume the error at the initial step is zero then the error at time is bounded and satisfies

 sup0≤k≤T/ΔtEk=O(Δt+hγϵ2).

The constant is given by the explicit formula where is described by equation Equation 72.

Furthermore for the Crank Nicholson scheme, , if we assume the solutions belong to , then the approximation error satisfies

 sup0≤k≤T/ΔtEk=O((Δt)2+hγϵ2). (31)

As before we assume that the error in the initial data is zero for ease of exposition. The proofs of Theorem 3 and Theorem 4 are given in the following sections.

### 3.2 Error analysis

Theorem 3 and Theorem 4 are proved along similar lines. In both cases we define the -projections of the actual solutions onto the space of piecewise constant functions defined over the cells . These are given as follows. Let be the average of the exact solution in the unit cell given by

 ~uki :=1hd∫Uiuk(x)dx (32) ~vki :=1hd∫Uivk(x)dx (33)

and the projection of the solution onto piecewise constant functions are given by

 ~uk(x) :=∑i,xi∈D~ukiχUi(x) (34) ~vk(x) :=∑i,xi∈D~vkiχUi(x) (35)

The error between with is now split into two parts. From the triangle inequality, we have

 ≤∥∥^uk−~uk∥∥L2(D;Rd)+∥∥~uk−uk∥∥L2(D;Rd) (36) ≤∥∥^vk−~vk∥∥L2(D;Rd)+∥∥~vk−vk∥∥L2(D;Rd) (37)

In subsection 3.3 and subsection 3.4 we will show that the error between the projections of the actual solution and the discrete approximation for both forward Euler and implicit one step methods decay according to

 =O(Δt+hγϵ2). (38)

In what follows we can estimate the terms

 (39)

and show they go to zero at a rate of uniformly in time. The estimates given by Equation 38 together with the estimates for Equation 39 establish Theorem 3 and Theorem 4. We now establish the estimates for the differences and .

We write

 =∑i,xi∈D∫Ui∣∣∣1hd∫Ui(uk(y)−uk(x))dy∣∣∣2dx =∑i,xi∈D∫Ui[1h2d∫Ui∫Ui(uk(y)−uk(x))⋅(uk(z)−uk(x))dydz]dx ≤∑i,xi∈D∫Ui[1hd∫Ui∣∣uk(y)−uk(x)∣∣2dy]dx (40)

where we used Cauchy’s inequality and Jensen’s inequality. For , , where for and for . Since we have

 ∣∣uk(x)−uk(y)∣∣ (41)

and substitution in Equation 40 gives

 ≤c2γh2γ∑i,xi∈D∫Uidx(supt∥u(t)∥C0,γ(D;Rd))2 ≤c2γ(|D|+hd−1|∂D|)h2γ(supt∥u(t)∥C0,γ(D;Rd))2. (42)

A similar estimate can be derived for and substitution of the estimates into Equation 39 gives

 supk(∥∥~uk−u(tk)∥∥L2(D;Rd)+∥∥~vk−v(tk)∥∥L2(D;Rd))=O(hγ). (43)

In the next section we establish the error estimate Equation 38 for both forward Euler and general one step schemes in subsection 3.3 and subsection 3.4.

### 3.3 Error analysis for approximation of L2 projection of the the exact solution

In this sub-section, we estimate the difference between approximate solution and the projection of the exact solution onto piece wise constant functions given by , see Equation 34 and Equation 35 . Let and . Let the time and spatial discretization error of the peridynamic equation be denoted by and respectively.

Substracting from Equation 23 to get

 ^uk+1i−^ukiΔt−~uk+1i−~ukiΔt =^vki−~uk+1i−~ukiΔt =^vki−~vki+(~vki−∂~uki∂t)+(∂~uki∂t−~uk+1i−~ukiΔt). (44)

Taking the average over unit cell of the exact peridynamic equation Equation 16 at time , we will get . Therefore, the equation for is given by

 ek+1i(u)=eki(u)+Δteki(v)+Δtτki(u), (45)

where we identify the discretization error as

 τki(u) :=∂~uki∂t−~uk+1i−~ukiΔt. (46)

Similarly, we substract from Equation 24 and add and subtract terms to get

 ^vk+1i−^vkiΔt−~vk+1i−~vkiΔt =−∇PDϵ(^uk)(xi)+bki−∂vki∂t+(∂vki∂t−~vk+1i−~vkiΔt) =−∇PDϵ(^uk)(xi)+bki−∂vki∂t +(∂~vki∂t−~vk+1i−~vkiΔt)+(∂vki∂t−∂~vki∂t), (47)

where we identify as follows

 τki(v) :=∂~vki∂t−~vk+1i−~vkiΔt. (48)

Note from the exact peridynamic equation, we have

 bki−∂vki∂t=∇PDϵ(uk)(xi). (49)

Combining subsection 3.3, Equation 48, and Equation 49, to get

 ek+1i(v) =eki(v)+Δtτki(v)+Δt(∂vki∂t−∂~vki∂t) +Δt(−∇PDϵ(^uk)(xi)+∇PDϵ(uk)(xi)) =eki(v)+Δtτki(v)+Δt(∂vki∂t−∂~vki∂t) +Δt(−∇PDϵ(^uk)(xi)+∇PDϵ(~uk)(xi)) +Δt(−∇PDϵ(~uk)(xi)+∇PDϵ(uk)(xi)). (50)

The spatial discretization error and is given by

 σki(u) :=(−∇PDϵ(~uk)(xi)+∇PDϵ(uk)(xi)) (51) σki(v) :=∂vki∂t−∂~vki∂t. (52)

We finally have

 ek+1i(v) =eki(v)+Δt(τki(v)+σki(u)+σki(v)) +Δt(−∇PDϵ(^uk)(xi)+∇PDϵ(~uk)(xi)). (53)

We now show the consistency and stability property of the numerical scheme.

#### 3.3.1 Consistency

We deal with time discretization error and spatial discretization error separately. Time discretization error follows easily using the Taylor’s series while spatial discretization error uses the property of nonlinear peridynamic force.

Time discretization: We first estimate the time discretization error. We do Taylor series expansion to estimate as follows

 τki(u) =1hd∫Ui(∂uk(x)∂t−uk+1(x)−uk(x)Δt)dx =1hd∫Ui(−12∂2uk(x)∂t2Δt+O((Δt)2))dx.

Computing the norm of and using Jensen’s inequality to get

 ∥∥τk(u)∥∥L2(D;Rd) ≤Δt∥∥∥∂2uk∂t2∥∥∥L2(D;Rd)+O((Δt)2) ≤Δtsupt∥∥∥∂2u(t)∂t2∥∥∥L2(D;Rd)+O((Δt)2). (54)

Therefore, we have

 supk∥∥τk(u)∥∥L2(D;Rd)=O(Δt). (55)

Similarly, we have

 supk∥∥τk(v)∥∥L2(D;Rd)=O(Δt). (56)

Equation 55 and Equation 55 show that the numerical scheme is consistent in the time discretization.

Spatial discretization: We now estimate the spatial discretization error. Substituting the definition of and following the similar steps employed in subsection 3.2, to get

 ∣∣σki(v)∣∣ =∣∣ ∣∣∂vki∂t−1hd∫Ui∂vk(x)∂tdx∣∣ ∣∣≤cγhγ∫Ui1|xi−x|γ∣∣∣∂vk(xi)∂t−∂vk(x)∂t∣∣∣dx ≤cγhγ∥∥∥∂vk∂t∥∥∥C0,γ(D;Rd)≤cγhγsupt∥∥∥∂v(t)∂t∥∥∥C0,γ(D;Rd). (57)

Taking the norm of error and substituting the estimate above, to get

 ∥∥σk(v)∥∥2L2(D;Rd) ≤h2γc2γ(|D|+hd−1|∂D|)(supt∥∥∥∂v(t)∂t∥∥∥C0,γ(D;Rd))2. (58)

From above we conclude that

 supk∥∥σk(v)∥∥L2(D;Rd)=O(hγ). (59)

To estimate , we make use of Equation 125. We have

 ∣∣σki(u)∣∣