A continuous/discontinuous deformation analysis (CDDA) method based on deformable blocks for fracture modeling

Abstract

In the framework of finite element meshes, a novel continuous/discontinuous deformation analysis (CDDA) method is proposed in this paper for modeling of crack problems. In the present CDDA, simple polynomial interpolations are defined at the deformable block elements, and a link element is employed to connect the adjacent block elements. The CDDA is particularly suitable for modeling the fracture propagation because the switch from continuous deformation analysis to discontinuous deformation analysis is natural and convenient without additional procedures. The SIFs (stress intensity factors) for various types of cracks, such as kinked cracks or curved cracks, can be easily computed in the CDDA by using the virtual crack extension technique (VCET). Both the formulation and implementation of the VCET in CDDA are simple and straightforward. Numerical examples indicate that the present CDDA can obtain high accuracy in SIF results with simple polynomial interpolations and insensitive to mesh sizes, and can automatically simulate the crack propagation without degrading accuracy.

Keywords

fracture / crack / propagation / deformable block / continuous/discontinuous deformation analysis (CDDA)

Introduction

The accurate analysis of cracked solid are of vital importance for the safety assessment and life prediction of cracked engineering structures and materials [1]. The traditional finite element method (FEM) with embedded-singularity elements by [2,3], and others, has been widely used for fracture modeling with quite good accuracy, because of its simplicity and popularity. However, the method finds difficulties in modeling crack propagation due to the requirement of remeshing during crack propagation. To overcome this difficulty, a wide range of new numerical methods, such as the MM (Meshless Methods) [4-7], the NMM (numerical manifold method) [8,9], and the XFEM (extended finite element method) [10-12], have been proposed. Despite clear general progress with these methods, there are still some technical issues in their application to fracture problems, for instance, the complexity in the construction of the discontinuous interpolation functions along the crack, the expensive cost to refine the nodal arrangement near the crack tip, and the complicated algorithm for the subdividing of integral cells around the crack.

Among the above-mentioned methods for crack simulation, the other difficulty is the accurate calculation of the SIFs of the crack tip. The embedded-singularity elements or the singular quarter-point elements in the FEM work well for static cracks, but it is not easy to be implemented when crack grows. The J-integral or contour integral scheme for the calculation of SIF has been widely welcomed by the engineers because it has the advantages of high accuracy and path independent, and it is suitable for both elastic and plastic analyses. However, due to the complexity in mathematic formulation, it is not an easy task to apply J-integral to modeling crack growth, especially for kinked cracks or curved cracks [13,14]. The VCET (virtual crack extension technique) [15] has the same merits as J-integral, and is suitable for the analyses of various cracks including kinked cracks and curved cracks. However, the requirement of two sets of analysis data and the separating of crack tip node along the virtual crack extension direction makes it is not so convenient to implement. The VCCT (virtual crack closure technique) [16-18], which is derived from the VCET, only requires one step results of FEM analysis, but it will cause some loss of precision and has certain demand of the element sizes around the crack tip. Despite all this, the VCET and VCCT get increased attentions in the field of fracture analysis, due to their simple formulation, high accuracy and insensitivity to mesh sizes.

On the basis of existing research and in the framework of finite element meshes, a continuous/discontinuous analysis (CDA) method is proposed for the modeling of crack propagation. In the present CDA, simple polynomial interpolations are defined at the deformable block elements, and a link element is employed to connect the adjacent block elements.

It should be noted that the proposed method is different from existing formulation of the NMM and DDA. For the NMM, its most appealing advantage is the possibility of a unified approach from continuum analysis to discontinuum analysis. However, for a problem containing large quantities of discontinuities, e.g., faults and joint when modeling rock, the generation of the cover system is quite complicated and also computationally expensive. On the other hand, DDA uses independent degrees of freedom for each block and can be easily applied for analyzing discrete system. DDA can be regarded as a special form of the NMM, however, applying DDA from continuum to discontinuum analyses is numerically cumbersome. The contribution of the present work is the use of thin element previously developed by Goodman and Desai in the context of the FEM. The use of thin element connecting between blocks enables the DDA formulation having both continuum and discontinuum capability. Here we coin the method as Continuous/Discontinuous Analysis (CDA) method. The CDA is particularly suitable for modeling the fracture propagation because the switch from continuous deformation analysis to discontinuous deformation analysis is natural and convenient without additional procedures. The SIFs for various types of cracks, such as kinked cracks or curved cracks, can be easily computed in the CDA by using the VCET, and the formulation and implementation of the VCET in the CDA are very simple and convenient. Numerical examples indicate that the present method can obtain the high accuracy SIF results with simple polynomial interpolations and shows insensitivity to mesh sizes, and can automatically simulate the crack propagation with fairly high accuracy.

A novel continuous/discontinuous method based on deformable block

Stiffness matrix of the block element

For an arbitrary geometry of analysis domain \mathrm{\Omega } as shown in Fig. 1, a triangular finite element mesh is first set up by using the mature mesh generation technique in the FEM. The finite element mesh can also be the arbitrary polygonal elements. If each triangular element is regarded as an independent deformable block, the displacement function in each block element e can be approximated as

(1)$$\begin{gathered} u(x,y)=p^T(x,y)\tilde{a}={\displaystyle \sum _{k=\mathrm{1}}^m{p}_k(x,y)a_k}, \\ v(x,y)=p^T(x,y)\tilde{b}={\displaystyle \sum _{k=\mathrm{1}}^m{p}_k(x,y)b_k,} \end{gathered}$$

where m is the number of terms in the basis, p^T(x,y)={[\mathrm{1},x_c,y_c,x_c{y}_c,\mathrm{...}]}^{}are monomial basis functions, x_c=x-x_{\mathrm{0}},y_c=y-y_{\mathrm{0}},x_{\mathrm{0}},y_{\mathrm{0}} is the center of the element e, \tilde{a}={[a_{\mathrm{1}},a_{\mathrm{2}},a_{\mathrm{3}},a_{\mathrm{4}},\mathrm{...}]}^T and \tilde{b}={[b_{\mathrm{1}},b_{\mathrm{2}},b_{\mathrm{3}},b_{\mathrm{4}},\mathrm{...}]}^T are the interpolation coefficients to be determined. The linear basis p^T(x,y)={[\mathrm{1},x_c,y_c]}^{} is used here as an example to introduce the CDA.

Equation (1) can be rewritten in the form

(2)u=\mathrm{\Phi }a,

where u={\left[\begin{array}{ll}u& v\end{array}\right]}^T, a={\left[\begin{array}{llllll}{a}_{\mathrm{1}}& a_{\mathrm{2}}& a_{\mathrm{3}}& b_{\mathrm{1}}& b_{\mathrm{2}}& b_{\mathrm{3}}\end{array}\right]}^T, and

(3)\mathrm{\Phi }=\left[\begin{array}{llllll}\mathrm{1}& x_c& y_c& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& x_c& y_c\end{array}\right].

From Eq. (2), we get the strain and stress of the block element e as

(4)ϵ=\mathrm{\Lambda }u=\mathrm{\Lambda }\mathrm{\Phi }a=Ba,

(5)\sigma =Dϵ=DBa=Sa,

where the strain matrix B=\mathrm{\Lambda }\mathrm{\Phi }, the stress matrix S=DB,\mathrm{\Lambda } is the differential operator of the plane problem, and

(6)\mathrm{\Lambda }=\left[\begin{array}{ll}\partial /\partial x& \mathrm{0}\\ \mathrm{0}& \partial /\partial y\\ \partial /\partial y& \partial /\partial x\end{array}\right].

D is the elastic matrix, and

(7)D={\displaystyle \frac{E}{\mathrm{1}-{\mu }^{\mathrm{2}}}}\left[\begin{array}{lll}\mathrm{1}& \mu & \mathrm{0}\\ \mu & \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\displaystyle \frac{\mathrm{1}-\mu }{\mathrm{2}}}\end{array}\right]for\; plane\; stress\; case,

(8)D={\displaystyle \frac{E}{(\mathrm{1}+\mu )(\mathrm{1}-\mathrm{2}\mu )}}\left[\begin{array}{lll}\mathrm{1}-\mu & \mu & \mathrm{0}\\ \mu & \mathrm{1}-\mu & \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\displaystyle \frac{\mathrm{1}-\mathrm{2}\mu }{\mathrm{2}}}\end{array}\right]for\; plane\; strain\; case,

where E is the Young’s modulus and \mu is the Poisson’s ration.

Hence, the stiffness matrix of the block element e can be expressed as

(9)k^e={{\iint }_{\mathrm{\Delta }}B}^{T}DBt\mathrm{d}x\mathrm{d}y,

where t is the thickness of the block element e.

Stiffness matrix of the link element

Let’s consider two adjacent block elements e_{\mathrm{1}} and e_{\mathrm{2}} in Fig. 1. From Eq. (2), the displacement functions of the elements e_{\mathrm{1}} and e_{\mathrm{2}} are expressed as

(10)u^{e_{\mathrm{1}}}={\mathrm{\Phi }}^{e_{\mathrm{1}}}{a}^{e_{\mathrm{1}}}

(11)u^{e_{\mathrm{2}}}={\mathrm{\Phi }}^{e_{\mathrm{2}}}{a}^{e_{\mathrm{2}}}

The elements e_{\mathrm{1}} and e_{\mathrm{2}} have their own independent freedoms and independent deformations. Suppose that e_{\mathrm{1}} and e_{\mathrm{2}} are connected by a thin-layer link element L_{ik} with the thickness b as shown in Fig. 1. The thickness b is taken as b_a/1000 in this study, where b_a is the average length of all element sides. For a given point x belonging to the link element L_{ik}, a displacement function at x can by defined by

(12){\overline{u}}^L=\lambda (u^{e_{\mathrm{1}}}-u^{e_{\mathrm{2}}})=\lambda ({\mathrm{\Phi }}^{e_{\mathrm{1}}}{a}^{e_{\mathrm{1}}}-{\mathrm{\Phi }}^{e_{\mathrm{2}}}{a}^{e_{\mathrm{2}}}),

where u^{e_{\mathrm{1}}} and u^{e_{\mathrm{2}}} are defined by Eq. (10) and Eq. (11), {\overline{u}}^L={\left[\begin{array}{ll}\overline{u}& \overline{v}\end{array}\right]}^T is the displacement function defined at the local coordinate of the link element L_{ik}, \lambda is the coordinate transformation matrix, and

(13)\lambda =\left[\begin{array}{ll}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right].

Equation (12) can be rewritten in the form

(14){\overline{u}}^L={\mathrm{\Phi }}^L{a}^L,

where,

(15){\mathrm{\Phi }}^L=\lambda \left[\begin{array}{ll}{\mathrm{\Phi }}^{e_{\mathrm{1}}}& -{\mathrm{\Phi }}^{e_{\mathrm{2}}}\end{array}\right],

(16)a^L=\left\{\begin{array}{l}{a}^{e_{\mathrm{1}}}\\ a^{e_{\mathrm{2}}}\end{array}\right\},

Due to the thickness b of the link element L_{ik} is much less than its side length ik, the element strain of L_{ik} can be approximated as

(17){ϵ}^L=\left\{\begin{array}{l}{\gamma }_s\\ {ϵ}_n\end{array}\right\}={\displaystyle \frac{{\overline{u}}^L}{b}}={\displaystyle \frac{{\mathrm{\Phi }}^L}{b}}{a}^L=B^L{a}^L,

where B^L={\mathrm{\Phi }}^L/b is the strain matrix. A similar formulation as Eq. (17) was first used by Zhu [19] to construct the weak Interlayer element in rock mechanics. But here it is used to connect the two separate block elements and to develop a new type of CDA method.

Accordingly, the element stress of L_{ik} is expressed as

(18){\sigma }^L=\left\{\begin{array}{l}{\tau }_s\\ {\sigma }_n\end{array}\right\}=D^L{ϵ}^L=S^L{a}^L,

where S^L=D^L{B}^L is the stress matrix, D^L is the elasticity matrix, and

(19)D^L=\left[\begin{array}{ll}G& \mathrm{0}\\ \mathrm{0}& E_{\mathrm{0}}\end{array}\right].

where the shear modulus G=E/\mathrm{2}(\mathrm{1}+v),E_{\mathrm{0}}=Efor plane stress case, and E_{\mathrm{0}}=E/(\mathrm{1}-v^{\mathrm{2}})for plane strain case.

Hence, the stiffness matrix of the link element L_{ik} can be expressed as

(20)k^L={\displaystyle {\int }_{L_{ik}}{(B^L)}^T{D}^L{B}^{L}tb\mathrm{d}\overline{x}}.

The summation of all stiffness matrices of block elements e and link elements L leads to the following discrete system equation

(21)KU=F,

where U is the freedom vector of the block elements to be solved, F is the right hand side load vector, and

(22)F={\displaystyle \sum _e{F}^e},

(23)K={\displaystyle \sum _e{k}^e}+{\displaystyle \sum _L{k}^L}.

Calculation of SIFs by the VCET

For the problem containing a crack with a tip m in Fig. 1, by using Eq. (21), the initial discrete system equation can be expressed as

(24)K_{\mathrm{0}}{U}_{\mathrm{0}}=F_{\mathrm{0}}.

We can get the freedom vector U_{\mathrm{0}} of all block elements by solving Eq. (24). Then, by using Eq. (18), we can obtain the element stress {\sigma }^{L_{mn}}={\left[\begin{array}{ll}{\tau }_s& {\sigma }_n\end{array}\right]}^T of the link element L_{mn} at the crack tip, which is shown in Fig. 2.

Suppose that the crack virtually propagates from point m to point n (Fig. 2). The discrete system equation after the crack virtual propagation can be easily obtained in the CDA as

(25)K_{\mathrm{1}}{U}_{\mathrm{1}}=F_{\mathrm{0}},

where K_{\mathrm{1}} is the global stiffness matrix after the crack virtual propagation, and

(26)K_{\mathrm{1}}=K_{\mathrm{0}}-k^{L_{mn}}.

It means that K_{\mathrm{1}} is a simple subtraction of the element stiffness matrix k^{L_{mn}} of the link element L_{mn} from the initial discrete system equation K_{\mathrm{0}}.

Then we can get the freedom vector U_{\mathrm{1}} of the block elements by solving Eq. (25). From Eq. (14), we can obtain the element displacement {\overline{u}}^{L_{mn}}={\left[\begin{array}{ll}\overline{u}& \overline{v}\end{array}\right]}^T after the crack virtual propagation.

Taken the length of the link element L_{mn} as \mathrm{\Delta }a, the energy release rate at the crack tip m can be approximately calculated by {\overline{u}}^{L_{mn}} and {\sigma }^{L_{mn}}, which is

(27)G_{\mathrm{1}}\cong {\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta }a}}{\displaystyle {\int }_{\mathrm{0}}^{\mathrm{\Delta }a}{\sigma }_n\overline{v}\mathrm{d}\overline{x}},

(28)G_{\mathrm{2}}\cong {\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta }a}}{\displaystyle {\int }_{\mathrm{0}}^{\mathrm{\Delta }a}{\tau }_s\overline{u}\mathrm{d}\overline{x}},

in which G_{\mathrm{1}} is the energy release rate of crack mode I, G_{\mathrm{2}} is the energy release rate of crack mode II. Then the corresponding SIFs K_I and K_{II} can be obtained from the relation between the energy release rate and the stress intensity factor [20].

It is seen that, the formulation and implementation of the VCET in the CDA are much more simple and convenient than they are in the commonly used J-integral method. The present formulation is also suitable for the analyses of various cracks including kinked cracks and curved cracks. The subsequent numerical examples in Section 5 indicate that the present method can obtain the high accuracy SIF results with simple polynomial interpolations and shows insensitivity to mesh sizes.

Simulation of the crack propagation

After getting the SIFs K_I and K_{II} of crack tip m by using Eq. (27) and Eq. (28), we can compute the equivalent SIF of the mixed crack and the crack propagation direction with the maximum hoop stress criterion. The formulation for the calculation of cracking angle {\theta }_0 is

(29)K_I\sin {\theta }_{\mathrm{0}}+K_{II}(\mathrm{3}\cos {\theta }_{\mathrm{0}}-\mathrm{1})=\mathrm{0.}

The crack propagation direction \overrightarrow{l} at the crack tip m can then be obtained via the cracking angle {\theta }_0 in Eq. (29), which is shown in Fig. 2(a). The automatic simulation of the crack propagation in the CDA can be implemented according to the following steps:

1) Divide all the block elements along the crack propagation direction \overrightarrow{l}. For example, as shown in Fig. 2, the block element \mathrm{\Delta }m-n-q is divided into the block element \mathrm{\Delta }m-n-m_{\mathrm{1}} and the block element \mathrm{\Delta }m-m_{\mathrm{1}}-q, the block element \mathrm{\Delta }n-p-q is divided into the block element \mathrm{\Delta }n-m_{\mathrm{2}}-m_{\mathrm{1}} and the block element ▱p-q-m_{\mathrm{1}}-m_{\mathrm{2}}.

2) Add new element freedoms for the new added block elements, and recalculate the element stiffness matrices of the new added block elements and the new added link elements.

3) Repeat the above steps until the crack stops propagation or the crack tip reaches the boundary.

As can be seen, the local geometry algorithm in the CDA is some similar as it is in the XFEM or the NMM, but the adding/removing of the freedoms and the construction of the discontinuity functions along the crack in the CDA are much more simple than they are in the XFEM or the NMM.

The second method to automatically simulate the crack propagation in the CDA is shown in Fig. 3. In general, the crack propagation direction \overrightarrow{l} from Eq. (29) doesn’t walk along one of the link elements around the crack tip m. To accurately simulate the crack propagation, we can employ the following algorithm:

1) Find out the link element L_{mn} which has the minimal clamp angle between the crack propagation direction \overrightarrow{l} and the link elements connecting to the crack tip m, as shown in Fig. 3(a).

2) Calculate the cross point n^{\prime } of the side n-q and direction \overrightarrow{l}(Fig. 3(a)).

3)Move point n to the cross point n^{\prime } as shown in Fig. 3(b).

4)For the purpose of the implementation of the virtual crack extension of the next step, follow the same steps of 1) to 3) and move point p to the cross point p^{\prime } as shown in Fig. 3(c).

5) Recalculate the element stiffness matrices of the block elements and the link elements which are related to the moving points.

6) Repeat the above steps until the crack stops propagation or the crack tip reaches the boundary.

These two methods for simulating the crack propagation in Figs. 2 and 3 have their own advantages and disadvantages. The method in Fig. 2 is more stable and robust, but the geometry algorithm and the implementation of the second method in Fig. 3 are much more simple and convenient than they are in Fig. 2.

Numerical examples

Central-cracked plate

A plate with a central crack of length 2a, as shown in Fig. 4 is first tested by the present method. The dimensions of the plate used in the test are L = 2 m and W = 1 m. The plate is subjected to uniform traction of s = 3 MPa in y direction. The elastic material parameters used are E = 3.0 × 10⁴ MPa and u = 0.3. The problem is solved under plane stress assumption. A number of tests have been performed by varying a between 0.2W, 0.4W and 0.6W. The problem is modeled using 402 and 1262 irregular nodes and solved for plane stress case. The analytical solution for K_I is available in Ref. [21] as

(30)K_I=\sigma \sqrt{\mathrm{\pi }a}\sqrt{\sec {\displaystyle \frac{\mathrm{\pi }a}{\mathrm{2}W}}}\left[\mathrm{1}-\mathrm{0.025}{\left({\displaystyle \frac{a}{W}}\right)}^{\mathrm{2}}+\mathrm{0.06}{\left({\displaystyle \frac{a}{W}}\right)}^{\mathrm{4}}\right].

The computational results of the SIFs by the present CDA are listed in Table 1. The relative errors for the SIFs in Table 1 (also the following tables) are defined as r_{ϵ}=(k^{num}-k^{ana})/k^{ana}, where k^{num} is the numerical solution of SIF, and k^{ana} is the analytical solution or reference solution of SIF. Table 1 indicates that the CDA solution is in a good agreement with the analytical solution even for sparsely distributed nodes.

Three-point bend specimen

The problem of a stationary crack in a three-point bend specimen is considered. The geometry is shown in Fig. 5. The crack is located at the midspan of the beam so that only mode I cracking develops. The dimensions of the specimen are S=\mathrm{12} and W=\mathrm{6}. The load is F=\mathrm{1} applied over unit length and unit depth. The discretisation with 1148 nodes is also shown in Fig. 5 and the computed SIFs are compared with the analytical solution obtained by John [22] in Table 2. As can be seen from the table, the CDA shows high solution accuracy for SIFs with the maximum error less than 0.74%.

Curved crack in an infinite plate

The problem of a curved crack in an infinite plate is considered. The dimensions of the model are shown in Fig. 6, with R=\mathrm{4.25} and \beta =\mathrm{28.0725}°. The analytical solutions for SIFs of the problem are available in Refs [23,24].

The discrete model with 2083 nodes is shown in Fig. 7. The comparison of the SIF results from the present CDA and the analytical solutions is listed in Table 3. It is seen that the CDA shows good accuracy for this problem.

Branched crack in an infinite plate

A branched crack in an infinite plate under uniform tension is examined as shown in Fig. 8. The dimensions of the model are fixed to H=\mathrm{4},W=\mathrm{5},a/W=\mathrm{0.05},b/a=\mathrm{0.9} and \theta =\mathrm{45}°. The material constants are the Young’s modulus E=\mathrm{1000} and the Poisson’s ration v=\mathrm{0.3}. 3688 discrete nodes were used for discretisation as shown in Fig. 9. The normalized stress intensity factors for tips A and B are defined as

(31)F_I^A=K_I^A/\sigma \sqrt{\mathrm{\pi }c},F_I^B=K_I^B/\sigma \sqrt{\mathrm{\pi }c},F_{II}^B=K_{II}^B/\sigma \sqrt{\mathrm{\pi }c},

where c=(a+b\cos \theta )/\mathrm{2}.

The SIF results from the present CDA and the reference solution [25] are listed in Table 4. It is seen that the relative error is less than 1.82% for all the cases.

Curvilinear crack propagation in PMMA beam with circle holes

The CDA is applied to model the propagation of crack from initial notch in polymethylmethacrylate (PMMA) beam with circle holes (Fig. 10). The laboratory experimental result for the PMMA beam and its digitized photograph of observed crack trajectory can be found in Refs [26,27]. Please refer to Refs [26,27]. for more details of the PMMA beam experiment. It is seen that the present CDA can accurately simulate the curvilinear crack propagation in PMMA beam for which experimental result is available (Figs. 11 and 12). The numerical studies for the crack growth of the PMMA beam also indicate that the crack simulating algorithm in Fig. 3 is a little more efficient than the algorithm in Fig. 2 because the algorithm in Fig. 3 avoids the subdividing of the existing block and the dealing with the increasing freedom of new added block.

Conclusions

1) It is seen that, if we employ the definition of link element stiffness matrices k^L in Eq. (20), the present CDA has the ability to do the continuous analysis like the most commonly used FEM. But if we replace the link element stiffness matrices k^L with the contact stiffness k_s and k_n at the normal and tangent directions, the CDA becomes the general DDA (Discontinuous Deformation Analysis) method. Different from other similar methods, the present CDA can switch from continuous deformation analysis to discontinuous deformation analysis naturally and conveniently by simply adding or removing the link elements.

2) Although the triangular block elements are used to introduce the implementation of the present CDA in the paper, arbitrary shape elements can actually be used in the CDA when it is required.

3) The stress intensity factors for various types of cracks, such as kinked cracks or curved cracks, can be easily computed in the CDA by using the VCET, and the formulation and implementation of the VCET in the CDA are very simple and convenient. Numerical examples indicate that the present method can obtain the high accuracy SIF results with simple polynomial interpolations and shows insensitivity to mesh sizes, and can automatically simulate the crack propagation with fairly high accuracy.

4) Only two-dimensional elastic fracture problems are investigated and discussed in this paper. Extensions of the present work to three dimensional analysis or elasto-plastic analysis are natural and feasible. They will also be the further work of the present authors.

References

[1] Atluri S N. Methods of Computer Modeling in Engineering and the Sciences. Tech Science Press, 2005

[2] Tong P, Pian T H H, Lasry S J. A hybrid-element approach to crack problems in plane elasticity. International Journal for Numerical Methods in Engineering, 1973, 7(3): 297-308

[3] Atluri S N, Kobayashi A S, Nakagaki M. An assumed displacement hybrid finite element model for linear fracture mechanics. International Journal of Fracture, 1975, 11(2): 257-271

[4] Xu Y, Saigal S. An element free Galerkin formulation for stable crack growth in an elastic solid. Computer Methods in Applied Mechanics and Engineering, 1998, 154(3-4): 331-343

[5] Belytschko T, Fleming M. Smoothing, enrichment and contact in the element-free Galerkin method. Computers & Structures, 1999, 71(2): 173-195

[6] Ching H K, Batra R C. Determination of crack tip fields in linear elastostatics by the Meshless Local Petrov-Galerkin (MLPG) method. Computer Modeling in Engineering & Sciences, 2001, 2(2): 273-290

[7] Gu Y T, Wang W, Zhang L C, Feng X Q. An enriched radial point interpolation method (e-RPIM) for analysis of crack tip fields. Engineering Fracture Mechanics, 2011, 78(1): 175-190

[8] Ma G W, An X M, Zhang H H, Li L X. Modeling complex crack problems using the numerical manifold method. International Journal of Fracture, 2009, 156(1): 21-35

[9] Wu Z J, Wong L N Y. Frictional crack initiation and propagation analysis using the numerical manifold method. Computers and Geotechnics, 2012, 39: 38-53

[10] Moes N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 1999, 46(1): 131-150

[11] Sukumar N, Chopp D L, Moes N, Belytschko T. Modeling holes and inclusions by level sets in the extended finite-element method. Computer Methods in Applied Mechanics and Engineering, 2001, 190(46-47): 6183-6200

[12] Sukumar N, Chopp D L, Béchet E, Moës N. Three-dimensional non-planar crack growth by a coupled extended finite element and fast marching method. International Journal for Numerical Methods in Engineering, 2008, 76(5): 727-748

[13] Fleming M, Chu Y A, Moran B, Belytschko T. Enriched element-free Galerkin methods for crack tip fields. International Journal for Numerical Methods in Engineering, 1997, 40(8): 1483-1504

[14] Gosz M, Dolbow J, Moran B. Domain integral formulation for stress intensity factor computation along curved three-dimensional interface cracks. International Journal of Solids and Structures, 1998, 35(15): 1763-1783

[15] Hellen T K. On the method of virtual crack extensions. International Journal for Numerical Methods in Engineering, 1975, 9(1): 187-207

[16] Rybicki E F, Kanninen M F. A finite element calculation of stress intensity factors by a modified crack closure integral. Engineering Fracture Mechanics, 1977, 9(4): 931-938

[17] Raju I S. Calculation of strain-energy release rates with high-order and singular finite-elements. Engineering Fracture Mechanics, 1987, 28(3): 251-274

[18] Xie D, Waas A M, Shahwan K W, Schroeder J A, Boeman R G. Computation of Energy Release Rates for Kinking Cracks based on Virtual Crack Closure Technique. CMES: Computer Modeling in Engineering & Sciences, 2004, 6(6): 515-524

[19] Zhu B F. The Finite Element Method Theory and Application. Beijing: China Water Conservancy and Hydropower Press, 1998

[20] Irwin G R. One set of fast crack propagation in high strength steel and aluminum alloys. Sagamore Resrach Conference Proceedings, 1956, 2: 289-305

[21] Anderson T L. Fracture Mechanics: Foundamentals and Application. 2nd ed. Boca Raton: CRC, 1995

[22] John E S. Wide range stress intensity factor expressions for ASTM E399 standard fracture toughness specimens. International Journal of Fracture, 1976, 12: 475-476

[23] Gdoutos E. Fracture Mechanics. Boston: Kluver Academics Publisher, 1979

[24] Budyn E R L. Mutiple Crack Growth by the Extended Finite Element Method. Northwestern University, 2004

[25] Chen Y, Hasebe N. New integration scheme for the branch crack problem. Engineering Fracture Mechanics, 1995, 52(5): 791-801

[26] Bittencourt T N, Wawrzynek P A, Ingraffea A R, Sousa J L. Quasiautomatic simulation of crack propagation for 2d lefm problems. Engineering Fracture Mechanics, 1996, 55(2): 321-334

[27] Ingraffea A R, Grigoriu M. Probabilistic fracture mechanics: A validation of predictive capability, Report 90-8. Department of Structural Engineering, Cornell University, 1990.

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg