On the computation of energy release rates by a peridynamic virtual crack extension method

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 363 (2020) 112883 www.elsevier.com/locate/cma

On the computation of energy release rates by a peridynamic virtual crack extension method Heng Zhanga , Pizhong Qiaoa,b ,∗ a

State Key Laboratory of Ocean Engineering, School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, PR China b Department of Civil and Environmental Engineering, Washington State University, Sloan Hall 117, Pullman, WA 99164-2910, USA Received 4 October 2018; received in revised form 20 January 2020; accepted 25 January 2020 Available online xxxx

Abstract A new peridynamic virtual crack extension (PVCE) method is developed to predict the energy release rate G in fracture mechanics, and the directionally and virtually broken bonds are uniquely applied to constitute the virtual crack for crack extension in the numerical implementation. An algorithm of the PVCE method considering serial crack extensions and quasistatic analysis is proposed, and the modified PVCE method is further introduced to reduce the computational error. As validation and application examples, the developed PVCE method is applied to compute the energy release rates of single edge-notched tension (SENT), double cantilever beam (DCB), and end loaded split (ELS) tests, and the results are compared with the existing analytical solutions and numerical finite element analysis. The new PVCE method builds a natural “bridge” for peridynamic fracture analysis between the peridynamics and classical fracture mechanics, and it is capable of accurately and effectively computing the energy release rates for both mode I and mode II fracture cases as well as predicting the critical load and displacement in fracture analysis. c 2020 Elsevier B.V. All rights reserved. ⃝ Keywords: Peridynamics; Virtual crack extension; Energy release rate; Fracture analysis

1. Introduction Virtual crack extension method [1,2] is a direct way to compute the energy release rate, G, in fracture mechanics. This method consists of calculating the strain energy change due to serial virtual crack propagation and employing numerical differentiation to determine the strain energy release rate. It has been commonly performed with finite element method (FEM) for prediction of crack behavior [2]. Peridynamics [3] is a nonlocal theory, and it reformulates the classical continuum mechanics, in which the integral–differential equations are established to replace the partial differential equations in continuum mechanics. Thus, peridynamics is naturally suitable for discontinuity and fracture analysis, and it has been successfully applied to fracture analysis of various materials [4–12]. In addition, peridynamics has also been extended to solve ∗ Correspondence to: Department of Engineering Mechanics, Shanghai Jiao Tong University, Shanghai 200240, PR China.

E-mail addresses: [email protected], [email protected] (P. Qiao). https://doi.org/10.1016/j.cma.2020.112883 c 2020 Elsevier B.V. All rights reserved. 0045-7825/⃝

2

H. Zhang and P. Qiao / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112883

typical discontinuous problems of delamination of bimaterial structures under thermal loading [13], axisymmetric indentation fracture [14], holed plate failure under tensile loading [15,16]. Several failure criteria have been proposed for peridynamic fracture analysis. Silling and Askari [17] first introduced the bond-based critical stretch model for bond broken criterion, and it was later extended to the statebased peridynamics one by the authors [18]. Recently, Zhang and Qiao [19] proposed a new bond failure (so-called critical skew) criterion for state-based peridynamic mode II fracture analysis. Foster et al. [20] presented an energybased failure criterion used in peridynamic state. Wang et al. [21] considered both the critical stretch and critical shear energy density criteria for the 3D conjugated bond-pair-based peridynamic formulation. Associated with classical fracture mechanics, Hu et al. [22] established the formulation of the J-integral in the framework of bondbased peridynamics, and Stenstr¨om and Eriksson [23] recently reformulated the J-integral in peridynamics by the displacement derivatives entirely. Panchadhara and Gordon [24] applied the peridynamic stress intensity factors to dynamic fracture initiation and propagation using the formulated J-integral, which can be utilized for mixed mode crack problems. Diehl et al. [25] estimated the critical traction for mode I crack opening by the bond-based peridynamic simulation. However, there is no direct method to compute the energy release rates in state-based peridynamics. Thus, the goal of this paper is to implement the virtual crack extension by the peridynamic model to computationally determine the strain energy release rate. In this paper, the measurement of energy release rates is performed by a new peridynamic model, which is called peridynamic virtual crack extension (PVCE) method. The peridynamic strain energies for serial but slightly different crack lengths are calculated, and the numerical differentiation is utilized for computation of energy release rate. First, the theory of virtual crack extension in peridynamics is introduced. An algorithm of the PVCE method considering serial crack extensions and quasi-static analysis is developed, and modified PVCE method to reduce computational error is further introduced. Then, the single edge-notched tension (SENT) test is performed for verification of the PVCE method in calculation of energy release rates, followed by a parametric study. Furthermore, the double cantilever beam (DCB) and end loaded split (ELS) tests are considered for applications of the PVCE method to mode I and II fracture problems, respectively. The predictions of critical load and displacement by the PVCE method are presented and compared with the available analytical solutions and FEM results. 2. Review of peridynamic model In this section, a review about the state-based peridynamic model is briefly presented, and the theory of the virtual crack extension by the peridynamics model is introduced. 2.1. General state-based peridynamic model In the state-based peridynamic model, the motion equation of the material point x can be expressed as [26]: ∫ ⟨ ⟩ [ ]⟨ ⟩} { (1) ρ (x) u¨ (x, t) = T [x, t] x ′ − x − T x ′ , t x − x ′ d Vx ′ + b (x, t) Hx

where ρ is the mass density of point x, u is the displacement of point x at time t, Hx is the neighborhood of x with a cutoff distance δ, which is called as ‘horizon’, x ′ is another material point in the neighborhood Hx , and [ b](x, t) is the body force density of point x. As shown in Fig. 1, ξ = x ′ − x is the bond vector, T [x, t] and T x ′ , t are the force vector states of points x and x ′ , respectively. In the ordinary state-based peridynamic model [26], the force state is parallel to the deformed position of connected bond, as shown in Fig. 1. For the linear elastic isotropic material, the general strain energy density of point x in the ordinary state-based peridynamic model takes the form of [26]: W =

kθ 2 α ( d) d + ωe •e 2 2

(2)

where k and α are the peridynamic parameters, respectively. θ is the volume dilatation, ed is the deviatoric extension state, and ω is the influence function. The symbol (•) is the dot product that is defined in [26]. And these parameters

H. Zhang and P. Qiao / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112883

3

Fig. 1. Ordinary state-based peridynamic model.

can be expressed for the three-dimensional (3D) and two-dimensional (2D) cases, respectively as [18]: ⎧ ωx • e d 15µ ⎪ ⎪ , e = e − θ x/3, k = k ′ , α = ; 3D ⎨θ = 3 q q ωx • e d 8µ ⎪ ⎪ ⎩θ = 2 , e = e − θ x/2, k = k ′ , α = ; 2D q q Thus, the scalar force state t takes the form of [18]: ⎧ ωx 15µ d ⎪ ⎪3k ′ θ + ωe 3D ⎨ q q t= 8µ d ωx ⎪ ⎪ ⎩2k ′ θ + ωe 2D q q

(3)

(4)

where e is the extension scalar state, x is the scalar state of bond ξ that is equal to the bond length |ξ |, q is the weighted volume defined as ωx • x, µ is the shear modulus, and k ′ is the bulk modulus that can be expressed in terms of the Young’s modulus E and the Poisson’s ratio v as: ⎧ E ⎪ ⎪ 3D ⎪ ⎪ 3 − 2v) (1 ⎪ ⎪ ⎨ E Plane stress k′ = (5) ⎪ 2 (1 − v) ⎪ ⎪ ⎪ ⎪ E ⎪ ⎩ Plane strain 2 (1 + v) (1 − 2v) 2.2. Theory of peridynamic virtual crack extension (PVCE) method For the brittle elastic solid material, the energy released from the body must be equal to the dissipated energy (surface energy) during the crack propagation. It means that, for the crack extension of da, the variations of the energies of the whole body must satisfy [27]: d We = dUe + dUk + d Ws

(6)

where We is the total potential of the external force, Ue is the total elastic energy of the cracked body, Uk is the total kinetic energy of the cracked body, and Ws is the total crack surface energy associated with the entire crack. As for quasi-static problems, Uk can be neglected, and the corresponding energy release rate during crack propagation can be expressed as: d Ws d We − dUe G= = (7) da da

4

H. Zhang and P. Qiao / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112883

Considering the fixed load or displacement boundary conditions, the energy release rate can be thus calculated by: ⎧ dUe ⎪ ⎨+ Fixed load da G= (8) ⎪ ⎩− dUe Fixed displacement da Thus, the energy release rate is related to the differentiation of the total strain energy of the cracked body to the crack extension length, and it can be numerically calculated with the total strain energy stored for two or more slightly different crack lengths. For the peridynamic model, after the total strain energies of the cracked body are calculated for several different crack lengths, the numerical differentiation can be utilized to obtain the corresponding energy release rate. For a normal two-step differential, the calculated energy release rate is originally expressed as: ⎧ Ue (a + da) − Ue (a) ⎪ ⎨+ Fixed load da G (a) = (9) ⎪ ⎩− Ue (a + da) − Ue (a) Fixed displacement da where G (a) is the numerical energy release rate with the crack length a, and Ue (a + da) and Ue (a) are the total strain energies with the pre-crack lengths a + da and a, respectively. For linear elastic materials, the total strain energies can be obtained with the formulation of the strain energy density given in Eq. (2) by the state-based peridynamics. 3. Numerical implementation of the PVCE method In the peridynamic numerical simulation, the continuum body is discretized into the finite nodes, and the peridynamic motion equation in Eq. (1) can then be rewritten as: ∑ ⟨ ⟩ ⟨ ⟩ ρi u¨ i (xi , t) = {T [xi , t] x j − xi − T [x j , t] xi − x j }V j + bi (xi , t) (10) j∈Hi

where i and j are the node numbers, and ⏐V j is the⏐ volume of node j. The uniform grid size ∆x is used in the following examples. All nodes j in Hi satisfying ⏐x j − xi ⏐ ≤ δ are summed, where δ is the horizon size and δ = m∆x. 3.1. Quasi-static analysis in peridynamics The quasi-static analysis is performed to obtain deformation response of the cracked body under static load condition. In the following examples, the adaptive dynamic relaxation (ADR) method [28] is used for the quasi-static analysis. Based on the ADR method using the central-explicit integration, the displacements and velocities for the next time step can be obtained as: ⎧ ( ) ˙ n−1/2 + 2∆tD−1 Fn ⎪ (2 − cn ∆t) U ⎨U n+1/2 ˙ = (11) (2 + cn ∆t) ⎪ ⎩Un+1 = Un + ∆t U n+1/2 ˙ where n is the nth iteration step, ∆t is the time step size, c and D are the damping coefficient and fictitious diagonal density matrix, respectively, which can be calculated as shown in [29], and F is the summation of internal and external forces, in which its ith component can be written as: ∑ ⟨ ⟩ ⟨ ⟩ Fi = {T [xi , t] x j − xi − T [x j , t] xi − x j }V j + bi (xi , t) (12) j∈Hi

˙ 0 = 0, the integration can be started by By assuming that U0 ̸= 0 and U ˙ 1/2 = ∆tD−1 F0 /2 U

(13)

H. Zhang and P. Qiao / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112883

5

Fig. 2. Virtual crack extension in peridynamics (The black thick solid line: crack path; the blue dashed line: virtual crack extension; the red thin solid line: applied broken bond.). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

With the application of the ADR, the uniform load applied to the cracked body converges to the stable static solution. In the case of continuously-simulated quasi-static elastic behavior, a non-equilibrium criterion is proposed to predict whether the system arrives at the balance condition: |∆Ue | ≤ βUe

(14)

where Ue is the total strain energy of the whole system, ∆Ue is the variation of Ue between the adjacent iteration steps, β is the convergence standard to describe the allowable ratio of total strain energy, which is chosen as β = 1.0×10−8 in this paper. It can also be written in detail in nth iteration as: ⏐( )⏐ ⏐ ⏐ ∑ ∑ ∑ ⏐ ⏐ n−1 n Win Vi (15) Wi Vi − Wi Vi ⏐ ≤ β ⏐ ⏐ ⏐ i∈H

i∈H

i∈H

where Win is the strain energy density of node i in nth iteration step, and it can be calculated with Eq. (2) for elastic materials. When the non-equilibrium in Eq. (15) is satisfied, the responding deformation of loaded body converges in the stable state and the virtual crack can then be applied. 3.2. Virtual crack in peridynamic model In peridynamic model, the nodes interact with each other with “peridynamic bonds”. The visible crack in body can be formulated by a series of directional broken bonds. The crack can be regarded as the accumulation of broken bonds. Thus, for virtual crack extension, the typical virtual broken bonds are applied to form the virtual crack. As shown in Fig. 2, the system (material) is discretized into nodes, and the nodes interact with each other by bonds. The virtual crack (the blue dashed line in Fig. 2) can be formulated by all the broken bonds across the crack line. When the virtual crack extension da is applied, the bonds across da are all artificially broken, such as the bond connecting nodes i and j (the red thin solid line shown in Fig. 2). Unlike the failure criterion model in [18], because bonds are artificially broken in virtual crack analysis, no bond broken criterion is needed. Thus, it is not a real crack propagation process. Furthermore, different from the virtual crack analysis by the FEM [1], no extra nodes or bonded nodes are created or needed in peridynamic model and the virtual crack is applied naturally and smoothly with the artificially-broken bonds by the proposed PVCE method. 3.3. An algorithm of the PVCE method A computational algorithm for implementing the PVCE method is flowcharted in Fig. 3, and the serial virtual crack extensions are considered to capture the function of energy release rates G (a) of different crack lengths. As the flowchart shown in Fig. 3, the cracked body with the initial crack length a0 is considered and unit crack extension da = ∆h is used in the initial step. The ADR method is utilized for the quasi-static analysis of

6

H. Zhang and P. Qiao / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112883

Fig. 3. Flowchart of the PVCE method for serial virtual crack extensions.

the cracked body under static load. With the ADR time integration, the total strain energy Ue (a) of a quasi-static solution is captured when Ue satisfies the non-equilibrium in Eq. (14). After this stable solution is reached, the unit crack extension N*∆h is serially applied along the pre-crack as shown in Fig. 2. Meanwhile, the ADR method is continuously utilized to obtain the stable solutions for different crack lengths. After the iterations, the total strain energy Ue (a) with different crack lengths of a = a 0 + N*∆h is acquired. This total strain energy Ue (a) is then substituted in Eq. (9) to obtain the function of energy release rate G (a). Typically, when the applied time of maximum virtual crack M = 2 is reached, the strain energies of Ue (a0 ) and Ue (a0 + ∆h) are acquired by the iterations. The single energy release rate G (a0 ) is then obtained by Eq. (9). 3.4. Numerical modification of the PVCE method The energy release rate is related to the differentiation of total strain energy with respect to the crack extension length as given in Eq. (8), and it can be calculated with the total strain energy stored for two or more slightly different crack lengths. The most commonly used numerical calculation way is the original two-step difference method as given in Eq. (9).

7

H. Zhang and P. Qiao / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112883

Under the assumption that the function of the strain energy Ue is continuously differentiable to the crack length a, the strain energy Ue is expanded near a using the first few terms of a Taylor series as (da)2 + O (da)3 2 Combining Eqs. (9) and (16) and employing the unit extension da = ∆h result in Ue (a + da) = Ue (a) + Ue′ (a) da + Ue′′ (a)

G (a) − G (a) = Ue′′ (a)

(16)

(∆h) + O (∆h)2 = O (∆h) 2

(17)

As given in Eq. (17), the relative difference of the numerical energy release rate G (a) to G(a) is O (∆h). To reduce the error of the numerical differentiation, the modified PVCE method for calculating the energy release rate G (a) is defined as: ⎧ 4Ue (a + ∆h) − Ue (a + 2∆h) − 3Ue (a) ⎪ ⎨+ Fixed load 2∆h G (a) = (18) ⎪ ⎩− 4Ue (a + ∆h) − Ue (a + 2∆h) − 3Ue (a) Fixed displacement 2∆h where Ue (a), Ue (a + ∆h) and Ue (a + 2∆h) are the total strain energy with the pre-crack lengths a, a + ∆h and a + 2∆h, respectively. Substituting da = ∆h and da = 2∆h into Eq. (16) and combining them in Eq. (18) lead to G (a) − G (a) = O (∆h)2

(19)

Thus, the relative difference of the modified numerical energy release rate G (a) to G(a) is O ∆h . Compared to the two-step numerical expression ( of) G (a) in Eq. (9), the modified model demonstrated in Eq. (19) efficiently reduces the error from O (∆h) to O ∆h 2 . Because Ue (a), Ue (a + ∆h) and Ue (a + 2∆h) are all needed in Eq. (19), the serial virtual crack extensions should be performed in the modified model. Overall, the numerical energy release rate in state-based peridynamic model can be computed in Eq. (8), in which the total peridynamic strain energies of loaded system with serial crack extensions are needed and the entire boundary problem needs to be solved multiple times. In the peridynamic J-integral method [22–24], the energy release rate can be directly obtained with only the present solved state. In those methods, the functions of displacement derivatives and force interactions in [22] and displacement derivatives in [23] associated with the contour integral path are required to be computed, while an approximate non-local peridynamic stress tensor from peridynamic bond forces are utilized in [24]. However, the proposed method is much simple and easy to be implemented using only the total strain energy changes, and it has fewer restrictions and assumptions. Thus, it can be potentially extended for much more general cases, such as heterogeneous materials and curved crack. (

2

)

4. Validation The single edge-notched tension (SENT) specimens under constant stress and displacement load (Fig. 4) which are called as SENT SL and SENT DL tests, respectively, are considered to validate accuracy of the proposed PVCE method in this study. The elastic isotropic material 1CrMoV steel [30] is utilized in this example, and the material properties are: E = 214 GPa, v = 0.27 and ρ = 7850 kg/m3 . Both the plane stress and plane strain conditions are considered, and the uniform thickness value of 1 mm is used. 4.1. Single edge-notched tension specimens under stress load (SENT_SL) In this section, the single edge-notched tension (SENT) plate under the fixed tensile stress of σ = 1 MPa is first considered (see Fig. 4). From the standard mode I fracture test, the stress intensity factor K I of this typical cracked specimens can be given as [27]: √ K I = σ πa F (b) , b = a/w (20) where F (b) = 1.122 − 0.231b + 10.55b2 − 21.71b3 + 30.38b4

(21)

8

H. Zhang and P. Qiao / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112883

Fig. 4. The single edge-notched tension (SENT) specimen.

As shown in Fig. 4, a and w are the crack length and width of the specimen, respectively. The energy release rate G can then be calculated from K I as: ⎧ K I2 ⎪ ⎪ ⎨ Plane stress E) G= ( (22) ⎪ 1 − v 2 K I2 ⎪ ⎩ Plane strain E Because the expression in Eq. (20) is accurate with 1% over the range of a/w < 0.6, it can be used as the standard (benchmark) values of energy release rate to validate the numerical peridynamic simulation results by the PVCE method. 4.1.1. Behavior of the SENT_SL with a/w = 0.5 The typical crack length a = 0.1 m is first considered, which leads to the value of the a/w = 0.5. Comparisons of the displacement and strain energy density by the peridynamic and FEM models are shown in Figs. 5 and 6, respectively, where the typical mesh of δ = 4 mm and m = 4 are utilized and the plane stress condition is considered. As shown in Figs. 5 and 6, the distributions of displacement and strain energy density by the peridynamic model show great consistence with those by the FEM. By summing the strain energy of cracked body in Fig. 6, the total strain energy of SENT specimen by the peridynamic model is 4.77 × 10−4 J, within 1.1% difference to the FEM result (i.e., 4.72 × 10−4 J). The serial virtual crack extensions are then applied as shown in Fig. 3. The total strain energy Ue and external energy We during application of the serial virtual crack extensions are presented in Fig. 7. As shown in Fig. 7, the stable value of Ue (a) is obtained by the finite iterations using the ADR method, and the unit crack ∆h is then serially applied. The following strain energies of Ue (a + ∆h) and Ue (a + 2∆h) are obtained at their stable solution steps. For the fixed load condition, the numerical energy release rate G (a) based on the modified PVCE method in Eq. (18) is 11.86 J/m2 , within the difference of 0.94% to the FEM prediction (i.e., 11.75 J/m2 ). In addition, in the stable steps, the values of total external energy are twice of the computed strain energy, and the extra portion of external energy is absorbed by the adaptive damping. 4.1.2. Parametric studies and modification The crack extension length ∆h and the convergence by the mesh parameters δ and m are discussed for the PVCE method, and the modified model is then performed. The numerical energy release rates G (a) of the SENT SL with these parameters and modification are calculated and compared to the analytical values G (a) by Eq. (22), which

H. Zhang and P. Qiao / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112883

9

Fig. 5. Comparisons of displacement distributions in x-direction: (a) FEM and (b) Peridynamics, and y-direction: (c) FEM and (d) peridynamics, of the SENT_SL with crack length a = 0.1 m.

Fig. 6. Comparison of strain energy density distributions: (a) FEM and (b) Peridynamics, of the SENT_SL with crack length a = 0.1 m.

10

H. Zhang and P. Qiao / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112883

Fig. 7. The energies during tension simulation of the SENT_SL using the serial virtual crack extensions.

Table 1 The numerical energy release rates (J/m2 ) by the PVCE method with different values of ∆h. Plane stress δ = 4 mm δ = 2 mm

Plane strain

∆h = 1/2∆x

∆h = ∆x

∆h = 2∆x

∆h = 1/2∆x

∆h = ∆x

∆h = 2∆x

13.41 13.22

12.13 11.96

12.39 12.05

12.54 12.36

11.24 11.00

11.49 11.16

are 11.75 J/m2 and 10.89 J/m2 for the plane stress and plane strain cases, respectively. The relative difference e is defined to examine accuracy of the PVCE method as: ⏐ ⏐ ⏐ G (a) − G (a) ⏐ ⏐ ⏐ e=⏐ (23) ⏐ ⏐ ⏐ G (a) The energy release rates of the SENT SL using the PVCE method with different sizes of virtual crack length ∆h are reported in Table 1. As shown in Table 1, the value of ∆h = ∆x is the most effective and accurate size when compared to the analytical values. Thus, the value of ∆h = ∆x is used in the following examples. Then, the convergence studies of δ-convergence and m-convergence [31] are performed with the horizons of δ = 8 mm, 4 mm and 2 mm, and m = 4, 5, 6, respectively. The numerical peridynamic results of energy release rate for different values of δ and m are presented in Table 2, and the relative differences are shown in Fig. 8. Overall, when the values of m are fixed, the predictions based on the PVCE method converge to the analytical values as δ decreases; while the same to m-convergence but with less efficiency. The relative difference reduces within 2% when the value of δ = 4 mm is used. Because the dense mesh size would lead to more expensive computational cost, the suitable values of m and δ should be chosen for different accuracy in the application of the PVCE method. In comparisons with the PVCE method of the original two-step difference in Eq. (9), the modified PVCE method in Eq. (18) is then performed. The numerical energy release rates and relative differences using the original (Eq. (9)) and modified (Eq. (18)) PVCE are presented in Table 3 and Fig. 9, where the fixed value of m = 4 and different values of δ are considered. As shown in Fig. 9, the relative differences of modified PVCE method are much smaller than the original PVCE method, especially for the sparse mesh size. Thus, the modified PVCE method improves the accuracy in computing the strain energy release rate, and it is much effective in converging to the analytical result, which is consistent to the relative differences defined in Eqs. (17) and (19).

H. Zhang and P. Qiao / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112883 Table 2 The numerical energy release rates by the PVCE method with different values of δ and m (J/m2 ). Plane stress

m=4 m=5 m=6

Plane strain

δ = 8 mm

δ = 4 mm

δ = 2 mm

δ = 8 mm

δ = 4 mm

δ = 2 mm

12.56 12.51 12.37

12.13 12.10 12.04

11.96 11.95 11.95

11.66 11.63 11.48

11.24 11.23 11.17

11.00 11.00 11.00

Fig. 8. Relative differences of PVCE results to analytical values in the cases of plane stress (a) and plane strain (b).

Fig. 9. Relative differences of original and modified PVCE predictions in comparison to analytical values. Table 3 The numerical energy release rates by the original and modified PVCE methods (J/m2 ). Plane stress

Original Modified

Plane strain

δ = 8 mm

δ = 4 mm

δ = 2 mm

δ = 8 mm

δ = 4 mm

δ = 2 mm

12.56 12.01

12.13 11.86

11.96 11.86

11.66 11.15

11.24 11.00

11.00 10.84

11

12

H. Zhang and P. Qiao / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112883

Fig. 10. Predictions of energy release rate of the SENT_SL by the modified PVCE method with different values of a/w in the cases of plane stress (a) and plane strain (b).

4.1.3. Behavior of the SENT_SL with different values of a/w The different values of crack length varying from a = 0.04 to 0.12 m are considered, which lead to a/w changing from 0.2 to 0.6. The fixed mesh sizes of δ = 4 mm and m = 4 are used, and only the modified PVCE method is considered. The energy release rate predictions of the SENT SL by the modified PVCE method with different values of a/w are presented in Fig. 10. As shown in Fig. 10, the energy release rate increases with the increasing value of initial crack length. In both the plane stress and plane strain cases, the predicted energy release rate values by the modified PVCE method greatly match those analytical values from Eq. (22), with the maximum relative difference of 3.15%. 4.2. Single edge-notched tension specimens under fixed displacement load (SENT_DL) In a similar fashion, the single edge-notched tension specimen under fixed displacement load (SENT DL) is evaluated in this section. The fixed displacement load is u y = 2.0 × 10−6 m. Considering the parametric study above, the typical mesh sizes of δ = 4 mm and m = 4 and the modified PVCE method in Eq. (18) are considered in the SENT DL analysis. The virtual crack closure technique (VCCT) is performed by the FEM to obtain the standard values of energy release rate. 4.2.1. Behavior of the SENT_DL with a/w = 0.5 The typical crack length a = 0.1 m is also first considered. Comparisons of the displacement and strain energy density by the numerical peridynamic and FEM models are shown in Figs. 11 to 12, respectively, where the plane stress condition is considered. Again, the distributions of displacement and strain energy density by the peridynamic model show great consistence with the FEM results. The total strain energy of specimen is 5.43 × 10−4 J, within 0.4% difference of the FEM result (i.e., 5.45 × 10−4 J). The computed total strain energy during tension simulation of the SENT DL by the modified PVCE method is presented in Fig. 13. As shown in Fig. 13, the stable solution under the fixed displacement load is much easier to be reached than that with the fixed stress load (see Fig. 7). The total strain energy Ue (a) decreases with the increasing crack length a. In addition, the typical numerical energy release rate calculated by Eq. (18) is given in Fig. 14. 4.2.2. Behavior of the SENT_DL with different values of a/w The different values of crack length varying from a = 0.04 to 0.14 m are also considered for the SENT DL specimen, which lead to a/w changing from 0.2 to 0.7. The predicted energy release rate values with different a/w of the SENT DL are presented in Fig. 14. As shown in Fig. 14, the energy release rate increases with the values of a/w when a/w < 0.5, whereas it decreases when a/w > 0.5. In addition, the predicted energy release rate values by the modified PVCE method greatly coincide with those predicted by the VCCT method of FEM with the maximum relative deviation of 1.0%.

H. Zhang and P. Qiao / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112883

13

Fig. 11. Comparisons of displacement distributions in x-direction: (a) FEM and (b) Peridynamics, and y-direction: (c) FEM and (d) peridynamics, of the SENT_DL with crack length a = 0.1 m.

Fig. 12. Comparison of strain energy density distributions: (a) FEM and (b) Peridynamics, of the SENT_DL with crack length a = 0.1 m.

14

H. Zhang and P. Qiao / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112883

Fig. 13. The total strain energy during tension simulation of the SENT_DL using the PVCE method.

Fig. 14. Values of numerical energy release rate of the SENT_DL by the PVCE and VCCT methods with different values of a/w in the cases of plane stress (a) and plane strain (b).

5. Applications The double cantilever beam (DCB) [32] and end loaded split (ELS) [33] tests are considered for applications of the proposed PVCE method for the mode I and II fracture problems, respectively. The critical load and displacement of the specimens are predicted by the modified PVCE method. The critical load Pc and the critical displacement uc can be calculated for each delamination length a as [34]: √ √ Gc Gc Pc = P , uc = u (24) G (a) G (a) where P and u are the applied load and displacement of the specimen, respectively. For the fixed displacement condition of u = u 0 , P is the function of crack length a and u0 (i.e., P = P(a, u0 )). G (a) is the corresponding energy release rate under the given load with crack length a, and Gc is the critical energy release rate of materials that is equal to G I c or G I I c for the DCB or ELS tests, respectively. Thus, for the DCB and ELS specimens with different crack length a, when the corresponding energy release rate G (a) under typical P and u is calculated by the modified PVCE method, the unique critical load Pc and uc are obtained.

H. Zhang and P. Qiao / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112883

15

Fig. 15. Double cantilever beam (DCB) (a) and end loaded split (ELS) specimens (b) under fixed displacement load.

As shown in Fig. 15, the special dimensions of DCB and ELS specimens are given as: L = 240 mm and h = 20 mm. In both the DCB and ELS tests, the elastic isotropic material is considered, and the material properties are E = 200 GPa, v = 0.3, G I c = 9627.8 J/m2 and G I I c = 100 J/m2 . The plane stress condition is considered, and the uniform thickness value of 10 mm is utilized. The fixed displacement load is u y = 1.0×10−4 m. In the peridynamic simulation model, the typical mesh sizes of δ = 1 mm and m = 5 are utilized. The modified PVCE method with the serial crack extensions is performed from the initial crack length a = 100 mm to 200 mm. For each crack length from a = 100 mm to 200 mm, the modified PVCE method is utilized to calculate the energy release rate. The responding load and energy release rate values of DCB and ELS tests are calculated by the modified PVCE method, and their predictions are presented in Figs. 16 and 17, respectively. By substituting those calculated values into Eq. (24), the critical load and displacement are obtained. The critical load Pc and displacement uc of DCB and ELS specimens of different crack lengths are presented in Figs. 18 and 19, where the VCCT of FEM and peridynamic fracture criterion (PFC) [18] are also considered for comparison with and verification of the present modified PVCE method. As shown in Figs. 18 and 19, in both the DCB and ELS tests, the proposed modified PVCE method accurately predicts the critical load and displacement values, within the maximum difference of 2.2% to the VCCT results of FEM. In summary, the proposed modified PVCE method is capable of accurately and effectively calculating the energy release rate for both the mode I and mode II fracture problems and greatly predicting the critical load and displacement of typical fracture specimens. 6. Conclusions In this paper, the virtual crack extension method is implemented to numerically calculate the energy release rates by the state-based peridynamic model, and the new peridynamic virtual crack extension (PVCE) method is proposed. The directional virtual broken bonds are uniquely applied to constitute the virtual crack for crack extension in the proposed peridynamic model. A non-equilibrium criterion is introduced to capture the stable solution of loaded system, and the modified PVCE method is developed to reduce the computational error.

16

H. Zhang and P. Qiao / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112883

Fig. 16. Responding load (a) and energy release rate (b) of DCB specimens under fixed displacement load.

Fig. 17. Responding load (a) and energy release rate (b) of ELS specimens under fixed displacement load.

Fig. 18. Critical load (a) and critical displacement (b) of DCB specimen of different crack lengths.

H. Zhang and P. Qiao / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112883

17

Fig. 19. Critical load (a) and critical displacement (b) of ELS specimen of different crack lengths.

Several fracture mechanics examples using the standard test specimens, i.e., single edge-notched tension (SENT), double cantilever beam (DCB) and end loaded split (ELS) tests, are analyzed using the proposed PVCE method. Overall, the energy release rates of the SENT specimens are accurately captured by the PVCE method under either the constant stress or displacement load condition. The modified PVCE method (see Eq. (18)) is much effective to converge to the analytical result when compared to the originally defined one (Eq. (9)). Meanwhile, in both the DCB and ELS tests, the modified PVCE method successfully predicts the critical load and displacement values, within the maximum difference of 2.2% to the VCCT results by the FEM. In conclusion, the proposed PVCE method is capable of accurately and effectively predicting the energy release rates for both mode I and mode II fracture problems, and it builds a natural “bridge” for fracture analysis between the peridynamics and classical fracture mechanics. With availability of the proposed PVCE method, the fracture criteria in fracture mechanics can be effectively applied in peridynamic theory, as demonstrated in the mode I and mode II fracture analysis using the respective DCB and ELS tests. Because only the total strain energy change during crack extension is needed for energy release rate computation in the proposed PVCE method, only few restrictions and assumptions exist. Thus, it can easily be extended for more general situations, such as orthotropic materials, curved crack, dynamic fracture, etc. Acknowledgment The partial financial support from the National Natural Science Foundation of China (NSFC Grant Nos. 11972224, 51679136 and 51678360) to this study is acknowledged. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]

V.B. Watwood, The finite element method for prediction of crack behavior, Nucl. Eng. Des. 11 (1969) 323–332. T.K. Hellen, On the method of virtual crack extensions, Int. J. Numer. Methods Engrg. 9 (1975) 187–207. S.A. Silling, Reformulation of elasticity theory for discontinuities and long-range forces, J. Mech. Phys. Solids 48 (2000) 175–209. Y.D. Ha, F. Bobaru, Studies of dynamic crack propagation and crack branching with peridynamics, Int. J. Fract. 162 (2010) 229–244. Y.D. Ha, F. Bobaru, Characteristics of dynamic brittle fracture captured with peridynamics, Eng. Fract. Mech. 78 (2011) 1156–1168. M. Ghajari, L. Iannucci, P. Curtis, A peridynamic material model for the analysis of dynamic crack propagation in orthotropic media, Comput. Methods Appl. Mech. Engrg. 276 (2014) 431–452. Y.L. Hu, E. Madenci, Bond-based peridynamic modeling of composite laminates with arbitrary fiber orientation and stacking sequence, Compos. Struct. 153 (2016) 139–175. Y.L. Hu, N.V. De Carvalho, E. Madenci, Peridynamic modeling of delamination growth in composite laminates, Compos. Struct. 132 (2015) 610–620. X.P. Zhou, X.B. Gu, Y.T. Wang, Numerical simulations of propagation, bifurcation and coalescence of cracks in rocks, Int. J. Rock Mech. Min. Sci. 80 (2015) 241–254. Y. Wang, X. Zhou, Y. Shou, The modeling of crack propagation and coalescence in rocks under uniaxial compression using the novel conjugated bond-based peridynamics, Int. J. Mech. Sci. 128–129 (2017) 614–643.

18

H. Zhang and P. Qiao / Computer Methods in Applied Mechanics and Engineering 363 (2020) 112883

[11] D. Huang, G. Lu, P. Qiao, An improved peridynamic approach for quasi-static elastic deformation and brittle fracture analysis, Eng. Fract. Mech. 94–95 (2015) 111–122. [12] H. Zhang, P. Qiao, A state-based peridynamic model for quantitative elastic and fracture analysis of orthotropic materials, Eng. Fract. Mech. 206 (2019) 147–171. [13] H. Zhang, P. Qiao, An extended state-based peridynamic model for damage growth prediction of bimaterial structures under thermomechanical loading, Eng. Fract. Mech. 189 (2018) 81–97. [14] Y. Zhang, P. Qiao, An axisymmetric ordinary state-based peridynamic model for linear elastic solids, Comput. Methods Appl. Mech. Engrg. 341 (2018) 517–550. [15] H. Zhang, P. Qiao, A coupled peridynamic strength and fracture criterion for open-hole failure analysis of plates under tensile load, Eng. Fract. Mech. 204 (2018) 103–118. [16] H. Zhang, P. Qiao, L. Lu, Failure analysis of plates with singular and non-singular stress raisers by a coupled peridynamic model, 158 (2019) 446–456. [17] S.A. Silling, E. Askari, A meshfree method based on the peridynamic model of solid mechanics, Comput. Struct. 83 (2005) 1526–1535. [18] H. Zhang, P. Qiao, A state-based peridynamic model for quantitative fracture analysis, Int. J. Fract. 211 (2018) 217–235. [19] Y. Zhang, P. Qiao, A new bond failure criterion for ordinary state-based peridynamic mode II fracture analysis, Int. J. Fract. 215 (2019) 105–128. [20] J.T. Foster, S.A. Silling, W. Chen, An energy based failure criterion for use with peridynamic states, Int. J. Multiscale Comput. Eng. 9 (2011) 675–687. [21] Y. Wang, X. Zhou, Y. Wang, Y. Shou, A 3-D conjugated bond-pair-based peridynamic formulation for initiation and propagation of cracks in brittle solids, Int. J. Solids Struct. 134 (2018) 89–115. [22] W. Hu, Y.D. Ha, F. Bobaru, S.A. Silling, The formulation and computation of the nonlocal J-integral in bond-based peridynamics, Int. J. Fract. 176 (2012) 195–206. [23] C. Stenström, K. Eriksson, The J-contour integral in peridynamics via displacements, Int. J. Fract. 216 (2019) 173–183. [24] R. Panchadhara, P.A. Gordon, R. Panchadhara, P.A. Gordon, Application of peridynamic stress intensity factors to dynamic fracture initiation and propagation, Int. J. Fract. 201 (2016) 81–96. [25] P. Diehl, F. Franzelin, D. Pflüger, G.C. Ganzenmüller, Bond-based peridynamics: a quantitative study of Mode I crack opening, Int. J. Fract. 201 (2016) 1–14. [26] S.A. Silling, M. Epton, O. Weckner, J. Xu, E. Askari, Peridynamic states and constitutive modeling, J. Elasticity 88 (2007) 151–184. [27] C.T. Sun, Z.H. Jin, Fracture Mechanics, Academic Press, Waltham, 2013. [28] B. Kilic, E. Madenci, An adaptive dynamic relaxation method for quasi-static simulations using the peridynamic theory, Theor. Appl. Fract. Mech. 53 (2010) 194–204. [29] E. Madenci, O. Erkan, Peridynamic Theory and its Applications, Springer, New York, 2014. [30] B.K. Neale, An investigation into the effect of thickness on the fracture behaviour of compact tension specimens, Int. J. Fract. 14 (1978) 203–212. [31] F. Bobaru, M. Yang, L.F. Alves, S.A. Silling, E. Askari, J. Xu, Convergence, adaptive refinement, and scaling in 1D peridynamics, Int. J. Numer. Methods Engrg. 77 (2009) 852–877. [32] ASTM E-399-12, Standard test method for linear elastic plane strain fracture toughness KIC of metallic materials, in: Annual Book of ASTM StandArds, Vol 03.01, American Society for Testing and Materials, West Conshohocken, PA, 2013. [33] H. Wang, Use of end-loaded-split (ELS) test to study stable fracture behaviour of composites under mode II loading, Compos. Struct. 36 (1997) 71–79. [34] K. Ronald, A summary of benchmark examples to assess the performance of quasi-static delamination propagation prediction capabilities in finite element codes, J. Compos. Mater. 49 (2005) 3297–3316.