# Two-Dimensional Fracture Mechanics Analysis Using a Single-Domain Boundary Element Method

Authors: Chien-Chung Ke; Cheng-Lung Kuo; Shih-Meng Hsu; Shang-Chia Liu; Chao-Shi Chen
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

This work calculates the stress intensity factors (SIFs) at the crack tips, predicts the crack initiation angles, and simulates the crack propagation path in the two-dimensional cracked anisotropic materials using the single-domain boundary element method (BEM) combined with maximum circumferential stress criterion.
The BEM formulation, based on the relative displacements of the crack tip, is used to determine the mixed-mode SIFs and simulate the crack propagation behavior. Numerical examples of the application of the formulation for different crack inclination angles, crack lengths, degree of material anisotropy, and crack types are presented. Furthermore, the propagation path in Cracked Straight Through Brazilian Disc (CSTBD) specimen is numerically predicted and the results of numerical and experimental data compared with the actual laboratory observations. Good agreement is found between the two approaches. The proposed BEM formulation is therefore suitable to simulate the process of crack propagation. Additionally, the anisotropic rock slope failure initiated by the tensile crack can also be analyzed by the proposed crack propagation simulation technique.


---

## Body

## 1. Introduction

Fracture mechanics theory has been developed to solve many geotechnical engineering problems, such as blasting, hydraulic fracturing, and slope stability. In two-dimensional fracture mechanics problems, SIFs are important parameters in analysis of cracked materials. The singularity of stresses near the crack tip is the challenge to numerical modelling methods, even to the BEM. Because the coincidence of the crack surfaces gives rise to a singular system of algebraic equations, the solution of cracked problem cannot be obtained with the direct formulation of the BEM. Several special methods within the scope of the BEM have been suggested for handling stress singularities, such as the Green’s function method [1], the subregional method [2–4], and the displacement discontinuity method (DDM) [5–7].The Green’s function method overcomes the crack modelling problem without considering any source point along the crack boundaries. This method has the advantage of avoiding crack surface modelling and gives excellent accuracy; it is, however, restricted to very simple crack geometries for which analytical Green’s function is available. The sub-regional method has the advantage of modelling cracks with any geometric shape. The method has the disadvantage of introducing the artificial boundaries of the original region into several subregions, thus resulting in a large system of equations. In crack propagation analysis, these artificial boundaries must be repeatedly introduced for each increment of the crack extension. Therefore, this method cannot be easily implemented as an automatic procedure in an incremental analysis of crack extension problems. The DDM overcomes the crack modelling by replacing each pair of coincident source points on crack boundaries by a single source point [8]. Instead of using the Green’s stresses and displacements from point forces, the DDM uses Green’s functions corresponding to point dislocations, that is, displacement discontinuities. This method is quite suitable for crack problems in infinite domains where there is no crack boundary. However, it alone may not be efficient for finite domain problems since the kernel functions in the DDM involve singularities with order higher than those in the traditional displacement BEM. Hence, this method is not suitable for problems involving finite domains.With the advances in single-domain BEM in recent decades, it involves two sets of boundary integral equations for the study of cracked media [9–21]. The single-domain analysis can eliminate remeshing problems, which are typical of the FEM and the subregional BEM. The single-domain BEM has received considerable attention and has been found to be a proper method for simulating crack propagation processes.The single-domain BEM formulation can be achieved by applying the displacement integral equation to the no-crack boundary only, and the traction integral equation on one side of the crack surface only. Since only one side of the crack surface is collocated, one needs to choose either the relative crack opening displacement (COD). This BEM formulation can be applied to the general fracture mechanics analysis in anisotropic media while keeping the single-domain merit.In this study, the BEM formulation combined with the maximum circumferential stress criterion is adopted to predict the crack initiation angles and to simulate the crack propagation paths. Crack propagation in an anisotropic homogeneous plate under mixed-mode I-II loading is simulated by an incremental crack growth with a piece-wise linear discretization. A new computer program, which can automatically generate a new mesh required for analyzing the changing boundary configuration sequentially, has been developed to simulate the fracture propagation process. To demonstrate the proposed BEM procedure for predicting crack propagation in anisotropic materials, the propagation path in a CSTBD is numerically predicted and compared with the actual laboratory observations.A geotechnical engineering problem, slope stability, was analyzed here by the proposed crack propagation simulation technique. Slope stability analysis in the past relied heavily on the use of the limit equilibrium method. An analysis is often carried out by first assuming a failure mode, then using an extensive search to determine the location of critical failure surface. From the type of analysis, estimates of the normal and shear stress distribution on the failure surface and a factor of safety are obtained. This type of study has proved very effective for most engineering applications. However, numerical models, such as continuum models and discontinuum models, tend to be general purpose in nature; that is, they are capable of solving a wide variety of problems. While it is often desirable to have a general-purpose tool available, it requires that each problem be constructed individually. The zones must be arranged by the user to fit the limits of geomechanical units and/or the slope geometry. Hence, numerical models often require more time to set up and run than special-purpose tools such as limit equilibrium method.In the past, the two most popular techniques in continuum mechanics, namely, finite element method (FEM) and finite difference method (FDM), are used in the analysis of the slope stability problems. Early numerical analyses of rock slopes were predominantly undertaken using continuum finite element codes. Kalkani and Piteau (1976) [22], for example, used this method to analyze toppling of rock slopes at Hells Gate in British Columbia, Canada, and Krahn and Morgenstern (1976) [23] undertook preliminary finite element modeling of the Frank Slide in Albetra, Canada. More recently, Stacey et al. in 2003 [24] used finite element analysis in an innovative analysis of extensile strain distributions associated with deep open pit mines. The use of finite difference codes has predominantly involved the use of the FLAC 2D and 3D codes [25]. Coggan et al. (2000) [26] demonstrated the use of both 2D and 3D finite difference analyses in the back analysis of highly kaolinised china clay slopes [27]. However, the tensile crack propagation simulations are scarce in these analyses. This study is interested in the modeling of the fracture propagation path resulting from an unstable slope. The approach is based on fracture mechanics in that fracture propagation is permitted only from the tips of existing cracks. An anisotropic rock slope with a tensile crack is modelled. Failure is triggered by a rise of water pressure within the previously existing crack. At present time, only the tensile failure is considered for the fracture propagation under mixed mode constraints.


## 2.Methodology

### 2.1.  Basic Equations for Anisotropic Elasticity

For the linear elastic, homogeneous, and anisotropic material, the stress and displacement fields can be formulated in terms of two analytical functions,ϕk(zk), of the complex variables zk=x+μky(k=1,2), where μk are the roots of the following characteristic equation [28]:
(2.1)a11μ4-2a16μ3+(2a12+a66)μ2-2a26μ+a22=0,
where the coefficients aij are the compliance components calculated in the x-y coordinate system. The detailed relationship of these components with the material elasticity can be found in Chen et al. (1998) [28]. If the roots μj of  (2.1) are assumed to be distinct, the general expression for the stress and displacements is [28, 29]. One has
(2.2)σx=2Re[μ12ϕ1′(z1)+μ22ϕ2′(z2)],σy=2Re[ϕ1′(z1)+ϕ2′(z2)],τxy=-2Re[μ1ϕ1′(z1)+μ2ϕ2′(z2)],u=2Re[A11ϕ1(z1)+A12ϕ2(z2)],ν=2Re[A21ϕ1(z1)+A22ϕ2(z2)],
where
(2.3)A1j=a11μj2+a12-a16μj,A2j=a12μj+a22μj-a26(j=1,2).
The traction components in the x and y directions are(2.4)Tx=2Re[μ1ϕ1(z1)+μ2ϕ2(z2)],Ty=-2Re[ϕ1(z1)+ϕ2(z2)].With the complex analytical functionsϕi(zi), one can, in general, express (2.2) and (2.4) as follows [28, 30, 31](2.5)ui=2Re[∑j=12Aijϕj(zj)],Ti=-2Re[∑j=12Bijϕj(zj)],σ2i=2Re[∑j=12Bijϕj′(zj)],σ1i=-2Re[∑j=12Bijμjϕj′(zj)](i=1,2),
where zj=x+μjy, the complex number μj, and the elements of the complex matrices A are defined in (2.3), and matrices B can be defined as(2.6)Bij=[-μ1-μ211].Considering the concentrated forces acting at the source point(x0,y0), the analytic functions (ϕk) with the complex variables can be expressed as [30](2.7)ϕk(zk)=-12π[Dklf1ln⁡(zk-zk0)+Dk2f2ln⁡(zk-zk0)],
where zk0=x0+μky0 and fk are the magnitudes of the point force in the k-direction, and(2.8)Dkl=U-1(V-1+V̅-1)-1,U=[A11A12A21A22],V=iUW-1,W=[-μ1-μ211],
where i=-1, overbar means the complex conjugate, and superscript −1 means matrix inverse.Green’s functions of the tractionsTij and displacements Uij can be obtained by substituting (2.7) into (2.5). Their complete expressions are as follows [32]:(2.9)Tij(zk,zk0)=2Re[Qi1(μ1nx-ny)Rj1/(z1-z10)+Qi2(μ2nx-ny)Rj2/(z2-z20)],Uij(zk,zk0)=2Re[Pi1Rj1ln⁡(z1-z10)+Pi2Rj2ln⁡(z2-z20)](i,j=1,2).
In (2.9),  nx and ny are the outward normal components of the field points, and(2.10)Qij=-W=[μ1μ2-1-1].
The complex coefficients Rjk are obtained from the requirements of unit loads at zk0 and from the displacement continuity for the fundamental solution. They are the solutions of the following equation:(2.11)[1-11-1μ1-μ¯1μ2-μ¯2A11-A¯11A12-A¯12A21-A¯21A22-A¯22]⋅[Rj1R¯j1Rj2R¯j2]=[δj2(2πi)-δj1(2πi)00],
where δjk is the Kronecker delta.


### 2.2.  Single-Domain Boundary Integral Equations

A single-domain boundary element method (BEM), based on the relative displacements at the crack tip, is used to determine the mixed-mode SIFs of anisotropic materials. The single-domain BEM formulation consists of the following displacement and traction integral equations (see Figure1).(1)
Displacement integral equation: We have(2.12)Cij(zk,B0)uj(zk,B0)+∫ΓBTij(zk,B,zk,B0)uj(zk,B)dΓ(zk,B)+∫ΓCTij(zk,C,zk,B0)[uj(zk,C+)-uj(zk,C-)]dΓ(zk,C)=∫ΓBUij(zk,B,zk,B0)tj(zk,B)dΓ(zk,B),
where zk and zk0 are the field and source points on the boundary. Here, the subscripts B and C denote the outer boundary and the crack surface, respectively. Cij are coefficients that depend only on the local geometry of the uncracked boundary at zk,B0; Tij and Uij are the Green’s traction and displacement given in (2.9); uj and tj are the boundary displacement and traction; ΓC has the same outward normal as ΓC+.(2)
Traction integral equation: One has(2.13)0.5tj(zk,C0)+nm(zk,C0)∫ΓBClmikTij,k(zk,B,zk,C0)uj(zk,B)dΓ(zk,B)+nm(zk,C0)∫ΓCClmikTij,k(zk,C,zk,C0)[uj(zk,C+)-uj(zk,C-)]dΓ(zk,C)=nm(zk,C0)∫ΓBClmikUij,k(zk,B,zk,C0)tj(zk,B)dΓ(zk,B),
where Clmik is the fourth-order stiffness tensor for anisotropic medium; nm is the unit outward normal to the contour path; the gradient tensors Tij,k and Uij,k denote differentiation with respect to zk0.Figure 1
Geometry of a two-dimensional cracked domain.The Cauchy singularity in (2.12) can be avoided by the rigid-body motion method. The integrand on the right-hand side of (2.12) has only integrable singularity, which can be resolved by the bicubic transformation method [33]. The hypersingularity in (2.13) is resolved by the Gauss quadrature formula, which is very similar to the traditional weighted Gauss quadrature but with a different weight. Therefore, (2.12) and (2.13) can be discretized and solved numerically for the unknown displacements (or the relative crack opening displacements (CODs) on the crack surface) and tractions. In the following section, we present an approach for the evaluation of mixed-mode SIFs.


### 2.3.  Evaluation of Mixed-Mode SIFs

The mixed-mode SIFs for anisotropic materials can be determined by using the extrapolation method of the relative COD, combined with a set of the shape functions. The relative COD is defined as [18](2.14)Δui=∑k=13ϕkΔuik,
where the subscript i(=1,2) denotes the components of the relative COD, and the superscript k(=1,2,3) denotes the relative COD at nodes s=-2/3,0,2/3, respectively; ϕk are the shape functions which can be expressed as follows:(2.15)ϕ1=338s+1[5-8(s+1)+3(s+1)2],ϕ2=14s+1[-5+18(s+1)-9(s+1)2],ϕ3=3385s+1[1-4(s+1)+3(s+1)2].It is well known that for a crack in a homogeneous and anisotropic solid, the relative COD is proportional tor, where r is the distance behined the crack tip. Therefore, the relation of the relative CODs and the SIFs can be found as [18, 33](2.16)Δu1=22rπ(H11KI+H12KII),Δu2=22rπ(H21KI+H22KII),
where(2.17)H11=Im⁡(μ2A11-μ1A12μ1-μ2),H12=Im⁡(A11-A12μ1-μ2),H21=Im⁡(μ2A21-μ1A22μ1-μ2),H22=Im⁡(A21-A22μ1-μ2).Substituting (2.14) into (2.16), a set of algebraic equations is obtained, and then the SIFs KI and KII can be solved. A sign convention for the corresponding SIFs is shown in Figure 2. Using the relative COD, the sign of SIFs can then be determined.Figure 2
Sign convention used for determining SIFs in mode I and mode II.


### 2.4.  Particular Solutions of Gravity and Far-Field Stresses

As mentioned above, if the particular solutions corresponding to the body force of gravity and far-field stresses can be derived in exact closed form, the single-domain BEM formulation presented in this study can then be applied to solve the body force and far-field problems. For the gravity force, the exact close-form solutions can be obtained in a similar way as for the corresponding half-space case [17]. Assuming that the gravity has the components gx and gy, respectively, in the x- and y-directions, the particular solution for the displacements can be found as [34](2.18)uxp=a1ρgxx2+b1ρgyy2,uyp=a2ρgxx2+b2ρgyy2,
where the coefficients ai and bi depend on the elastic stiffness and their expressions are(2.19)a1=0.5(c55c66-c56c56)Δ1,a2=0.5(c15c56-c55c16)Δ1,b1=0.5(c24c46-c26c44)Δ2,b2=0.5(c44c66-c46c46)Δ2,
where(2.20)Δ1=c11c55c66+2c16c15c56-c55c162-c11c562-c66c152,Δ2=c22c44c66+2c26c24c46-c44c262-c22c462-c66c242,
and cij are the elastic stiffness coefficients.Similarly, the particular stresses can be expressed as(2.21)[σ11pσ22pσ23pσ13pσ12p]=[d11d12d21d22d41d42d51d52d61d62][ρgxxρgyy].
Again, dij depend on the elastic stiffness coefficients, and their expressions are(2.22)d11=1,d21=2(a1c12+a2c26),d41=2(a1c14+a2c46),d51=0,d61=0,d12=2(b2c12+b1c16),d22=1,d42=0,d52=2(b2c25+b1c56),d62=0.


### 2.5.  Crack Initiation and Propagation

Three criteria are commonly utilized to predict the crack initiation angle in fracture mechanics problems: the maximum circumferential stress criterion, orσ-criterion [35]; the maximum energy release rate criterion, or G-criterion [36]; the minimum strain energy density criterion, or S-criterion [37]. Among them, the σ-criterion has been found to predict well the directions of crack initiation compared to the experimental results for polymethylmethacrylate [38, 39] and brittle clay [40]. Because of its simplicity, the σ-criterion seems to be the most popular criterion in mixed mode I-II fracture studies [41]. Therefore, the σ-criterion is utilized in this paper to determine the crack initiation angle.For anisotropic materials, the general form of the elastic stress field near the crack tip in the local Cartesian coordinatesx"-y", as shown in Figure 3, can be expressed in terms of the two SIFs KI and KII as follows [29]:(2.23)σx′′=KI2πrRe[μ1μ2μ1-μ2(μ2cos⁡θ+μ2sin⁡θ-μ1cos⁡θ+μ1sin⁡θ)]+KII2πrRe[1μ1-μ2(μ22cos⁡θ+μ2sin⁡θ-μ12cos⁡θ+μ1sin⁡θ)],σy′′=KI2πrRe[1μ1-μ2(μ1cos⁡θ+μ2sin⁡θ-μ2cos⁡θ+μ1sin⁡θ)]+KII2πrRe[1μ1-μ2(1cos⁡θ+μ2sin⁡θ-1cos⁡θ+μ1sin⁡θ)],τx′′y′′=KI2πrRe[μ1μ2μ1-μ2(1cos⁡θ+μ1sin⁡θ-1cos⁡θ+μ2sin⁡θ)]+KII2πrRe[1μ1-μ2(μ1cos⁡θ+μ1sin⁡θ-μ2cos⁡θ+μ2sin⁡θ)],
where μ1 and μ2 are the roots of (2.1). We have(2.24)σθ=σx′′+σy′′2-σx′′-σy′′2cos⁡2θ-τx′′y′′sin⁡2θ,τrθ=-σx"-σy"2sin⁡2θ+τx"y"cos⁡2θ.
For the σ-criterion, the angle of crack initiation, θ0, must satisfy the following:(2.25)∂σθ∂θ=0(orτrθ=0),∂2σθ∂θ2<0.Figure 3
Crack tip coordinate system and stress components.A numerical procedure is applied to find the angleθ0 when σθ is maximum for known values of the material elastic constants, the anisotropic orientation angle ψ, and the crack geometry.Since the proposed BEM formulation is simple and can be used for any kind of crack geometry, it is straightforward to extend to analyze the crack propagation in anisotropic materials. The process of crack propagation in anisotropic homogeneous material under mixed mode I-II loading is simulated by incremental crack extension with a piecewise linear discretization. For each incremental analysis, the crack extension is conveniently modeled by a new boundary element. A computer program has been developed to automatically generate new data required for sequentially analyzing the changing boundary configuration. Based on the calculation of the SIFs and crack initiation angle for each increment, the procedure of crack propagation can be simulated. The steps in the crack propagation process are summarized as follows (Figure4):(1)
compute the SIFs using the proposed BEM formulation;(2)
determine the angle of crack initiation based on the maximum circumferential stress criterion;(3)
extend the crack by a linear element (of length selected by the user) along the direction determined in step (2);(4)
automatically generate the new BEM mesh;(5)
repeat all of the above steps until the crack is near the outer boundary.Figure 4
Process of crack propagation by increasing the number of linear elements.


## 2.1.  Basic Equations for Anisotropic Elasticity

For the linear elastic, homogeneous, and anisotropic material, the stress and displacement fields can be formulated in terms of two analytical functions,ϕk(zk), of the complex variables zk=x+μky(k=1,2), where μk are the roots of the following characteristic equation [28]:
(2.1)a11μ4-2a16μ3+(2a12+a66)μ2-2a26μ+a22=0,
where the coefficients aij are the compliance components calculated in the x-y coordinate system. The detailed relationship of these components with the material elasticity can be found in Chen et al. (1998) [28]. If the roots μj of  (2.1) are assumed to be distinct, the general expression for the stress and displacements is [28, 29]. One has
(2.2)σx=2Re[μ12ϕ1′(z1)+μ22ϕ2′(z2)],σy=2Re[ϕ1′(z1)+ϕ2′(z2)],τxy=-2Re[μ1ϕ1′(z1)+μ2ϕ2′(z2)],u=2Re[A11ϕ1(z1)+A12ϕ2(z2)],ν=2Re[A21ϕ1(z1)+A22ϕ2(z2)],
where
(2.3)A1j=a11μj2+a12-a16μj,A2j=a12μj+a22μj-a26(j=1,2).
The traction components in the x and y directions are(2.4)Tx=2Re[μ1ϕ1(z1)+μ2ϕ2(z2)],Ty=-2Re[ϕ1(z1)+ϕ2(z2)].With the complex analytical functionsϕi(zi), one can, in general, express (2.2) and (2.4) as follows [28, 30, 31](2.5)ui=2Re[∑j=12Aijϕj(zj)],Ti=-2Re[∑j=12Bijϕj(zj)],σ2i=2Re[∑j=12Bijϕj′(zj)],σ1i=-2Re[∑j=12Bijμjϕj′(zj)](i=1,2),
where zj=x+μjy, the complex number μj, and the elements of the complex matrices A are defined in (2.3), and matrices B can be defined as(2.6)Bij=[-μ1-μ211].Considering the concentrated forces acting at the source point(x0,y0), the analytic functions (ϕk) with the complex variables can be expressed as [30](2.7)ϕk(zk)=-12π[Dklf1ln⁡(zk-zk0)+Dk2f2ln⁡(zk-zk0)],
where zk0=x0+μky0 and fk are the magnitudes of the point force in the k-direction, and(2.8)Dkl=U-1(V-1+V̅-1)-1,U=[A11A12A21A22],V=iUW-1,W=[-μ1-μ211],
where i=-1, overbar means the complex conjugate, and superscript −1 means matrix inverse.Green’s functions of the tractionsTij and displacements Uij can be obtained by substituting (2.7) into (2.5). Their complete expressions are as follows [32]:(2.9)Tij(zk,zk0)=2Re[Qi1(μ1nx-ny)Rj1/(z1-z10)+Qi2(μ2nx-ny)Rj2/(z2-z20)],Uij(zk,zk0)=2Re[Pi1Rj1ln⁡(z1-z10)+Pi2Rj2ln⁡(z2-z20)](i,j=1,2).
In (2.9),  nx and ny are the outward normal components of the field points, and(2.10)Qij=-W=[μ1μ2-1-1].
The complex coefficients Rjk are obtained from the requirements of unit loads at zk0 and from the displacement continuity for the fundamental solution. They are the solutions of the following equation:(2.11)[1-11-1μ1-μ¯1μ2-μ¯2A11-A¯11A12-A¯12A21-A¯21A22-A¯22]⋅[Rj1R¯j1Rj2R¯j2]=[δj2(2πi)-δj1(2πi)00],
where δjk is the Kronecker delta.


## 2.2.  Single-Domain Boundary Integral Equations

A single-domain boundary element method (BEM), based on the relative displacements at the crack tip, is used to determine the mixed-mode SIFs of anisotropic materials. The single-domain BEM formulation consists of the following displacement and traction integral equations (see Figure1).(1)
Displacement integral equation: We have(2.12)Cij(zk,B0)uj(zk,B0)+∫ΓBTij(zk,B,zk,B0)uj(zk,B)dΓ(zk,B)+∫ΓCTij(zk,C,zk,B0)[uj(zk,C+)-uj(zk,C-)]dΓ(zk,C)=∫ΓBUij(zk,B,zk,B0)tj(zk,B)dΓ(zk,B),
where zk and zk0 are the field and source points on the boundary. Here, the subscripts B and C denote the outer boundary and the crack surface, respectively. Cij are coefficients that depend only on the local geometry of the uncracked boundary at zk,B0; Tij and Uij are the Green’s traction and displacement given in (2.9); uj and tj are the boundary displacement and traction; ΓC has the same outward normal as ΓC+.(2)
Traction integral equation: One has(2.13)0.5tj(zk,C0)+nm(zk,C0)∫ΓBClmikTij,k(zk,B,zk,C0)uj(zk,B)dΓ(zk,B)+nm(zk,C0)∫ΓCClmikTij,k(zk,C,zk,C0)[uj(zk,C+)-uj(zk,C-)]dΓ(zk,C)=nm(zk,C0)∫ΓBClmikUij,k(zk,B,zk,C0)tj(zk,B)dΓ(zk,B),
where Clmik is the fourth-order stiffness tensor for anisotropic medium; nm is the unit outward normal to the contour path; the gradient tensors Tij,k and Uij,k denote differentiation with respect to zk0.Figure 1
Geometry of a two-dimensional cracked domain.The Cauchy singularity in (2.12) can be avoided by the rigid-body motion method. The integrand on the right-hand side of (2.12) has only integrable singularity, which can be resolved by the bicubic transformation method [33]. The hypersingularity in (2.13) is resolved by the Gauss quadrature formula, which is very similar to the traditional weighted Gauss quadrature but with a different weight. Therefore, (2.12) and (2.13) can be discretized and solved numerically for the unknown displacements (or the relative crack opening displacements (CODs) on the crack surface) and tractions. In the following section, we present an approach for the evaluation of mixed-mode SIFs.


## 2.3.  Evaluation of Mixed-Mode SIFs

The mixed-mode SIFs for anisotropic materials can be determined by using the extrapolation method of the relative COD, combined with a set of the shape functions. The relative COD is defined as [18](2.14)Δui=∑k=13ϕkΔuik,
where the subscript i(=1,2) denotes the components of the relative COD, and the superscript k(=1,2,3) denotes the relative COD at nodes s=-2/3,0,2/3, respectively; ϕk are the shape functions which can be expressed as follows:(2.15)ϕ1=338s+1[5-8(s+1)+3(s+1)2],ϕ2=14s+1[-5+18(s+1)-9(s+1)2],ϕ3=3385s+1[1-4(s+1)+3(s+1)2].It is well known that for a crack in a homogeneous and anisotropic solid, the relative COD is proportional tor, where r is the distance behined the crack tip. Therefore, the relation of the relative CODs and the SIFs can be found as [18, 33](2.16)Δu1=22rπ(H11KI+H12KII),Δu2=22rπ(H21KI+H22KII),
where(2.17)H11=Im⁡(μ2A11-μ1A12μ1-μ2),H12=Im⁡(A11-A12μ1-μ2),H21=Im⁡(μ2A21-μ1A22μ1-μ2),H22=Im⁡(A21-A22μ1-μ2).Substituting (2.14) into (2.16), a set of algebraic equations is obtained, and then the SIFs KI and KII can be solved. A sign convention for the corresponding SIFs is shown in Figure 2. Using the relative COD, the sign of SIFs can then be determined.Figure 2
Sign convention used for determining SIFs in mode I and mode II.


## 2.4.  Particular Solutions of Gravity and Far-Field Stresses

As mentioned above, if the particular solutions corresponding to the body force of gravity and far-field stresses can be derived in exact closed form, the single-domain BEM formulation presented in this study can then be applied to solve the body force and far-field problems. For the gravity force, the exact close-form solutions can be obtained in a similar way as for the corresponding half-space case [17]. Assuming that the gravity has the components gx and gy, respectively, in the x- and y-directions, the particular solution for the displacements can be found as [34](2.18)uxp=a1ρgxx2+b1ρgyy2,uyp=a2ρgxx2+b2ρgyy2,
where the coefficients ai and bi depend on the elastic stiffness and their expressions are(2.19)a1=0.5(c55c66-c56c56)Δ1,a2=0.5(c15c56-c55c16)Δ1,b1=0.5(c24c46-c26c44)Δ2,b2=0.5(c44c66-c46c46)Δ2,
where(2.20)Δ1=c11c55c66+2c16c15c56-c55c162-c11c562-c66c152,Δ2=c22c44c66+2c26c24c46-c44c262-c22c462-c66c242,
and cij are the elastic stiffness coefficients.Similarly, the particular stresses can be expressed as(2.21)[σ11pσ22pσ23pσ13pσ12p]=[d11d12d21d22d41d42d51d52d61d62][ρgxxρgyy].
Again, dij depend on the elastic stiffness coefficients, and their expressions are(2.22)d11=1,d21=2(a1c12+a2c26),d41=2(a1c14+a2c46),d51=0,d61=0,d12=2(b2c12+b1c16),d22=1,d42=0,d52=2(b2c25+b1c56),d62=0.


## 2.5.  Crack Initiation and Propagation

Three criteria are commonly utilized to predict the crack initiation angle in fracture mechanics problems: the maximum circumferential stress criterion, orσ-criterion [35]; the maximum energy release rate criterion, or G-criterion [36]; the minimum strain energy density criterion, or S-criterion [37]. Among them, the σ-criterion has been found to predict well the directions of crack initiation compared to the experimental results for polymethylmethacrylate [38, 39] and brittle clay [40]. Because of its simplicity, the σ-criterion seems to be the most popular criterion in mixed mode I-II fracture studies [41]. Therefore, the σ-criterion is utilized in this paper to determine the crack initiation angle.For anisotropic materials, the general form of the elastic stress field near the crack tip in the local Cartesian coordinatesx"-y", as shown in Figure 3, can be expressed in terms of the two SIFs KI and KII as follows [29]:(2.23)σx′′=KI2πrRe[μ1μ2μ1-μ2(μ2cos⁡θ+μ2sin⁡θ-μ1cos⁡θ+μ1sin⁡θ)]+KII2πrRe[1μ1-μ2(μ22cos⁡θ+μ2sin⁡θ-μ12cos⁡θ+μ1sin⁡θ)],σy′′=KI2πrRe[1μ1-μ2(μ1cos⁡θ+μ2sin⁡θ-μ2cos⁡θ+μ1sin⁡θ)]+KII2πrRe[1μ1-μ2(1cos⁡θ+μ2sin⁡θ-1cos⁡θ+μ1sin⁡θ)],τx′′y′′=KI2πrRe[μ1μ2μ1-μ2(1cos⁡θ+μ1sin⁡θ-1cos⁡θ+μ2sin⁡θ)]+KII2πrRe[1μ1-μ2(μ1cos⁡θ+μ1sin⁡θ-μ2cos⁡θ+μ2sin⁡θ)],
where μ1 and μ2 are the roots of (2.1). We have(2.24)σθ=σx′′+σy′′2-σx′′-σy′′2cos⁡2θ-τx′′y′′sin⁡2θ,τrθ=-σx"-σy"2sin⁡2θ+τx"y"cos⁡2θ.
For the σ-criterion, the angle of crack initiation, θ0, must satisfy the following:(2.25)∂σθ∂θ=0(orτrθ=0),∂2σθ∂θ2<0.Figure 3
Crack tip coordinate system and stress components.A numerical procedure is applied to find the angleθ0 when σθ is maximum for known values of the material elastic constants, the anisotropic orientation angle ψ, and the crack geometry.Since the proposed BEM formulation is simple and can be used for any kind of crack geometry, it is straightforward to extend to analyze the crack propagation in anisotropic materials. The process of crack propagation in anisotropic homogeneous material under mixed mode I-II loading is simulated by incremental crack extension with a piecewise linear discretization. For each incremental analysis, the crack extension is conveniently modeled by a new boundary element. A computer program has been developed to automatically generate new data required for sequentially analyzing the changing boundary configuration. Based on the calculation of the SIFs and crack initiation angle for each increment, the procedure of crack propagation can be simulated. The steps in the crack propagation process are summarized as follows (Figure4):(1)
compute the SIFs using the proposed BEM formulation;(2)
determine the angle of crack initiation based on the maximum circumferential stress criterion;(3)
extend the crack by a linear element (of length selected by the user) along the direction determined in step (2);(4)
automatically generate the new BEM mesh;(5)
repeat all of the above steps until the crack is near the outer boundary.Figure 4
Process of crack propagation by increasing the number of linear elements.


## 3. Results of Numerical Analysis

A computer program, which is based on the aforementioned BEM formulation, has been accomplished to analyze the SIFs for isotropic and anisotropic materials. The SIFs for both isotropic and anisotropic materials with different crack angles, crack lengths, and anisotropic orientations are analyzed by the program.The stresses in a linear elastic anisotropic medium depend merely on ratios of elastic constants and ratios of geometric dimensions. Thus, units of the problems can be eliminated by the method of normalization. Accordingly, the SIFs of the following examples are normalized with respect to the applied load and to the square root of half crack length.Totally 9 numerical examples including isotropic and anisotropic materials are presented to illustrate the accuracy and versatility of the proposed BEM program for determining the SIFs, predicting the crack initiation angle, and simulating crack propagation path. The examples include cases for finite/infinite domains, curved/edge cracks, and isotropic/anisotropic conditions. A generalized plane stress is assumed in all the examples except for crack propagation simulation of slope failure.


### 3.1.  Stress Intensity Factors Determination

Example 1: Isotropic Cracked Brazilian Disc
In order to compare our results with the existing published results, an isotropic and cracked Brazilian disc with a central slant crack is considered. The geometry of the problem is that of a thin circular disc of radiusR and thickness t with a central crack of length 2a, loaded with a pair of concentrated and diametral compressive load W, as shown in Figure 5. The outer boundary and crack surface are discretized with 28 continuous and 10 discontinuous quadratic elements, respectively. Two cases are analyzed: (1) a/R=0.5, the crack angle β varies between 0 and π/2, and (2) β=45°, a/R varies between 0.1 and 0.7. The two normalized SIFs, FI(=KI/σπa) and  FII(=KII/σπa), calculated with the BEM program for these two cases, are compared with those obtained numerically by Atkinson et al. [42] and Chen et al. [28]. The results are shown in Tables 1 and 2. In general, a good agreement is found among these three methods.Table 1
Normalized SIFs for a central slant in an isotropic Brazilian disc subjected to a concentrated load (a/R=0.5).
β (rad.)Atkinson et al. (1982) [42]Chen et al. (1998) [28]This studyKI/σπaKII/σπaKI/σπaKII/σπaKI/σπaKII/σπa01.38701.33901.3430π/160.970−1.3400.960−1.2750.952−1.2812π/160.030−2.1130.074−2.0610.056−2.0503π/16−0.946−2.300−0.903−2.275−0.915−2.262π/4−1.784−2.132−1.737−2.103−1.749−2.0985π/16−2.446−1.728−2.377−1.711−2.395−1.7146π/16−2.885−1.188−2.826−1.197−2.851−1.2027π/16−3.127−0.604−3.092−0.614−3.123−0.617π/2−3.2080−3.1800−3.2130Table 2
Normalized SIFs for a central slant in an isotropic Brazilian disc subjected to a concentrated load (β=45°).
a/RAtkinson et al. (1982) [42]Chen et al. (1998) [28]This studyKI/σπaKII/σπaKI/σπaKII/σπaKI/σπaKII/σπa0.1−1.035−2.010−1.020−1.968−1.018−1.9650.2−1.139−2.035−1.116−1.995−1.116−1.9920.3−1.306−2.069−1.272−2.036−1.277−2.0290.4−1.528−2.100−1.484−2.069−1.492−2.0650.5−1.784−2.132−1.737−2.103−1.749−2.0980.6−2.048−2.200−2.020−2.148−2.039−2.1390.7N.A.N.A.−2.337−2.213−2.364−2.224Figure 5
The geometry of cracked Brazilian disc.Example 2: A Curved Crack in an Infinite Domain
Consider a circular-arc crack of a radiusR embedded in an infinite domain under a far-field tensile stress σ and out-of-plane shear stress τ as shown in Figure 6. The center of the circular arc is taken at the origin of the coordinate system, the midpoint of the crack is located on the x-axis, and the angle subtended by the crack is 2α. In this example only 20 discontinuous quadratic elements are used to discretize the curved crack surface. For α=30° and α=45°, the numerical solutions of SIFs calculated by this study as well as the analytic ones by Tada et al. [43] are shown in Table 3.Table 3
Normalized SIFs for a circular arc crack in an infinite domain.
αTada et al. (1985) [43]This studyKI/σπaKII/σπaKI/σπaKII/σπa30°0.55010.33040.55280.328845°0.45740.51120.46430.5126Figure 6
A curved crack under a far-field tensile stressσ and out-of-plane shear stress τ.Example 3: Anisotropic Rectangular Plate under a Uniform Tension
In order to evaluate the influence of material anisotropy on the SIFs, consider an anisotropic rectangular plate of width 2w, and height 2h with a central crack inclined 45° to the x-axis, as shown in Figure 7. The plate is loaded with a uniform tensile stress in the y direction. The ratios of crack length and of height to width are a/W=0.2 and h/W=2.0, respectively. The material is glass-epoxy with elastic properties E=48.26GPa, E′=17.24GPa, ν′=0.29, and G′=6.89GPa. The direction of the fibers is rotated from ψ=0° to ψ=180°. The outer boundary and crack surface are discretized with 32 continuous and 10 discontinuous quadratic elements, respectively. Table 4 shows the results obtained by the proposed method as well as those by Sollero and Aliabadi [44], Gandhi [45], and Chen et al. [28]. Again, an excellent agreement is obtained.Table 4
Normalized SIFs for a central slant in anisotropic rectangular plate subjected to a uniform tension (β=45°,a/R=0.5).
ψ (Deg.)Gandhi (1972) [45]Sollero and Aliabadi (1995) [15]Chen et al. (1998) [28]This studyKI/σπaKII/σπaKI/σπaKII/σπaKI/σπaKII/σπaKI/σπaKII/σπa00.5220.5070.5170.5060.5190.5040.5240.514450.5150.5050.5130.5020.5160.5050.5180.513900.5130.5090.5150.5100.5370.5320.5330.5271050.5170.5100.5180.5120.5070.5020.5400.5281200.5240.5120.5260.5130.5200.5080.5410.5261350.5320.5110.5350.5140.5320.5110.5360.5221800.5220.5070.5170.5060.5190.5040.5240.514Figure 7
An anisotropic rectangular plate with a central crack inclined 45° under uniform tensile stressσ.Example 4: Cracked Body Hanging under Its Own Weight
Consider a body with a crack of length 2a under the gravitational loading as shown in Figure 8. The ratios of distance from crack to free end to height  are h/l=0.5. The ratio of crack length to width is varied from 0.05 to 0.5. The outer boundary and crack surface are discretized with 44 continuous and 30 discontinuous quadratic elements, respectively. Plane strain condition is considered. Table 5 shows the results obtained by the proposed method as well as the analytic ones by Shi’s handbook [46], and the numerical ones by Ostanin et al. (2011) [47]. Again, an excellent agreement is obtained.Table 5
Normalized SIFs of mode I for a cracked body under its own weight.
a/bExact Sol.Numerical Sol.Sih (1973) [46]Ostanin et al. (2011) [47]This studyKI/σπaKI/σπaDiff. (%)KI/σπaDiff. (%)0.050.7930.7950.250.7950.250.101.1211.1250.351.1240.270.151.3731.3810.581.3790.440.201.5851.6000.941.5910.380.251.7721.7981.51.7790.400.301.9421.9802.01.9560.720.352.0972.5142.72.0850.570.402.2422.3213.52.2300.540.452.3782.4844.92.3680.800.502.5072.6445.52.4502.27Figure 8
A cracked body hanging under its own weight.Example 5: Edge Crack in an Elastic Body Hanging under Its Own Weight
Figure9 shows a square block l×l containing an edge crack of length a is hanging under its own weight. The distance of crack tip to the free end is h. Poisson’s ratio is assumed to be zero. In this example, the outer boundary and crack surface are discretized with 44 continuous and 30 discontinuous quadratic elements, respectively. Table 6 shows the normalized SIFs of mode I. This problem was solved previously by Ostanin et al. (2011) [47] using a complex variable boundary element method (CVBEM). As shown in Table 6, the present numerical results are in excellent agreement with those obtained by Ostanin et al. (2011) [47].Table 6
Normalized SIFs of mode I for an edge-cracked body hanging under its own weight.
h/la/l0.10.20.30.40.5Ostanin et al. (2011) [47]This studyDiff. (%)Ostanin et al. (2011) [47]This studyDiff. (%)Ostanin et al. (2011) [47]This studyDiff. (%)Ostanin et al. (2011) [47]This studyDiff. (%)Ostanin et al. (2011) [47]This studyDiff. (%)0.10.0490.0502.040.1360.1370.740.2580.25800.4150.41500.6100.6090.160.20.0770.0761.300.1670.1680.600.2820.28200.4270.4260.230.6120.6100.330.30.1080.1090.930.2070.2060.480.3300.3290.300.4880.48800.7050.7060.140.40.1410.1400.710.2520.25200.3900.3890.260.5760.5770.170.8440.8430.120.50.1750.1760.570.3000.2990.330.4590.4600.220.6780.6790.151.0091.0110.20Figure 9
An edge crack elastic body hanging under its own weight.


### 3.2.  Crack Initiation Angle Prediction and Propagation Path Simulation

The proposed BEM formulation combined with the maximum circumferential stress criterion is developed to predict the angle of crack initiation and to simulate the path of crack propagation under mixed-mode loading. The crack propagation process in the cracked materials is numerically estimated by two-dimensional stress and displacement analysis. In order to understand the behavior of cracks under mixed-mode loading, the BEM program is applied.Example 1: Crack Initiation Angle Prediction
The proposed BEM formulation is also used to predict the initial growth of cracks in anisotropic materials. To examine the validity of our crack initiation prediction procedure, the tests of Erdogan and Sih [35] and Vallejo [40] are reproduced numerically with our BEM program. Erdogan and Sih [35] conducted uniaxial tension test on isotropic Plexiglass sheets 229 × 457 × 4.8 mm in size containing a 50.8 mm central crack. The crack inclination angle β between the crack plane and the tensile stress is varied. Figure 10 shows the variation of the crack initiation angle θ0 with the crack angle β determined numerically and experimentally. A good agreement is found between the experimental results of Erdogan and Sih [35] and our numerical predictions.
Other verification is done using the experimental results of Vallejo [40]. The uniaxial compression tests were conducted on cracked prismatic specimens of kaolinite clay 76.2 × 76.2 × 25.4 mm in size containing a central crack 24.9 mm in length by Vallejo [40]. Several tests are carried out by varying the crack angle β between the crack plane and the compressive stress. Figure 11 shows a comparison between the crack initiation angles measured experimentally and those predicted numerically. Again, a good agreement is found with the experimental results.Figure 10
Variation of crack initiation angleθ0 with the crack angle β. Plexiglass plate subjected to uniaxial tension.Figure 11
Variation of crack initiation angleθ0 with the crack angle β. Prismatic sample of kaolinite clay subjected to uniaxial compression.Example 2: Isotropic Cracked Plate under Pure Mode I Loading
The behavior of an existing crack under pure mode I loading is studied for the application. The first problem, as shown in Figure12, is a square plate with a horizontal edge crack subjected to uniaxial tension. The width of the square plate is w. The initial crack length, a, is equal to w/3. The no-crack boundary and crack surface are discretized with 55 continuous and 10 discontinuous quadratic elements, respectively. According to the experimental results from Erdogan and Sih [35], it is known that the crack propagation angle (θ) is zero when the crack inclination angle (β) is 90 degree with respect to y-axis, which means that the crack will propagate along the horizontal direction. Figure 13 shows the path of crack propagation under pure mode I loading. It is shown that the path of crack propagation is a horizontal line, which is in full agreement with the experimental results made by Erdogan and Sih [35].Figure 12
Square plate with a straight edge crack under pure mode I loading.Figure 13
Crack propagation path and normalized SIFs of the cracked plate initially under pure mode I loading.Example 3: Crack Propagation Path Simulation on Anisotropic Material
To demonstrate the proposed BEM procedure when predicting crack propagation in the anisotropic materials under mixed-mode I-II loading, the propagation path in a CSTBD specimen is numerically predicted and compared with the actual laboratory observations. In these experiments, a crack initially inclined with respect to the applied stress is allowed to grow under concentrated diametrical loading. The Brazilian tests on CSTBD specimens with a diameter of 7.4 cm, a thickness of 1 cm, and a crack length of 2.2 cm are conducted to observe the actual propagation paths and are compared with the numerical predictions. Details of the experimental procedure can be found in the paper by Ke et al. [48]. The five elastic constants of anisotropic marble are E=78.302 GPa, E′=67.681 GPa, ν=0.267, ν′=0.185, G=30.735 GPa, and G′=25.336 GPa, respectively. The ratios of  E/E′ and E/G′ are 1.156 and 3.091, respectively. Since the value of  E/E′=1.156, this marble can be classified as a slightly anisotropic rock.
Two specimens with the material inclination angleψ=45°, defined as the AM-4, and DM-4, have crack angles β=0° and β=45° respectively. After Brazilian tests with cracked specimens, the paths of crack propagation for AM-4 and DM-4 are shown in Figures 14 and 16, respectively. It can be observed that the crack propagates nearly perpendicular to the crack surface in the initial stage and then rapidly approaches toward the loading point. The proposed BEM procedure is also used to simulate crack propagation in the CSTBD specimens. The outer boundary and crack surface are discretized with 28 continuous and 20 discontinuous quadratic elements, respectively. Figures 15 and 17 are the comparisons of crack propagation paths between experimental observations and numerical predictions in AM4 and DM-4, respectively. Again, the proposed BEM procedure accurately simulates the crack propagation in these anisotropic specimens. According to the simulations of foregoing examples, it can be concluded that the proposed BEM is capable of predicting the crack propagation in anisotropic rocks.Figure 14
Photograph of specimen AM-4 after failure (ψ=0°and β=45°).Figure 15
Propagation of a crack at the center of a CSTBD specimen withψ=0°  and β=45°. Comparison between experimental observations and numerical predictions for specimen AM-4.Figure 16
Photograph of specimen DM-4 after failure (ψ=45°  and β=45°).Figure 17
Propagation of a crack at the center of a CSTBD specimen withψ=45°  and β=45°. Comparison between experimental observations and numerical predictions for specimen DM-4.


### 3.3.  Crack Propagation Simulation of Slope Failure

A slope 15 m high with a 2 m deep initial tension crack located on the top of the slope is used for the analysis. 10 discontinuous quadratic elements and 46 quadratic elements were used to discretize the crack surface and outer boundary. Plane strain condition is considered. The top of the slope has surface recharge due to heavy rain. Figure18 gives a description of the problem layout that includes the slope geometry and the crack. The crack propagation process under body force of gravity is considered. The crack increment length is fixed at a size 1/6 times the crack tip element. An anisotropic rock slope was simulated with different values of rock anisotropy (ψ=0°,15°,30°,45°,60°,75°,90°). The dimensionless elastic constant (E/E′=3) is considered. Figure 18 gives a description of the problem layout that includes the slope geometry, the crack, and anisotropy orientation. Figure 19 shows the slope failure surfaces simulated by the BEM. It can be found that the slope failure surfaces depend on the different anisotropy inclination angles. The anisotropy orientation angle ψ has a strong influence on the surface of slope failure.Figure 18
The description of an anisotropic rock slope.Figure 19
The anisotropic rock slope failure surface simulated by the proposed BEM.


## 3.1.  Stress Intensity Factors Determination

Example 1: Isotropic Cracked Brazilian Disc
In order to compare our results with the existing published results, an isotropic and cracked Brazilian disc with a central slant crack is considered. The geometry of the problem is that of a thin circular disc of radiusR and thickness t with a central crack of length 2a, loaded with a pair of concentrated and diametral compressive load W, as shown in Figure 5. The outer boundary and crack surface are discretized with 28 continuous and 10 discontinuous quadratic elements, respectively. Two cases are analyzed: (1) a/R=0.5, the crack angle β varies between 0 and π/2, and (2) β=45°, a/R varies between 0.1 and 0.7. The two normalized SIFs, FI(=KI/σπa) and  FII(=KII/σπa), calculated with the BEM program for these two cases, are compared with those obtained numerically by Atkinson et al. [42] and Chen et al. [28]. The results are shown in Tables 1 and 2. In general, a good agreement is found among these three methods.Table 1
Normalized SIFs for a central slant in an isotropic Brazilian disc subjected to a concentrated load (a/R=0.5).
β (rad.)Atkinson et al. (1982) [42]Chen et al. (1998) [28]This studyKI/σπaKII/σπaKI/σπaKII/σπaKI/σπaKII/σπa01.38701.33901.3430π/160.970−1.3400.960−1.2750.952−1.2812π/160.030−2.1130.074−2.0610.056−2.0503π/16−0.946−2.300−0.903−2.275−0.915−2.262π/4−1.784−2.132−1.737−2.103−1.749−2.0985π/16−2.446−1.728−2.377−1.711−2.395−1.7146π/16−2.885−1.188−2.826−1.197−2.851−1.2027π/16−3.127−0.604−3.092−0.614−3.123−0.617π/2−3.2080−3.1800−3.2130Table 2
Normalized SIFs for a central slant in an isotropic Brazilian disc subjected to a concentrated load (β=45°).
a/RAtkinson et al. (1982) [42]Chen et al. (1998) [28]This studyKI/σπaKII/σπaKI/σπaKII/σπaKI/σπaKII/σπa0.1−1.035−2.010−1.020−1.968−1.018−1.9650.2−1.139−2.035−1.116−1.995−1.116−1.9920.3−1.306−2.069−1.272−2.036−1.277−2.0290.4−1.528−2.100−1.484−2.069−1.492−2.0650.5−1.784−2.132−1.737−2.103−1.749−2.0980.6−2.048−2.200−2.020−2.148−2.039−2.1390.7N.A.N.A.−2.337−2.213−2.364−2.224Figure 5
The geometry of cracked Brazilian disc.Example 2: A Curved Crack in an Infinite Domain
Consider a circular-arc crack of a radiusR embedded in an infinite domain under a far-field tensile stress σ and out-of-plane shear stress τ as shown in Figure 6. The center of the circular arc is taken at the origin of the coordinate system, the midpoint of the crack is located on the x-axis, and the angle subtended by the crack is 2α. In this example only 20 discontinuous quadratic elements are used to discretize the curved crack surface. For α=30° and α=45°, the numerical solutions of SIFs calculated by this study as well as the analytic ones by Tada et al. [43] are shown in Table 3.Table 3
Normalized SIFs for a circular arc crack in an infinite domain.
αTada et al. (1985) [43]This studyKI/σπaKII/σπaKI/σπaKII/σπa30°0.55010.33040.55280.328845°0.45740.51120.46430.5126Figure 6
A curved crack under a far-field tensile stressσ and out-of-plane shear stress τ.Example 3: Anisotropic Rectangular Plate under a Uniform Tension
In order to evaluate the influence of material anisotropy on the SIFs, consider an anisotropic rectangular plate of width 2w, and height 2h with a central crack inclined 45° to the x-axis, as shown in Figure 7. The plate is loaded with a uniform tensile stress in the y direction. The ratios of crack length and of height to width are a/W=0.2 and h/W=2.0, respectively. The material is glass-epoxy with elastic properties E=48.26GPa, E′=17.24GPa, ν′=0.29, and G′=6.89GPa. The direction of the fibers is rotated from ψ=0° to ψ=180°. The outer boundary and crack surface are discretized with 32 continuous and 10 discontinuous quadratic elements, respectively. Table 4 shows the results obtained by the proposed method as well as those by Sollero and Aliabadi [44], Gandhi [45], and Chen et al. [28]. Again, an excellent agreement is obtained.Table 4
Normalized SIFs for a central slant in anisotropic rectangular plate subjected to a uniform tension (β=45°,a/R=0.5).
ψ (Deg.)Gandhi (1972) [45]Sollero and Aliabadi (1995) [15]Chen et al. (1998) [28]This studyKI/σπaKII/σπaKI/σπaKII/σπaKI/σπaKII/σπaKI/σπaKII/σπa00.5220.5070.5170.5060.5190.5040.5240.514450.5150.5050.5130.5020.5160.5050.5180.513900.5130.5090.5150.5100.5370.5320.5330.5271050.5170.5100.5180.5120.5070.5020.5400.5281200.5240.5120.5260.5130.5200.5080.5410.5261350.5320.5110.5350.5140.5320.5110.5360.5221800.5220.5070.5170.5060.5190.5040.5240.514Figure 7
An anisotropic rectangular plate with a central crack inclined 45° under uniform tensile stressσ.Example 4: Cracked Body Hanging under Its Own Weight
Consider a body with a crack of length 2a under the gravitational loading as shown in Figure 8. The ratios of distance from crack to free end to height  are h/l=0.5. The ratio of crack length to width is varied from 0.05 to 0.5. The outer boundary and crack surface are discretized with 44 continuous and 30 discontinuous quadratic elements, respectively. Plane strain condition is considered. Table 5 shows the results obtained by the proposed method as well as the analytic ones by Shi’s handbook [46], and the numerical ones by Ostanin et al. (2011) [47]. Again, an excellent agreement is obtained.Table 5
Normalized SIFs of mode I for a cracked body under its own weight.
a/bExact Sol.Numerical Sol.Sih (1973) [46]Ostanin et al. (2011) [47]This studyKI/σπaKI/σπaDiff. (%)KI/σπaDiff. (%)0.050.7930.7950.250.7950.250.101.1211.1250.351.1240.270.151.3731.3810.581.3790.440.201.5851.6000.941.5910.380.251.7721.7981.51.7790.400.301.9421.9802.01.9560.720.352.0972.5142.72.0850.570.402.2422.3213.52.2300.540.452.3782.4844.92.3680.800.502.5072.6445.52.4502.27Figure 8
A cracked body hanging under its own weight.Example 5: Edge Crack in an Elastic Body Hanging under Its Own Weight
Figure9 shows a square block l×l containing an edge crack of length a is hanging under its own weight. The distance of crack tip to the free end is h. Poisson’s ratio is assumed to be zero. In this example, the outer boundary and crack surface are discretized with 44 continuous and 30 discontinuous quadratic elements, respectively. Table 6 shows the normalized SIFs of mode I. This problem was solved previously by Ostanin et al. (2011) [47] using a complex variable boundary element method (CVBEM). As shown in Table 6, the present numerical results are in excellent agreement with those obtained by Ostanin et al. (2011) [47].Table 6
Normalized SIFs of mode I for an edge-cracked body hanging under its own weight.
h/la/l0.10.20.30.40.5Ostanin et al. (2011) [47]This studyDiff. (%)Ostanin et al. (2011) [47]This studyDiff. (%)Ostanin et al. (2011) [47]This studyDiff. (%)Ostanin et al. (2011) [47]This studyDiff. (%)Ostanin et al. (2011) [47]This studyDiff. (%)0.10.0490.0502.040.1360.1370.740.2580.25800.4150.41500.6100.6090.160.20.0770.0761.300.1670.1680.600.2820.28200.4270.4260.230.6120.6100.330.30.1080.1090.930.2070.2060.480.3300.3290.300.4880.48800.7050.7060.140.40.1410.1400.710.2520.25200.3900.3890.260.5760.5770.170.8440.8430.120.50.1750.1760.570.3000.2990.330.4590.4600.220.6780.6790.151.0091.0110.20Figure 9
An edge crack elastic body hanging under its own weight.


## 3.2.  Crack Initiation Angle Prediction and Propagation Path Simulation

The proposed BEM formulation combined with the maximum circumferential stress criterion is developed to predict the angle of crack initiation and to simulate the path of crack propagation under mixed-mode loading. The crack propagation process in the cracked materials is numerically estimated by two-dimensional stress and displacement analysis. In order to understand the behavior of cracks under mixed-mode loading, the BEM program is applied.Example 1: Crack Initiation Angle Prediction
The proposed BEM formulation is also used to predict the initial growth of cracks in anisotropic materials. To examine the validity of our crack initiation prediction procedure, the tests of Erdogan and Sih [35] and Vallejo [40] are reproduced numerically with our BEM program. Erdogan and Sih [35] conducted uniaxial tension test on isotropic Plexiglass sheets 229 × 457 × 4.8 mm in size containing a 50.8 mm central crack. The crack inclination angle β between the crack plane and the tensile stress is varied. Figure 10 shows the variation of the crack initiation angle θ0 with the crack angle β determined numerically and experimentally. A good agreement is found between the experimental results of Erdogan and Sih [35] and our numerical predictions.
Other verification is done using the experimental results of Vallejo [40]. The uniaxial compression tests were conducted on cracked prismatic specimens of kaolinite clay 76.2 × 76.2 × 25.4 mm in size containing a central crack 24.9 mm in length by Vallejo [40]. Several tests are carried out by varying the crack angle β between the crack plane and the compressive stress. Figure 11 shows a comparison between the crack initiation angles measured experimentally and those predicted numerically. Again, a good agreement is found with the experimental results.Figure 10
Variation of crack initiation angleθ0 with the crack angle β. Plexiglass plate subjected to uniaxial tension.Figure 11
Variation of crack initiation angleθ0 with the crack angle β. Prismatic sample of kaolinite clay subjected to uniaxial compression.Example 2: Isotropic Cracked Plate under Pure Mode I Loading
The behavior of an existing crack under pure mode I loading is studied for the application. The first problem, as shown in Figure12, is a square plate with a horizontal edge crack subjected to uniaxial tension. The width of the square plate is w. The initial crack length, a, is equal to w/3. The no-crack boundary and crack surface are discretized with 55 continuous and 10 discontinuous quadratic elements, respectively. According to the experimental results from Erdogan and Sih [35], it is known that the crack propagation angle (θ) is zero when the crack inclination angle (β) is 90 degree with respect to y-axis, which means that the crack will propagate along the horizontal direction. Figure 13 shows the path of crack propagation under pure mode I loading. It is shown that the path of crack propagation is a horizontal line, which is in full agreement with the experimental results made by Erdogan and Sih [35].Figure 12
Square plate with a straight edge crack under pure mode I loading.Figure 13
Crack propagation path and normalized SIFs of the cracked plate initially under pure mode I loading.Example 3: Crack Propagation Path Simulation on Anisotropic Material
To demonstrate the proposed BEM procedure when predicting crack propagation in the anisotropic materials under mixed-mode I-II loading, the propagation path in a CSTBD specimen is numerically predicted and compared with the actual laboratory observations. In these experiments, a crack initially inclined with respect to the applied stress is allowed to grow under concentrated diametrical loading. The Brazilian tests on CSTBD specimens with a diameter of 7.4 cm, a thickness of 1 cm, and a crack length of 2.2 cm are conducted to observe the actual propagation paths and are compared with the numerical predictions. Details of the experimental procedure can be found in the paper by Ke et al. [48]. The five elastic constants of anisotropic marble are E=78.302 GPa, E′=67.681 GPa, ν=0.267, ν′=0.185, G=30.735 GPa, and G′=25.336 GPa, respectively. The ratios of  E/E′ and E/G′ are 1.156 and 3.091, respectively. Since the value of  E/E′=1.156, this marble can be classified as a slightly anisotropic rock.
Two specimens with the material inclination angleψ=45°, defined as the AM-4, and DM-4, have crack angles β=0° and β=45° respectively. After Brazilian tests with cracked specimens, the paths of crack propagation for AM-4 and DM-4 are shown in Figures 14 and 16, respectively. It can be observed that the crack propagates nearly perpendicular to the crack surface in the initial stage and then rapidly approaches toward the loading point. The proposed BEM procedure is also used to simulate crack propagation in the CSTBD specimens. The outer boundary and crack surface are discretized with 28 continuous and 20 discontinuous quadratic elements, respectively. Figures 15 and 17 are the comparisons of crack propagation paths between experimental observations and numerical predictions in AM4 and DM-4, respectively. Again, the proposed BEM procedure accurately simulates the crack propagation in these anisotropic specimens. According to the simulations of foregoing examples, it can be concluded that the proposed BEM is capable of predicting the crack propagation in anisotropic rocks.Figure 14
Photograph of specimen AM-4 after failure (ψ=0°and β=45°).Figure 15
Propagation of a crack at the center of a CSTBD specimen withψ=0°  and β=45°. Comparison between experimental observations and numerical predictions for specimen AM-4.Figure 16
Photograph of specimen DM-4 after failure (ψ=45°  and β=45°).Figure 17
Propagation of a crack at the center of a CSTBD specimen withψ=45°  and β=45°. Comparison between experimental observations and numerical predictions for specimen DM-4.


## 3.3.  Crack Propagation Simulation of Slope Failure

A slope 15 m high with a 2 m deep initial tension crack located on the top of the slope is used for the analysis. 10 discontinuous quadratic elements and 46 quadratic elements were used to discretize the crack surface and outer boundary. Plane strain condition is considered. The top of the slope has surface recharge due to heavy rain. Figure18 gives a description of the problem layout that includes the slope geometry and the crack. The crack propagation process under body force of gravity is considered. The crack increment length is fixed at a size 1/6 times the crack tip element. An anisotropic rock slope was simulated with different values of rock anisotropy (ψ=0°,15°,30°,45°,60°,75°,90°). The dimensionless elastic constant (E/E′=3) is considered. Figure 18 gives a description of the problem layout that includes the slope geometry, the crack, and anisotropy orientation. Figure 19 shows the slope failure surfaces simulated by the BEM. It can be found that the slope failure surfaces depend on the different anisotropy inclination angles. The anisotropy orientation angle ψ has a strong influence on the surface of slope failure.Figure 18
The description of an anisotropic rock slope.Figure 19
The anisotropic rock slope failure surface simulated by the proposed BEM.


## 4. Conclusions

A formulation of the BEM, based on the relative displacements near the crack tip, is utilized to determine the mixed-mode SIFs of anisotropic rocks. Numerical examples for the determination of the mixed-mode SIFs for a CSTBD specimen are presented for isotropic and anisotropic media. The numerical results obtained by the proposed method are in good agreement with those reported by previously published results. In addition, the SIFs for various crack geometry and loading type such as a curved crack under far-field tensile stress and a cracked body hanging by its own weight are also determined by the proposed BEM formulation. The numerical results obtained by this study are in agreement with those reported by previously published results.This paper presents the development of BEM procedure based on the maximum circumferential stress criterion for predicting the crack initiation directions and propagation paths in isotropic and anisotropic materials under mixed-mode loading. Good agreements are found between crack initiation and propagation predicted with the BEM and experimental observations reported by previous researchers of isotropic materials. Numerical simulations of crack initiation and propagation in CSTBD specimens of the anisotropic rock are also found to compare well with experimental results. Additionally, the crack propagation simulation technique is used to apply for analyzing the rock slope within a preexisting tensile crack. We can find that the slope failure surfaces strongly depend on the different anisotropy inclination angles.



---

*Source: 581493-2012-06-12.xml*