A new fracture criterion for peridynamic and dual-horizon peridynamics

Abstract

A new fracture criterion based on the crack opening displacement for peridynamic (PD) and dual-horizon peridynamics (DH-PD) is proposed. When the relative deformation of the PD bond between the particles reaches the critical crack tip opening displacement of the fracture mechanics, we assume that the bond force vanishes. A new damage rule similar to the local damage rule in conventional PD is introduced to simulate fracture. The new formulation is developed for a linear elastic solid though the extension to nonlinear materials is straightforward. The performance of the new fracture criterion is demonstrated by four examples, i.e. a bilateral crack problem, double parallel crack, monoclinic crack and the double inclined crack. The results are compared to experimental data and the results obtained by other computational methods.

Keywords

Castigliano’s theorem / breaking energy / critical extension / XFEM / COD / PD-COD

Introduction

The prediction of fracture is a long-standing problem in the field of computational solid mechanics. The inherent difficulty arises from the basic incompatibility of cracks with the partial differential equations that are used in the classical theory of solid mechanics. Many fracture and multiscale modiling of fracture and numerical methods have been proposed such as Fracture modeling [¹–⁴] and Multiscale modeling [⁵–¹¹] including fourth order phase-field model, a novel two-stage discrete crack method, a higher-order stress-based gradient-enhanced damage model, an adaptive multiscale method, efficient coarse graining in multiscale modeling and so on, Fracture Finite Element methods (FEM) including efficient remeshing techniques [¹²–¹⁷], extended finite element method (XFEM) [¹⁸,¹⁹] or XIGA [²⁰,²¹], the numerical manifold method (NMM) [²²,²³], element-free Galerkin (EFG) methods [²⁴],the reproducing kernel particle method (RKPM) [²⁵], and many other meshless methods [²⁶–²⁸] and enriched meshfree method (MM), see e.g. the contributions by Rabczuk [²⁹–³⁷]. The cracking particle methods (CPM) [³⁸–⁴¹] is a method that can easily deal with complex crack patterns as fracture is a natural outcome of the simulation. For all other methods, mixed mode fracture is usually studied either theoretical based on different failure criteria [⁴²–⁴⁵] or using test methods. Researchers prefer to conduct their experiments, and the specimens must be designed since they will be able to provide the same states, e.g. the centrally cracked Brazilian disk specimen [⁴⁶–⁴⁹], the compact tension shear specimen [⁵⁰–⁵²], but the experimental fracture studies on real components are very expensive and difficult. Peridynamics is a method that reformulates the fundamental equations of continuum mechanics, so that fracture is a natural outcome of the simulation [⁵³,⁵⁴]. Therefore, the object is discretized with particles and the material is modelled by the interaction of those particles. The first version is the so-called bond-based peridynamics formulation which is limited to a very specific material behaviour. The state-based peridynamics formulation overcomes this difficulty but both formulations (bond-based as well as state-based) require a uniform particle spacing. Dual-horizon peridynamics is a formulation which allows for adaptive refinement with non-uniform spaced particles [⁵⁵,⁵⁶].

Silling et al. [⁵⁷], investigated the deformation of an infinite bar; Weckner et al. [⁵⁸], used the Laplace and Fourier transforms in three-dimensional PD; Yu et al. [⁵⁹], proposed an adaptive trapezoidal integration scheme. Kilic [⁶⁰] described an efficient load distribution scheme. Peridynamics [⁶¹,⁶²] has attracted great attention due to its flexibility in modeling complex fracture patterns. Silling and Bobaru proposed a weighted local function of the particle weight method to determine the particle damage problem [⁶³]. Silling and Askari [⁶⁴] derived the critical energy release rate for bond-based in integral form. Foster et al. [²⁶] proposed the critical energy density as a failure criterion in ate-dependent situations. Silling and Bobaru [⁶⁶] used the function \mu to modify the force density vector. Ayatollahi and Aliha [⁶⁷] demonstrated the failure parameter of a critical stretch by experiments. Feng and Zhang simulated the cracking process of concrete [⁶⁸]. The article of Zhou [⁶⁹] examined rock-like materials. While Ren et al. [⁷⁰] proposed a new criteria for damage determination of shear deformation.

In this paper, we present a new fracture criterion in PD. Though the formulation is devised for linear elastic solids, it can easily be extended to non-linear materials. This new idea is comparable to the COD criterion employed in fracture mechanics theory. When the relative deformation of two adjacent particles in PD reaches the critical COD value, the PD forces between the two particles vanish and a crack is formed. The contents of this article are summarized below: Section 2 introduces the PD theory. In section 3, the COD criterion is derived in the context of PD. Also, the critical value of PD-COD which governs the crack propagation is provided. In section 4, four examples are presented which verify the new formulation. The results of the first example is compared to experimental data and results obtained by other methods included the 'classical' PD method and XFEM. The subsequent three examples include the center double parallel crack problem, the center monoclinic crack problem and the center double oblique crack problem. Finally, the manuscript closes with conclusions in section 5.

The peridynamic ( PD ) theory

The PD theory discretizes the objects into many particles. The force between the particles changes with the distance between the particles. As seen in Fig. 1: {\mathit{u}}_{(k)},{\mathit{u}}_{(j)} are the displacements of particle k and j, respectively. The position vector of those particles in the initial configuration is denoted by {\mathit{x}}_{(k)} and {\mathit{x}}_{(j)}, respectively while {\mathit{y}}_{(k)},{\mathit{y}}_{(j)} indicate the position of particles k and j in the deformed configuration. The initial relative position vector ({\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}) prior deformation becomes ({\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}) after deformation. The relative position vector ({\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}) and the stretch between material points{\mathit{x}}_{(k)}and{\mathit{x}}_{(j)}can be defined as:

(1)({\mathit{y}}_{(j)}-{\mathit{y}}_{(k)})=\underline{Y}(x_{(k)},t)\langle {\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\rangle

and

(2)s_{(k)(j)}={\displaystyle \frac{(\left|{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}\right|-\left|{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\right|)}{\left|{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\right|}}

The force density vectort_{(k)(j)} that the material point at location{\mathit{x}}_{(j)} exerts on the material point at location {\mathit{x}}_{(k)} can be expressed as:

(3)t_{(k)(j)}({\mathit{u}}_{(j)}-{\mathit{u}}_{(k)},{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)},t)=\underset{\_}{T}({\mathit{x}}_{(k)},t)\langle {\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\rangle

The forcet_{(k)(j)} that a material point {\mathit{x}}_{(j)} exerts on the material point {\mathit{x}}_{(k)} in turn is given by

(4)t_{(k)(j)}=\mathrm{2}\delta \left\{ad{\displaystyle \frac{{\Lambda }_{(k)(j)}}{\left|{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\right|}}{\theta }_{(k)}+b({\mathit{x}}_{(j)}-{\mathit{x}}_{(k)},t)s_{(k)(j)}\right\}{\displaystyle \frac{{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}}{\left|{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}\right|}}

where {\Lambda }_{(k)(j)} is defined as:

(5){\Lambda }_{(k)(j)}=\left({\displaystyle \frac{{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}}{\left|{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}\right|}}\right)\cdot \left({\displaystyle \frac{{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}}{\left|{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\right|}}\right)

The micropotential w_{(k)(j)} between material points {\mathit{x}}_{(k)}and {\mathit{x}}_{(j)}and the strain energy density W_{(k)} of material point {\mathit{x}}_{(k)}can be expressed as:

(6)w_{(k)(j)}=w_{(k)(j)}({\mathit{y}}_{({\mathrm{1}}^k)}-{\mathit{y}}_{(k)},{\mathit{y}}_{({\mathrm{2}}^k)}-{\mathit{y}}_{(k)},\cdots )

and

(7)\begin{array}{l}{W}_{(k)}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \sum _{j=\mathrm{1}}^{\infty }({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{w}_{(k)(j)}({\mathit{y}}_{({\mathrm{1}}^k)}-{\mathit{y}}_{(k)},{\mathit{y}}_{({\mathrm{2}}^k)}-{\mathit{y}}_{(k)},\cdots )}\\ +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{w}_{(j)(k)}({\mathit{y}}_{({\mathrm{1}}^j)}-{\mathit{y}}_{(j)},{\mathit{y}}_{({\mathrm{2}}^j)}-{\mathit{y}}_{(j)},\cdots )\cdot V_{(j)}\end{array}

It was shown in [⁷¹] that the dynamic equation of PD can be recast as:

(8)\rho (x)\ddot{\mathit{u}}(x,t)={\displaystyle {\int }_{H_x}\mathit{f}({\mathit{u}}_{(j)}-{\mathit{u}}_{(k)},{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)})d{H}_x+b(\mathit{x},t)}

in which H_xis the neighborhood of \mathit{x} and b(\mathit{x},t) is the body force density field.

PD-COD crack fracture criterion

Consider a crack and a crack extension as illustrated in Fig. 2 ( SZW: the length of extension zone, which is the amount of crack forward extension). The opening displacement of the crack tip is twice the height of the extension zone ( COD= 2×SZD ) and the stretch height is equal to the length of extension of the crack front ( SZD= SZW ). The actual crack tip opening displacement measurement is shown in Fig. 3 exploiting the crack tip symmetrical primitive crack at right angle to intersect up and down the crack surface points a and a^{\prime }. This distance of two points is the value of the opening displacement COD at the crack tip.

As the load increases, the opening displacement of the crack tip \delta increases, when a critical value{\delta }_{cr}is reached. Therefore, the crack propagation criterion established by COD is written as:

(9)\delta \ge {\delta }_{cr}

where {\delta }_{cr} is the critical COD value, determining the beginning of the crack extension.

There are two main methods to calculated the COD value, one is the D-M model derived from the BCS formula [⁷²,⁷³], the other is the Wells formula [⁷⁴,⁷⁵]. In this paper, the BCS formula is used to solve the crack propagation problem.

Previous studies on tensile tests of a large sheets have revealed a flat plastic zone as depicted in Fig. 4 (a). The plastic zone is simplified into a triangular curve. Let us assume the material of the plastic zone is ideal plasticity, and the material surrounded by the plastic zones and cracks is the elastic zone. The plastic zone of the crack tip is excavated, formatted a center through crack which length is \mathrm{2}a=\mathrm{2}c+\mathrm{2}R in the elastic infinite plate, as shown in Fig. 4 (b) above. So, the above problem is a simplified D-M elasticity model where the displacement of the crack tip \delta can be obtained by the Paris displacement formula.

As shown in Fig. 5 by the Castigliano`s theorem, the relative displacement between the two particles is obtained by the derivative of the elastic strain energy as:

(10){\delta }_i={\displaystyle \frac{\partial U}{\partial P_i}}

where U is the strain energy of the elastic body, P_i is the force and, {\delta }_i is the displacement.

When F\to \mathrm{0}, the actual displacement between any two points is recast as:

(11){\delta }_i=\underset{F\to \mathrm{0}}{\lim }{\displaystyle \frac{\partial U}{\partial F}}

The relationship between the crack propagation force G_c and U can be obtained by

(12)G_c={\left({\displaystyle \frac{\partial U}{\partial A}}\right)}_P

Integrated Eq.(12) over A yields

(13)U=U_{\mathrm{0}}+{\displaystyle {\int }_{\mathrm{0}}^A{G}_c}\mathrm{d}A

(14)G_c=G_{\mathrm{{\rm I}}}+G_{\mathrm{{\rm I}}\mathrm{{\rm I}}}={\displaystyle \frac{K_c^{\mathrm{2}}}{E}}={\displaystyle \frac{K_{\mathrm{{\rm I}}}^{\mathrm{2}}}{E}}+{\displaystyle \frac{K_{\mathrm{{\rm I}}\mathrm{{\rm I}}}^{\mathrm{2}}}{E}}

With

(15)K_{\mathrm{{\rm I}}}=K_{\mathrm{{\rm I}}P}+K_{\mathrm{{\rm I}}F}\mathrm{and}{K}_{\mathrm{{\rm I}}\mathrm{{\rm I}}}=K_{\mathrm{{\rm I}}\mathrm{{\rm I}}P}+K_{\mathrm{{\rm I}}\mathrm{{\rm I}}F}

where U_{\mathrm{0}} is the strain energy of the elastic body with crack length \mathrm{2}a=\mathrm{0} and U is the strain energy of the elastomer with crack length \mathrm{2}a\ne \mathrm{0}.

Substituting Eqs.(13), (14) and (15) into Eq.(11) leads to

(16)\begin{array}{l}\delta =\underset{F\to \mathrm{0}}{\lim }\left[{\displaystyle \frac{\partial U_{\mathrm{0}}}{\partial F}}+{\displaystyle \frac{\partial }{\partial F}}{\displaystyle {\int }_{\mathrm{0}}^A{G}_c\mathrm{d}A}\right]\\ =\underset{F\to \mathrm{0}}{\lim }\left[{\displaystyle \frac{\partial U_{\mathrm{0}}}{\partial F}}+{\displaystyle \frac{\partial }{\partial F}}{\displaystyle {\int }_{\mathrm{0}}^A\left({\displaystyle \frac{K_{\mathrm{{\rm I}}}^{\mathrm{2}}}{E}}+{\displaystyle \frac{K_{\mathrm{{\rm I}}\mathrm{{\rm I}}}^{\mathrm{2}}}{E}}\right)\mathrm{d}A}\right]\\ =\underset{F\to \mathrm{0}}{\lim }\left[{\displaystyle \frac{\partial U_{\mathrm{0}}}{\partial F}}+{\displaystyle \frac{\mathrm{1}}{E}}\cdot {\displaystyle \frac{\partial }{\partial F}}{\displaystyle {\int }_{\mathrm{0}}^A(K_{\mathrm{{\rm I}}}^{\mathrm{2}}+K_{\mathrm{{\rm I}}\mathrm{{\rm I}}}^{\mathrm{2}})\mathrm{d}A}\right]\\ =\underset{F\to \mathrm{0}}{\lim }\left[{\displaystyle \frac{\partial U_{\mathrm{0}}}{\partial F}}+{\displaystyle \frac{\mathrm{1}}{E}}\cdot {\displaystyle \frac{\partial }{\partial F}}\left({\displaystyle {\int }_{\mathrm{0}}^A{(K_{\mathrm{{\rm I}}P}+K_{\mathrm{{\rm I}}F})}^{\mathrm{2}}\mathrm{d}A}+{\displaystyle {\int }_{\mathrm{0}}^A{(K_{\mathrm{{\rm I}}\mathrm{{\rm I}}P}+K_{\mathrm{{\rm I}}\mathrm{{\rm I}}F})}^{\mathrm{2}}\mathrm{d}A}\right)\right]\\ =\underset{F\to \mathrm{0}}{\lim }\left\{{\displaystyle \frac{\partial U_{\mathrm{0}}}{\partial F}}+{\displaystyle \frac{\mathrm{2}}{E}}\cdot \left({\displaystyle {\int }_{\mathrm{0}}^A(K_{\mathrm{{\rm I}}P}+K_{\mathrm{{\rm I}}F}){\displaystyle \frac{\partial K_{\mathrm{{\rm I}}F}}{\partial F}}\mathrm{d}A}+{\displaystyle {\int }_{\mathrm{0}}^A(K_{\mathrm{{\rm I}}\mathrm{{\rm I}}P}+K_{\mathrm{{\rm I}}\mathrm{{\rm I}}F}){\displaystyle \frac{\partial K_{\mathrm{{\rm I}}\mathrm{{\rm I}}F}}{\partial F}}\mathrm{d}A}\right)\right\}\end{array}

Since K_{\mathrm{{\rm I}}F} is proportional to F, when F\to \mathrm{0}, Eq.(16) can be rewritten as:

(17)\delta ={\left({\displaystyle \frac{\partial U_{\mathrm{0}}}{\partial F}}\right)}_{F=\mathrm{0}}+{\displaystyle \frac{\mathrm{2}}{E}}{\displaystyle {\int }_{\mathrm{0}}^A{K}_{\mathrm{{\rm I}}P}{\displaystyle \frac{\partial K_{{\rm I}F}}{\partial F}}\mathrm{d}A}+{\displaystyle \frac{\mathrm{2}}{E}}{\displaystyle {\int }_{\mathrm{0}}^A{K}_{\mathrm{{\rm I}}\mathrm{{\rm I}}P}{\displaystyle \frac{\partial K_{\mathrm{{\rm I}}\mathrm{{\rm I}}F}}{\partial F}}\mathrm{d}A}

where K_{{\rm I}P},K_{{\rm I}{\rm I}P} and K_{{\rm I}F},K_{{\rm I}{\rm I}F} are the stress intensity factor under force P and the virtual equilibrium force F, respectively.

From the relationship between the critical energy release rate G_c and the strain energy U, the crack opening displacement {\delta }_i is obtained as:

(18)\delta ={\displaystyle \frac{\mathrm{2}}{E}}{\displaystyle {\int }_{\mathrm{0}}^A{K}_{\mathrm{{\rm I}}P}{\displaystyle \frac{\partial K_{\mathrm{{\rm I}}F}}{\partial F}}}\mathrm{d}A+{\displaystyle \frac{\mathrm{2}}{E}}{\displaystyle {\int }_{\mathrm{0}}^A{K}_{\mathrm{{\rm I}}\mathrm{{\rm I}}P}{\displaystyle \frac{\partial K_{\mathrm{{\rm I}}\mathrm{{\rm I}}F}}{\partial F}}}\mathrm{d}A

where K_{{\rm I}P}=\sigma \sqrt{\pi \xi },K_{{\rm I}F}=F/\sqrt{\pi \xi } and K_{{\rm I}{\rm I}P}=\tau \sqrt{\pi \xi },K_{{\rm I}{\rm I}F}=F/\sqrt{\pi \xi }are the stress intensity factor of the force and the virtual equilibrium force at the crack tip, \xi is the instantaneous crack length.

Then when the crack length A=\mathrm{2}a, the opening displacement of the crack tip can be expressed as:

(19)\begin{array}{l}\delta ={\displaystyle \frac{\mathrm{2}}{E}}{\displaystyle {\int }_{\mathrm{0}}^A\sigma \sqrt{\pi \xi }{\displaystyle \frac{\mathrm{1}}{\sqrt{\pi \xi }}}}\mathrm{d}A+{\displaystyle \frac{\mathrm{2}}{E}}{\displaystyle {\int }_{\mathrm{0}}^A\tau \sqrt{\pi \xi }{\displaystyle \frac{\mathrm{1}}{\sqrt{\pi \xi }}}}\mathrm{d}A\\ \mathrm{=2}{\displaystyle \frac{\sigma }{E}}A+2{\displaystyle \frac{\tau }{E}}A\\ ={\displaystyle \frac{\mathrm{4}}{E}}a(\tau +\sigma )\end{array}

Combining PD with the COD method as shown in Fig. 6, the opening displacement between adjacent material points can be expressed as:

(20){\mathrm{\Delta }}_{cr}=\left|{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}\right|-\left|{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\right|={\delta }_{cr}

where {\mathrm{\Delta }}_{cr} is the critical opening displacement of PD adjacent material points and {\delta }_{cr} is the critical crack opening displacement of the fracture mechanics.

In order to include damage initiation in the material response, a history-dependent scalar-valued functional {\varphi }_{\mathrm{\Delta }} can be introduced

(21)

and the force density vector t_{(k)(j)}can be modified as

(22)t_{(k)(j)}=\mathrm{2}\delta \left\{ad{\displaystyle \frac{{\Lambda }_{(k)(j)}}{\left|{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\right|}}{\theta }_{(k)}+b\cdot \mathrm{\Delta }({\mathit{x}}_{(j)}-{\mathit{x}}_{(k)})s_{(k)(j)}\right\}{\displaystyle \frac{{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}}{\left|{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}\right|}}

where {\theta }_{(k)} is the dilatation term

(23){\theta }_{(k)}=d\delta {\displaystyle \sum _i^N{\Lambda }_{(k)(i)}}{\varphi }_{\mathrm{\Delta }}({\mathit{x}}_{(j)}-{\mathit{x}}_{(k)},t)s_{(k)(i)}{V}_i

When the object under the external load with the time changes, we continue to calculate the open displacement of the crack tip. When the displacement between the two particles satisfies \mathrm{\Delta }\ge {\mathrm{\Delta }}_{cr}, the PD force t_{(k)(j)} is set to zero and the crack propagates (based on the history value {\varphi }_{\mathrm{\Delta }}=\mathrm{0}). When the crack tip opened displacement , the force t_{(k)(j)} is unequal to zero and {\varphi }_{\mathrm{\Delta }}=\mathrm{1}.

When the crack does not propagate, the local damage function \psi (\mathit{x},t) is introduced into the PD-COD model in order to express the relationship between the crack opening displacement and damage of the particles

(24)\psi (\mathit{x},t)=\mathrm{1}-{\displaystyle \frac{{\displaystyle {\int }_H{\varphi }_{\mathrm{\Delta }}({\mathit{x}}_{(j)}-{\mathit{x}}_{(k)},t)\mathrm{d}V}}{{\displaystyle {\int }_H\mathrm{d}V}}}

The local damage of the PD-COD model ranges from 0 to 1 as shown in Fig. 7. When the crack opening displacement satisfies \mathrm{\Delta }\ge {\mathrm{\Delta }}_{cr}, then the local damage is one and all the interactions initially with the point have been eliminated. A local damage value of \psi (\mathit{x},t)=\mathrm{0} means that all interactions are intact. However, the creation of a crack terminates half of the interactions with its horizon, resulting in a local damage value of one-half as shown in Fig. 7 (b).

Numerical example

In this paper, four numerical examples as illustrated in Figures 8 to 11 are studied in order to demonstrate the performance of the new fracture criterion.

Examples 1 Double-notched specimen made of Q345 steel under uniaxial tension

Q345 material is a low alloy high strength structural steel, with elasticity modulusE=\mathrm{203}\mathrm{GPa}, poisson’s ratio\nu =\mathrm{0.3}, elongation\delta =\mathrm{27.96}\% and density \rho =\mathrm{7850}{\mathrm{kg/m}}^3. The length and width of the specimen are 70mm and 40mm, respectively as shown in Fig. 8. The specimen is loaded under uniaxial tension with a constant loading rate of \mathrm{2.217}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{m/s}ensuring quasi-static conditions. The crack length is 10 mm and three different crack size distances in loading directions are tested according to Table 1. The fractured specimens of these experiments are illustrated in Fig. 12 and compared to the fractured specimen from FEM simulations in Fig. 13. Furthermore, Figure 14 shows results obtained by he C3D8R XFEM element in ABAQUS. The calculated results by XFEM agree fairly well with experimental results in Fig. 12 and the FEM in Fig. 13 though both numerical simulations are not able to capture the curvature of the crack in the second specimen.

Fig.12 is made by experiment, and Fig.13 is made by FEM. We used the maximum strain energy to model the crack propagation in FEM. When the strain energy of crack tip reached the critical strain energy, we scattered the units at crack tip. So, the crack will propagation forward for some distance.

As shown in Fig.12 and Fig.13. if the longitudinal initial crack spacing is zero, the two propagating cracks join in a horizontal line. For a longitudinal crack spacing of 10 mm, the crack initially propagates horizontally. When the horizontal and longitudinal distances of the two crack tips are equal, the cracks extend along a 45°direction into a slanted crack until they join. No crack coalescence is observed for the specimen with a 20 mm longitudinal crack spacing.

The crack propagation results are shown in Fig. 15, using the PD elongation s_{ij} crack fracture criterion. The results of three different specimens are consistent with the results of the experiment, FEM and XFEM. It indicates the PD elongation criterion can simulated the multi-crack propagation and fusion.

As shown in Fig. 16, the results of the crack propagation of the three specimens simulated by PD-COD fracture criterion are consistent with the results of FEM, XFEM and PD elongation. The variety of simulation results are very similar. Therefore, we can infer that the PD-COD fracture criterion can simulate the crack propagation and fusion.

Although both of the fracture criteria can simulate the fracture process, there are nuances between two methods at the simulation the multi-crack propagation and interaction process. It can be seen from Fig. 13 (c) and Fig. 14 (c) the two cracks have a longitudinal distance of 20mm. Though the cracks propagation path is along the horizontal straight line, two cracks however, have mutual influence during the expansion process.

The stress distribution in Fig. 13 (c) and Fig. 14 (c) are consistent. The stress in the middle of two cracks is obviously higher than other parts. There are local damage in the middle test, but the damage value is small and no cracks are formed.

It can be seen that the two crack propagation does not affect the stress state at the middle of the model in Fig. 15 (c). However from Fig. 16 (c), the crack propagation is not only able to show the crack propagation path, but also shows the damage distribution in the middle part of the model. The most important result is that the damage distribution is consistent with the results of FEM and XFEM.

Examples 2 Double center crack problem

We consider a specimen with two double-center cracks under uniaxial tension with a constant loading rate of \mathrm{20}\mathrm{m/s}. The model is shown in Fig. 10. The material parameters from the previous example are adopted. The thickness of the plate is 0.05m, the initial lengths are 0.01m for both cracks. The longitudinal spacing of crack is 4mm (case 1), 5mm (case 2) and 10mm (case 3), respectively. A tensile load is applied at both ends of the plate in vertical direction. The particle spacing in all PD simulations is 0.5mm. The influence of the longitudinal crack spacing on damage rate and crack propagation path is analyzed subsequently.

Fig. 17 and Fig. 18 present the fracture patterns for the PD extension s_{ij} crack propagation criterion and the PD-COD crack propagation criterion. Comparing Fig. 17 with Fig. 18, both methods yield similar results. An increaseing distance of the crack leads to a higher curvature of the crack path.

Examples 3 Slanted center crack problem

The third example is a specimen with a slanted center crack. Three different specimen with different crack angle versus the horizontal axis are tested. All specimen are loaded under uniaxial tension with a constant loading rate of \mathrm{20}\mathrm{m/s}. The model is illustrated in Fig. 10.The material parameters from the first example are adopted. The thickness of the plate is 0.05m, the length of one crack is 0.01m. The inclination angle a of the slanted crack is 30° (case 1), 45° (case 2) and 60° (case 3), respectively. A tensile load is applied at both ends of the plate in a vertical direction. The particle spacing in all PD simulations is 0.5mm. The influence of the longitudinal crack spacing on the damage rate and crack propagation path was analyzed by the PD elongations_{ij}criterion and the new COD criteria.

Fig. 19 and Fig. 20 show the final fracture pattern from the PD simulations based on the extension s_{ij} crack propagation criterion and the COD criterion, respectively. As expected, the final crack path is perpendicular to the loading direction in this mode I dominated fracture problem. The results of both cracking criterion agree well.

Examples 4 Specimen with two slanted initial cracks

The fourth example is a specimen with two slanted initial cracks. The specimen is again loaded under uniaxial tension with a constant loading rate of \mathrm{20}\mathrm{m/s}and the material parameters are the same as in all other examples. The model is presented in Fig. 11. The thickness of the plate is 0.05m, the length of one crack is 0.01m. The inclination angle a of both slanted cracks are 30° (case 1), 45° (case 2) and 60° (case 3), respectively. The particle spacing in all PD simulations is 0.5mm.The influence of the longitudinal crack spacing on damage rate and crack propagation path was analyzed by PD elongations_{ij}criterion and PD-COD criterion and compared with two kinds of crack propagation path.

Fig. 21 and Fig. 22 illustrates the crack propagation path for the three cases considered here. The results from the PD simulations based on the extension s_{ij} crack propagation criterion and our novel criterion agree well. In all cases crack shielding occurs and the crack propagates towards the boundary of the plates. The crack shielding is more pronounced with increasing inclination angle for our new cracking criterion. Such a tendency is not apparent for the standard cracking criterion in PD. For case 3, the standard cracking criterion provides a short slightly inclined crack path which does not seem reasonable due to the mode I dominated fracture mode. Such an artifact does not occur for our new fracture criterion.

Conclusions

In this paper, a new PD crack propagation criterion is proposed. It is based on the opening displacement COD method. The key idea is when the opened displacement of PD adjacent material points equal to the critical crack opened displacement at crack tip (the critical COD), the bonds between particles are broken and a crack is formed. A local damage formulation\psi (x,t) is also introduced in analogy to the original PD damage formulation to simulate the local damage.

Four examples are studied to verify the correctness of the novel cracking criterion. The first example is compared to experimental data and results of other numerical methods and shows excellent agreement. The following three examples are the double parallel crack, the monoclinic crack and the double oblique crack. They further verify that PD-COD can accurately simulate the crack propagation. From the four examples, we can conclude that the crack propagation path will be affected by both the crack longitudinal spacing and crack oblique angle. In the future, the COD criterion can be extended to nonlinear materials and SB-peridynamics which should be straightforward.

References

[1] Amiri F, Millan D, Arroyo M, Silani M, Rabczuk T. Fourth order phase-field model for local max-ent approximants applied to crack propagation. Computer Methods in Applied Mechanics and Engineering, 2016, 312(C): 254–275

[2] Amiri F, Millán D, Shen Y, Rabczuk T, Arroyo M. Phase-field modeling of fracture in linear thin shells. Theoretical and Applied Fracture Mechanics, 2014, 69: 102–109

[3] Areias P, Rabczuk T, de Sá J C. A novel two-stage discrete crack method based on the screened Poisson equation and local mesh re_nement. Computational Mechanics, 2016, 58(6): 1003–1018

[4] Thai T Q, Rabczuk T, Bazilevs Y, Meschke G. A higher-order stress-based gradient-enhanced damage model based on isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 2016, 304: 584–604

[5] Budarapu P, Gracie R, Bordas S, Rabczuk T. An adaptive multiscale method for quasi-static crack growth. Computational Mechanics, 2014, 53(6): 1129–1148

[6] Budarapu P R, Gracie R, Yang S W, Zhuang X, Rabczuk T. Efficient Coarse Graining in Multiscale Modeling of Fracture. Theoretical and Applied Fracture Mechanics, 2014, 69: 126–143

[7] Silani M, Talebi H, Hamouda A S, Rabczuk T. Nonlocal damage modelling in clay/epoxy nanocomposites using a multiscale approach. Journal of Computational Science, 2016, 15: 18–23

[8] Talebi H, Silani M, Rabczuk T. Concurrent Multiscale Modelling of Three Dimensional Crack and Dislocation Propagation. Advances in Engineering Software, 2015, 80: 82–92

[9] Silani M, Ziaei-Rad S, Talebi H, Rabczuk T. A Semi-Concurrent Multiscale Approach for Modeling Damage in Nanocomposites. Theoretical and Applied Fracture Mechanics, 2014, 74: 30–38

[10] Talebi H, Silani M, Bordas S, Kerfriden P, Rabczuk T. A Computational Library for Multiscale Modelling of Material Failure. Computational Mechanics, 2014, 53(5): 1047–1071

[11] Talebi H, Silani M, Bordas S P A, Kerfriden P, Rabczuk T. Molecular Dynamics/XFEM Coupling by a Three-Dimensional Extended Bridging Domain with Applications to Dynamic Brittle Fracture. International Journal for Multiscale Computational Engineering, 2013, 11(6): 527–541

[12] Areias P, Rabczuk T, Msekh M. Phase-field analysis of finite-strain plates and shells including element subdivision. Computer Methods in Applied Mechanics and Engineering, 2016, 312(C): 322–350

[13] Areias P, Msekh M A, Rabczuk T. Damage and fracture algorithm using the screened Poisson equation and local remeshing. Engineering Fracture Mechanics, 2016, 158: 116–143

[14] Areias P M A, Rabczuk T, Camanho P P. Finite strain fracture of 2D problems with injected anisotropic softening elements. Theoretical and Applied Fracture Mechanics, 2014, 72: 50–63

[15] Areias P, Rabczuk T, Dias-da-Costa D. Element-wise fracture algorithm based on rotation of edges. Engineering Fracture Mechanics, 2013, 110: 113–137

[16] Areias P, Rabczuk T, Camanho P P. Initially rigid cohesive laws and fracture based on edge rotations. Computational Mechanics, 2013, 52(4): 931–947

[17] Areias P, Rabczuk T. Finite strain fracture of plates and shells with configurational forces and edge rotation. International Journal for Numerical Methods in Engineering, 2013, 94(12): 1099–1122

[18] Dolbow J E. An extended finite element method with discontinuous enrichment for applied mechanics. Northwestern university, 1999.

[19] Fries T P, Belytschko T. The extended/generalized finite element method: an overview of the method and its applications. International Journal for Numerical Methods in Engineering, 2010, 84(3): 253–304

[20] Ghorashi S, Valizadeh N, Mohammadi S, Rabczuk T. T-spline based XIGA for Fracture Analysis of Orthotropic Media. Computers & Structures, 2015, 147: 138–146

[21] Nguyen-Thanh N, Valizadeh N, Nguyen M N, Nguyen-Xuan H, Zhuang X, Areias P, Zi G, Bazilevs Y, De Lorenzis L, Rabczuk T. An extended isogeometric thin shell analysis based on Kirchhoff-Love theory. Computer Methods in Applied Mechanics and Engineering, 2015, 284: 265–291

[22] Shi G H. Numerical manifold method and discontinuous deformation analysis. 1997.

[23] Cai Y, Zhuang X, Zhu H. A generalized and efficient method for finite cover generation in the numerical manifold method. International Journal of Computational Methods, 2013, 10(05): 1350028

[24] Belytschko T, Lu Y Y, Gu L. Element‐free Galerkin methods. International Journal for Numerical Methods in Engineering, 1994, 37(2): 229–256

[25] Liu W K, Jun S, Zhang Y F. Reproducing kernel particle methods. International Journal for Numerical Methods in Fluids, 1995, 20(8‐9): 1081–1106

[26] Li S, Liu W K. Meshfree and particle methods and their applications. Applied Mechanics Reviews, 2002, 55(1): 1–34

[27] Nguyen V P, Rabczuk T, Bordas S, Duflot M. Meshless methods: a review and computer implementation aspects. Mathematics and Computers in Simulation, 2008, 79(3): 763–813

[28] Zhuang X, Augarde C E, Mathisen K M. Fracture modeling using meshless methods and level sets in 3D: framework and modeling. International Journal for Numerical Methods in Engineering, 2012, 92(11): 969–998

[29] Amiri F, Anitescu C, Arroyo M, Bordas S, Rabczuk T. XLME interpolants, a seamless bridge between XFEM and enriched meshless methods. Computational Mechanics, 2014, 53(1): 45–57

[30] Rabczuk T, Areias P M A. A meshfree thin shell for arbitrary evolving cracks based on an external enrichment. CMES-Computer Modeling in Engineering and Sciences, 2006, 16(2): 115–130

[31] Rabczuk T, Bordas S, Zi G. A three-dimensional meshfree method for continuous multiple-crack initiation, propagation and junction in statics and dynamics. Computational Mechanics, 2007, 40(3): 473–495

[32] Rabczuk T, Areias P M A, Belytschko T. A meshfree thin shell method for nonlinear dynamic fracture. International Journal for Numerical Methods in Engineering, 2007, 72(5): 524–548

[33] Rabczuk T, Zi G. A meshfree method based on the local partition of unity for cohesive cracks. Computational Mechanics, 2007, 39(6): 743–760

[34] Rabczuk T, Bordas S, Zi G. On three-dimensional modelling of crack growth using partition of unity methods. Computers & Structures, 2010, 88(23-24): 1391–1411

[35] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A geometrically non-linear three dimensional cohesive crack method for reinforced concrete structures. Engineering Fracture Mechanics, 2008, 75(16): 4740–4758

[36] Rabczuk T, Gracie R, Song J H, Belytschko T. Immersed particle method for fluid-structure interaction. International Journal for Numerical Methods in Engineering, 2010, 81(1): 48–71

[37] Bordas S, Rabczuk T, Zi G. Three-dimensional crack initiation, propagation, branching and junction in non-linear materials by extrinsic discontinuous enrichment of meshfree methods without asymptotic enrichment. Engineering Fracture Mechanics, 2008, 75(5): 943–960

[38] Rabczuk T, Song J H, Belytschko T. Simulations of instability in dynamic fracture by the cracking particles method. Engineering Fracture Mechanics, 2009, 76(6): 730–741

[39] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A simple and robust three-dimensional cracking-particle method without enrichment. Computer Methods in Applied Mechanics and Engineering, 2010, 199(37-40): 2437–2455

[40] Rabczuk T, Belytschko T. Cracking particles: a simpli_ed meshfree method for arbitrary evolving cracks. International Journal for Numerical Methods in Engineering, 2004, 61(13): 2316–2343

[41] Rabczuk T, Belytschko T. A three dimensional large deformation meshfree method for arbitrary evolving cracks. Computer Methods in Applied Mechanics and Engineering, 2007, 196(29-30): 2777–2799

[42] Erdogan F, Sih G C. On the crack extension in plates under plane loading and transverse shear. Journal of Basic Engineering, 1963, 85(4): 519–527

[43] Sih G C. Strain-energy-density factor applied to mixed mode crack problems. International Journal of Fracture, 1974, 10(3): 305–321

[44] Hussain M A, Pu S L, Underwood J. Strain energy release rate for a crack under combined mode I and mode II. In: Proceedings of the 1973 National Symposium on Fracture Mechanics, Part II. ASTM International, 1974.

[45] Ayatollahi M R, Abbasi H. Prediction of fracture using a strain based mechanism of crack growth. Building Research Journal, 2001, 49: 167–180

[46] Awaji H, Sato S. Combined Mode Fracture Toigliess measurement by the Disk Test. Journal of Engineering Materials and Technology, 1978, 100(2): 175

[47] Shetty D K, Rosenfield A R, Duckworth W H. Mixed-mode fracture in biaxial stress state: application of the diametral-compression (Brazilian disk) test. Engineering Fracture Mechanics, 1987, 26(6): 825–840

[48] Chang S H, Lee C I, Jeon S. Measurement of rock fracture toughness under modes I and II and mixed-mode conditions by using disc-type specimens. Engineering Geology, 2002, 66(1-2): 79–97

[49] Aliha M R M, Ashtari R, Ayatollahi M R. Mode I and mode II fracture toughness testing for a coarse grain marble. Applied Mechanics and Materials, 2006, 5: 181–188.

[50] Richard H A, Benitz K.A loading device for the creation of mixed mode in fracture mechanics. international Journal of Fracture, 1983, 22(2): R55–R58.

[51] Arcan M, Hashin Z, Voloshin A. A method to produce uniform plane-stress states with applications to fiber-reinforced materials. Experimental Mechanics, 1978, 18(4): 141–146

[52] Zipf R K Jr, Bieniawski Z T. Mixed mode testing for fracture toughness of coal based on critical-energy-density. In: The 27th US Symposium on Rock Mechanics (USRMS). American Rock Mechanics Association, 1986.

[53] Silling S A. Reformulation of elasticity theory for discontinuities and long-range forces. Journal of the Mechanics and Physics of Solids, 2000, 48(1): 175–209

[54] Silling S A, Zimmermann M, Abeyaratne R. Deformation of a peridynamic bar. Journal of Elasticity, 2003, 73(1-3): 173–190

[55] Ren H, Zhuang X, Cai Y, Rabczuk T. Dual-Horizon Peridynamics. International Journal for Numerical Methods in Engineering, 2016, 108(12): 1451–1476

[56] Ren H, Zhuang X, Rabczuk T. Dual-horizon peridynamics: A stable solution to varying horizons. Computer Methods in Applied Mechanics and Engineering, 2017, 318: 762–782

[57] Silling S A. Dynamic fracture modeling with a meshfree peridynamic code. Computational Fluid & Solid Mechamcs,. 2003, 641–644

[58] Weckner O, Brunk G, Epton M A, Silling S A. Green’s functions in non-local three-dimensional linear elasticity In: Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. The Royal Society, 2009

[59] Yu K, Xin X J, Lease K B. A new adaptive integration method for the peridynamic theory. Modelling and Simulation in Materials Science and Engineering, 2011, 19(4): 045003

[60] Kilic B. Peridynamic theory for progressive failure prediction in homogeneous and heterogeneous materials. ProQuest, 2008.

[61] Silling S A. Reformulation of elasticity theory for discontinuities and long-range forces. Journal of the Mechanics and Physics of Solids, 2000, 48(1): 175–209

[62] Silling S A, Askari E. A meshfree method based on the peridynamic model of solid mechanics. Computers & Structures, 2005, 83(17-18): 1526–1535

[63] Silling S A, Bobaru F. Peridynamic modeling of membranes and fibers. International Journal of Non-linear Mechanics, 2005, 40(2-3): 395–409

[64] Silling S A, Askari E. A meshfree method based on the peridynamic model of solid mechanics. Computers & Structures, 2005, 83(17-18): 1526–1535

[65] Foster J T, Silling S A, Chen W. An energy based failure criterion for use with peridynamic states. International Journal for Multiscale Computational Engineering, 2011, 9(6): 675–688

[66] Silling S A, Bobaru F. Peridynamic modeling of membranes and fibers. International Journal of Non-linear Mechanics, 2005, 40(2-3): 395–409

[67] Ayatollahi M R, Aliha M R M. Analysis of a new specimen for mixed mode fracture tests on brittle materials. Engineering Fracture Mechanics, 2009, 76(11): 1563–1573

[68] Shen F, Zhang Q, Huang D. Damage and failure process of concrete structure under uni-axial tension based on peridynamics modeling. Chinese Journal of Computational Mechanics, 2013, 30: 79–83

[69] Zhou X P, Shou Y D. Numerical Simulation of Failure of Rock-Like Material Subjected to Compressive Loads Using Improved Peridynamic Method. International Journal of Geomechanics, 2016: 04016086

[70] Ren H, Zhuang X, Rabczuk T. A new peridynamic formulation with shear deformation for elastic solid. Journal of Micromechanics and Molecular Physics, 2016, 1(02): 1650009

[71] Oterkus E. Peridynamic theory for modeling three-dimensional damage growth in metallic and composite structures. University of Arizona, 2010.

[72] Bilby B A, Cottrell A H, Swinden K H. The spread of plastic yield from a notch. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. The Royal Society, 1963, 272(1350): 304–314

[73] Bilby B A, Cottrell A H, Smith E, Swinden K H. Plastic yielding from sharp notches. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. The Royal Society, 1964, 279(1376): 1–9

[74] Dugdale D S. Yielding of steel sheets containing slits. Journal of the Mechanics and Physics of Solids, 1960, 8(2): 100–104

[75] Barenblatt G I. The mathematical theory of equilibrium cracks in brittle fracture. Advances in Applied Mechanics, 1962, 7: 55–129

[76] Zeng S. Bilateral crack growth behavior of Q345 steel under uniaxial tensile load. Guangxi University, 2014

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg