1. Introduction

Small cracks can have significant effect on the dynamic propagation behaviour of macrocracks if they are within the range of the macrocrack's influence domain. The location of small cracks and their orientation with respect to the tip of a macrocrack are important, since the stress intensity of the macrocrack tip changes accordingly. Numerical tools can be useful to investigate such an interesting and important problem.

Numerical prediction of crack growth in computational mechanics has been and it still is considerable problem that can't be easily solved by conventional numerical methods such as cohesive elements ^[1,2] and extended finite element method (XFEM) ^[3,4]. XFEM require damage criterion and careful stress tracking around the crack tip in order to decide if crack is going to branch or not. On the other hand, mesh dependency is a problem for cohesive elements. All of the above problems are making it difficult to correctly simulate crack propagation and especially multiple crack propagation/interaction ^[5]. Meshfree methods can be a good alternative to finite element method and has been applied to dynamic fracture and fragmentation ^[6,7]. Moreover, Cracking Particle Method (CPM) ^[8] was introduced for complex fracture patterns such as crack branching and coalescence. Another promising approach for fracture modelling is lattice methods which represent a medium as a connection of interacting nodes or particles ^{[9,10,11,12,13]}.

In this study, it is shown that using a new continuum mechanics formulation, peridynamics, as an alternative method, it is possible to correctly model and simulate dynamic fracture, in particular multiple small cracks interacting with a macrocrack in brittle materials. Peridynamics ^[14,15,16] doesn't need fracture criteria for crack propagation as in XFEM or cohesive FEM methods. Moreover, multiple cracks can easily be analysed since cracks are not treated as special objects in the formulation. Peridynamics has been successfully used to analyze different material systems and geometrical configurations ^{[17,18,19,20,21]}. An extensive literature survey on peridynamics is given in Madenci and Oterkus ^[22].

---

2. Peridynamic Formulation

By considering all of the shortcomings of classical theories such as Finite Element Analysis (FEM), the peridynamic theory was introduced by Silling ^[14] as a new continuum mechanics formulation in order to solve problems with discontinuities such as cracks. Peridynamic theory uses displacements rather than displacement derivatives in its formulation. Hence, as opposed to classical continuum mechanics formulation which uses partial differential equations, peridynamic theory is based on integral equations which are defined at fracture surfaces. Furthermore, material damage is part of the peridynamic constitutive laws ^[22].

In peridynamic theory, material points inside the body interact with each other through response functions. Material points (there is infinitely many) constitute the continuum and each response function contains all the information associated with the material ^[22].

Moreover, peridynamic theory can be considered as continuum version of molecular dynamics, where acceleration of any material point at position x in reference state at time t can be found from the equation of motion as follows:

$\rho \mathbf{\ddot{u}}\left( {{\mathbf{x}}_{k}},t \right)=\int\limits_{{{H}_{x}}}{\mathbf{f}\left( \mathbf{u}\left( {{\mathbf{x}}_{j}},t \right)-\mathbf{u}\left( {{\mathbf{x}}_{k}},t \right),{{\mathbf{x}}_{j}}-{{\mathbf{x}}_{k}} \right)}\,d{{V}_{x}}+\mathbf{b}\left( {{\mathbf{x}}_{k}},t \right)$

(1)

where H_x represents horizon, x_j are family members of the material point x_k inside its horizon, u is the displacement vector, b is the body load and ρ is the mass density. In Eq. (1), ρ is a pairwise force function which is the force that the material point x_j exerts on the material point x_k. The relative position of the material points x_j and x_k is denoted as (see Figure 1):

$\xi ={{\mathbf{x}}_{j}}-{{\mathbf{x}}_{k}}$

(2)

and their relative displacement is defined as:

$\eta =\mathbf{u}\left( {{\mathbf{x}}_{j}},t \right)-\mathbf{u}\left( {{\mathbf{x}}_{k}},t \right)$

(3)

In the original peridynamic formulation, i.e., bond based peridynamics, we can define pairwise force function as:

$\mathbf{f}\left( \eta ,\xi \right)=\frac{\xi +\eta }{\left| \xi +\eta \right|}f\left( \left| \xi +\eta \right|,\xi \right)$

(4)

where $f\left( \left| \xi +\eta \right|,\xi \right)$ is a scalar-valued function and it depends on the bond stretch, s and the bond constant, c which can be written as ^[15]

$f\left( \left| \xi +\eta \right|,\xi \right)=c\,s$

(5)

In Eq. (5), the stretch can be expressed as

$s=\frac{\left| \xi +\eta \right|-\left| \xi \right|}{\left| \xi \right|}$

(6)

The bond constant c, can be specified in terms of elastic modulus, E, thickness, δ and horizon size as

$c=\frac{9E}{\pi h{{\delta }^{3}}}\text{for}\,2D$

(7)

Eq. (7) can be obtained by comparing strain energy density of a material point inside a body subjected to a simple loading condition obtained from classical continuum mechanics and peridynamics. Derivation procedure is given in Madenci and Oterkus ^[22].

Note that stretch is positive when bond (interaction) is in tension and negative when bond is in compression. Combining equations given above and using numerical discretization based on meshless scheme, we can perform the time integration of the peridynamic equation of motion by using explicit time integration ^[15]. If the displacement of the material point at ${{n}^{th}}$ time step is represented as $\mathbf{u}_{k}^{n}=\mathbf{u}_{k}^{{}}\left( t=n\,\Delta t \right)$, then the peridynamic form of the equation of motion can be written as

${{\rho }_{k}}\mathbf{\ddot{u}}_{k}^{n}=\sum\limits_{j=1}^{N}{\frac{\xi +{{\eta }^{n}}}{\left| \xi +{{\eta }^{n}} \right|}}\cdot c\cdot s\cdot {{V}_{j}}\cdot \nu {{c}_{j}}\cdot s{{c}_{j}}+\mathbf{b}_{k}^{n}$

(8)

or

${{\rho }_{k}}\mathbf{\ddot{u}}_{k}^{n}=\sum\limits_{j=1}^{N}{\frac{\xi +{{\eta }^{n}}}{\left| \xi +{{\eta }^{n}} \right|}}\cdot c\cdot \frac{\left| \xi +{{\eta }^{n}} \right|-\left| \xi \right|}{\left| \xi \right|}\cdot {{V}_{j}}\cdot \nu {{c}_{j}}\cdot s{{c}_{j}}+\mathbf{b}_{k}^{n}$

(9)

where N is the number of the family members inside the horizon of the material point x_k, V_j is the volume of the material point x_j, vc_j and sc_j represent volume correction factor and surface correction factor, respectively ^[22].

Introducing failure in peridynamics is straightforward. Failure can be introduced by eliminating the interactions between the material points. A history dependent scalar valued function can be defined for each interaction as

$f\left( \left| \xi +\eta \right|,\xi \right)=c\,s\mu \left( t,\xi \right)$

(10)

where

$\mu \left( t,\xi \right)=\left\{ \begin{align} & 1\,\text{if}\,s\left( {t}',\xi \right)＜{{s}_{0}} \\ & 0\,\,\,\,\,\,\,\,\text{otherwise} \\ \end{align} \right.$

(11)

and s₀ is called critical stretch. The local damage at a material point can be defined as ^[15]

$\varphi \left( \mathbf{x},t \right)=1-\frac{\int\limits_{{{H}_{x}}}{\mu \left( \mathbf{x},t,\xi \right)d{{V}_{\xi }}}}{\int\limits_{{{H}_{x}}}{d{{V}_{\xi }}}}$

(12)

---

3. Validation Study

In order to validate the peridynamic formulation, a Diagonally Loaded Square Plate (DLSP) problem is considered as shown in Figure 2.

The material is chosen as Poly(methyl methacrylate) (PMMA), also known as acrylic or acrylic glass. The material properties are summarized in Table 1. The edge length of the plate is specified as $2w=0.15m$. The length of the pre-crack length is $2a=0.045m$ and the pre-crack inclination angle is defined by angle, α. Different values of α indicates different types of loading conditions ranging from mode-I to mode-Ⅱ. When α is zero degree, then the initial crack becomes perpendicular to the loading direction which induces crack surfaces to open in loading direction and represents mode-I loading. Mode-Ⅱ type of loading, which forces crack surfaces to slide against each other, can be obtained when α is equal to 62.50° and if crack-edge length ratio ($a/w$) is equal to 0.3. For the other α values between mode-I and mode-Ⅱ angles, mixed-mode loading conditions can be induced. In this study, the following α values are considered: α=0°, 15°, 30°, 45°, 62.5°. For the spatial discretization, meshless scheme is used and the distance between material points is specified as $\Delta x={{10}^{-3}}m$. Moreover, the horizon size is chosen as $\delta =3.015\Delta x$. The plate is loaded at the upper and lower edges through a region with a length of $b=0.025m$. A velocity boundary constraint of $V=1\times {{10}^{-8}}m/s$ is considered to achieve a quasi-static loading condition. To achieve such a condition in peridynamics, the problem is solved by using Adaptive Dynamic Relaxation method ^[23]. Hence, a time step size of $\Delta t=1\sec .$ is used. The total number of time steps is 40000. In order to prevent unrealistic failure around the loading area, a no-fail zone is introduced with a length of $n=0.0375m$. The critical stretch is specified as ${{s}_{c}}=9.735\times {{10}^{-3}}$.

As can be seen in Figure 3, peridynamic fracture pattern results agree very well with experimental results obtained by Ayatollahi and Aliha ^[24] for different pre-crack inclination angles, α.

Moreover, for all cases, load variation (P) with respect to displacement (δ) results are plotted as shown in Figure 4. Since PMMA material has a linear elastic behavior, P-δ curves follow a straight line prior to the fracture load. After crack starts to propagate, the curves smoothly drop down and suddenly reduce to zero value when the plate fractures completely.

Evaluated fracture loads from Figure 4, are also compared with experimental results (see Table 1 of ^[24]) as shown in Figure 5 and comparable results are obtained for all crack inclination angles, α.

Moreover, crack (propagation) initiation angles are also compared with experimental results (see Figure 11 of ^[24]) as demonstrated in Figure 6 and again a good agreement is obtained for all crack inclination angles, α.

Hence, it can be concluded that peridynamics can make accurate fracture predictions and can be used for more sophisticated problem cases as considered in Section 5.

---

4. Problem Definition

A rectangular plate with dimensions of 0.05m by 0.05m is considered as shown in Figure 7. As in the validation case, the material is chosen as PMMA. Plate contains a macrocrack that interacts with multiple small cracks. In the case of bond based peridynamics, there is a constraint on Poisson's ratio as 1/3 for 2-Dimensional problems due to the assumption of pairwise forces between material points which is slightly lower than the actual Poisson's ratio of PMMA. However, for dynamic fracture problems, Poisson's ratio has insignificant influence on speed and direction of crack propagation ^[5].

In this study, the plate is considered under tension loading and it is subjected to a high velocity boundary condition of 5 m/s. Several types of small crack configurations are considered; single small crack collinear to the main crack, two symmetrical small cracks, horizontal and transverse array of small cracks, to investigate the effect of small cracks on the macrocrack behaviour by calculating the crack propagation speed.

---

5. Numerical Results

Peridynamic model used for this study is defined with a fixed horizon size of δ = 3 × Δx, where Δx is the spacing between material points and it is specified as 0.0001m. Therefore, a total number of 250000 material points exist in the model. Boundary region is equal to horizon size, δ and thickness of the plate is specified as h = Δx (Figure 8). Critical stretch is set to s₀ = 0.030857. The time step is chosen as Δt = 4 × 10⁻⁸ s and the number of time steps is 2000.

Finally, the length of a small crack is defined as

${{l}_{small\,crack}}＜＜{{l}_{macrocrack}}\to {{l}_{small\,crack}}={{l}_{macrocrack}}/50$

(13)

---

5.1. Macrocrack Propagation

In the first example, the crack propagation of a macrocrack is investigated without considering small cracks in the model in order to compare with those cases including small cracks. It is observed that the macrocrack started to propagate around a time step of 500 and reached the end of the plate after 1700 time steps. The crack speed is calculated by comparing the crack length at 4 × 10⁻⁵ s with the initial crack length (Figure 9a) as

${{v}_{cp}}=\frac{\Delta l}{dt}=480\,\text{m/s}$

(14)

---

5.2. One Small Crack Collinear to a Macrocrack

In the second example, a single small crack is aligned with the macrocrack as shown in Figure 10 ^[25,26].

The macrocrack propagation speed is calculated for different ratios of a/b as given in Table 2. The minimum crack propogation speed is obtained for the highest a/b ratio which corresponds to furthest small crack. By comparing with the only macrocrack case, the collinear small crack causes an increase in macrocrack propagation speed. However, no effect on initiation time is observed, i.e., all the cases have the same initiation time as benchmark case—500 time steps. The effect of a/b ratios on the crack shapes can be seen in Figure 11.

---

5.3. Two Parallel Small Cracks

In the third example, two symmetric small cracks are positioned on both sides of the macrocrack (Figure 12).

The macrocrack propagation speed is calculated for varying h/l and s/l values. It is found that for s/l = −2, the influence of the small cracks on the speed of macrocrack propagation is insignificant and the shape of the crack path is very similar to that of a single macrocrack (Figure 13).

Moreover, h/l ratio has the least influence on the overall macrocrack propagation as shown in Table 3. The minimum crack speed is obtained for h/l = 1.25 and s/l = 0 (Figure 14). It is worth noting that for s/l = 0 and s/l = 2, the crack propagation is observed in small cracks only. This behaviour is also observed by Wang at al. ^[25]. They observed that for h/l = 0.75 and 1, the influence of the stress intensity factor of small cracks is greater than the macrocrack. A similar behaviour is also observed for h/l = 1.25.

---

5.4. Multiple Small Cracks Interacting with the Macrocrack

For the last example, interactions between macrocrack and multiple small cracks are investigated. First, a set of horizontal cracks (see Figure 15) for h/l = 1 and H/l = 2.5 are investigated. It is observed that H/l should be greater than 2 in order to exclude the effect between adjacent columns of small cracks ^[25].

The crack propagation speed is calculated for different ratios of s/l as given in Table 4. It is found that the present small crack configuration stopped the macrocrack propagation. However, the present configuration has no influence on crack initiation time and a very low influence on the overall crack propagation length (Table 4 and Figure 16).

Moreover, the effect of multiple small cracks in the vertical direction is investigated (Figure 17) for h/l = 1 and s/l = −2, 0 and 2.

The crack propagation speed is calculated for different ratios of s/l as given in Table 5. It is found that the crack propagation speed is reduced by 34% for s/l = 0 compared to a single macrocrack case (Table 5). The results show good correlation with conclusions from Wang at al. ^[25] where they observed relatively high amplification decrease for h/l = 1 (Figure 18).

Finally, the superior capability of peridynamics for fracture prediction is demonstrated by considering a significantly large number of small cracks. In the first case, thirty-two small cracks are considered as shown in Figure 19. The second case includes eighty small cracks as depicted in Figure 20.

Corresponding fracture patterns for these two cases are presented in Figures 21 and 22, respectively. In both cases, it can be seen that small cracks behind the macrocrack tip does not effect crack propagation. Furthermore, only small cracks ahead of the macrocrack tip close to the centerline influences the crack propagation pattern and causes crack branching behavior.

---

6. Conclusion

In this study, quantitative comparison of effects of small cracks on dynamic macrocrack propagation is obtained by using periydnamics. Several small crack configurations are considered including single small crack collinear to the main crack, two symmetrical small cracks, and horizontal and transverse array of small cracks. The results show expected behaviour when compared to similar cases found in the literature. Finally, two cases including significantly large number of small cracks are demonstrated. These cases show the superior capability of peridynamics in capturing sophisticated fracture patterns where significantly large number of cracks present. Moreover, presented results also show that this analysis can be useful for the design of new microstructurally toughened materials.

---

Conflict of Interest

The authors declare that there is no conflict of interest regarding the publication of this manuscript.

---