Non-local modeling for understanding the fracture process and the creation of interface

We formulate a nonlocal model for calculating the evolution of deformation inside a cracking body. We work within the small deformation setting and the strain $S(y,x)$ is expressed as a difference quotient of the displacement $u$ resolved along the direction $\mathbf{e}=\frac{y - x}{|y - x|}$, here \begin{equation} \label{strain} S=S(y,x) = \frac{u(y) - u(x)}{|y - x|} \cdot \mathbf{e}. \end{equation}

The acceleration of a material point is now due to the sum of the forces acting on the point from near by points. The constitutive relation relating the strain to force along the direction $\mathbf{e}$ for fixed $x$ and $y$ is initially elastic for small gradients and softens after a critical strain is reached. The nonlocal model treated here and given in (Lipton, 2014) and (Lipton, 2016) is a smooth version of the prototypical microelastic bond (PMB) model introduced by Silling in 2000 (Silling, 2000). The force versus strain relation is illustrated in Figure 1 below.

Stress vs Force
Fig.1 - Force as a function of strain.

The length scale of non-local interaction of a point $x$ with its environment is $\epsilon$ and $x$ interacts with all points other points $y$ inside a ball of radius $\epsilon$. The ball centered at $x$ with radius $\epsilon$ is denoted by $\mathcal{H}_\epsilon(x)$.

The force strain relation is derived from the elastic energy density at a material point $x$ and is assumed to be given by \begin{equation} W(x)=\frac{1}{V_\epsilon}\int_{\mathcal{H}_\epsilon(x)} |y-x|\mathcal{W}^\epsilon\big(\mathcal{S},x,y\big)dy \label{eqn-W} \end{equation} where $\mathcal{W}^\epsilon(\mathcal{S},x,y)$ is the force potential per unit length between $x$ and $y$ and $V_\epsilon$ is the area (in 2D) or the volume (in 3D) of the neighborhood $\mathcal{H}_\epsilon(x)$.

By Hamilton's principle applied to a bounded body $D\subset\mathbb{R}^d$, $d=2$, $3$, the equation of motion describing the displacement field $u(x,t)$ is \begin{equation} \rho\, u_{tt}(x,t)=\frac{2}{V_\epsilon}\int_{\mathcal{H}_\epsilon(x)\cap D}\, \big(\partial_{\mathcal{S}} \mathcal{W}^\epsilon(\mathcal{S}(y,x)\big)\mathbf{e}\,dy+b(x,t). \label{eqofmotion} \end{equation}

We assume the general form \begin{equation} \mathcal{W}^\epsilon(\mathcal{S},x,y)=\frac{J^\epsilon(|y-x|)}{\epsilon|y-x|} \Psi(|y-x|\mathcal{S}^2) \label{eqn-WJPsi} \end{equation} where $J^\epsilon(|y-x|)=J(|y-x|/\epsilon)>0$ is a weight function and $\Psi:[0,\infty)\rightarrow\mathbb{R}^+$ is a continuously differentiable function such that $\Psi(0)=0$, $\Psi'(0)>0$, and $ \Psi_\infty:=\lim_{r\rightarrow \infty}\Psi(r) < \infty$. The pairwise force density is then given by \begin{equation} \partial_\mathcal{S}\mathcal{W}^\epsilon(\mathcal{S},x,y)= \frac{2J^\epsilon(|y-x|)}{\epsilon}\Psi'\big(|y-x|\{\mathcal{S}^2\big)\mathcal{S}. \label{force} \end{equation} For fixed $x$ and $y$, there is a unique maximum in the curve of force versus strain, see Figure 1. The location of this maximum can depend on the distance between $x$ and $y$ and occurs at the strain $\mathcal{S}_c$ such that $\partial^2\mathcal{W}^\epsilon/\partial\mathcal{S}^2(\mathcal{S}_c,y-x)=0$. This value is $\mathcal{S}_c=\sqrt{\overline{r}/|y-x|}$, where $\overline{r}$ is the unique number such that $\Psi'(\overline{r})+2\overline{r}\Psi''(\overline{r})=0$. The evolution equation for the displacement inside the body is complemented by bounded initial conditions for the displacement $u_0\in L^2(D;\mathbb{R}^d)$ and velocity ${v}_0\in L^2(D;\mathbb{R}^d)$. One easily finds that the problem admits a unique evolution $ u(x,t)$ in $C^2([0,T];L^2(D;\mathbb{R}^d))$. Direct calculation shows that energy balance holds at each instant of the evolution, see (Lipton, 2014), (Lipton, 2016).

For finite horizon $\epsilon>0$ the elastic moduli and critical energy release rate are recovered directly from the strain potential $\mathcal{W}^\epsilon(\mathcal{S},x,y)$. First suppose the displacement inside $\mathcal{H}_\epsilon(x)$ is affine, that is, $u(x)=Fx$ where $F$ is a constant matrix. For small strains, i.e., $\mathcal{S}=Fe\cdot e\ll\mathcal{S}_c$, the strain potential is linear elastic to leading order and characterized by elastic moduli $\mu$ and $\lambda$ associated with a linear elastic isotropic material \begin{eqnarray} W(x)&=&\frac{1}{V_\epsilon}\int_{H_\epsilon(x)}|y-x|\mathcal{W}^\epsilon(\mathcal{S}(x,y)\,dy\nonumber\\ &=&2\mu |F|^2+\lambda |Tr\{F\}|^2+O(\epsilon|F|^4). \label{LEFMequality} \end{eqnarray} The elastic moduli $\lambda$ and $\mu$ are calculated directly from the strain energy density and are given by \begin{equation}\label{calibrate1} \mu=\lambda=M\frac{1}{d+2} \Psi'(0) \,, \end{equation} where the constant $M=\int_{0}^1r^{d}J(r)dr$ for dimensions $d=2,3$. In regions of discontinuity the same strain potential is used to calculate the amount of energy consumed by a crack per unit area of crack growth, i.e., the critical energy release rate $\mathcal{G}$. Calculation shows that $\mathcal{G}_c$ equals the work necessary to eliminate force interaction on either side of a fracture surface per unit fracture area and is given in three dimensions by \begin{equation} \mathcal{G}_c=\frac{4\pi}{V_d}\int_0^\epsilon\int_z^\epsilon\int_0^{\cos^{-1}(z/\zeta)} \mathcal{W}^\epsilon(\infty,\zeta)\zeta^2\sin{\phi}\,d\phi\,d\zeta\,dz \label{calibrate2formula} \end{equation} where $\zeta=|y-x|$ and $V_d$ is the volume of the $d$ dimensional unit ball, see Figure 2. In $d$ dimensions, the result is \begin{equation} \mathcal{G}_c=M\frac{2\omega_{d-1}}{\omega_d}\, \Psi_\infty\,, \label{calibrate2formulad} \end{equation} where $\omega_{d}$ is the volume of the $d$ dimensional unit ball, $\omega_1=2,\omega_2=\pi,\omega_3=4\pi/3$.

Spherical Cap
Fig.2 - Evaluation of the critical energy release rate $\mathcal{G}_c$. For each point $x$ along the dashed line, $0\leq z\leq \epsilon$, the work required to break the interaction between $x$ and $y$ in the spherical cap is summed up usinsing spherical coordinates centered at $x$, which depends on $z$.

It is crucial to note that $\mathcal{G}_c$ is independent of the length scale of nonlocal interaction. This is due to our choice of scaling by $\epsilon^{-1}$ in our formulation of the energy.

We prescribe a bounded initial displacement field $u_0(x)$ and bounded initial velocity field ${v}_0(x)$ with bounded Griffith fracture energy given by \begin{eqnarray} \int_{D}\,2\mu |\mathcal{E} u_0|^2+\lambda |{\rm div}\, u_0|^2\,dx+\mathcal{G}|\mathcal{J}_{u_0}|\leq C \label{LEFMboundinitil} \end{eqnarray} for some $C<\infty$. Here $\mathcal{J}_{u_0}$ is the initial crack set across which the displacement $u_0$ has a jump discontinuity. This jump set need not be geometrically simple; it can be a complex network of cracks. $|\mathcal{J}_{u_0}|=H^{d-1}(J_{u_0})$ is the $d-1$ dimensional Hausdorff measure of the jump set. This agrees with the total surface area (length) of the crack network for sufficently regular cracks for $d=3(2)$. The strain tensor associated with the initial displacement $u_0$ is denoted by $\mathcal{E} u_0$. Consider the sequence of solutions $u^{\epsilon}$ of the initial value problem associated with progressively smaller peridynamic horizons $\epsilon$. The peridynamic evolutions $u^{\epsilon}$ converge in mean square uniformly in time to a limit evolution $u^0(x,t)$ in $C([0,T];L_0^2(D,\mathbb{R}^d )$ and $u_t^0(x,t)$ in $L^2([0,T]\times D;\mathbb{R}^{d})$ with the same initial data, i.e., \begin{eqnarray} \lim_{\epsilon\rightarrow 0}\max_{0\leq t\leq T}\int_D|u^{\epsilon}(x,t)-u^0(x,t)|^2\,dx=0, \label{unifconvg} \end{eqnarray} see (Lipton, 2014) and (Lipton, 2016). It is found that the limit evolution $u^0(x,t)$ has bounded Griffith surface energy and elastic energy given by \begin{eqnarray} \int_{D}\,2\mu |\mathcal{E} u^0(t)|^2+\lambda |{\rm div}\,u^0(t)|^2\,dx+\mathcal{G}|\mathcal{J}_{u^0(t)}|\leq C, \label{LEFMbound} \end{eqnarray} for $0\leq t\leq T$, where $\mathcal{J}_{u^0(t)}$ denotes the evolving fracture surface inside the domain $D$, (Lipton, 2014), (Lipton, 2016). The limit evolution $u^0(t)$ is found to lie in the space of special functions of bounded deformation SBD, see (Lipton, 2016). For functions in SBD the strain $\mathcal{S}(x,y)$ is related to the strain tensor $\mathcal{E}u$ by \begin{equation} \lim_{\epsilon\rightarrow 0}\frac{1}{V_\epsilon}\int_{\mathcal{H}_\epsilon(x)}|\mathcal{S}-\mathcal{E}u(x) e\cdot e|\,dy, \label{equatesandE} \end{equation} for almost every $x$ in $D$. The jump set $\mathcal{J}_{u^0(t)}$ is the countable union of rectifiable surfaces (arcs) for $d=3(2)$, see (Ambrosio, Coscia, and Dal Maso, 1997).

In domains away from the crack set the limit evolution satisfies local linear elastodynamics (the PDEs of the standard theory of solid mechanics). Fix a tolerance $\tau>0$. If for subdomains $D'\subset D$ and for times $0\leq t\leq T $ the associated strains $\mathcal{S}_{u^\epsilon}$ satisfy $|\mathcal{S}_{u^\epsilon}|\leq S_c$ for every $\epsilon<\tau$ then it is found that the limit evolution $u^0(t,x)$ is governed by the PDE \begin{eqnarray} \rho u^0_{tt}(t,x)= {\rm div}\sigma(t,x)+b(t,x), \hbox{on $[0,T]\times D'$}, \label{waveequationn} \end{eqnarray} where the stress tensor $\sigma$ is given by \begin{eqnarray} \sigma =\lambda I_d Tr(\mathcal{E}\,u^0)+2\mu \mathcal{E}u^0, \label{stress} \end{eqnarray} $I_d$ is the identity on $\mathbb{R}^d$, and $Tr(\mathcal{E}\,u^0)$ is the trace of the strain. To summarize in the limit of small horizon $\epsilon\rightarrow 0$ peridynamic solutions converge in mean square to limit solutions that are linear elastodynamic away from the crack set, that is, the PDEs of the local theory hold at points off of the crack. The elastodynamic balance laws are characterized by elastic moduli $\mu$, $\lambda$. The evolving crack set possesses bounded Griffith surface free energy associated with the critical energy release rate $\mathcal{G}_c$. These convergence and localization results are established in the papers (Lipton, 2014), (Lipton, 2016).

The strategy for proving (\ref{unifconvg}) and (\ref{LEFMbound}) consists of five steps:

  1. Constructing upper bounds on the kinetic and potential energy of the evolutions that hold uniformly for $0\leq t\leq T$ and $0<\epsilon$.
  2. Showing compactness of the evolution $u^\epsilon(t)$ in $L^2(D)$ for each time $0\leq t\leq T$.
  3. Showing uniform Lipschitz continuity in time with respect to $\epsilon$ via the Gronwall inequality.
  4. Showing limit points of the sequence $\{u^\epsilon(t)\}$ belong to $SBD$ for each time $0\leq t\leq T$.
  5. Point wise convergence of the nonlocal potential energies to the Griffith fracture energy using the SBD structure theorem (Ambrosio, Coscia, and Dal Maso, 1997).
The strategy for proving \eqref{waveequationn} consists of composition of weak limit of the strains with a elastic-softening force characterized by a monotonic decreasing potential function $\Psi'(x)$. This allows passage to the the weak limit and identifying it as linear function of the absolutely continuous part of the symmetric strain as an element of SBD, see (Lipton, 2016). As it stands we have yet to derive kinetic relations for evolution of the crack tip in the limit of vanishing non locality. Moreover the interaction between elastic waves and crack growth should be articulated. Future work seeks to obtain kinetic relations for the crack front from passing to the limit of vanishing non-locality in the non-local model.

Current work focuses on establishing apriori convergence rates of finite difference and finite element approximations to nonlocal models. Initial work on apriori convergence rates have appered in (Jha and Lipton, 2018a) and (Jha and Lipton, 2018b)

Future work also will seek to validate the non-local fracture models by comparison of simulation with bench mark experiments.

The work discussed here has appeared in Dynamic Brittle Fracture as a Small Horizon Limit of Peridynamics. Journal of Elasticity, 117, Issue 1, 2014 pp 21-50, Cohesive Dynamics and Brittle Fracture. Journal of Elasticity, 124, Issue 2, 2016, pp. 143-191, Numerical analysis of nonlocal fracture models in Holder space. SIAM J. Numer. Anal. 56 (2018), pp. 906-941, and Numerical convergence of nonlinear nonlocal continuum models to local elastodynamics. International Journal of Numerical Methods in Engineering. 114 (2018), pp. 1389-1410.
This work is supported by ARO Grant W911NF-16-1-0456 and NSF Grant DMS-0807265


  1. L. Ambrosio, A. Coscia, and G. Dal Maso. Fine properties of functions with bounded deformation. Arch. Rational Mech. Anal., 139:201-238, 1997.
  2. P.K. Jha and R. Lipton. Numerical analysis of nonlocal fracture models in Holder space. SIAM J. Numer. Anal. 56 (2018a), pp. 906-941.
  3. P.K. Jha and R. Lipton. Numerical convergence of nonlinear nonlocal continuum models to local elastodynamics. International Journal of Numerical Methods in Engineering. 114 (2018b), pp. 1389-1410.
  4. R. Lipton. Cohesive Dynamics and Brittle Fracture. Journal of Elasticity, 124, Issue 2, 2016, pp. 143-191.
  5. R. Lipton. Dynamic Brittle Fracture as a Small Horizon Limit of Peridynamics. Journal of Elasticity, 117, Issue 1, 2014 pp 21-50.
  6. S.A. Silling. Reformulation of elasticity theory for discontinuities and long-range forces. Journal of the Mechanics and Physics of Solids, 48:175-209, 2000.