Peridynamics versus XFEM: a comparative study for quasi-static crack problems

Abstract

Peridynamics (PD) is a nonlocal continuum theory based on integro-differential equations without spatial derivatives. The fracture criterion is implicitly incorporated in the PD theory and fracture is a natural outcome of the simulation. However, capturing of complex mixed-mode crack patterns has been proven to be difficult with PD. On the other hand, the extended finite element method (XFEM) is one of the most popular methods for fracture which allows crack propagation with minimal remeshing. It requires a fracture criterion which is independent of the underlying discretization though a certain refinement is needed in order to obtain suitable results. This article presents a comparative study between XFEM and PD. Therefore, two examples are studied. The first example is crack propagation in a double notched specimen under uniaxial tension with different crack spacings in loading direction. The second example is the specimens with two center cracks. The results show that PD as well as XFEM are well suited to capture this type of behaviour.

Keywords

XFEM / peridynamic / bilateral crack / parallel double cracks / nonlocal theory

Introduction

For crack propagation problems in the engineering structure, it is difficult to obtain an analytical solution for the stress, strain and displacement field and to determine the direction of crack propagation. With the advancements in computer science and computational science, the finite element method has become a powerful tool which was also applied to fracture problems. Analyzing the extension characteristics of cracks and predicting the service life of structures are of high importance in engineering. The traditional finite element method (FEM) is not well suited for crack propagation problems as the crack can only propagating along existing element edges by node splitting. The results will highly depend on the mesh topology and therefore remeshing techniques are commonly exploited. Besides some recent advancements in remeshing techniques, see for instance the work by Areias and co-workers [¹–⁶], their implementation to 3D and complex fracture patterns remains still a challenge.

The extended finite element method (XFEM) [⁷] allows crack propagation without remeshing. Though XFEM can principally capture crack propagation without remeshing, a certain mesh refinement is still required in order to obtain results with satisfactory accuracy. For problems in linear elastic fracture mechanics (LEFM), it was shown that its accuracy is related to the position of the crack tip [⁸,⁹]. It is also complicated to deal with complex fracture patterns such as crack coalescence or crack branching. Daux et al. developed strategies to deal with such type of problems [⁸]. Nevertheless, reliable fracture criteria for such applications are still missing. While the original XFEM was developed for problems in LEFM, it was meanwhile applied to non-linear problems including cohesive cracks; see for instance the work by Fang et al. [¹⁰,¹¹]. The phantom node method [¹²,¹³] is similar to XFEM and it was shown in Ref. [¹⁴] that the crack kinematics of the phantom node method can be derived from XFEM. However, it has some advantages: 1) It models the crack by overlapping elements (not by enrichment) and therefore, existing element formulations can be readily exploited. In contrast, efficient FEM formulations for problems with constraints for instance need to be redeveloped and tested for the ‘standard’ XFEM. 2) No ‘mixed terms’ are needed in the phantom node method and 3) the implementation is easier. Nonetheless, the phantom node method requires that the crack propagates through an entire element though some special crack tip elements have been developed which allow the crack tip inside an element [¹⁵,¹⁶]. Both, the phantom node method as well as XFEM has been implemented in the commercial software package ABAQUS.

Meshless methods (MM) eliminate the mesh and can therefore deal well with fracture problems as well, see e.g., the contributions by Rabczuk [¹⁷–²⁵] or Zhuang and co-workers [²⁶–²⁸]. However, as XFEM, they still require fracture criteria. An interesting and powerful meshfree method for fracture is the cracking particles method (CPM) [²⁹–³¹] which can—similarly to PD—also complex fracture pattern such as crack coalescence, crack nucleation and crack coalescence without any special criteria. In the CPM, the crack is modelled by a set of crack segments and therefore the crack path is not continuous. It has been shown by [³²,³³] that certain meshless methods are nearly identical to the so-called state-based PD formulation which will be explained in the next section.

Peridynamics (PD) is a new numerical method for characterizing materials based on nonlocal models. It was first proposed in 2000 by Silling [³⁴,³⁵]. For a linear elastic solid, the original bond-based PD formulation [³⁴] modeled fracture by equating the critical energy release rate of the continuum with the one from PD theory which leads finally to a critical stretch between the particles as fracture criterion which can be implemented with ease. However, the extension to nonlinear materials and complex mixed-mode fracture cannot easily be treated with the bond-based PD. Therefore, the state-based PD was developed which can handle more complex constitutive models. However, it was shown in [³²,³³] that such models are quite similar to meshfree methods. Nonetheless, the theory has been further developed and applied to some numerical simulation of discontinuous damage problems in recent years [³⁶]. The third generation of PD is called dual-horizon peridynamics (DH-PD) and allows for adaptive refinement [³⁷]; it can be applied to bond-based PD as well as state-based PD.

Silling and Bobaru proposed a weighted local function of the particle weight method to determine the particle damage problem [³⁸]. In multi-scale problems, there are some big challenges. Some researchers have done some work on it [³⁹–⁴⁷], see also the contribution from Costa et al. on a multi-scale PD finite element model [⁴⁸]. Jung and Seok modelled fatigue cracks and computed stress intensity factors of stratified inhomogeneous materials [⁴⁹]. Shen et al. simulated the cracking process of concrete [⁵⁰]. Gu and Zhou captured the fracture process of a plate with a round hole and provided new ideas and new methods for material failure [⁵¹]. Cheng and Liu focused on dynamic fracture of functionally graded materials and analyzed the dynamic response of those materials [⁵²]. The article of Zhou and Shou [⁵³] deals with rock-like materials while Ren et al. proposed a new criteria for damage determination of shear deformation [⁵⁴].

In this paper, a comparative study of PD and XFEM is done for problems involving crack propagation including crack interactions for quasi-static problems. Two examples are considered: crack propagation in a double-notched specimen under uniaxial tension and a specimen with a double center crack. Both specimen are made of Q345 steel and experimental data is available. We use the XFEM capabilities in ABAQUS to solve this problem while our own PD program is employed to solve these problems. Subsequently, the PD formulation is reviewed. Section 3 presents the numerical examples before the manuscript closes with conclusions in section 4.

The peridynamic (PD) theory

Basic model

PD discretizes objects with many particles. In this article, we employ the state-based PD formulation which is subsequently summarized. The key in state-based PD is how to compute the force state \underline{\mathit{T}}depending on the deformation state \underline{\mathit{Y}}between the particles in the current 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{\mathit{Y}}({\mathit{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 state for material point {\mathit{x}}_{(k)}depends on other material points within its horizon. It can also be expressed as

(3)\underline{\mathit{T}}({\mathit{x}}_{(k)},t)=\underline{\mathit{T}}\left({\underline{\mathit{Y}}}_{(k)},t\right).

Figure 1 shows the force density vector {\mathit{t}}_{(k)(j)} that the material point {\mathit{x}}_{(j)}exerts on the material point {\mathit{x}}_{(k)}. It is given as

(4){\mathit{t}}_{(k)(j)}({\mathit{u}}_{(j)}-{\mathit{u}}_{(k)},{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)},t)=\underline{\mathit{T}}({\mathit{x}}_{(k)},t)\langle {\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\rangle .

The interaction between material points {\mathit{x}}_{(k)}and {\mathit{x}}_{(j)} can be derived from a scalar-valued micropotential, w_{(k)(j)}

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

where {\mathit{y}}_{(k)}is the positon vector of point {\mathit{x}}_{(k)}in the deformed configuration and {\mathit{y}}_{({\mathrm{1}}^k)}is the positon vector of the first material point that interacts with point {\mathit{x}}_{(k)}.

Then strain energy density, W_{(k)}, of material point {\mathit{x}}_{(k)}can be expressed as the summation of micro-potentials, w_{(k)(j)}within its horizon of material points, {\mathit{x}}_{(k)} and is given by

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 )}

(6)+{\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)},\mathrm{...})\cdot V_{(j)},

in which w_{(k)(j)}=\mathrm{0}for k=j.

Applying the principle of virtual work, the PD equations of motion at material point {\mathit{x}}_{(k)}can be derived as

(7)\delta {\displaystyle {\int }_{t_{\mathrm{0}}}^{t_{\mathrm{1}}}(T-U)}\mathrm{d}t=\mathrm{0},

where Tand Urepresent the total kinetic and potential energies in the body. Solving for the Lagrange’s equation, this principle is satisfied as

(8){\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{\partial L}{\partial {\dot{\mathit{u}}}_{(k)}}}\right)-{\displaystyle \frac{\partial L}{\partial {\mathit{u}}_{(k)}}}=\mathrm{0},

where the Lagrangian Lis defined as

(9)L=T-U.

By summing the kinetic and potential energies of all material points, the total kinetic and potential energies in the body is obtained by

(10)T={\displaystyle \sum _{i=\mathrm{1}}^{\infty }\frac{\mathrm{1}}{\mathrm{2}}}{\rho }_{(i)}{\dot{\mathit{u}}}_{(i)}\cdot {\dot{\mathit{u}}}_{(i)}{V}_{(i)},

and

(11)U={\displaystyle \sum _{i=\mathrm{1}}^{\infty }{W}_{(i)}{V}_{(i)}-}{\displaystyle \sum _{i=\mathrm{1}}^{\infty }({\mathit{b}}_{(i)}\cdot {\mathit{u}}_{(i)})}{V}_{(i)}.

Through the Eqs. (9), (10) and (11), Lagrange’s equation of the material point {\mathit{x}}_{(k)}can be recast as

{\rho }_{(k)}{\mathit{u}}_{(k)}={\displaystyle \sum _{j=\mathrm{1}}^{\infty }[{\mathit{t}}_{(k)(j)}({\mathit{u}}_{(j)}-{\mathit{u}}_{(k)},{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)},t)}

(12)-{\mathit{t}}_{(j)(k)}({\mathit{u}}_{(k)}-{\mathit{u}}_{(j)},{\mathit{x}}_{(k)}-{\mathit{x}}_{(j)},t)]V_{(j)}+{\mathit{b}}_{(k)},

where

(13){\mathit{t}}_{(k)(j)}({\mathit{u}}_{(j)}-{\mathit{u}}_{(k)},{\mathit{x}}_{(j)}-{\mathit{x}}_{(k)},t)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{1}}{V_{(k)}}}\left({\displaystyle \sum _{i=\mathrm{1}}^{\infty }{\displaystyle \frac{\partial w_{(k)(i)}}{\partial ({\mathit{y}}_{(i)}-{\mathit{y}}_{(k)})}}{V}_{(i)}}\right),

and

(14){\mathit{t}}_{(j)(k)}({\mathit{u}}_{(k)}-{\mathit{u}}_{(j)},{\mathit{x}}_{(k)}-{\mathit{x}}_{(j)},t)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{1}}{V_{(j)}}}\left({\displaystyle \sum _{i=\mathrm{1}}^{\infty }{\displaystyle \frac{\partial w_{(i)(k)}}{\partial ({\mathit{y}}_{(k)}-{\mathit{y}}_{(i)})}}{V}_{(i)}}\right).

Finally, the PD equations of motion at material point {\mathit{x}}_{(k)}can be rewritten as

(15){\rho }_{(k)}{\ddot{\mathit{u}}}_{(k)}={\displaystyle \sum _{j=\mathrm{1}}^{\infty }\left(\underline{\mathit{T}}({\mathit{x}}_{(k)},t)\langle {\mathit{x}}_{(j)}-{\mathit{x}}_{(k)}\rangle -\underline{\mathit{T}}({\mathit{x}}_{(j)},t)\langle {\mathit{x}}_{(k)}-{\mathit{x}}_{(j)}\rangle \right)V_{(j)}+{\mathit{b}}_{(k)}}.

Because the volume of each material point V_{(j)}is infinitesimally small, the infinite summation can be expressed as integration while considering only the material points within the horizon,

(16){\displaystyle \sum _{j=\mathrm{1}}^{\infty }(\cdot )V_{(j)}}\to {\displaystyle {\int }_H(\cdot )\mathrm{d}{V}^{\prime }\to }{\displaystyle {\int }_H(\cdot )\mathrm{d}H}.

Using the expression of Eq. (16), with the domain of integral H and the thickness of the plate, the PD dynamic equation is expressed as

\rho (\mathit{x}){\ddot{\mathit{u}}}_{(k)}={\displaystyle {\int }_H(\mathit{t}({\mathit{u}}^{\prime }-\mathit{u},{\mathit{x}}^{\prime }-\mathit{x},t)-{\mathit{t}}^{\prime }(\mathit{u}-{\mathit{u}}^{\prime },\mathit{x}-\mathit{x},t))\mathrm{d}H}

(17){\displaystyle +\mathit{b}(\mathit{x},t)},

in which \mathit{b}(\mathit{x},t)denotes the body force.

Local damage model

Fracture is modelled by simply removing the interaction of two particles when a certain stretch is exceed as illustrated in Fig. 2. All of the micropotentials between the material points {\mathit{x}}_{(k^{+})} and {\mathit{x}}_{(j^{-})} are located above and below the new crack surface can be obtained as

(18)w_{(k^{+})(j^{-})}=\left(\mathrm{4}a{d}^{\mathrm{2}}{\delta }^{\mathrm{2}}\left({\displaystyle \sum _{i=\mathrm{1}}^{K^{-}}{s}_c^{\mathrm{2}}{V}_{(i)}}\right)+\mathrm{4}\delta b{s}_c^{\mathrm{2}}\left|{\mathit{x}}_{(j^{-})}-{\mathit{x}}_{(k^{+})}\right|\right),

and

(19)w_{(j^{-})(k^{+})}=\left(\mathrm{4}a{d}^{\mathrm{2}}{\delta }^{\mathrm{2}}\left({\displaystyle \sum _{i=\mathrm{1}}^{J^{+}}{s}_c^{\mathrm{2}}{V}_{(i)}}\right)+\mathrm{4}\delta b{s}_c^{\mathrm{2}}\left|{\mathit{x}}_{(k^{+})}-{\mathit{x}}_{(j^{-})}\right|\right),

in which, a and d are the parameters of PD.

Then, total strain energy required to remove all of the interactions across the newly created crack surface can be recast as

(20)W={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \sum _{k=\mathrm{1}}^{K^{+}}\frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \sum _{j=\mathrm{1}}^{J^{-}}{w}_{(k^{+})(j^{-})}{V}_{(k^{+})}{V}_{(j^{-})}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \sum _{k=\mathrm{1}}^{K^{+}}\frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \sum _{j=\mathrm{1}}^{J^{-}}{w}_{(j^{-})(k^{+})}{V}_{(j^{-})}{V}_{(k^{+})}},

where K^{+}, J^{-}indicate the number of material points {\mathit{x}}_{(k^{+})}and {\mathit{x}}_{(j^{-})} above and below the crack surface.

The corresponding strain energy release rate G_cis given by

(21)G_c={\displaystyle \frac{s_c^{\mathrm{2}}{\displaystyle \sum _{k=\mathrm{1}}^{K^{+}}{\displaystyle \sum _{j=\mathrm{1}}^{J^{-}}\left(\mathrm{2}\delta b\left|{\mathit{x}}_{(j^{-})}-{\mathit{x}}_{(k^{+})}\right|+a{d}^{\mathrm{2}}{\delta }^{\mathrm{2}}\left({\displaystyle \sum _{i=\mathrm{1}}^{K^{-}}{V}_{(i)}}+{\displaystyle \sum _{i=\mathrm{1}}^{J^{+}}{V}_{(i)}}\right)\right)}}{V}_{(k^{+})}{V}_{j^{-}}}{A}},

where s_cindicates the critical stretch.

For a 2D analysis, the expression for the critical energy release rate becomes

(22)G_c=\mathrm{2}h{\displaystyle {\int }_{\mathrm{0}}^{\delta }\left\{{\displaystyle {\int }_z^{\delta }{\displaystyle {\int }_{\mathrm{0}}^{{\cos }^{-\mathrm{1}}z/\xi }\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}c\xi s_c^{\mathrm{2}}\xi \right)\mathrm{d}\phi \mathrm{d}\xi }}\right\}\mathrm{d}z={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}}c{s}_c^{\mathrm{2}}\left({\displaystyle \frac{h{\delta }^{\mathrm{4}}}{\mathrm{2}}}\right),

in which h represents the thickness of the material. Hence, the critical stretch can be expressed as

(23)s_c=\sqrt{{\displaystyle \frac{G_c}{\left({\displaystyle \frac{\mathrm{6}}{\mathrm{\pi }}}\mu +{\displaystyle \frac{\mathrm{16}}{\mathrm{9}{\mathrm{\pi }}^{\mathrm{2}}}}(\mathrm{\kappa }-\mathrm{2}\mu )\right)\delta }}}.

In order to include damage initiation in the material response, a history-dependent scalar-valued function \muis defined as

(24)

Note that when the stretch between these material points exceeds its critical stretch, failure occurs and \mu =\mathrm{1}, otherwise \mu =\mathrm{0}. Then the force density vector can be modified through \mu

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

(25)\cdot {\displaystyle \frac{{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}}{\left|{\mathit{y}}_{(j)}-{\mathit{y}}_{(k)}\right|}}.

Local damage at a point is defined through the weight ratio function. The local damage at a point can be quantified as

(26)\varphi ({\mathit{x}}_{(k)},t)=\mathrm{1}-{\displaystyle \frac{{\displaystyle {\int }_H\mu ({\mathit{x}}_{(j)}-{\mathit{x}}_{(k)},t)\mathrm{d}{V}_{(j)}}}{{\displaystyle {\int }_H\mathrm{d}{V}_{(j)}}}}.

The local damage ranges from 0 to 1. When the damage is zero, it means that all interactions are intact and when the local damage is one, all the interactions initially associated with the point have been eliminated. The crack has then be completed formed and propagated. Then, half of the interaction within its horizon, resulting in a local damage value of one half.

Numerical example

Double-notched specimen made of Q345 steel under uniaxial tension

Q345 material is a low alloy high strength structural steel, with elasticity modulus E=203 GPa, Poisson’s ratio v=0.3, yield strength f=410 MPa, ultimate strength {\sigma }_s=572 MPa, elongation \delta=27.96% and density \rho=7850 kg/m³. The length and width of the specimen are 70 mm and 40 mm, respectively as shown in Fig. 3. The specimen is loaded under uniaxial tension with a constant loading rate of 2.217×10^–5 m/s. The crack length are 10 mm and three different crack size distances in loading directions are tested according to Table 1. The initial and fracture specimen from the experiments are illustrated in Fig. 4 and Fig. 5.

The C3D8R XFEM element in ABAQUS is employed to simulate the crack propagation of the three specimens. The XFEM region of the specimen is divided into two parts. Each of them contains only one crack. Three kinds of specimen crack growth results are given in Fig. 6. The calculated results by XFEM are in good agreement with the experimental results in Fig. 5.

The Q345 steel specimen is simplified model in 2D in PD. A particle distance of \Delta x=0.5 mm is employed in PD and the load is applied displacement-controlled over a width d=\mathrm{3}\Delta x. To better visualize the crack, only the area near the crack-- i.e., 40 mm in longitudinal direction-- is shown.

The fracture patterns obtained by the PD simulations for the 3 specimen are shown in Fig. 7. When the crack tip damage value is close to 1, the PD force of two particles disappears, and the crack propagates forward. With increasing loading, the damage site expands and eventually runs through the entire specimen as shown in Fig. 8.

For a vertical crack tip distance of 0 mm, two cracks propagate expectedly along in the horizontal direction and eventually join into one single crack as depicted in Fig. 5(a) and Fig. 8 (a). When the vertical crack tip distance is increased to 10 mm, the two cracks propagate initially in the horizontal direction as well extended along the horizontal direction before they turn in a 45° angle versus the horizontal axis and finally join as shown in Fig. 5(b) and Fig. 8(b). Some artificial damage zone can be observed in the PD simulation which do not occur in XFEM but the overall agreement is good. The specimen with a crack tip distance of 20 mm has a quite distinguished fracture pattern. Instead of crack coalescence, both cracks propagate through the entire specimen. No crack shielding occurs which is captured by both XFEM and PD as illustrated in Fig. 5(c) and Fig. 8(c).

Double center crack problem

Let us consider a specimen with two double-center cracks under uniaxial tension with a constant loading rate of 20 m/s. The material parameters from the previous example are adopted. The thickness of the plate is 0.05 m, the length of one crack is 0.01 m. The longitudinal spacing of the crack is varied from 5 mm (case 1) to 10 mm (case 2). 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 was analyzed.

Figures 9 and 10 depict the displacements y-direction at different load steps. Note that the maximum displacement for case 2 is less than for case 1 after the crack has grown through the entire specimen. In both cases the maximum displacement position change from the ends of the plate to near the crack tip.

Figures 11 and 12 show that the crack distance has a slight effect on the crack propagation path. The greater the crack distance, the earlier the crack extension begins. At 27.2 ms, the crack length increases with respect to the increase of the crack length of 10 mm, and the crack propagation direction changes as the crack tip expands, while the crack with the spacing of 5 mm spreads along the vertical load direction. At 34 ms, cracks with 5 mm distance propagation are changed to the opposite direction and the crack with 10 mm distance propagation direction is perpendicular to the load direction. At 40.8 ms, the crack propagates through the entire length of the plate and the spacing between the crack tips in both plates is much larger than the initial crack spacing.

Figure 13 shows the damage value of the material in the plates. The abscissa represents the ordinal value of the particle and the ordinate represents the damage value. The larger the damage value, the faster the crack will propagate. So, the crack propagation rate of big crack spacing is much higher than that of the small crack spacing, and this phenomenon would become weaker with the crack propagation. The damage rate of the two cracks with the crack spacing of 10 mm at 27.2 ms is significantly higher than that at the crack spacing of 5 mm, the gap between the two plates decreases gradually at 34 ms, and the maximum of the two plates is nearly the same at 40.8 ms.

Conclusions

In this paper, the influence of bilateral crack and central longitudinal crack on the crack propagation path is simulated by PD theory. At the same time, the experiment and XFEM are used to validate the simulation of the bilateral crack propagation. The accuracy of the crack propagation law of PD theory and the advantage the XFEM method are obtained. A comparison between the results of XFEM and PD demonstrates that PD is a good competitor to XFEM.

For the center double crack spacing, the bigger the spacing, the earlier the crack initiation time and the faster the initial expansion rate. As the crack propagates, this phenomenon is progressively weakened. The results show that PD theory can accurately simulate the relationship between crack propagation and crack propagation with time.

References

[1] 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: 322–350

[2] 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

[3] 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

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

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

[6] 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

[7] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering, 1999, 45(5): 601–620

[8] Daux C, Moës N, Dolbow J, Sukumar N, Belytschko T. Arbitrary branched and intersecting cracks with the extended finite element method. International Journal for Numerical Methods in Engineering, 2000, 48(12): 1741–1760

[9] Nagashima T, Omoto Y, Tani S. Stress intensity factor analysis of interface cracks using X-FEM. International Journal for Numerical Methods in Engineering, 2003, 56(8): 1151–1173

[10] Fang X, Jin F, Wang J. Cohesive crack model based on extended finite element method. Journal of Tsinghua University, 2007, 47(3): 344–347

[11] Fang X, Jin F, Wang J. Seismic fracture simulation of the Koyna gravity dam using an extended finite element method. Journal of Tsinghua University, 2008, 48(12): 2065–2069

[12] Song J H, Areias P M A, Belytschko T. A method for dynamic crack and shear band propagation with phantom nodes. International Journal for Numerical Methods in Engineering, 2006, 67(6): 868–893

[13] Hansbo A, Hansbo P. A finite element method for the simulation of strong and weak discontinuities in solid mechanics. Computer Methods in Applied Mechanics and Engineering, 2004, 193(33–35): 3523–3540

[14] Areias P M A, Belytschko T.A comment on the article “A finite element method for simulation of strong and weak discontinuities in solid mechanics” by A. Hansbo and P. Hansbo [Comput. Methods Appl. Mech. Engrg. 193 (2004) 3523–3540]. Computer Methods in Applied Mechanics & Engineering, 2006, 195: 1275–1276

[15] Rabczuk T, Zi G, Gerstenberger A, Wall W A. A new crack tip element for the phantom node method with arbitrary cohesive cracks. International Journal for Numerical Methods in Engineering, 2008, 75(5): 577–599

[16] Chau-Dinh T, Zi G, Lee P S, Rabczuk T, Song J H. Phantom-node method for shell models with arbitrary cracks. Computers & Structures, 2012, 92: 242–256

[17] 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

[18] Rabczuk T, Areias P M A. A meshfree thin shell for arbitrary evolving cracks based on an extrinsic basis. Cmes Computer Modeling in Engineering & Ences, 2006, 16(2): 115–130

[19] 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

[20] 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

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

[22] 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

[23] 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

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

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

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

[27] Zhuang X, Zhu H, Augarde C. An improved meshless Shepard and least square method possessing the delta property and requiring no singular weight function. Computational Mechanics, 2014, 53(2): 343–357

[28] Zhuang X, Cai Y, Augarde C. A meshless sub-region radial point interpolation method for accurate calculation of crack tip fields. Theoretical and Applied Fracture Mechanics, 2014, 69: 118–125

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

[30] 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

[31] 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

[32] Ganzenmüller G C, Hiermaier S, May M. On the similarity of meshless discretizations of peridynamics and Smoothed Particle Hydrodynamics. Computers & Structures, 2015, 150: 71–78

[33] Bessa M A, Foster J T, Belytschko T, Liu W K. A meshfree unification: reproducing kernel peridynamics. Computational Mechanics, 2014, 53(6): 1251–1264

[34] 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

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

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

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

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

[39] 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

[40] Talebi H, Silani M, Rabczuk T. Concurrent multiscale modelling of three dimensional crack and dislocation propagation. Advances in Engineering Software, 2015, 80: 82–92

[41] Silani M, Talebi H, Ziaei-Rad S, Hamouda A M, Zi G, Rabczuk T. A three dimensional extended arlequin method for dynamic fracture. Computational Materials Science, 2015, 96: 425–431

[42] 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

[43] Talebi H, Silani M, Bordas S P A, Kerfriden P, Rabczuk T. A computational library for multiscale modelling of material failure. Computational Mechanics, 2014, 53(5): 1047–1071

[44] 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

[45] Yang S W, Budarapu P R, Mahapatra D R, Bordas S P A, Zi G, Rabczuk T. A meshless adaptive multiscale method for fracture. Computational Materials Science, 2015, 96: 382–395

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

[47] 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

[48] Costa T, Bond S, Littlewood D, Moore S. Peridynamic Multiscale Finite Element Methods. Sandia National Laboratories (SNL-NM), Albuquerque, NM (United States), 2015, SAND2015–10472

[49] Jung J, Seok J. Fatigue crack growth analysis in layered heterogeneous material systems using peridynamic approach. Composite Structures, 2016, 152: 403–407

[50] Shen F, Zhang Q, Huang D. Damage and failure process of concrete structure under uni-axial tension based on peridynamics modeling. Chin J Comput Mech, 2013, 30: 79–83

[51] Gu X, Zhou X P. The numerical simulation of tensile plate with circular hole using peridynamic theory. Chin J Comput Mech, 2015, 36(5): 376–383

[52] Cheng Z, Liu J. Fracture analysis of functionally graded materials under impact loading based on peridynamics. Chinese Journal of Applied Mechanics, 2016, (4): 634–639

[53] Zhou X, Shou Y. Numerical simulation of failure of rock-like material subjected to compressive loads using improved peridynamic method. International Journal of Geomechanics, 2017, 17(3): 04016086

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

[55] Zeng S. Bilateral crack growth behavior of Q345 steel under uniaxial tensile load. Dissertation for master degree. Guangxi University, 2014

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg