# Studies of dynamic crack propagation and crack branching

## Transcript Of Studies of dynamic crack propagation and crack branching

University of Nebraska - Lincoln

[email protected] of Nebraska - Lincoln

Faculty Publications from the Department of Engineering Mechanics

Mechanical & Materials Engineering, Department of

2010

Studies of dynamic crack propagation and crack branching with peridynamics

Youn Doh Ha Ph.D. University of Nebraska at Lincoln

Florin Bobaru Ph.D. University of Nebraska at Lincoln, [email protected]

Follow this and additional works at: https://digitalcommons.unl.edu/engineeringmechanicsfacpub Part of the Applied Mechanics Commons, Ceramic Materials Commons, Computational Engineering

Commons, Engineering Mechanics Commons, Mechanics of Materials Commons, Structural Engineering Commons, Structural Materials Commons, and the Structures and Materials Commons

Ha, Youn Doh Ph.D. and Bobaru, Florin Ph.D., "Studies of dynamic crack propagation and crack branching with peridynamics" (2010). Faculty Publications from the Department of Engineering Mechanics. 71. https://digitalcommons.unl.edu/engineeringmechanicsfacpub/71

This Article is brought to you for free and open access by the Mechanical & Materials Engineering, Department of at [email protected] of Nebraska - Lincoln. It has been accepted for inclusion in Faculty Publications from the Department of Engineering Mechanics by an authorized administrator of [email protected] of Nebraska - Lincoln.

Published in International Journal of Fracture 162 (2010), pp. 229–244; doi: 10.1007/s10704-010-9442-4 Copyright © 2010 Springer Science+Business Media B.V. Used by permission. Submitted August 31, 2009; accepted January 5, 2010; published online January 20, 2010.

Studies of dynamic crack propagation and crack branching with peridynamics

Youn Doh Ha and Florin Bobaru

Department of Engineering Mechanics, University of Nebraska–Lincoln, Lincoln, NE 68588-0526, USA Corresponding author — F. Bobaru, email [email protected]

Abstract In this paper we discuss the peridynamic analysis of dynamic crack branching in brittle materials and show results of convergence studies under uniform grid refinement (m-convergence) and under decreasing the peridynamic horizon (δ-convergence). Comparisons with experimentally obtained values are made for the crack-tip propagation speed with three different peridynamic horizons. We also analyze the influence of the particular shape of the micro-modulus function and of different materials (Duran 50 glass and soda-lime glass) on the crack propagation behavior. We show that the peridynamic solution for this problem captures all the main features, observed experimentally, of dynamic crack propagation and branching, as well as it obtains crack propagation speeds that compare well, qualitatively and quantitatively, with experimental results published in the literature. The branching patterns also correlate remarkably well with tests published in the literature that show several branching levels at higher stress levels reached when the initial notch starts propagating. We notice the strong influence reflecting stress waves from the boundaries have on the shape and structure of the crack paths in dynamic fracture. All these computational solutions are obtained by using the minimum amount of input information: density, elastic stiffness, and constant fracture energy. No special criteria for crack propagation, crack curving, or crack branching are used: dynamic crack propagation is ob-

tained here as part of the solution. We conclude that peridynamics is a reliable formulation for modeling dynamic crack propagation.

Keywords: dynamic fracture, crack branching, brittle fracture, peridynamics, nonlocal methods, meshfree methods

1 Introduction

1.1 Literature review of dynamic crack propagation

In a brittle material, a propagating crack can depart from its original straight trajectory and curve or split into two or more branches. Under very high states of stress, the propagating crack will divide into a river-delta crack pattern (Bowden et al. 1967; Ramulu and Kobayashi 1985). This fragmentation of highly loaded, brittle materials is often a succession of multiple branching of what was initially a single crack. Increases in the roughness of the fracture surface prior to branching were consistently observed in all reported investigations (Ramulu and Kobayashi 1985; Döll 1975). In crack branching of edge notch specimens of brittle materials it has also been observed that the crack tip velocity drops by no more than 5–10% in the branching region (Döll 1975).

In atomistic models, under conditions that lead to instability of the crack path, cracks can branch without a specific criterion (see Zhou et al. 1996).

229

230

Y. D. Ha & F. Bobaru in International Journal of Fracture 162 (2010)

Particle-type models (see Bolander and Saito 1998) are also capable of simulating crack branching. However, none of these methods is able to capture the crack propagation speed or the angle of crack branching correctly. For instance, MD simulations show instabilities that lead, shortly after the bifurcation of a crack, to the propagation of only one of the two branches, the other being arrested. Moreover, the branching angle computed with MD (see Zhou et al. 1996) is greater than 90°, whereas experiments show much smaller crack branching angles (Ramulu and Kobayashi 1985). One may ask whether quantum mechanical calculations are needed to predict the phenomenon of dynamic fracture in brittle materials (see Cox et al. 2005) which is one of the great challenges in dynamic fracture. One likely reason for MD simulations’ failure to correctly predict dynamic fracture is that, for example, crack branching events are controlled by the interaction and wave reflections from the boundaries (Ravi-Chandar 1998). Because of this, one would have to either model the entire structure with MD (not a viable option) or use a multiscale model that is capable of transferring the waves between the scales correctly (still an open problem). Numerical simulations based on continuum methods of dynamic crack propagation behavior have been very difficult to develop and, to this date, a reliable method for simulating this complex problem has not been found in spite of considerable efforts in this direction (e.g. Xu and Needleman 1994; Camacho and Ortiz 1996; Ortiz and Pandolfi 1999; Belytschko et al. 2003; Rabczuk and Belytschko 2004; Song et al. 2006). All these methods use some version of cohesive-zone models. As such, they all modify the local continuum mechanics equations and introduce a nonlocal effect given by the parameters and length scales in the cohesive-zone model. To reduce mesh dependency when the grids are refined special methodologies have to be used (Zhou and Molinari 2004). For the existing approaches, the difficulties in modeling dynamic fracture processes like crack branching are many. For example, continuum-type methods using the cohesive FEM or the XFEM require a damage criterion and a tracking of the stresses around the crack tip to decide when to branch the crack. Decisions also have to be taken in terms of the angle of propagation of the branches and about how many branches will be allowed to form. In methods in which the crack advances along the element sides by separating elements from one another,

the crack path becomes non-smooth (see Xu and Needleman 1994; Camacho and Ortiz 1996; Ortiz and Pandolfi 1999). Since the correct path (which minimizes the strain energy) of the crack propagation is not computed correctly, there are significant departures from the true energy released during the crack propagation event. In such cases, reliable prediction of strength of brittle ceramics under impact, for example, becomes difficult. Mesh dependency is an additional problem in cohesive-zone FEM-based methods. Important progress has been recently made by using the XFEM method which allows cracks to pass through the finite elements (see e.g. Belytschko et al. 2003). Subdivision of the cut elements for numerical integration purposes increases the complexity and the cost of the method. This method requires phenomenological damage models and branching criteria, as well as tracking of the crack path using level sets, for example. It is not yet clear if the method is applicable to problems that involve fragmentation and/or multiple crack interactions, branching, and coalescence. The method does not predict the experimentally observed crack propagation speeds (see Song et al. 2008). Cohesivezone based models need to modify the experimental values of the fracture energy by several factors in order to get propagation velocities in the range of measured ones (Song et al. 2008).

In the present contribution we try to answer whether quantum, atomistic, or multiscale models are needed in dynamic fracture in order to correctly simulate the observed crack propagation velocities and crack paths (Cox et al. 2005; Song et al. 2008). We will show that peridynamics is able to correctly model and simulate dynamic fracture, in particular crack branching in brittle materials. Peridynamics, which is a reformulation of continuum mechanics (Silling 2000; Silling et al. 2007b), does not require criteria for crack propagation or crack branching: these happen spontaneously in this method and are autonomously generated by a simple bond-failure criterion that is correlated to the material’s energy release rate. The name “peridynamics” comes from the Greek “peri” which means “nearby,” and dynamics. Peridynamics is a nonlocal method in which material points interact not only with their nearest neighbors but also with points nearby, inside a horizon. This is what physically happens at the atomic scale, for example, but peridynamics extends this idea to the continuum scale. We will show convergence in terms of the crack path and the crack propagation speed under

S t u d i e s o f dy n a m i c c r a c k p r o p a g a t i o n a n d c r a c k b r a n c h i n g w i t h p e r i dy n a m i c s

231

grid refinement (the so-called δ-convergence, Bobaru et al. 2009) and under decreasing peridynamic horizon (the so-called δ-convergence, Bobaru et al. 2009). The crack branching patterns obtained using peridynamics follow remarkably close the experimental results which show secondary branching taking place when higher stress levels are reached at the tip of the pre-notch prior to crack propagation. Moreover, the only input parameters in the model are the Young’s modulus, the density, and the fracture energy (which is kept constant, and not a function of the propagation velocity or of the incurred damage, in this first study).

The paper is organized as follows: in the next section we describe the sample problem setup. In Section 2 we briefly review the peridynamic formulation and the connections between the parameters in the formulation and the material properties like the energy release rate. In Section 3 we present the numerical results for the convergence studies. We look at both the crack path and the propagation speed of the crack, as measures of convergence. In Section 4 we analyze the influence of the micromodulus function on crack branching results as well as the solutions for two different brittle materials under higher loading conditions that lead to cascading branching. We also comment on the roughening zones that take place in the branching regions and on the effect of the reflection waves on the propagation paths of the dynamic cracks. The conclusions are given in Section 5.

1.2 Problem setup

We consider the following setup as a benchmark problem for analyzing crack branching phenomena: a prenotched thin rectangular plate with 0.1m by 0.04m as shown in Figure 1. All simulations in this paper are 2D simulations. For some 3D results we refer the readers to Ha et al. (2010). The materials chosen for this study are selected because for these materials there are experimental results published on the crack propagation velocity in the region of branching or the maximum propagation velocity measured. The two materials used here are a Duran 50 glass (taken form Döll 1975) and a soda-lime glass (taken from Bowden et al. 1967). The material properties are summarized in Table 1. Please note that in the bond-based peridynamic implementation used here, the numerical models will be limited to using a fixed Poisson ratio of 1/3 (for 2D plane stress problems). If other Poisson ratios are desired, then the state-based

Figure 1. Description of the problem setup for the crack branching study.

peridynamics formulation should be used (see Silling et al. 2007b). For dynamic fracture problems, the Poisson ratio value does not have a significant influence on the propagation speed or the crack path shapes (Silling et al. 2007a).

In the experimental settings the loading of the sample may take tens of seconds or more. In explicit dynamic simulations that would be too expensive to compute. Instead, we choose to apply, along the upper and lower edges (see Figure 1), traction loadings σ suddenly at the initial time step and maintain this loading constant after that. The theoretical background for the peridynamics analysis is based on Silling’s original peridynamics paper (Silling 2000), the imposition of the traction boundary conditions is as in Ha and Bobaru (2009), and the numerical implementation of failure is like in Silling and Askari (2005). The same geometrical setup for studying crack branching simulations has been used in other studies (Belytschko et al. 2003; Rabczuk and Belytschko 2004; Song et al. 2006).

While there is no analytical solution for the crack branching problem, we can compare our simulation results with experiments. Unfortunately, the experimental papers we found do not provide a complete description of the conducted experiment on crack branching: some papers show the crack paths but do not provide crack propagation speed data, others give the propagation speed but do not show the crack paths, and most do not describe in detail the loading conditions. We decided to perform the peridynamic simulations using a setup similar to that used in a few recent simulation papers (Belytschko et al. 2003; Rabczuk and Belytschko 2004; Song et al. 2006). The material parameters, however, are like those used in the experiments (Bowden et al. 1967; Döll 1975). The maximum crack propagation speed, or the

232

Y. D. Ha & F. Bobaru in International Journal of Fracture 162 (2010)

Table 1. Material properties for Duran 50 and soda-lime glasses

Density (ρ) (kg/m3) Young’s modulus (E) (GPa)

Duran 50 glass

2,235

65

Soda-lime glass

2,440

72

Poisson ratio (υ) Fracture energy (G0) (J/m2)

0.2

204

0.22

135

crack propagation speed in the region of branching, is data that is fairly reproducible in experiments and this is reported in Bowden et al. (1967) and Döll (1975), for example. We are not aware of any numerical method that can reproduce the experimentally measured dynamic crack propagation velocity. Note that in previous studies, for certain methods, the fracture energy has to be significantly modified (by several factors) in order to bring the dynamic crack propagation speed closer to the measured values (see, e.g. Belytschko et al. 2003; Rabczuk and Belytschko 2004).

2 The peridynamic formulation

The peridynamic formulation (Silling 2000) relies on integration of forces acting on a material point and thus it does not face any of the mathematical inconsistencies seen in the classical continuum mechanics equations. The integration takes place over a “horizon” (which, in principle, extends to infinity but, for convenience is finite) within which the material points are interacting with each other. In certain problems, the size of the horizon can be correlated to an intrinsic material length-scale. However, in many cases a material length scale is not “visible” either because the micro-structure and the specific loading and boundary conditions do not lead to a measurable effect of the length-scale. In such cases, the horizon is selected by the user according to convenience (see Bobaru et al. 2009). Allowing a variable horizon (with a correspondingly scaled micromodulus parameter) defines away of introducing adaptive refinement for this nonlocal method. It is important to notice that peridynamics is a continuum theory, not a particle-type method. This allows the convergence results of the peridynamic solution to the classical elasticity solutions in the limit of the horizon going to zero (Bobaru et al. 2009; Silling and Lehoucq 2008).

An important advantage of peridynamics is the way damage is introduced: material points are connected within the horizon via elastic (linear or nonlinear) bonds that have a critical relative elongation, s0, at which they break (Silling 2000).

The critical relative elongation for brittle materials is computed from the experimentally measured value of the fracture energy for a specific material (Silling and Askari 2005). Damage is implemented as the fraction between the number of broken bonds and the number of initial bonds (Silling and Askari 2005). Cracks in peridynamics form as surfaces between material points form, as a consequence of sequential breaking of bonds. Thus, there is no need to track the cracks like in other continuum methods, or to impose criteria for when cracks should branch, change direction, turn, coalesce, etc. Moreover, peridynamics allows for spontaneous generation of cracks where no flaws were present before. This is shown, for example, in Silling et al. (2009) for the crack nucleation and in simulation of spallation (see Xie 2005) where spallation is treated as real fracture and not modeled by void-growth formulations as in existing literature results.

We now briefly review the peridynamic formulation based on Silling’s original peridynamics paper (Silling 2000). Also, we consider the summary of the numerical implementation of the traction boundary conditions in peridynamics (Ha and Bobaru 2009) and the formulation for the damage model in peridynamics (Silling and Askari 2005).

The peridynamic equations of motion are given by:

∫ ρü (x, t) = f(u (xˆ , t) − u(x, t), xˆ − x) dxˆ + b(x, t) (1) H

where f is the pairwise force function in the peridynamic bond that connects node xˆ to x and u is the displacement vector field. ρ is the density and b (x, t) is the body force. The integral is defined over a region H called the “horizon,” which is the compact supported domain of the pairwise force function around point x.

A micro-elastic material (Silling 2000) is defined as one for which the pairwise force derives from a potential ω:

∂ω (η, ξ )

f (η, ξ ) = ∂η

(2)

where ξ = xˆ − x is the relative position in the refer-

ence configuration and η = uˆ − u is the relative dis-

S t u d i e s o f dy n a m i c c r a c k p r o p a g a t i o n a n d c r a c k b r a n c h i n g w i t h p e r i dy n a m i c s

233

placement. A linear micro-elastic potential is obtained if we take

c (ξ ) s2ξ

ω (η, ξ ) = 2

(3)

where ξ = ξ and the relative elongation of a bond is

ζ − ξ

s= ξ

(4)

where ζ = ξ + η . The function c (ξ ) is called micro-modulus and has the meaning of the bond elastic stiffness. There are various formulations for the micromodulus function and in Section 4.1 we perform tests for dynamic crack propagation to assess the influence of the particular shape of the micro-modulus function on the crack path structure. We will observe that the crack propagation speed is not influenced by the shape of the micro-modulus, once the horizon is reasonably small compared to the dimensions of the structure analyzed. The pairwise force corresponding to a linear microelastic potential has the following form:

{ ξ+η

f (η, ξ ) = ξ + η c (ξ ) s,

ξ ≤ δ

0,

ξ>δ

(5)

where δ is the radius of the horizon region (which we will also refer to as the horizon). In this paper, we use the constant and conical 2D micro-modulus functions (see Figure 2). Following the same procedure performed to calculate the micro-modulus function in 1D (Bobaru et al. 2009), we obtain the constant micro-modulus function in 2D, plane stress conditions:

6E

c (ξ ) = c0 = πδ3 (1 − υ)

(6)

Similarly, the conical micro-modulus function is obtained as

c

(ξ ) = c1 ( 1 −

ξ δ

)

=

π

δ

3

24E (1 −

υ)

(

1

−

ξ δ

)

(7)

The shapes of the constant and conical micromodulus functions are illustrated in Figure 2.

In the bond-based peridynamics, any particle inside the horizon of another particle interacts only through a central potential. This assumption results (for an isotropic, linear, micro-elastic material) in an effective Poisson ratio of 1/3 in 2D (and 1/4 in 3D), but this limitation is readily eliminated by using the state-based peridynamics (Silling et al. 2007a,b). In this paper, we utilize the bond-

Figure 2. Constant (left) and conical (right) micromodulus functions.

based peridynamics, thus, in all the reported simulations here the effective Poisson ratio is 1/3.

In order to introduce failure into the peridynamic model, we consider that the peridynamic bonds can be broken when they are stretched beyond a predefined limit. We call this limit the “critical relative elongation, s0.” According to Silling and Askari (2005), there is no force sustained by the bond after its failure. Also, once a bond fails, it is failed forever; this makes the model history dependent. To completely separate a body into two halves across a fracture plane requires breaking all the bonds that initially connected points in the opposite halves (see Silling and Askari 2005). The energy per unit fracture length (in 2D, fracture area in 3D) for complete separation of the two halves of the body is called fracture energy, G0. In 3D, Silling and Askari relate the critical elongation, s0, with this measurable quantity (G0) (Silling and Askari 2005). In 2D with plane stress conditions, the fracture energy can be derived as

δ δ cos −1(z/ξ )

∫ ∫ ∫ G0 = 2

0 z

0

[c

(ξ

) 2

s02

ξ

]

ξ

dθ

dξ

dz

(8)

(See Figure 3 for an explanation of this computation.) Substituting the constant micro-modulus function (Equation 6) and rearranging for s0, we can rewrite this equation to obtain s0

s0

=

√

4 π G0 9Eδ

(9)

In similar way, the critical relative elongation for the conical micro-modulus function (Equation 7) is

√s0 = 5πG0 (10) 9Eδ

The critical relative elongation depends on the material properties and the horizon δ. Note that as the horizon goes to zero, the critical relative elongation goes to infinity, thus breaking such bonds

234

Y. D. Ha & F. Bobaru in International Journal of Fracture 162 (2010)

Figure 3. Evaluation of fracture energy. For each point A along the dashed line, 0 ≤ z ≤ δ, the work required to break the bonds connecting A to each point B in the circular cap is given by Equation (8).

requires larger and larger forces. This agrees with the physical experience at atomic and subatomic scales where, in order to break apart smaller and smaller sized bonds, one needs higher and higher forces. The values for the fracture energy used in this paper are the ones given in Döll (1975) for Duran 50 glass and soda-lime glass materials and measured at the instant of crack branching of a dynamically running crack. Note that there is evidence that the fracture energy varies with the crack propagation speed (see, e.g. Döll 1975). However, other authors point out the fact that the apparent increase in the fracture energy with crack speed may be due to the presence of microcracks (see, e.g. Cox et al. 2005; Ravi-Chandar 1998). For simplicity, in this work, we keep the fracture energy constant and equal to that measured at crack branching. Nevertheless, it is very easy to introduce the velocity-dependent fracture energy (Döll 1975) in our model and we plan to do so in the future. Moreover, it has been observed that peridynamics generates fragment size distributions closer to experimentally measured ones if the critical relative elongation (and therefore the fracture energy) depends on the damage index (ratio of number of broken bonds and number of initial intact bonds) in the following way: if the damage index is larger than some fraction (say 0.2) then the critical relative elongation value increases with the damage index (see Silling 2005). Such damage-dependent model is used in Ha and Bobaru (2010) for dynamic fracture problems. Note that other models have been used in the past to explain dynamic instabilities (see Buehler et al. 2003) in dynamic crack propagation using MD models, but this does not appear to be needed in peridynamics to trigger instabilities or crack branching.

Additionally, since we apply the load abruptly along the upper and bottom boundaries, relatively high tensile stress act along these boundar-

ies (see the first figure (a) in Figure 4) at the early stages of the simulation. In order to prevent tearing of the first few layers of nodes from the rest of the plate we set the boundary nodes as no-fail zones. In a no-fail zone the damage index is always zero. The traction boundary conditions are applied to a single layer of nodes at the surface in peridynamics, which is similar to how one imposes these conditions in the FEM, for example. The numerical implementation of traction boundary conditions and the convergence studies are shown in Ha and Bobaru (2009). These loads are applied suddenly and a shock wave propagates. In Figure 4, we show a few snapshots of the strain energy and close-ups around the tip of the pre-notch of the damage for dynamic crack propagation in the setup shown in Figure 1, for Duran glass. The model has a uniform grid spacing with Δx = 0.125 mm, the horizon δ = 0.5 mm, and a uniform tensile stress σ = 12 MPa is applied. In Figure 4a, the colors denote the magnitude of the elastic strain energy. In Figure 4a note the ripples behind the wave-front caused by the wave dispersion which, in peridynamics, is due to the size of the horizon and the size of the discretization (see discussion of the 1D case in Xie 2005, pp. 40–44). The colors in the right-hand side plots of Figure 4 represent damage levels. Right after the shock wave reaches the center line, the crack starts propagating (compare Figure 4a, c at 6 and 9 μs).

3 Numerical studies of convergence in dynamic crack branching

In peridynamics, we can talk about three types of convergence (see Bobaru et al. 2009 and Figure 5):

• The δ-convergence: δ → 0 and m (= δ/Δx) is fixed or increases with decreasing δ but at a slower rate. In this case the numerical peridynamic approximation converges to an approximation of the classical solution (if this exists), almost everywhere. The larger m is, the closer this approximation becomes. (The convergence is not guaranteed to be uniform in problems with singularities.)

• The m-convergence: δ is fixed and m → ∞. The numerical peridynamic approximation converges to the exact nonlocal peridynamic solution for the given δ.

• The (δm)-convergence: δ → 0 and m increases with decreasing δ, with m increasing faster than δ decreases. In this case numerical peridynamic approximation converges to the an-

S t u d i e s o f dy n a m i c c r a c k p r o p a g a t i o n a n d c r a c k b r a n c h i n g w i t h p e r i dy n a m i c s

235

Figure 4. Elastic strain energy and damage map around the precrack tip, at the initial stages of crack propagation (δ = 0.0005m, Δx = 0.000125 m, applied load σ = 12MPa). a) Elastic strain energy; b) Time; c Crack-tip near-view.

Figure 5. Graphical descriptions for the a) m-convergence and b) δ -convergence.

alytical peridynamic solution and converges uniformly to the local classical solution (if this exists), almost everywhere.

Here we are studying the m-convergence and we make some observations related to the δ-convergence for dynamic crack propagation problems. The problem to be analyzed is shown in Figure 1: an edge-notch plate. One way to introduce a pre-crack in peridynamic model is to break all bonds that cross the pre-crack line. Another way is to erase nodes that are along the precrack line in addition to breaking all bonds crossing the lines. The first option, under m-convergence, will result in the same “effective” pre-crack. The second option, under m-convergence, can also maintain the pre-crack but only if the total volume of the nodes removed remains the same for any grid spacing used. Note, however, that under δ-convergence things are more delicate. The reason is that if we change the horizon the damage area along the pre-

crack and in front of the pre-crack tip changes independent of the way in which we introduce the pre-crack in the peridynamic model.

In Figure 6, the initial damage areas of two models with different grid spacings are compared. The color-bar represents damage levels. The thick black line denotes the pre-notch. The triangles and the squares are nodes of the coarse and fine models, respectively. In Figure 6a, the damage area of the coarse model matches with the area of the fine model for the same horizon size of δ = 2 mm. However, the damage area of the fine model with δ = 1 mm is slightly smaller than the area of the coarse model with δ = 2 mm, as shown in Figure 6b.

3.1 The m-convergence study

For time integration we use an explicit method, the Velocity–Verlet algorithm. The Velocity–Verlet algorithm (Hairer et al. 2003) is:

236

Y. D. Ha & F. Bobaru in International Journal of Fracture 162 (2010)

Figure 6. Relations between the horizon, grid spacing, and the damage area. (triangles = nodes of the coarse model, squares = nodes of the fine model). a) Damage areas on two grids with same horizon (δ = 0.002 m for both models); b) Damage areas on two grids with different horizons (δ = 0.002 m for coarse model, 0.001 m for fine model).

u˙ n + ½ = u˙ n + Δ2t ün

un + 1 = un + Δtu˙ n + ½

u˙ n + 1 = u˙ n+ ½ + Δ2t ün + 1

(11a) (11b) (11c)

where u, u˙ , and ü denote the displacement, velocity, and acceleration vectors, respectively.

We perform the m-convergence tests for two different horizon sizes: δ = 3 mm and δ = 2 mm. Please note that these horizon sizes are relatively large compared to the structural dimensions. All models have the uniform grid spacing. A uniform time step size of 25 ns is used and this is a stable time step for the finest model among all tests performed in this paper, with δ = 0.5 mm and m = 4. A uniform tensile stress σ = 12 MPa is applied (as described in the previous section) for all the tests in this section. All computations in this section use the constant micro-modulus function (Equation 6) and the Duran 50 glass material parameters.

We first perform the m-convergence study for the fixed horizon δ = 3 mm. The peridynamic models used for this study are the ones with Δx = 1 mm (4,326 nodes), Δx = 0.5 mm (16,646 nodes), and Δx = 0.25 mm (65,448 nodes). For this test, the model with Δx = 1 mm has maximum 29 nodes in the horizon (m = 3), the model with Δx = 0.5 mm has a maximum of 113 nodes in the horizon (m = 6), and the one with Δx = 0.25 mm has maximum 441 nodes in the horizon (m = 12). The crack paths for these cases at 46 μs are compared in Figure 7. In all damage map plots we use the same range for the color-bar of the damage index as in Figure 6. The results in Figure 7 show that m-convergence of the crack path is obtained even for m-values as

small as 3. In all cases in Figure 7, the crack starts propagating around 7 μs, and the crack branches around 25 μs nearby 0.071 m measured from the left-side of the plate. In these results we notice that a thicker damage zone is produced before the crack branches. This may be an indication of the fracture surface roughness prior to branching that is observed consistently in all reported experimental investigations of crack branching in brittle materials (see, e.g. Ramulu and Kobayashi 1985 and references therein).

For the fixed horizon δ = 2 mm, the m-convergence study is performed using the same three numerical grids as above. Since the horizon is smaller, the corresponding m values will be smaller than previously: the coarsest model with Δx = 1 mm has maximum 13 nodes in the horizon (m = 2), the one with Δx = 0.5 mm has maximum 49 nodes in the horizon (m = 4), and the one with Δx = 0.25 mm has maximum 197 nodes in the horizon (m = 8). The crack branching path at 40 μs for our peridynamic simulations are shown in Figure 8. The crack path given by the coarsest model with m = 2 is slightly different from the others, but the paths with m = 4 and 8 are very similar to each other. In Figure 8, crack branching takes place around 24 μs and around 0.068 m from the left edge.

These peridynamic results for two different horizon sizes indicate that m-convergence takes place for the dynamic crack branching problem in terms of the crack path and the crack propagation speed, the latter because the crack tip locations at the same times are similar among the different solutions. This is expected to hold for any given δ. It appears that m = 4 is a good choice because the

S t u d i e s o f dy n a m i c c r a c k p r o p a g a t i o n a n d c r a c k b r a n c h i n g w i t h p e r i dy n a m i c s

237

Figure 7. Crack path computed with different grids for δ = 0.003 m at 46 μs. a) m = 3, Δx = 0.001 m; b) m = 6, Δx = 0.0005 m; c) m = 12, Δx = 0.00025 m.

Figure 8. Crack path computed with different grids for δ = 0.002 m at 46 μs. a) m = 2, Δx = 0.001 m; b) m = 4, Δx = 0.0005 m; c) m = 8, Δx = 0.00025 m.

number of nodes inside horizon allows a sufficiently large number of directions along which the true crack path can develop. Using a larger value of m requires higher computational cost, while the results are not affected. In all remaining tests we will use this value of m.

3.2 Crack path under changing horizon δ

The δ-convergence has to be treated carefully for problems with initial cracks or notches due to the changing size of the initial damage area as discussed in the beginning of Section 3 and Figure 6. The tests in this section are also for the Duran 50 glass material and the constant micro-modulus function (Equation 6).

For a fixed m = 4, the impact of a changing δ is investigated by using the four different kinds of horizons, and therefore, four different grids: the

coarsest model has δ = 4 mm with the uniform grid spacing of Δx = 1 mm (4,326 nodes), the subsequent models have half the horizon size of the previous model and half the grid spacing. Thus, the other three models have, respectively, δ = 2 mm and Δx = 0.5 mm (16,646 nodes), δ = 1 mm and Δx = 0.25 mm (65,448 nodes), and δ = 0.5 mm and Δx = 0.125 mm (258,566 nodes). A uniform time step size of 25 ns (which is a stable time step for the finest model) and the uniform tensile stress σ = 12 MPa are applied for all peridynamic models. In all models, the maximum number of nodes in each horizon is 49. The critical relative elongation s0 also changes with a decreasing horizon, see Equation (9). Also, the damage area becomes smaller as the horizon deceases, and this is especially important for the area in front of the pre-crack tip. The δ-convergence, here, has to be understood within this context. Note that there exists the possibility

[email protected] of Nebraska - Lincoln

Faculty Publications from the Department of Engineering Mechanics

Mechanical & Materials Engineering, Department of

2010

Studies of dynamic crack propagation and crack branching with peridynamics

Youn Doh Ha Ph.D. University of Nebraska at Lincoln

Florin Bobaru Ph.D. University of Nebraska at Lincoln, [email protected]

Follow this and additional works at: https://digitalcommons.unl.edu/engineeringmechanicsfacpub Part of the Applied Mechanics Commons, Ceramic Materials Commons, Computational Engineering

Commons, Engineering Mechanics Commons, Mechanics of Materials Commons, Structural Engineering Commons, Structural Materials Commons, and the Structures and Materials Commons

Ha, Youn Doh Ph.D. and Bobaru, Florin Ph.D., "Studies of dynamic crack propagation and crack branching with peridynamics" (2010). Faculty Publications from the Department of Engineering Mechanics. 71. https://digitalcommons.unl.edu/engineeringmechanicsfacpub/71

This Article is brought to you for free and open access by the Mechanical & Materials Engineering, Department of at [email protected] of Nebraska - Lincoln. It has been accepted for inclusion in Faculty Publications from the Department of Engineering Mechanics by an authorized administrator of [email protected] of Nebraska - Lincoln.

Published in International Journal of Fracture 162 (2010), pp. 229–244; doi: 10.1007/s10704-010-9442-4 Copyright © 2010 Springer Science+Business Media B.V. Used by permission. Submitted August 31, 2009; accepted January 5, 2010; published online January 20, 2010.

Studies of dynamic crack propagation and crack branching with peridynamics

Youn Doh Ha and Florin Bobaru

Department of Engineering Mechanics, University of Nebraska–Lincoln, Lincoln, NE 68588-0526, USA Corresponding author — F. Bobaru, email [email protected]

Abstract In this paper we discuss the peridynamic analysis of dynamic crack branching in brittle materials and show results of convergence studies under uniform grid refinement (m-convergence) and under decreasing the peridynamic horizon (δ-convergence). Comparisons with experimentally obtained values are made for the crack-tip propagation speed with three different peridynamic horizons. We also analyze the influence of the particular shape of the micro-modulus function and of different materials (Duran 50 glass and soda-lime glass) on the crack propagation behavior. We show that the peridynamic solution for this problem captures all the main features, observed experimentally, of dynamic crack propagation and branching, as well as it obtains crack propagation speeds that compare well, qualitatively and quantitatively, with experimental results published in the literature. The branching patterns also correlate remarkably well with tests published in the literature that show several branching levels at higher stress levels reached when the initial notch starts propagating. We notice the strong influence reflecting stress waves from the boundaries have on the shape and structure of the crack paths in dynamic fracture. All these computational solutions are obtained by using the minimum amount of input information: density, elastic stiffness, and constant fracture energy. No special criteria for crack propagation, crack curving, or crack branching are used: dynamic crack propagation is ob-

tained here as part of the solution. We conclude that peridynamics is a reliable formulation for modeling dynamic crack propagation.

Keywords: dynamic fracture, crack branching, brittle fracture, peridynamics, nonlocal methods, meshfree methods

1 Introduction

1.1 Literature review of dynamic crack propagation

In a brittle material, a propagating crack can depart from its original straight trajectory and curve or split into two or more branches. Under very high states of stress, the propagating crack will divide into a river-delta crack pattern (Bowden et al. 1967; Ramulu and Kobayashi 1985). This fragmentation of highly loaded, brittle materials is often a succession of multiple branching of what was initially a single crack. Increases in the roughness of the fracture surface prior to branching were consistently observed in all reported investigations (Ramulu and Kobayashi 1985; Döll 1975). In crack branching of edge notch specimens of brittle materials it has also been observed that the crack tip velocity drops by no more than 5–10% in the branching region (Döll 1975).

In atomistic models, under conditions that lead to instability of the crack path, cracks can branch without a specific criterion (see Zhou et al. 1996).

229

230

Y. D. Ha & F. Bobaru in International Journal of Fracture 162 (2010)

Particle-type models (see Bolander and Saito 1998) are also capable of simulating crack branching. However, none of these methods is able to capture the crack propagation speed or the angle of crack branching correctly. For instance, MD simulations show instabilities that lead, shortly after the bifurcation of a crack, to the propagation of only one of the two branches, the other being arrested. Moreover, the branching angle computed with MD (see Zhou et al. 1996) is greater than 90°, whereas experiments show much smaller crack branching angles (Ramulu and Kobayashi 1985). One may ask whether quantum mechanical calculations are needed to predict the phenomenon of dynamic fracture in brittle materials (see Cox et al. 2005) which is one of the great challenges in dynamic fracture. One likely reason for MD simulations’ failure to correctly predict dynamic fracture is that, for example, crack branching events are controlled by the interaction and wave reflections from the boundaries (Ravi-Chandar 1998). Because of this, one would have to either model the entire structure with MD (not a viable option) or use a multiscale model that is capable of transferring the waves between the scales correctly (still an open problem). Numerical simulations based on continuum methods of dynamic crack propagation behavior have been very difficult to develop and, to this date, a reliable method for simulating this complex problem has not been found in spite of considerable efforts in this direction (e.g. Xu and Needleman 1994; Camacho and Ortiz 1996; Ortiz and Pandolfi 1999; Belytschko et al. 2003; Rabczuk and Belytschko 2004; Song et al. 2006). All these methods use some version of cohesive-zone models. As such, they all modify the local continuum mechanics equations and introduce a nonlocal effect given by the parameters and length scales in the cohesive-zone model. To reduce mesh dependency when the grids are refined special methodologies have to be used (Zhou and Molinari 2004). For the existing approaches, the difficulties in modeling dynamic fracture processes like crack branching are many. For example, continuum-type methods using the cohesive FEM or the XFEM require a damage criterion and a tracking of the stresses around the crack tip to decide when to branch the crack. Decisions also have to be taken in terms of the angle of propagation of the branches and about how many branches will be allowed to form. In methods in which the crack advances along the element sides by separating elements from one another,

the crack path becomes non-smooth (see Xu and Needleman 1994; Camacho and Ortiz 1996; Ortiz and Pandolfi 1999). Since the correct path (which minimizes the strain energy) of the crack propagation is not computed correctly, there are significant departures from the true energy released during the crack propagation event. In such cases, reliable prediction of strength of brittle ceramics under impact, for example, becomes difficult. Mesh dependency is an additional problem in cohesive-zone FEM-based methods. Important progress has been recently made by using the XFEM method which allows cracks to pass through the finite elements (see e.g. Belytschko et al. 2003). Subdivision of the cut elements for numerical integration purposes increases the complexity and the cost of the method. This method requires phenomenological damage models and branching criteria, as well as tracking of the crack path using level sets, for example. It is not yet clear if the method is applicable to problems that involve fragmentation and/or multiple crack interactions, branching, and coalescence. The method does not predict the experimentally observed crack propagation speeds (see Song et al. 2008). Cohesivezone based models need to modify the experimental values of the fracture energy by several factors in order to get propagation velocities in the range of measured ones (Song et al. 2008).

In the present contribution we try to answer whether quantum, atomistic, or multiscale models are needed in dynamic fracture in order to correctly simulate the observed crack propagation velocities and crack paths (Cox et al. 2005; Song et al. 2008). We will show that peridynamics is able to correctly model and simulate dynamic fracture, in particular crack branching in brittle materials. Peridynamics, which is a reformulation of continuum mechanics (Silling 2000; Silling et al. 2007b), does not require criteria for crack propagation or crack branching: these happen spontaneously in this method and are autonomously generated by a simple bond-failure criterion that is correlated to the material’s energy release rate. The name “peridynamics” comes from the Greek “peri” which means “nearby,” and dynamics. Peridynamics is a nonlocal method in which material points interact not only with their nearest neighbors but also with points nearby, inside a horizon. This is what physically happens at the atomic scale, for example, but peridynamics extends this idea to the continuum scale. We will show convergence in terms of the crack path and the crack propagation speed under

S t u d i e s o f dy n a m i c c r a c k p r o p a g a t i o n a n d c r a c k b r a n c h i n g w i t h p e r i dy n a m i c s

231

grid refinement (the so-called δ-convergence, Bobaru et al. 2009) and under decreasing peridynamic horizon (the so-called δ-convergence, Bobaru et al. 2009). The crack branching patterns obtained using peridynamics follow remarkably close the experimental results which show secondary branching taking place when higher stress levels are reached at the tip of the pre-notch prior to crack propagation. Moreover, the only input parameters in the model are the Young’s modulus, the density, and the fracture energy (which is kept constant, and not a function of the propagation velocity or of the incurred damage, in this first study).

The paper is organized as follows: in the next section we describe the sample problem setup. In Section 2 we briefly review the peridynamic formulation and the connections between the parameters in the formulation and the material properties like the energy release rate. In Section 3 we present the numerical results for the convergence studies. We look at both the crack path and the propagation speed of the crack, as measures of convergence. In Section 4 we analyze the influence of the micromodulus function on crack branching results as well as the solutions for two different brittle materials under higher loading conditions that lead to cascading branching. We also comment on the roughening zones that take place in the branching regions and on the effect of the reflection waves on the propagation paths of the dynamic cracks. The conclusions are given in Section 5.

1.2 Problem setup

We consider the following setup as a benchmark problem for analyzing crack branching phenomena: a prenotched thin rectangular plate with 0.1m by 0.04m as shown in Figure 1. All simulations in this paper are 2D simulations. For some 3D results we refer the readers to Ha et al. (2010). The materials chosen for this study are selected because for these materials there are experimental results published on the crack propagation velocity in the region of branching or the maximum propagation velocity measured. The two materials used here are a Duran 50 glass (taken form Döll 1975) and a soda-lime glass (taken from Bowden et al. 1967). The material properties are summarized in Table 1. Please note that in the bond-based peridynamic implementation used here, the numerical models will be limited to using a fixed Poisson ratio of 1/3 (for 2D plane stress problems). If other Poisson ratios are desired, then the state-based

Figure 1. Description of the problem setup for the crack branching study.

peridynamics formulation should be used (see Silling et al. 2007b). For dynamic fracture problems, the Poisson ratio value does not have a significant influence on the propagation speed or the crack path shapes (Silling et al. 2007a).

In the experimental settings the loading of the sample may take tens of seconds or more. In explicit dynamic simulations that would be too expensive to compute. Instead, we choose to apply, along the upper and lower edges (see Figure 1), traction loadings σ suddenly at the initial time step and maintain this loading constant after that. The theoretical background for the peridynamics analysis is based on Silling’s original peridynamics paper (Silling 2000), the imposition of the traction boundary conditions is as in Ha and Bobaru (2009), and the numerical implementation of failure is like in Silling and Askari (2005). The same geometrical setup for studying crack branching simulations has been used in other studies (Belytschko et al. 2003; Rabczuk and Belytschko 2004; Song et al. 2006).

While there is no analytical solution for the crack branching problem, we can compare our simulation results with experiments. Unfortunately, the experimental papers we found do not provide a complete description of the conducted experiment on crack branching: some papers show the crack paths but do not provide crack propagation speed data, others give the propagation speed but do not show the crack paths, and most do not describe in detail the loading conditions. We decided to perform the peridynamic simulations using a setup similar to that used in a few recent simulation papers (Belytschko et al. 2003; Rabczuk and Belytschko 2004; Song et al. 2006). The material parameters, however, are like those used in the experiments (Bowden et al. 1967; Döll 1975). The maximum crack propagation speed, or the

232

Y. D. Ha & F. Bobaru in International Journal of Fracture 162 (2010)

Table 1. Material properties for Duran 50 and soda-lime glasses

Density (ρ) (kg/m3) Young’s modulus (E) (GPa)

Duran 50 glass

2,235

65

Soda-lime glass

2,440

72

Poisson ratio (υ) Fracture energy (G0) (J/m2)

0.2

204

0.22

135

crack propagation speed in the region of branching, is data that is fairly reproducible in experiments and this is reported in Bowden et al. (1967) and Döll (1975), for example. We are not aware of any numerical method that can reproduce the experimentally measured dynamic crack propagation velocity. Note that in previous studies, for certain methods, the fracture energy has to be significantly modified (by several factors) in order to bring the dynamic crack propagation speed closer to the measured values (see, e.g. Belytschko et al. 2003; Rabczuk and Belytschko 2004).

2 The peridynamic formulation

The peridynamic formulation (Silling 2000) relies on integration of forces acting on a material point and thus it does not face any of the mathematical inconsistencies seen in the classical continuum mechanics equations. The integration takes place over a “horizon” (which, in principle, extends to infinity but, for convenience is finite) within which the material points are interacting with each other. In certain problems, the size of the horizon can be correlated to an intrinsic material length-scale. However, in many cases a material length scale is not “visible” either because the micro-structure and the specific loading and boundary conditions do not lead to a measurable effect of the length-scale. In such cases, the horizon is selected by the user according to convenience (see Bobaru et al. 2009). Allowing a variable horizon (with a correspondingly scaled micromodulus parameter) defines away of introducing adaptive refinement for this nonlocal method. It is important to notice that peridynamics is a continuum theory, not a particle-type method. This allows the convergence results of the peridynamic solution to the classical elasticity solutions in the limit of the horizon going to zero (Bobaru et al. 2009; Silling and Lehoucq 2008).

An important advantage of peridynamics is the way damage is introduced: material points are connected within the horizon via elastic (linear or nonlinear) bonds that have a critical relative elongation, s0, at which they break (Silling 2000).

The critical relative elongation for brittle materials is computed from the experimentally measured value of the fracture energy for a specific material (Silling and Askari 2005). Damage is implemented as the fraction between the number of broken bonds and the number of initial bonds (Silling and Askari 2005). Cracks in peridynamics form as surfaces between material points form, as a consequence of sequential breaking of bonds. Thus, there is no need to track the cracks like in other continuum methods, or to impose criteria for when cracks should branch, change direction, turn, coalesce, etc. Moreover, peridynamics allows for spontaneous generation of cracks where no flaws were present before. This is shown, for example, in Silling et al. (2009) for the crack nucleation and in simulation of spallation (see Xie 2005) where spallation is treated as real fracture and not modeled by void-growth formulations as in existing literature results.

We now briefly review the peridynamic formulation based on Silling’s original peridynamics paper (Silling 2000). Also, we consider the summary of the numerical implementation of the traction boundary conditions in peridynamics (Ha and Bobaru 2009) and the formulation for the damage model in peridynamics (Silling and Askari 2005).

The peridynamic equations of motion are given by:

∫ ρü (x, t) = f(u (xˆ , t) − u(x, t), xˆ − x) dxˆ + b(x, t) (1) H

where f is the pairwise force function in the peridynamic bond that connects node xˆ to x and u is the displacement vector field. ρ is the density and b (x, t) is the body force. The integral is defined over a region H called the “horizon,” which is the compact supported domain of the pairwise force function around point x.

A micro-elastic material (Silling 2000) is defined as one for which the pairwise force derives from a potential ω:

∂ω (η, ξ )

f (η, ξ ) = ∂η

(2)

where ξ = xˆ − x is the relative position in the refer-

ence configuration and η = uˆ − u is the relative dis-

S t u d i e s o f dy n a m i c c r a c k p r o p a g a t i o n a n d c r a c k b r a n c h i n g w i t h p e r i dy n a m i c s

233

placement. A linear micro-elastic potential is obtained if we take

c (ξ ) s2ξ

ω (η, ξ ) = 2

(3)

where ξ = ξ and the relative elongation of a bond is

ζ − ξ

s= ξ

(4)

where ζ = ξ + η . The function c (ξ ) is called micro-modulus and has the meaning of the bond elastic stiffness. There are various formulations for the micromodulus function and in Section 4.1 we perform tests for dynamic crack propagation to assess the influence of the particular shape of the micro-modulus function on the crack path structure. We will observe that the crack propagation speed is not influenced by the shape of the micro-modulus, once the horizon is reasonably small compared to the dimensions of the structure analyzed. The pairwise force corresponding to a linear microelastic potential has the following form:

{ ξ+η

f (η, ξ ) = ξ + η c (ξ ) s,

ξ ≤ δ

0,

ξ>δ

(5)

where δ is the radius of the horizon region (which we will also refer to as the horizon). In this paper, we use the constant and conical 2D micro-modulus functions (see Figure 2). Following the same procedure performed to calculate the micro-modulus function in 1D (Bobaru et al. 2009), we obtain the constant micro-modulus function in 2D, plane stress conditions:

6E

c (ξ ) = c0 = πδ3 (1 − υ)

(6)

Similarly, the conical micro-modulus function is obtained as

c

(ξ ) = c1 ( 1 −

ξ δ

)

=

π

δ

3

24E (1 −

υ)

(

1

−

ξ δ

)

(7)

The shapes of the constant and conical micromodulus functions are illustrated in Figure 2.

In the bond-based peridynamics, any particle inside the horizon of another particle interacts only through a central potential. This assumption results (for an isotropic, linear, micro-elastic material) in an effective Poisson ratio of 1/3 in 2D (and 1/4 in 3D), but this limitation is readily eliminated by using the state-based peridynamics (Silling et al. 2007a,b). In this paper, we utilize the bond-

Figure 2. Constant (left) and conical (right) micromodulus functions.

based peridynamics, thus, in all the reported simulations here the effective Poisson ratio is 1/3.

In order to introduce failure into the peridynamic model, we consider that the peridynamic bonds can be broken when they are stretched beyond a predefined limit. We call this limit the “critical relative elongation, s0.” According to Silling and Askari (2005), there is no force sustained by the bond after its failure. Also, once a bond fails, it is failed forever; this makes the model history dependent. To completely separate a body into two halves across a fracture plane requires breaking all the bonds that initially connected points in the opposite halves (see Silling and Askari 2005). The energy per unit fracture length (in 2D, fracture area in 3D) for complete separation of the two halves of the body is called fracture energy, G0. In 3D, Silling and Askari relate the critical elongation, s0, with this measurable quantity (G0) (Silling and Askari 2005). In 2D with plane stress conditions, the fracture energy can be derived as

δ δ cos −1(z/ξ )

∫ ∫ ∫ G0 = 2

0 z

0

[c

(ξ

) 2

s02

ξ

]

ξ

dθ

dξ

dz

(8)

(See Figure 3 for an explanation of this computation.) Substituting the constant micro-modulus function (Equation 6) and rearranging for s0, we can rewrite this equation to obtain s0

s0

=

√

4 π G0 9Eδ

(9)

In similar way, the critical relative elongation for the conical micro-modulus function (Equation 7) is

√s0 = 5πG0 (10) 9Eδ

The critical relative elongation depends on the material properties and the horizon δ. Note that as the horizon goes to zero, the critical relative elongation goes to infinity, thus breaking such bonds

234

Y. D. Ha & F. Bobaru in International Journal of Fracture 162 (2010)

Figure 3. Evaluation of fracture energy. For each point A along the dashed line, 0 ≤ z ≤ δ, the work required to break the bonds connecting A to each point B in the circular cap is given by Equation (8).

requires larger and larger forces. This agrees with the physical experience at atomic and subatomic scales where, in order to break apart smaller and smaller sized bonds, one needs higher and higher forces. The values for the fracture energy used in this paper are the ones given in Döll (1975) for Duran 50 glass and soda-lime glass materials and measured at the instant of crack branching of a dynamically running crack. Note that there is evidence that the fracture energy varies with the crack propagation speed (see, e.g. Döll 1975). However, other authors point out the fact that the apparent increase in the fracture energy with crack speed may be due to the presence of microcracks (see, e.g. Cox et al. 2005; Ravi-Chandar 1998). For simplicity, in this work, we keep the fracture energy constant and equal to that measured at crack branching. Nevertheless, it is very easy to introduce the velocity-dependent fracture energy (Döll 1975) in our model and we plan to do so in the future. Moreover, it has been observed that peridynamics generates fragment size distributions closer to experimentally measured ones if the critical relative elongation (and therefore the fracture energy) depends on the damage index (ratio of number of broken bonds and number of initial intact bonds) in the following way: if the damage index is larger than some fraction (say 0.2) then the critical relative elongation value increases with the damage index (see Silling 2005). Such damage-dependent model is used in Ha and Bobaru (2010) for dynamic fracture problems. Note that other models have been used in the past to explain dynamic instabilities (see Buehler et al. 2003) in dynamic crack propagation using MD models, but this does not appear to be needed in peridynamics to trigger instabilities or crack branching.

Additionally, since we apply the load abruptly along the upper and bottom boundaries, relatively high tensile stress act along these boundar-

ies (see the first figure (a) in Figure 4) at the early stages of the simulation. In order to prevent tearing of the first few layers of nodes from the rest of the plate we set the boundary nodes as no-fail zones. In a no-fail zone the damage index is always zero. The traction boundary conditions are applied to a single layer of nodes at the surface in peridynamics, which is similar to how one imposes these conditions in the FEM, for example. The numerical implementation of traction boundary conditions and the convergence studies are shown in Ha and Bobaru (2009). These loads are applied suddenly and a shock wave propagates. In Figure 4, we show a few snapshots of the strain energy and close-ups around the tip of the pre-notch of the damage for dynamic crack propagation in the setup shown in Figure 1, for Duran glass. The model has a uniform grid spacing with Δx = 0.125 mm, the horizon δ = 0.5 mm, and a uniform tensile stress σ = 12 MPa is applied. In Figure 4a, the colors denote the magnitude of the elastic strain energy. In Figure 4a note the ripples behind the wave-front caused by the wave dispersion which, in peridynamics, is due to the size of the horizon and the size of the discretization (see discussion of the 1D case in Xie 2005, pp. 40–44). The colors in the right-hand side plots of Figure 4 represent damage levels. Right after the shock wave reaches the center line, the crack starts propagating (compare Figure 4a, c at 6 and 9 μs).

3 Numerical studies of convergence in dynamic crack branching

In peridynamics, we can talk about three types of convergence (see Bobaru et al. 2009 and Figure 5):

• The δ-convergence: δ → 0 and m (= δ/Δx) is fixed or increases with decreasing δ but at a slower rate. In this case the numerical peridynamic approximation converges to an approximation of the classical solution (if this exists), almost everywhere. The larger m is, the closer this approximation becomes. (The convergence is not guaranteed to be uniform in problems with singularities.)

• The m-convergence: δ is fixed and m → ∞. The numerical peridynamic approximation converges to the exact nonlocal peridynamic solution for the given δ.

• The (δm)-convergence: δ → 0 and m increases with decreasing δ, with m increasing faster than δ decreases. In this case numerical peridynamic approximation converges to the an-

S t u d i e s o f dy n a m i c c r a c k p r o p a g a t i o n a n d c r a c k b r a n c h i n g w i t h p e r i dy n a m i c s

235

Figure 4. Elastic strain energy and damage map around the precrack tip, at the initial stages of crack propagation (δ = 0.0005m, Δx = 0.000125 m, applied load σ = 12MPa). a) Elastic strain energy; b) Time; c Crack-tip near-view.

Figure 5. Graphical descriptions for the a) m-convergence and b) δ -convergence.

alytical peridynamic solution and converges uniformly to the local classical solution (if this exists), almost everywhere.

Here we are studying the m-convergence and we make some observations related to the δ-convergence for dynamic crack propagation problems. The problem to be analyzed is shown in Figure 1: an edge-notch plate. One way to introduce a pre-crack in peridynamic model is to break all bonds that cross the pre-crack line. Another way is to erase nodes that are along the precrack line in addition to breaking all bonds crossing the lines. The first option, under m-convergence, will result in the same “effective” pre-crack. The second option, under m-convergence, can also maintain the pre-crack but only if the total volume of the nodes removed remains the same for any grid spacing used. Note, however, that under δ-convergence things are more delicate. The reason is that if we change the horizon the damage area along the pre-

crack and in front of the pre-crack tip changes independent of the way in which we introduce the pre-crack in the peridynamic model.

In Figure 6, the initial damage areas of two models with different grid spacings are compared. The color-bar represents damage levels. The thick black line denotes the pre-notch. The triangles and the squares are nodes of the coarse and fine models, respectively. In Figure 6a, the damage area of the coarse model matches with the area of the fine model for the same horizon size of δ = 2 mm. However, the damage area of the fine model with δ = 1 mm is slightly smaller than the area of the coarse model with δ = 2 mm, as shown in Figure 6b.

3.1 The m-convergence study

For time integration we use an explicit method, the Velocity–Verlet algorithm. The Velocity–Verlet algorithm (Hairer et al. 2003) is:

236

Y. D. Ha & F. Bobaru in International Journal of Fracture 162 (2010)

Figure 6. Relations between the horizon, grid spacing, and the damage area. (triangles = nodes of the coarse model, squares = nodes of the fine model). a) Damage areas on two grids with same horizon (δ = 0.002 m for both models); b) Damage areas on two grids with different horizons (δ = 0.002 m for coarse model, 0.001 m for fine model).

u˙ n + ½ = u˙ n + Δ2t ün

un + 1 = un + Δtu˙ n + ½

u˙ n + 1 = u˙ n+ ½ + Δ2t ün + 1

(11a) (11b) (11c)

where u, u˙ , and ü denote the displacement, velocity, and acceleration vectors, respectively.

We perform the m-convergence tests for two different horizon sizes: δ = 3 mm and δ = 2 mm. Please note that these horizon sizes are relatively large compared to the structural dimensions. All models have the uniform grid spacing. A uniform time step size of 25 ns is used and this is a stable time step for the finest model among all tests performed in this paper, with δ = 0.5 mm and m = 4. A uniform tensile stress σ = 12 MPa is applied (as described in the previous section) for all the tests in this section. All computations in this section use the constant micro-modulus function (Equation 6) and the Duran 50 glass material parameters.

We first perform the m-convergence study for the fixed horizon δ = 3 mm. The peridynamic models used for this study are the ones with Δx = 1 mm (4,326 nodes), Δx = 0.5 mm (16,646 nodes), and Δx = 0.25 mm (65,448 nodes). For this test, the model with Δx = 1 mm has maximum 29 nodes in the horizon (m = 3), the model with Δx = 0.5 mm has a maximum of 113 nodes in the horizon (m = 6), and the one with Δx = 0.25 mm has maximum 441 nodes in the horizon (m = 12). The crack paths for these cases at 46 μs are compared in Figure 7. In all damage map plots we use the same range for the color-bar of the damage index as in Figure 6. The results in Figure 7 show that m-convergence of the crack path is obtained even for m-values as

small as 3. In all cases in Figure 7, the crack starts propagating around 7 μs, and the crack branches around 25 μs nearby 0.071 m measured from the left-side of the plate. In these results we notice that a thicker damage zone is produced before the crack branches. This may be an indication of the fracture surface roughness prior to branching that is observed consistently in all reported experimental investigations of crack branching in brittle materials (see, e.g. Ramulu and Kobayashi 1985 and references therein).

For the fixed horizon δ = 2 mm, the m-convergence study is performed using the same three numerical grids as above. Since the horizon is smaller, the corresponding m values will be smaller than previously: the coarsest model with Δx = 1 mm has maximum 13 nodes in the horizon (m = 2), the one with Δx = 0.5 mm has maximum 49 nodes in the horizon (m = 4), and the one with Δx = 0.25 mm has maximum 197 nodes in the horizon (m = 8). The crack branching path at 40 μs for our peridynamic simulations are shown in Figure 8. The crack path given by the coarsest model with m = 2 is slightly different from the others, but the paths with m = 4 and 8 are very similar to each other. In Figure 8, crack branching takes place around 24 μs and around 0.068 m from the left edge.

These peridynamic results for two different horizon sizes indicate that m-convergence takes place for the dynamic crack branching problem in terms of the crack path and the crack propagation speed, the latter because the crack tip locations at the same times are similar among the different solutions. This is expected to hold for any given δ. It appears that m = 4 is a good choice because the

S t u d i e s o f dy n a m i c c r a c k p r o p a g a t i o n a n d c r a c k b r a n c h i n g w i t h p e r i dy n a m i c s

237

Figure 7. Crack path computed with different grids for δ = 0.003 m at 46 μs. a) m = 3, Δx = 0.001 m; b) m = 6, Δx = 0.0005 m; c) m = 12, Δx = 0.00025 m.

Figure 8. Crack path computed with different grids for δ = 0.002 m at 46 μs. a) m = 2, Δx = 0.001 m; b) m = 4, Δx = 0.0005 m; c) m = 8, Δx = 0.00025 m.

number of nodes inside horizon allows a sufficiently large number of directions along which the true crack path can develop. Using a larger value of m requires higher computational cost, while the results are not affected. In all remaining tests we will use this value of m.

3.2 Crack path under changing horizon δ

The δ-convergence has to be treated carefully for problems with initial cracks or notches due to the changing size of the initial damage area as discussed in the beginning of Section 3 and Figure 6. The tests in this section are also for the Duran 50 glass material and the constant micro-modulus function (Equation 6).

For a fixed m = 4, the impact of a changing δ is investigated by using the four different kinds of horizons, and therefore, four different grids: the

coarsest model has δ = 4 mm with the uniform grid spacing of Δx = 1 mm (4,326 nodes), the subsequent models have half the horizon size of the previous model and half the grid spacing. Thus, the other three models have, respectively, δ = 2 mm and Δx = 0.5 mm (16,646 nodes), δ = 1 mm and Δx = 0.25 mm (65,448 nodes), and δ = 0.5 mm and Δx = 0.125 mm (258,566 nodes). A uniform time step size of 25 ns (which is a stable time step for the finest model) and the uniform tensile stress σ = 12 MPa are applied for all peridynamic models. In all models, the maximum number of nodes in each horizon is 49. The critical relative elongation s0 also changes with a decreasing horizon, see Equation (9). Also, the damage area becomes smaller as the horizon deceases, and this is especially important for the area in front of the pre-crack tip. The δ-convergence, here, has to be understood within this context. Note that there exists the possibility