# Analysis of Crack Propagation Path on the Anisotropic Bi-Material Rock

Authors: Chao-Shi Chen; Chia-Huei Tu; Chen-Cheng Yang
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

This paper presents a single-domain boundary element method (SDBEM) for linear elastic fracture mechanics analysis in the 2D anisotropic bimaterial. In this formulation, the displacement integral equation is collocated on the uncracked boundary only, and the traction integral equation is collocated on one side of the crack surface only. The complete fundamental solution (Green's function) for anisotropic bi-materials was also derived and implemented into the boundary integral formulation so the discretization along the interface can be avoided except for the interfacial crack part. A special crack-tip element was introduced to capture exactly the crack-tip behavior. A computer program with the FORTRAN code has been developed to effectively calculate the stress intensity factors, crack initiation angle, and propagation path of an anisotropic bi-material. This SDBEM program has been verified having a good accuracy with the previous researches. In addition, a rock of type (1)/(2) disk specimen with a central crack was made to conduct the Brazilian test under diametrical loading. The result shows that the numerical analysis can predict relatively well the direction of crack initiation and the path of crack propagation.


---

## Body

## 1. Introduction

In rock masses, the interbed construction suffers from cracking which is caused by various factors. Of greater concern are those cracks that develop as a result of initiation and propagation path, leading to significant change in the failure resistance of the structure. Many relative researches and discussions are started like wildfire and never stopped.Because of the discontinuities of weak interbed on laccoliths (rock mass), there are many damages that occurred which are due to earth stress efforts or under geotechnical engineering, such as slope slip, tunnel collapse, deep excavation collapse, and crack openings which are due to deep well drilling. Crack is one of the fracture models to cause those damages because the crack openings and propagation on the field affect the rock mass structure stability. This paper will discuss the behavior of crack propagation on the basis of the theory of fracture mechanics. We defined the rock sample on interbed as “Bimaterial rock”.Theoretically, interfacial crack problems in isotropic Bimaterials were studied [1–3] where the authors showed that the stresses possess the singularity of γ-1/2±iε. Rice [4] re-examined the elastic fracture mechanics concepts for the isotropic interfacial crack and introduced an intrinsic material length scale so that the definition of the stress intensity factor (SIFs) possesses the same physical significance as those for the homogeneous cracks. Clements [5] and Willis  [6] studied interfacial crack problems in dissimilar anisotropic materials. They showed that the oscillatory behavior of the stresses and the phenomenon of interpenetrating of the crack faces were also present near the crack tip for anisotropic interface cracks. Recent studies on interfacial cracks in anisotropic materials have been conducted by many authors [7–21], and different definitions for the stress intensity factor exist. By introducing a characteristic length, however, the definition given by Gao et al. [15], and Wu [12, 13] is consistent with Rice’s general definition [4] and appears to be more explicit than other definitions.The study of fracture initiation and propagation in anisotropic rocks is subjected to Brazilian loads. A numerical procedure based on the SDBEM and maximum tensile stress criterion has been developed to predict the angle of crack initiation and the path of crack propagation in anisotropic rocks. Crack propagation in an anisotropic homogeneous rock disc under mixed mode I-II loading is simulated by an incremental crack extension with a piecewise linear discretization. A computer program, which can automatically generate a new mesh, has been developed to simulate the crack propagation process. Some experimental observations of crack initiation angles and crack propagation were obtained by conducting diametrical loading of initially cracked discs of a gypsum/cement. It was found that the numerical analysis could predict relatively well the direction of crack initiation and the crack propagation path.


## 2. Theoretical Background

### 2.1. Green’s Function in Anisotropic Bimaterials

With three complex analytical functionsϕi(zi), one can, in general, express displacements, stresses, and tractions as follows [11, 22, 23]:
(2.1)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, Re denotes the real part of a complex variable or function, a prime denotes the derivative, the three complex numbers μj(j=1,2,3) and the elements of the complex matrices B and Aare functions of the elastic properties [11, 22, 23].Assume that the medium is composed of two joined dissimilar anisotropic and elastic half-planes. We let the interface be along the x-axis and the upper(y>0) and lower (y<0) half-planes occupied by materials (1) and (2), respectively (Figure 1).Figure 1
Definition of the coordinate systems within an anisotropic Bimaterial.For concentrated force acting at the point (x0,y0) in material (2) (y0<0), we express the complex vector function as [11]
(2.2)ϕ(z)={ϕU(z),z∈1,ϕL(z)+ϕ(2)0(z),z∈2,
where the vector function
(2.3)ϕ(z)=[ϕ1(z),ϕ2(z),ϕ3(z)]T,
with the argument having the generic form z=x+μy.In (2.2), ϕ(2)0 is a singular solution corresponding to a point force acting at the point (x0,y0) in an anisotropic infinite plane with the elastic properties of material (2). This singular solution can be expressed as [11, 23]
(2.4)ϕk(zk)=-12π[Hk1P1ln(zk-sk)+Hk2P2ln(zk-sk)],
where sk=x0+μky0, Pk is the point force vector, and H is given by
(2.5)H=A-1(Y-1+Y̅-1)-1;Y=iAB-1.There are two unknown vector functions to be solved in (2.2), that is, ϕU(z) and ϕL(z). While the former is analytic in the upper (material (1)) half-plane, the latter is analytic in the lower (material (2)) half-plane. These expressions can be found by requiring continuity of the resultant traction and displacement across the interface, along with the standard analytic continuation arguments. Following this approach and after some complex algebraic manipulation, the complex vector functions in materials (1) and (2) are obtained as
(2.6)ϕ(z)={B(1)-1(Y(1)+Y̅(2))-1(Y̅(2)+Y(2))B(2)ϕ(2)0(z),z∈1,B(2)-1(Y̅(1)+Y(2))-1(Y̅(2)-Y̅(1))B̅(2)ϕ̅(2)0(z)+ϕ(2)0(z),z∈2.In (2.6), the special subscripts (1) and (2) are used exclusively to denote that the corresponding matrix or vector is in material (1) (y>0) and material (2) (y<0), respectively.Similarly, for a point force in material (2.1) (y0>0), these complex functions can be found as
(2.7)ϕ(z)={B(1)-1(Y̅(2)+Y(1))-1(Y̅(1)-Y̅(2))B̅(1)ϕ̅(1)0(z)+ϕ(1)0(z),z∈1,B(2)-1(Y(2)+Y̅(1))-1(Y̅(1)+Y(1))B(1)ϕ(1)0(z),z∈2,
where the vector functions ϕ(1)0 are the infinite plane solution given in (2.4) but with the elastic properties of material (1).With the complex function given in (2.6) and (2.7), Green’s functions of the displacement and traction can be obtained by substituting these complex functions into (2.1). These Green’s functions have four different forms depending on the relative location of the field and source points.It is noteworthy that these Green’s function can be used to solve both plane stress and plane strain problems in anisotropic Bimaterials. Although the isotropic solution cannot be analytically reduced from these Green’s functions one can numerically approximate it by selecting a very weak anisotropic (or nearly isotropic) medium [24, 25].


### 2.2. BEM Formulation for 2D Cracked Anisotropic Bimaterials

In this section, we present an SDBEM formulation in which neither the artificial boundary nor the discretization along the un-cracked interface is necessary. This SDBEM formulation was widely used recently by Chen et al. [26], for homogeneous materials and is now extended to anisotropic Bimaterials.The displacement integral equation applied to the outer boundary results in the following form (sk,B0∈ΓB only, Figure 2):
(2.8)bij(sk,b0)+uj(sk,b0)+∫BTij*(zk,B,sk,B0)uj(zk,B)dΓ(zk,B)+∫CTij*(zk,C,sk,B0)[uj(zk,C+)-uj(zk,C-)]dΓ(zk,C)=∫BUij*(zk,B,sk,B0)tj(zk,B)dΓ(zk,B),
where i,j,k = 1,2, Tij* and Uij* are Green’s tractions and displacements, uj and tj are the boundary displacement and tractions, respectively, bij are quantities that depend on the geometry of the boundary and are equal to δij/2 for a smooth boundary, and zk and sk0 are the field points and the source points on the boundary Γ of the domain. ΓC has the same outward normal as ΓC+. Here, the subscripts B and C denote the outer boundary and the crack surface, respectively.Figure 2
Geometry of a two-dimensional cracked domain.The traction integral equation (forsk0 being a smooth point on the crack) applied to one side of the crack surfaces is (sk,C0∈ΓC+ only)
(2.9)0.5tj(sk,C0)+nm(sk,C0)∫BClmikTij,k*(sk,C0,zk,B)uj(zk,B)dΓ(zk,B)+nm(sk,C0)∫CClmikTij,k*(sk,C0,zk,C)[uj(zk,C+)-uj(zk,C-)]dΓ(zk,C)=nm(sk,C0)∫BClmikUij,k*×(sk,C0,zk,B)tj(zk,B)dΓ(zk,B),
where nm is the outward normal at the crack surface sk,C0 and Clmik is the fourth-order stiffness tensor.Equations (2.8) and (2.9) form a pair of boundary integral equations [24, 27, 28] and can be used for the calculation of SIFs in anisotropic Bimaterials. The main feature of the BEM formulation is that it is a single-domain formulation with displacement integral (2.8) being collocated on the un-cracked boundary only and traction integral (2.9) on one side of the crack surface only. For problems without cracks, one needs (2.8) only, with the integral on the crack surface being discarded. Equation (2.8) then reduces to the well-known displacement integral on the un-cracked boundary being discarded. We emphasize here that since bimaterial Green’s functions are included in (2.8), discretization along the interface can be avoided, with the exception of the interfacial crack part which will be treated by the traction integral equation presented by (2.9).It is well known that a cracked domain poses certain difficulties for BEM modeling (Cruse, 1988 [29]). Previously, fracture mechanics problems in isotropic or anisotropic bimaterials were mostly handled by the multidomain method in which each side of the crack surface is put into different domains and artificial boundaries are introduced to connect the crack surface to the un-cracked boundary. For the bimaterial case, discretization along the interface is also required if one uses the Kelvin-type (infinite domain) Green’s functions.


### 2.3. Crack-Tip Modeling

In fracture mechanics analysis, especially in the calculation of the SIFs, one needs to know the asymptotic behavior of the displacements and stresses near the crack-tip. In our BEM analysis of the SIFs, we propose to use the extrapolation method of the crack tip displacements. We therefore need to know the exact asymptotic behavior of the relative crack displacement behind the cracktip. This asymptotic expression has different forms depending on the location of the cracktip. In this paper, two cases will be discussed, that is, a crack-tip within the homogeneous material and an interfacial crack-tip. Inclined cracks terminating at the interface will be discussed in a future paper.


#### 2.3.1. A Crack Tip within a Homogeneous Material

The mixed mode stress intensity factors (SIFs) for anisotropic media can be determined by using the extrapolation method of the relative crack displacement (RCD), combined with a set of the shape functions. The RCD is defined as [24]
(2.10)Δui=∑k=13ϕkΔuik.For this case, the relation of the RCDs at a distancer behind the crack tip and the SIFs can be found as [24, 25, 30]
(2.11)Δu1=22rπ(H11KI+H12KII),Δu2=22rπ(H21KI+H22KII),
where
(2.12)H11=Im(μ2P11-μ1P12μ1-μ2);H12=Im(P11-P12μ1-μ2),H21=Im(μ2P21-μ1P22μ1-μ2);H22=Im(P21-P22μ1-μ2).Substituting the RCDs into (2.10) and (2.11), we obtain a set of algebraic equations in which the SIFs KI and KII can be solved.


#### 2.3.2. An Interfacial Crack Tip

For this case, the relative crack displacements at a distancer behind the interfacial crack tip can be expressed, in terms of the three SIFs, as [15]
(2.13)Δu(r)=(∑j=13cjDQje-πδjr1/2+iδj)K,
where cj,δj,Qj, and D are the relative parameters in material (1) and material (2). Utilizing (2.5), we defined the matrix of material as
(2.14)Y(1)+Y̅(2)=D-iV,
where D and V are two real matrices, and then utilizing these two matrices, we define matrix P as
(2.15)P=-D-1V,
and the characteristic β relative to material
(2.16)β=-12tr(P2).
Then, we used the characteristic β obtain to define oscillation index ε as
(2.17)ε=12πln1+β1-β=1πtanh-1β,δ1=0,δ2=ε,δ3=-ε,Q1=P2+β2I,Q2=P(P-iβI),Q3=P(P+iβI),
where I is a 3×3 identity matrix.Utilizing the relationship between the characteristicβ and oscillation index ε, we define constant cj as
(2.18)c1=22πβ2,c2=-e-πεdiε2π(1+2iε)β2cosh(πε),c3=-e-πεdiε2π(1-2iε)β2cosh(πε),
where d is the characteristic distance along material interface to crack tip.Comparing (2.10) with (2.13), we noticed that while the relative crack displacement behaves as a square root for a crack tip within a homogeneous medium, for an interfacial crack tip, its behavior is r1/2+iδ, a square-root feature multiplied by weak oscillatory behaviors.Equation (2.13) can be recast into the following form, which is more convenient for the current numerical applications:
(2.19)Δu(r)=2rπM(rd)K,
where d is the characteristic length and M is a matrix function with its expression given by
(2.20)M(x)=Dβ2{(P2+β2I)-[cos(εlnx)+2εsin(εlnx)]P2+β[sin(εlnx)-2εcos(εlnx)]P(1+4ε2)cosh(πε)}.
Again, in order to capture the square-root and the weak oscillatory behavior, we construct a crack-tip element with tip at ξ=-1 in terms of which the relative crack displacement is expressed as
(2.21)Δu(r)=M(rd)[f1Δu11f2Δu12f3Δu13f1Δu21f2Δu22f3Δu23f1Δu31f2Δu32f3Δu33].In meshing the 2D anisotropic bimaterial problem (as shown in Figure3), we assume that the interface is along the X-axis and the upper (y>0) and lower (y<0) half-planes are occupied by materials (1) and (2), respectively. The corner of outer boundary is processed by the discontinuous elements Type I and Type III; the continuous element Type V is to deal with all outer smooth boundary; internal crack surface is processed by crack surface elements Type II; and crack tip element Type VI is to process crack tip problem. In order to avoid the oscillatory behavior of the interface, we mesh an anisotropic problem (as shown in Figure 3(6)) of bimaterial in which neither the interfacial elements nor the discretization along the un-cracked interface is necessary, with the exception of the interfacial crack part which will be treated by traction integral Equation (2.9) and (2.10) of RCD.Figure 3
Mesh 2D bimaterial problem with the five types of quadratic elements.


### 2.4. Crack Initiation and Fracture Propagation

In fracture mechanics, there are three criteria commonly used to predict the crack initiation angle: the maximum tensile stress criterion, orσ-criterion [31], the maximum energy release rate criterion, or G-criterion [32], and the minimum strain energy density criterion, or S-criterion [33]. Among them, the σ-criterion has been found to predict well the directions of crack initiation compared to the experimental results for polymethylmethacrylate [34, 35] and brittle clay [36]. Because of its simplicity, the σ-criterion seems to be the most popular criterion in mixed mode I-II fracture studies [37]. Therefore, the σ-criterion was used in this paper to determine the crack initiation angle for anisotropic plates.For anisotropic materials, the general form of the elastic stress field near the crack tip in the local Cartesian coordinatesx′′,y′′ of Figure 4 can be expressed in terms of the two stress intensity factors KI and KII a follows [30]:
(2.22)σx′′=KI2πrRe[μ1μ2μ1-μ2(μ2cosθ+μ2sinθ-μ1cosθ+μ1sinθ)]+KII2πr×Re[1μ1-μ2(μ22cosθ+μ2sinθ-μ12cosθ+μ1sinθ)],σy′′=KI2πrRe[1μ1-μ2(μ1cosθ+μ2sinθ-μ2cosθ+μ1sinθ)]+KII2πr×Re[1μ1-μ2(1cosθ+μ2sinθ-1cosθ+μ1sinθ)],τx′′y′′=KI2πrRe[μ1μ2μ1-μ2(1cosθ+μ1sinθ-1cosθ+μ2sinθ)]+KII2πr×Re[1μ1-μ2(μ1cosθ+μ1sinθ-μ2cosθ+μ2sinθ)].Figure 4
Crack tip coordinate system and stress components.Using coordinate transformation, the stress fields near the crack tip in the polar coordinates (r,θ) of Figure 4 are
(2.23)σθ=σx′′+σy′′2-σx′′+σy′′2cos2θ-τx′′y′′sin2θ,τθ=-σx′′+σy′′2sin2θ+τx′′y′′cos2θ.If the maximumσ-criterion is used, the angle of crack initiation, θ0, must satisfy
(2.24)∂σθ∂θ=0(orτrθ=0),∂2σθ∂θ2<0.A numerical procedure was applied to find the angleθ0 when σθ is a maximum for known values of the material elastic constants, the anisotropic orientation angle ψ, and the crack geometry.In this paper, the process of crack propagation in an anisotropic homogeneous plate under mixed mode I-II loading is simulated by incremental crack extension with a piecewise linear discretization. For each incremental analysis, crack extension is conveniently modeled by a new boundary element. A computer program has been developed to automatically generate new data required for analyzing sequentially 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 (Figure5):(1)
compute the SIFs using the proposed BEM;(2)
determine the angle of crack initiation based on the maximum tensile stress criterion;(3)
extend the crack by a linear element (of length selected by the user) along the direction determined in step2;(4)
automatically generate the new BEM mesh;(5)
repeat all the above steps until the crack is near the outer boundary.Figure 5
Process of crack propagation by increasing the number of linear elements.


## 2.1. Green’s Function in Anisotropic Bimaterials

With three complex analytical functionsϕi(zi), one can, in general, express displacements, stresses, and tractions as follows [11, 22, 23]:
(2.1)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, Re denotes the real part of a complex variable or function, a prime denotes the derivative, the three complex numbers μj(j=1,2,3) and the elements of the complex matrices B and Aare functions of the elastic properties [11, 22, 23].Assume that the medium is composed of two joined dissimilar anisotropic and elastic half-planes. We let the interface be along the x-axis and the upper(y>0) and lower (y<0) half-planes occupied by materials (1) and (2), respectively (Figure 1).Figure 1
Definition of the coordinate systems within an anisotropic Bimaterial.For concentrated force acting at the point (x0,y0) in material (2) (y0<0), we express the complex vector function as [11]
(2.2)ϕ(z)={ϕU(z),z∈1,ϕL(z)+ϕ(2)0(z),z∈2,
where the vector function
(2.3)ϕ(z)=[ϕ1(z),ϕ2(z),ϕ3(z)]T,
with the argument having the generic form z=x+μy.In (2.2), ϕ(2)0 is a singular solution corresponding to a point force acting at the point (x0,y0) in an anisotropic infinite plane with the elastic properties of material (2). This singular solution can be expressed as [11, 23]
(2.4)ϕk(zk)=-12π[Hk1P1ln(zk-sk)+Hk2P2ln(zk-sk)],
where sk=x0+μky0, Pk is the point force vector, and H is given by
(2.5)H=A-1(Y-1+Y̅-1)-1;Y=iAB-1.There are two unknown vector functions to be solved in (2.2), that is, ϕU(z) and ϕL(z). While the former is analytic in the upper (material (1)) half-plane, the latter is analytic in the lower (material (2)) half-plane. These expressions can be found by requiring continuity of the resultant traction and displacement across the interface, along with the standard analytic continuation arguments. Following this approach and after some complex algebraic manipulation, the complex vector functions in materials (1) and (2) are obtained as
(2.6)ϕ(z)={B(1)-1(Y(1)+Y̅(2))-1(Y̅(2)+Y(2))B(2)ϕ(2)0(z),z∈1,B(2)-1(Y̅(1)+Y(2))-1(Y̅(2)-Y̅(1))B̅(2)ϕ̅(2)0(z)+ϕ(2)0(z),z∈2.In (2.6), the special subscripts (1) and (2) are used exclusively to denote that the corresponding matrix or vector is in material (1) (y>0) and material (2) (y<0), respectively.Similarly, for a point force in material (2.1) (y0>0), these complex functions can be found as
(2.7)ϕ(z)={B(1)-1(Y̅(2)+Y(1))-1(Y̅(1)-Y̅(2))B̅(1)ϕ̅(1)0(z)+ϕ(1)0(z),z∈1,B(2)-1(Y(2)+Y̅(1))-1(Y̅(1)+Y(1))B(1)ϕ(1)0(z),z∈2,
where the vector functions ϕ(1)0 are the infinite plane solution given in (2.4) but with the elastic properties of material (1).With the complex function given in (2.6) and (2.7), Green’s functions of the displacement and traction can be obtained by substituting these complex functions into (2.1). These Green’s functions have four different forms depending on the relative location of the field and source points.It is noteworthy that these Green’s function can be used to solve both plane stress and plane strain problems in anisotropic Bimaterials. Although the isotropic solution cannot be analytically reduced from these Green’s functions one can numerically approximate it by selecting a very weak anisotropic (or nearly isotropic) medium [24, 25].


## 2.2. BEM Formulation for 2D Cracked Anisotropic Bimaterials

In this section, we present an SDBEM formulation in which neither the artificial boundary nor the discretization along the un-cracked interface is necessary. This SDBEM formulation was widely used recently by Chen et al. [26], for homogeneous materials and is now extended to anisotropic Bimaterials.The displacement integral equation applied to the outer boundary results in the following form (sk,B0∈ΓB only, Figure 2):
(2.8)bij(sk,b0)+uj(sk,b0)+∫BTij*(zk,B,sk,B0)uj(zk,B)dΓ(zk,B)+∫CTij*(zk,C,sk,B0)[uj(zk,C+)-uj(zk,C-)]dΓ(zk,C)=∫BUij*(zk,B,sk,B0)tj(zk,B)dΓ(zk,B),
where i,j,k = 1,2, Tij* and Uij* are Green’s tractions and displacements, uj and tj are the boundary displacement and tractions, respectively, bij are quantities that depend on the geometry of the boundary and are equal to δij/2 for a smooth boundary, and zk and sk0 are the field points and the source points on the boundary Γ of the domain. ΓC has the same outward normal as ΓC+. Here, the subscripts B and C denote the outer boundary and the crack surface, respectively.Figure 2
Geometry of a two-dimensional cracked domain.The traction integral equation (forsk0 being a smooth point on the crack) applied to one side of the crack surfaces is (sk,C0∈ΓC+ only)
(2.9)0.5tj(sk,C0)+nm(sk,C0)∫BClmikTij,k*(sk,C0,zk,B)uj(zk,B)dΓ(zk,B)+nm(sk,C0)∫CClmikTij,k*(sk,C0,zk,C)[uj(zk,C+)-uj(zk,C-)]dΓ(zk,C)=nm(sk,C0)∫BClmikUij,k*×(sk,C0,zk,B)tj(zk,B)dΓ(zk,B),
where nm is the outward normal at the crack surface sk,C0 and Clmik is the fourth-order stiffness tensor.Equations (2.8) and (2.9) form a pair of boundary integral equations [24, 27, 28] and can be used for the calculation of SIFs in anisotropic Bimaterials. The main feature of the BEM formulation is that it is a single-domain formulation with displacement integral (2.8) being collocated on the un-cracked boundary only and traction integral (2.9) on one side of the crack surface only. For problems without cracks, one needs (2.8) only, with the integral on the crack surface being discarded. Equation (2.8) then reduces to the well-known displacement integral on the un-cracked boundary being discarded. We emphasize here that since bimaterial Green’s functions are included in (2.8), discretization along the interface can be avoided, with the exception of the interfacial crack part which will be treated by the traction integral equation presented by (2.9).It is well known that a cracked domain poses certain difficulties for BEM modeling (Cruse, 1988 [29]). Previously, fracture mechanics problems in isotropic or anisotropic bimaterials were mostly handled by the multidomain method in which each side of the crack surface is put into different domains and artificial boundaries are introduced to connect the crack surface to the un-cracked boundary. For the bimaterial case, discretization along the interface is also required if one uses the Kelvin-type (infinite domain) Green’s functions.


## 2.3. Crack-Tip Modeling

In fracture mechanics analysis, especially in the calculation of the SIFs, one needs to know the asymptotic behavior of the displacements and stresses near the crack-tip. In our BEM analysis of the SIFs, we propose to use the extrapolation method of the crack tip displacements. We therefore need to know the exact asymptotic behavior of the relative crack displacement behind the cracktip. This asymptotic expression has different forms depending on the location of the cracktip. In this paper, two cases will be discussed, that is, a crack-tip within the homogeneous material and an interfacial crack-tip. Inclined cracks terminating at the interface will be discussed in a future paper.


### 2.3.1. A Crack Tip within a Homogeneous Material

The mixed mode stress intensity factors (SIFs) for anisotropic media can be determined by using the extrapolation method of the relative crack displacement (RCD), combined with a set of the shape functions. The RCD is defined as [24]
(2.10)Δui=∑k=13ϕkΔuik.For this case, the relation of the RCDs at a distancer behind the crack tip and the SIFs can be found as [24, 25, 30]
(2.11)Δu1=22rπ(H11KI+H12KII),Δu2=22rπ(H21KI+H22KII),
where
(2.12)H11=Im(μ2P11-μ1P12μ1-μ2);H12=Im(P11-P12μ1-μ2),H21=Im(μ2P21-μ1P22μ1-μ2);H22=Im(P21-P22μ1-μ2).Substituting the RCDs into (2.10) and (2.11), we obtain a set of algebraic equations in which the SIFs KI and KII can be solved.


### 2.3.2. An Interfacial Crack Tip

For this case, the relative crack displacements at a distancer behind the interfacial crack tip can be expressed, in terms of the three SIFs, as [15]
(2.13)Δu(r)=(∑j=13cjDQje-πδjr1/2+iδj)K,
where cj,δj,Qj, and D are the relative parameters in material (1) and material (2). Utilizing (2.5), we defined the matrix of material as
(2.14)Y(1)+Y̅(2)=D-iV,
where D and V are two real matrices, and then utilizing these two matrices, we define matrix P as
(2.15)P=-D-1V,
and the characteristic β relative to material
(2.16)β=-12tr(P2).
Then, we used the characteristic β obtain to define oscillation index ε as
(2.17)ε=12πln1+β1-β=1πtanh-1β,δ1=0,δ2=ε,δ3=-ε,Q1=P2+β2I,Q2=P(P-iβI),Q3=P(P+iβI),
where I is a 3×3 identity matrix.Utilizing the relationship between the characteristicβ and oscillation index ε, we define constant cj as
(2.18)c1=22πβ2,c2=-e-πεdiε2π(1+2iε)β2cosh(πε),c3=-e-πεdiε2π(1-2iε)β2cosh(πε),
where d is the characteristic distance along material interface to crack tip.Comparing (2.10) with (2.13), we noticed that while the relative crack displacement behaves as a square root for a crack tip within a homogeneous medium, for an interfacial crack tip, its behavior is r1/2+iδ, a square-root feature multiplied by weak oscillatory behaviors.Equation (2.13) can be recast into the following form, which is more convenient for the current numerical applications:
(2.19)Δu(r)=2rπM(rd)K,
where d is the characteristic length and M is a matrix function with its expression given by
(2.20)M(x)=Dβ2{(P2+β2I)-[cos(εlnx)+2εsin(εlnx)]P2+β[sin(εlnx)-2εcos(εlnx)]P(1+4ε2)cosh(πε)}.
Again, in order to capture the square-root and the weak oscillatory behavior, we construct a crack-tip element with tip at ξ=-1 in terms of which the relative crack displacement is expressed as
(2.21)Δu(r)=M(rd)[f1Δu11f2Δu12f3Δu13f1Δu21f2Δu22f3Δu23f1Δu31f2Δu32f3Δu33].In meshing the 2D anisotropic bimaterial problem (as shown in Figure3), we assume that the interface is along the X-axis and the upper (y>0) and lower (y<0) half-planes are occupied by materials (1) and (2), respectively. The corner of outer boundary is processed by the discontinuous elements Type I and Type III; the continuous element Type V is to deal with all outer smooth boundary; internal crack surface is processed by crack surface elements Type II; and crack tip element Type VI is to process crack tip problem. In order to avoid the oscillatory behavior of the interface, we mesh an anisotropic problem (as shown in Figure 3(6)) of bimaterial in which neither the interfacial elements nor the discretization along the un-cracked interface is necessary, with the exception of the interfacial crack part which will be treated by traction integral Equation (2.9) and (2.10) of RCD.Figure 3
Mesh 2D bimaterial problem with the five types of quadratic elements.


## 2.3.1. A Crack Tip within a Homogeneous Material

The mixed mode stress intensity factors (SIFs) for anisotropic media can be determined by using the extrapolation method of the relative crack displacement (RCD), combined with a set of the shape functions. The RCD is defined as [24]
(2.10)Δui=∑k=13ϕkΔuik.For this case, the relation of the RCDs at a distancer behind the crack tip and the SIFs can be found as [24, 25, 30]
(2.11)Δu1=22rπ(H11KI+H12KII),Δu2=22rπ(H21KI+H22KII),
where
(2.12)H11=Im(μ2P11-μ1P12μ1-μ2);H12=Im(P11-P12μ1-μ2),H21=Im(μ2P21-μ1P22μ1-μ2);H22=Im(P21-P22μ1-μ2).Substituting the RCDs into (2.10) and (2.11), we obtain a set of algebraic equations in which the SIFs KI and KII can be solved.


## 2.3.2. An Interfacial Crack Tip

For this case, the relative crack displacements at a distancer behind the interfacial crack tip can be expressed, in terms of the three SIFs, as [15]
(2.13)Δu(r)=(∑j=13cjDQje-πδjr1/2+iδj)K,
where cj,δj,Qj, and D are the relative parameters in material (1) and material (2). Utilizing (2.5), we defined the matrix of material as
(2.14)Y(1)+Y̅(2)=D-iV,
where D and V are two real matrices, and then utilizing these two matrices, we define matrix P as
(2.15)P=-D-1V,
and the characteristic β relative to material
(2.16)β=-12tr(P2).
Then, we used the characteristic β obtain to define oscillation index ε as
(2.17)ε=12πln1+β1-β=1πtanh-1β,δ1=0,δ2=ε,δ3=-ε,Q1=P2+β2I,Q2=P(P-iβI),Q3=P(P+iβI),
where I is a 3×3 identity matrix.Utilizing the relationship between the characteristicβ and oscillation index ε, we define constant cj as
(2.18)c1=22πβ2,c2=-e-πεdiε2π(1+2iε)β2cosh(πε),c3=-e-πεdiε2π(1-2iε)β2cosh(πε),
where d is the characteristic distance along material interface to crack tip.Comparing (2.10) with (2.13), we noticed that while the relative crack displacement behaves as a square root for a crack tip within a homogeneous medium, for an interfacial crack tip, its behavior is r1/2+iδ, a square-root feature multiplied by weak oscillatory behaviors.Equation (2.13) can be recast into the following form, which is more convenient for the current numerical applications:
(2.19)Δu(r)=2rπM(rd)K,
where d is the characteristic length and M is a matrix function with its expression given by
(2.20)M(x)=Dβ2{(P2+β2I)-[cos(εlnx)+2εsin(εlnx)]P2+β[sin(εlnx)-2εcos(εlnx)]P(1+4ε2)cosh(πε)}.
Again, in order to capture the square-root and the weak oscillatory behavior, we construct a crack-tip element with tip at ξ=-1 in terms of which the relative crack displacement is expressed as
(2.21)Δu(r)=M(rd)[f1Δu11f2Δu12f3Δu13f1Δu21f2Δu22f3Δu23f1Δu31f2Δu32f3Δu33].In meshing the 2D anisotropic bimaterial problem (as shown in Figure3), we assume that the interface is along the X-axis and the upper (y>0) and lower (y<0) half-planes are occupied by materials (1) and (2), respectively. The corner of outer boundary is processed by the discontinuous elements Type I and Type III; the continuous element Type V is to deal with all outer smooth boundary; internal crack surface is processed by crack surface elements Type II; and crack tip element Type VI is to process crack tip problem. In order to avoid the oscillatory behavior of the interface, we mesh an anisotropic problem (as shown in Figure 3(6)) of bimaterial in which neither the interfacial elements nor the discretization along the un-cracked interface is necessary, with the exception of the interfacial crack part which will be treated by traction integral Equation (2.9) and (2.10) of RCD.Figure 3
Mesh 2D bimaterial problem with the five types of quadratic elements.


## 2.4. Crack Initiation and Fracture Propagation

In fracture mechanics, there are three criteria commonly used to predict the crack initiation angle: the maximum tensile stress criterion, orσ-criterion [31], the maximum energy release rate criterion, or G-criterion [32], and the minimum strain energy density criterion, or S-criterion [33]. Among them, the σ-criterion has been found to predict well the directions of crack initiation compared to the experimental results for polymethylmethacrylate [34, 35] and brittle clay [36]. Because of its simplicity, the σ-criterion seems to be the most popular criterion in mixed mode I-II fracture studies [37]. Therefore, the σ-criterion was used in this paper to determine the crack initiation angle for anisotropic plates.For anisotropic materials, the general form of the elastic stress field near the crack tip in the local Cartesian coordinatesx′′,y′′ of Figure 4 can be expressed in terms of the two stress intensity factors KI and KII a follows [30]:
(2.22)σx′′=KI2πrRe[μ1μ2μ1-μ2(μ2cosθ+μ2sinθ-μ1cosθ+μ1sinθ)]+KII2πr×Re[1μ1-μ2(μ22cosθ+μ2sinθ-μ12cosθ+μ1sinθ)],σy′′=KI2πrRe[1μ1-μ2(μ1cosθ+μ2sinθ-μ2cosθ+μ1sinθ)]+KII2πr×Re[1μ1-μ2(1cosθ+μ2sinθ-1cosθ+μ1sinθ)],τx′′y′′=KI2πrRe[μ1μ2μ1-μ2(1cosθ+μ1sinθ-1cosθ+μ2sinθ)]+KII2πr×Re[1μ1-μ2(μ1cosθ+μ1sinθ-μ2cosθ+μ2sinθ)].Figure 4
Crack tip coordinate system and stress components.Using coordinate transformation, the stress fields near the crack tip in the polar coordinates (r,θ) of Figure 4 are
(2.23)σθ=σx′′+σy′′2-σx′′+σy′′2cos2θ-τx′′y′′sin2θ,τθ=-σx′′+σy′′2sin2θ+τx′′y′′cos2θ.If the maximumσ-criterion is used, the angle of crack initiation, θ0, must satisfy
(2.24)∂σθ∂θ=0(orτrθ=0),∂2σθ∂θ2<0.A numerical procedure was applied to find the angleθ0 when σθ is a maximum for known values of the material elastic constants, the anisotropic orientation angle ψ, and the crack geometry.In this paper, the process of crack propagation in an anisotropic homogeneous plate under mixed mode I-II loading is simulated by incremental crack extension with a piecewise linear discretization. For each incremental analysis, crack extension is conveniently modeled by a new boundary element. A computer program has been developed to automatically generate new data required for analyzing sequentially 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 (Figure5):(1)
compute the SIFs using the proposed BEM;(2)
determine the angle of crack initiation based on the maximum tensile stress criterion;(3)
extend the crack by a linear element (of length selected by the user) along the direction determined in step2;(4)
automatically generate the new BEM mesh;(5)
repeat all the above steps until the crack is near the outer boundary.Figure 5
Process of crack propagation by increasing the number of linear elements.


## 3. Numerical Examples

In this section, a lot of numerical examples are presented to verify the formulation and to show the efficiency and versatility of the present SDBEM for problems related to fracture in anisotropic Bimaterial.


### 3.1. Stress Intensity Factors (SIFs)

#### 3.1.1. Horizontal Crack in Material (1)

A horizontal crack under a uniform pressureP is shown in Figure 6. The crack has a length 2a, and is located at a distance d to the interface. The Poisson ratios for both materials (1) and (2) the same, that is, ν1 = ν2 = 0.3, while the ratio of the shear module ratio G2/G1 varies. A plane stress condition is assumed. In order to calculate the SIFs at crack tip (A) or (B), 20 quadratic elements were used to discretize the crack surface. The results are given in Table 1 for various values of the shear module ratio. They are compared to the results given by Isida and Noguchi [38], using a body force integral equation method and those by Yuuki and Cho [39], using a multidomain BEM formulation. As can be observed from this table, the results compare quite well.Table 1
Comparison of the SIFs (horizontal crack).
G2/G1d/2aPresentIsida et al. [40]Diff. (%)Yuuki et al. [39]Diff. (%)KI/(pπa) of Tip (A) or (B)0.250.051.4761.468-0.571.468-0.540.250.51.1981.197-0.091.197-0.122.00.050.8710.8720.140.869-0.172.00.50.9360.935-0.060.934-0.16KII/(pπa) of Tip (A) or (B)0.250.050.2850.2860.350.2922.500.250.50.0710.0710.700.0721.672.00.05-0.088-0.087-1.38-0.085-4.012.00.5-0.023-0.0242.50-0.023-3.54Figure 6
A horizontal crack under uniform pressure within material (1) of an infinite Bimaterial.


#### 3.1.2. Vertical Crack in Material (1)

Consider a vertical crack in material (1) subjected to far-field horizontal stresses as shown in Figure 7. The crack has a length 2a and is located at a distance d to the interface. The Poisson ratios for both materials (1) and (2) are the same, that is, ν1 = 0.35, ν2 = 0.3, and the ratio of the shear module is the same, that is, G1/G2 = 23.077, while the ratio of the crack length located at a distance (d/a) varies. A plane stress condition is assumed. In order to calculate the SIFs at crack tip (A) or (B), 20 quadratic elements were used to discretize the crack surface. The results are given in Table 2 for various values of the crack length located at a distance. They are compared to the results given by Isida and Noguchi [38] using a body force integral equation method and those by Cook and Erdogan [41] using a Wiener-Hopf technique and an asymptotic analysis. As can be observed from this table, the results compare quite well.Table 2
Comparison of the SIFs (vertical crack).
d/aPresentIsida et al. [40]Diff. (%)Cook et al. [41]Diff. (%)KI/(σπa) Of Tip (A)1.000.8820.8830.110.8830.142.000.9610.9620.100.9620.105.000.9950.993-0.170.993-0.1710.000.9960.9980.240.9980.24KI/(σπa) Of Tip (B)1.00—————2.000.9350.935-0.040.935-0.045.000.9960.991-0.430.991-0.4410.001.0030.998-0.510.998-0.51Figure 7
A vertical crack in material (1): an interface under far-field stresses.


#### 3.1.3. Vertical Crack Intersecting an Interface

Consider a vertical crack intersecting an interface and subjected to far-field horizontal stresses as shown in Figure8. The horizontal far field stresses applied in materials (1) and (2) are, respectively, σ1 and σ2(=σ1G2/G1). The Poisson ratio ν1 and ν2 are assumed to be equal to 0.3 and the shear modulus ratio G2/G1 is assumed to vary. The distance of the crack tip (A) and (B) to the interface are the same, that is, d1=d2= a, the half-length of the crack. Again, a plane stress condition is assumed and 20 quadratic elements were used to discretize the crack surface. The SIFs at the crack tips (A) and (B) are listed in Table 3 for several values of the shear modulus ratio, and are compare to those given by Isida and Noguchi [38]. Again, the results between the two numerical analyses compare quite well.Table 3
Comparison of the SIFs (vertical crack intersecting an interface).
G2/G1KI/(σπa) of Tip (A)KI/(σπa) of Tip (B)Isida and Noguchi  [38]PresentDiff. (%)Isida and Noguchi  [38]PresentDiff. (%)0.11.0621.0629-0.081.1531.1539-0.080.31.0151.0157-0.071.0641.06390.010.51.0001.0007-0.071.0281.02730.070.80.9970.9975-0.051.0061.00470.13Figure 8
A vertical crack intersecting an interface under far-field stresses.


#### 3.1.4. An Interface Kinked Crack in Infinite Bimaterials

Consider an interface kinked (45°) crack in infinite Bimaterials subjected to far-field tensile stresses as shown in Figure 9. The far-field tensile stresses applied in materials (1) and (2) are, respectively, σ1 and σ2(=σ1G2/G1). The kink crack length is 2a; the main crack length is 2c. Kink crack’s tip is (A). The Poisson ratios ν1 and ν2 are assumed to be equal to 0.3. The Young module in material (1) is 1 GPa, and according to the shear module ratio with material (1) and material (2) (G1/G2 = 0.25), we can get the Young module of material (2). Then according to equation (G=(E/2(1+ν))), we can get the relation between kink length and stress intensity factors. After numerical analysis, we compared with the results of Isida and Noguchi [38], using a body integral equation force method, as shown in Table 4, the results the two numerical analyses and found that they compare quite well.Table 4
Comparison of the SIFs (interfacial kinked crack).
a/cKI/(σπa) of Tip (A)KII/(σπa) of Tip (A)Isida and Noguchi  [38]PresentDiff. (%)Isida and Noguchi [38]PresentDiff. (%)0.20.7330.759-3.540.6310.6241.100.50.7080.730-3.100.6230.6121.761.00.6830.701-2.630.6100.6031.141.50.6690.685-2.390.6010.5941.16Figure 9
Interfacial kinked crack within infinite Bimaterials.


#### 3.1.5. Interfacial Horizontal Crack in Infinite Bimaterials

An interfacial crack along the x-axis of length2a is shown in Figure 10. The crack surface is under a uniform pressure P and the materials can be either isotropic or anisotropic. Twenty quadratic elements were used to discretize the crack, and the characteristic length is assumed as 2a.Figure 10
Interfacial crack within infinite Bimaterials.The SIFs at the crack tip of an interfacial crack were also calculated for the anisotropic Bimaterial case. The anisotropic elastic properties in material (1) were assumed to be those of glass/epoxy with E1 = 48.26 GPa, E2 = 17.24 GPa, G12 = 6.89 GPa, and υ12 = 0.29. For material (2), a graphite/epoxy with E1 = 114.8 GPa, E2 = 11.7 GPa, G12= 9.66 GPa, and υ12 = 0.21 was selected [42]. The material axis E1 in material (1) and material (2) makes angles ψ1 and ψ2, respectively, with respect to the horizontal direction (Figure 10). While the material axis E1 in material (2) was assumed to be along the horizontal direction  (i.e., ψ(2) = 0), the E1-axis in material (1) makes an angle ψ(1) with respect to the horizontal direction. The interfacial SIFs at crack tip (B) obtained by the present method are listed in Table 5 and compared to the exact solutions proposed by Wu [12]. A very good agreement is found between the numerical analysis and the exact solution.Table 5
Comparison of the SIFs for infinite anisotropic problemψ2=0° (interfacial crack).
ψ1KI/(pπa) of  Tip (A) or (B)KII/(pπa) of  Tip (A) or (B)Wu [12]PresentDiff. (%)Wu [12]PresentDiff. (%)01.00001.0053-0.53-0.0382-0.03810.26300.99681.0006-0.38-0.0349-0.0350-0.29450.99651.0001-0.36-0.0318-0.0319-0.31600.99711.0010-0.39-0.0290-0.0292-0.69901.00001.0054-0.54-0.0264-0.0265-0.38


#### 3.1.6. Interfacial Horizontal Crack in Finite Bimaterials

The example is included as a comparison with the literature in order to demonstrate the accuracy of SDBEM approach for an interfacial crack in anisotropic bimaterial plate. The geometry is that of rectangular plate and is shown in Figure11. For the comparison, crack length is taken as 2a = 2, h = 2w, and a/w = 0.4, and static tensile loading σ0 is applied on the upper and the lower boundary of the plate. Plane stress condition is assumed. The anisotropic elastic properties for materials (1) and (2) are given in Table 6. The normalized complex stress intensity factors at crack tip (A) or (B) are listed in Table 7 together with those from the work of Cho et al. [43], who used a multidomain BEM formulation and the results from Wünsche et al. [21] for a finite body. The outer boundary and interfacial crack surface were discretized with 40 continuous and 20 discontinuous quadratic elements, respectively. It is obvious from Table 7 that these are very close to those obtained by the other researchers.Table 6
Elastic properties for materials (1) and (2).
MaterialsE1 (MPa)E2 (MPa)ν12G12Material (1)100500.310.009Material (2)(i)100450.39.525(ii)100400.39.010(iii)100300.37.860(iv)100100.34.630Table 7
Comparison of the normalized complex SIFs for finite anisotropic problem (interfacial crack).
Material (2) E2/E1|K|/σ0πa of Tip (A) or (B)Cho et al. [43]Wünsche et al. [21]Diff. (%)PresentDiff. (%)(i) 0.451.3171.3120.381.31320.29(ii) 0.401.3371.3330.301.33510.14(iii) 0.301.3921.3860.431.3922-0.01(iv) 0.101.6971.6890.471.6968-0.06*|K|=KI2+KII2Figure 11
An interfacial crack within finite rectangular plate of Bimaterial.


### 3.2. Crack Initiation Angles

#### 3.2.1. Comparison of Numerical Predictions of Crack Initiation Angles with Experimental Results

In this section, we compared the numerical predictions of crack initiation angles with experimental results to verify the formulation. Al-Shayea [44] conducted uniaxial pressure test on limestone rock discs 98 mm and 84 mm in diameter and 22 mm in thickness with 30 mm notch.  The crack orientation angle β between the crack plane and the tensile stress was varied. Figure 12 shows the variation of the crack initiation angle θ with the crack angle β determined numerically and experimentally. A good agreement is found between the experimental results of Al-Shayea [44] and our numerical predictions.Figure 12
A variation of crack initiation angleθ with the crack angle β.


## 3.1. Stress Intensity Factors (SIFs)

### 3.1.1. Horizontal Crack in Material (1)

A horizontal crack under a uniform pressureP is shown in Figure 6. The crack has a length 2a, and is located at a distance d to the interface. The Poisson ratios for both materials (1) and (2) the same, that is, ν1 = ν2 = 0.3, while the ratio of the shear module ratio G2/G1 varies. A plane stress condition is assumed. In order to calculate the SIFs at crack tip (A) or (B), 20 quadratic elements were used to discretize the crack surface. The results are given in Table 1 for various values of the shear module ratio. They are compared to the results given by Isida and Noguchi [38], using a body force integral equation method and those by Yuuki and Cho [39], using a multidomain BEM formulation. As can be observed from this table, the results compare quite well.Table 1
Comparison of the SIFs (horizontal crack).
G2/G1d/2aPresentIsida et al. [40]Diff. (%)Yuuki et al. [39]Diff. (%)KI/(pπa) of Tip (A) or (B)0.250.051.4761.468-0.571.468-0.540.250.51.1981.197-0.091.197-0.122.00.050.8710.8720.140.869-0.172.00.50.9360.935-0.060.934-0.16KII/(pπa) of Tip (A) or (B)0.250.050.2850.2860.350.2922.500.250.50.0710.0710.700.0721.672.00.05-0.088-0.087-1.38-0.085-4.012.00.5-0.023-0.0242.50-0.023-3.54Figure 6
A horizontal crack under uniform pressure within material (1) of an infinite Bimaterial.


### 3.1.2. Vertical Crack in Material (1)

Consider a vertical crack in material (1) subjected to far-field horizontal stresses as shown in Figure 7. The crack has a length 2a and is located at a distance d to the interface. The Poisson ratios for both materials (1) and (2) are the same, that is, ν1 = 0.35, ν2 = 0.3, and the ratio of the shear module is the same, that is, G1/G2 = 23.077, while the ratio of the crack length located at a distance (d/a) varies. A plane stress condition is assumed. In order to calculate the SIFs at crack tip (A) or (B), 20 quadratic elements were used to discretize the crack surface. The results are given in Table 2 for various values of the crack length located at a distance. They are compared to the results given by Isida and Noguchi [38] using a body force integral equation method and those by Cook and Erdogan [41] using a Wiener-Hopf technique and an asymptotic analysis. As can be observed from this table, the results compare quite well.Table 2
Comparison of the SIFs (vertical crack).
d/aPresentIsida et al. [40]Diff. (%)Cook et al. [41]Diff. (%)KI/(σπa) Of Tip (A)1.000.8820.8830.110.8830.142.000.9610.9620.100.9620.105.000.9950.993-0.170.993-0.1710.000.9960.9980.240.9980.24KI/(σπa) Of Tip (B)1.00—————2.000.9350.935-0.040.935-0.045.000.9960.991-0.430.991-0.4410.001.0030.998-0.510.998-0.51Figure 7
A vertical crack in material (1): an interface under far-field stresses.


### 3.1.3. Vertical Crack Intersecting an Interface

Consider a vertical crack intersecting an interface and subjected to far-field horizontal stresses as shown in Figure8. The horizontal far field stresses applied in materials (1) and (2) are, respectively, σ1 and σ2(=σ1G2/G1). The Poisson ratio ν1 and ν2 are assumed to be equal to 0.3 and the shear modulus ratio G2/G1 is assumed to vary. The distance of the crack tip (A) and (B) to the interface are the same, that is, d1=d2= a, the half-length of the crack. Again, a plane stress condition is assumed and 20 quadratic elements were used to discretize the crack surface. The SIFs at the crack tips (A) and (B) are listed in Table 3 for several values of the shear modulus ratio, and are compare to those given by Isida and Noguchi [38]. Again, the results between the two numerical analyses compare quite well.Table 3
Comparison of the SIFs (vertical crack intersecting an interface).
G2/G1KI/(σπa) of Tip (A)KI/(σπa) of Tip (B)Isida and Noguchi  [38]PresentDiff. (%)Isida and Noguchi  [38]PresentDiff. (%)0.11.0621.0629-0.081.1531.1539-0.080.31.0151.0157-0.071.0641.06390.010.51.0001.0007-0.071.0281.02730.070.80.9970.9975-0.051.0061.00470.13Figure 8
A vertical crack intersecting an interface under far-field stresses.


### 3.1.4. An Interface Kinked Crack in Infinite Bimaterials

Consider an interface kinked (45°) crack in infinite Bimaterials subjected to far-field tensile stresses as shown in Figure 9. The far-field tensile stresses applied in materials (1) and (2) are, respectively, σ1 and σ2(=σ1G2/G1). The kink crack length is 2a; the main crack length is 2c. Kink crack’s tip is (A). The Poisson ratios ν1 and ν2 are assumed to be equal to 0.3. The Young module in material (1) is 1 GPa, and according to the shear module ratio with material (1) and material (2) (G1/G2 = 0.25), we can get the Young module of material (2). Then according to equation (G=(E/2(1+ν))), we can get the relation between kink length and stress intensity factors. After numerical analysis, we compared with the results of Isida and Noguchi [38], using a body integral equation force method, as shown in Table 4, the results the two numerical analyses and found that they compare quite well.Table 4
Comparison of the SIFs (interfacial kinked crack).
a/cKI/(σπa) of Tip (A)KII/(σπa) of Tip (A)Isida and Noguchi  [38]PresentDiff. (%)Isida and Noguchi [38]PresentDiff. (%)0.20.7330.759-3.540.6310.6241.100.50.7080.730-3.100.6230.6121.761.00.6830.701-2.630.6100.6031.141.50.6690.685-2.390.6010.5941.16Figure 9
Interfacial kinked crack within infinite Bimaterials.


### 3.1.5. Interfacial Horizontal Crack in Infinite Bimaterials

An interfacial crack along the x-axis of length2a is shown in Figure 10. The crack surface is under a uniform pressure P and the materials can be either isotropic or anisotropic. Twenty quadratic elements were used to discretize the crack, and the characteristic length is assumed as 2a.Figure 10
Interfacial crack within infinite Bimaterials.The SIFs at the crack tip of an interfacial crack were also calculated for the anisotropic Bimaterial case. The anisotropic elastic properties in material (1) were assumed to be those of glass/epoxy with E1 = 48.26 GPa, E2 = 17.24 GPa, G12 = 6.89 GPa, and υ12 = 0.29. For material (2), a graphite/epoxy with E1 = 114.8 GPa, E2 = 11.7 GPa, G12= 9.66 GPa, and υ12 = 0.21 was selected [42]. The material axis E1 in material (1) and material (2) makes angles ψ1 and ψ2, respectively, with respect to the horizontal direction (Figure 10). While the material axis E1 in material (2) was assumed to be along the horizontal direction  (i.e., ψ(2) = 0), the E1-axis in material (1) makes an angle ψ(1) with respect to the horizontal direction. The interfacial SIFs at crack tip (B) obtained by the present method are listed in Table 5 and compared to the exact solutions proposed by Wu [12]. A very good agreement is found between the numerical analysis and the exact solution.Table 5
Comparison of the SIFs for infinite anisotropic problemψ2=0° (interfacial crack).
ψ1KI/(pπa) of  Tip (A) or (B)KII/(pπa) of  Tip (A) or (B)Wu [12]PresentDiff. (%)Wu [12]PresentDiff. (%)01.00001.0053-0.53-0.0382-0.03810.26300.99681.0006-0.38-0.0349-0.0350-0.29450.99651.0001-0.36-0.0318-0.0319-0.31600.99711.0010-0.39-0.0290-0.0292-0.69901.00001.0054-0.54-0.0264-0.0265-0.38


### 3.1.6. Interfacial Horizontal Crack in Finite Bimaterials

The example is included as a comparison with the literature in order to demonstrate the accuracy of SDBEM approach for an interfacial crack in anisotropic bimaterial plate. The geometry is that of rectangular plate and is shown in Figure11. For the comparison, crack length is taken as 2a = 2, h = 2w, and a/w = 0.4, and static tensile loading σ0 is applied on the upper and the lower boundary of the plate. Plane stress condition is assumed. The anisotropic elastic properties for materials (1) and (2) are given in Table 6. The normalized complex stress intensity factors at crack tip (A) or (B) are listed in Table 7 together with those from the work of Cho et al. [43], who used a multidomain BEM formulation and the results from Wünsche et al. [21] for a finite body. The outer boundary and interfacial crack surface were discretized with 40 continuous and 20 discontinuous quadratic elements, respectively. It is obvious from Table 7 that these are very close to those obtained by the other researchers.Table 6
Elastic properties for materials (1) and (2).
MaterialsE1 (MPa)E2 (MPa)ν12G12Material (1)100500.310.009Material (2)(i)100450.39.525(ii)100400.39.010(iii)100300.37.860(iv)100100.34.630Table 7
Comparison of the normalized complex SIFs for finite anisotropic problem (interfacial crack).
Material (2) E2/E1|K|/σ0πa of Tip (A) or (B)Cho et al. [43]Wünsche et al. [21]Diff. (%)PresentDiff. (%)(i) 0.451.3171.3120.381.31320.29(ii) 0.401.3371.3330.301.33510.14(iii) 0.301.3921.3860.431.3922-0.01(iv) 0.101.6971.6890.471.6968-0.06*|K|=KI2+KII2Figure 11
An interfacial crack within finite rectangular plate of Bimaterial.


## 3.1.1. Horizontal Crack in Material (1)

A horizontal crack under a uniform pressureP is shown in Figure 6. The crack has a length 2a, and is located at a distance d to the interface. The Poisson ratios for both materials (1) and (2) the same, that is, ν1 = ν2 = 0.3, while the ratio of the shear module ratio G2/G1 varies. A plane stress condition is assumed. In order to calculate the SIFs at crack tip (A) or (B), 20 quadratic elements were used to discretize the crack surface. The results are given in Table 1 for various values of the shear module ratio. They are compared to the results given by Isida and Noguchi [38], using a body force integral equation method and those by Yuuki and Cho [39], using a multidomain BEM formulation. As can be observed from this table, the results compare quite well.Table 1
Comparison of the SIFs (horizontal crack).
G2/G1d/2aPresentIsida et al. [40]Diff. (%)Yuuki et al. [39]Diff. (%)KI/(pπa) of Tip (A) or (B)0.250.051.4761.468-0.571.468-0.540.250.51.1981.197-0.091.197-0.122.00.050.8710.8720.140.869-0.172.00.50.9360.935-0.060.934-0.16KII/(pπa) of Tip (A) or (B)0.250.050.2850.2860.350.2922.500.250.50.0710.0710.700.0721.672.00.05-0.088-0.087-1.38-0.085-4.012.00.5-0.023-0.0242.50-0.023-3.54Figure 6
A horizontal crack under uniform pressure within material (1) of an infinite Bimaterial.


## 3.1.2. Vertical Crack in Material (1)

Consider a vertical crack in material (1) subjected to far-field horizontal stresses as shown in Figure 7. The crack has a length 2a and is located at a distance d to the interface. The Poisson ratios for both materials (1) and (2) are the same, that is, ν1 = 0.35, ν2 = 0.3, and the ratio of the shear module is the same, that is, G1/G2 = 23.077, while the ratio of the crack length located at a distance (d/a) varies. A plane stress condition is assumed. In order to calculate the SIFs at crack tip (A) or (B), 20 quadratic elements were used to discretize the crack surface. The results are given in Table 2 for various values of the crack length located at a distance. They are compared to the results given by Isida and Noguchi [38] using a body force integral equation method and those by Cook and Erdogan [41] using a Wiener-Hopf technique and an asymptotic analysis. As can be observed from this table, the results compare quite well.Table 2
Comparison of the SIFs (vertical crack).
d/aPresentIsida et al. [40]Diff. (%)Cook et al. [41]Diff. (%)KI/(σπa) Of Tip (A)1.000.8820.8830.110.8830.142.000.9610.9620.100.9620.105.000.9950.993-0.170.993-0.1710.000.9960.9980.240.9980.24KI/(σπa) Of Tip (B)1.00—————2.000.9350.935-0.040.935-0.045.000.9960.991-0.430.991-0.4410.001.0030.998-0.510.998-0.51Figure 7
A vertical crack in material (1): an interface under far-field stresses.


## 3.1.3. Vertical Crack Intersecting an Interface

Consider a vertical crack intersecting an interface and subjected to far-field horizontal stresses as shown in Figure8. The horizontal far field stresses applied in materials (1) and (2) are, respectively, σ1 and σ2(=σ1G2/G1). The Poisson ratio ν1 and ν2 are assumed to be equal to 0.3 and the shear modulus ratio G2/G1 is assumed to vary. The distance of the crack tip (A) and (B) to the interface are the same, that is, d1=d2= a, the half-length of the crack. Again, a plane stress condition is assumed and 20 quadratic elements were used to discretize the crack surface. The SIFs at the crack tips (A) and (B) are listed in Table 3 for several values of the shear modulus ratio, and are compare to those given by Isida and Noguchi [38]. Again, the results between the two numerical analyses compare quite well.Table 3
Comparison of the SIFs (vertical crack intersecting an interface).
G2/G1KI/(σπa) of Tip (A)KI/(σπa) of Tip (B)Isida and Noguchi  [38]PresentDiff. (%)Isida and Noguchi  [38]PresentDiff. (%)0.11.0621.0629-0.081.1531.1539-0.080.31.0151.0157-0.071.0641.06390.010.51.0001.0007-0.071.0281.02730.070.80.9970.9975-0.051.0061.00470.13Figure 8
A vertical crack intersecting an interface under far-field stresses.


## 3.1.4. An Interface Kinked Crack in Infinite Bimaterials

Consider an interface kinked (45°) crack in infinite Bimaterials subjected to far-field tensile stresses as shown in Figure 9. The far-field tensile stresses applied in materials (1) and (2) are, respectively, σ1 and σ2(=σ1G2/G1). The kink crack length is 2a; the main crack length is 2c. Kink crack’s tip is (A). The Poisson ratios ν1 and ν2 are assumed to be equal to 0.3. The Young module in material (1) is 1 GPa, and according to the shear module ratio with material (1) and material (2) (G1/G2 = 0.25), we can get the Young module of material (2). Then according to equation (G=(E/2(1+ν))), we can get the relation between kink length and stress intensity factors. After numerical analysis, we compared with the results of Isida and Noguchi [38], using a body integral equation force method, as shown in Table 4, the results the two numerical analyses and found that they compare quite well.Table 4
Comparison of the SIFs (interfacial kinked crack).
a/cKI/(σπa) of Tip (A)KII/(σπa) of Tip (A)Isida and Noguchi  [38]PresentDiff. (%)Isida and Noguchi [38]PresentDiff. (%)0.20.7330.759-3.540.6310.6241.100.50.7080.730-3.100.6230.6121.761.00.6830.701-2.630.6100.6031.141.50.6690.685-2.390.6010.5941.16Figure 9
Interfacial kinked crack within infinite Bimaterials.


## 3.1.5. Interfacial Horizontal Crack in Infinite Bimaterials

An interfacial crack along the x-axis of length2a is shown in Figure 10. The crack surface is under a uniform pressure P and the materials can be either isotropic or anisotropic. Twenty quadratic elements were used to discretize the crack, and the characteristic length is assumed as 2a.Figure 10
Interfacial crack within infinite Bimaterials.The SIFs at the crack tip of an interfacial crack were also calculated for the anisotropic Bimaterial case. The anisotropic elastic properties in material (1) were assumed to be those of glass/epoxy with E1 = 48.26 GPa, E2 = 17.24 GPa, G12 = 6.89 GPa, and υ12 = 0.29. For material (2), a graphite/epoxy with E1 = 114.8 GPa, E2 = 11.7 GPa, G12= 9.66 GPa, and υ12 = 0.21 was selected [42]. The material axis E1 in material (1) and material (2) makes angles ψ1 and ψ2, respectively, with respect to the horizontal direction (Figure 10). While the material axis E1 in material (2) was assumed to be along the horizontal direction  (i.e., ψ(2) = 0), the E1-axis in material (1) makes an angle ψ(1) with respect to the horizontal direction. The interfacial SIFs at crack tip (B) obtained by the present method are listed in Table 5 and compared to the exact solutions proposed by Wu [12]. A very good agreement is found between the numerical analysis and the exact solution.Table 5
Comparison of the SIFs for infinite anisotropic problemψ2=0° (interfacial crack).
ψ1KI/(pπa) of  Tip (A) or (B)KII/(pπa) of  Tip (A) or (B)Wu [12]PresentDiff. (%)Wu [12]PresentDiff. (%)01.00001.0053-0.53-0.0382-0.03810.26300.99681.0006-0.38-0.0349-0.0350-0.29450.99651.0001-0.36-0.0318-0.0319-0.31600.99711.0010-0.39-0.0290-0.0292-0.69901.00001.0054-0.54-0.0264-0.0265-0.38


## 3.1.6. Interfacial Horizontal Crack in Finite Bimaterials

The example is included as a comparison with the literature in order to demonstrate the accuracy of SDBEM approach for an interfacial crack in anisotropic bimaterial plate. The geometry is that of rectangular plate and is shown in Figure11. For the comparison, crack length is taken as 2a = 2, h = 2w, and a/w = 0.4, and static tensile loading σ0 is applied on the upper and the lower boundary of the plate. Plane stress condition is assumed. The anisotropic elastic properties for materials (1) and (2) are given in Table 6. The normalized complex stress intensity factors at crack tip (A) or (B) are listed in Table 7 together with those from the work of Cho et al. [43], who used a multidomain BEM formulation and the results from Wünsche et al. [21] for a finite body. The outer boundary and interfacial crack surface were discretized with 40 continuous and 20 discontinuous quadratic elements, respectively. It is obvious from Table 7 that these are very close to those obtained by the other researchers.Table 6
Elastic properties for materials (1) and (2).
MaterialsE1 (MPa)E2 (MPa)ν12G12Material (1)100500.310.009Material (2)(i)100450.39.525(ii)100400.39.010(iii)100300.37.860(iv)100100.34.630Table 7
Comparison of the normalized complex SIFs for finite anisotropic problem (interfacial crack).
Material (2) E2/E1|K|/σ0πa of Tip (A) or (B)Cho et al. [43]Wünsche et al. [21]Diff. (%)PresentDiff. (%)(i) 0.451.3171.3120.381.31320.29(ii) 0.401.3371.3330.301.33510.14(iii) 0.301.3921.3860.431.3922-0.01(iv) 0.101.6971.6890.471.6968-0.06*|K|=KI2+KII2Figure 11
An interfacial crack within finite rectangular plate of Bimaterial.


## 3.2. Crack Initiation Angles

### 3.2.1. Comparison of Numerical Predictions of Crack Initiation Angles with Experimental Results

In this section, we compared the numerical predictions of crack initiation angles with experimental results to verify the formulation. Al-Shayea [44] conducted uniaxial pressure test on limestone rock discs 98 mm and 84 mm in diameter and 22 mm in thickness with 30 mm notch.  The crack orientation angle β between the crack plane and the tensile stress was varied. Figure 12 shows the variation of the crack initiation angle θ with the crack angle β determined numerically and experimentally. A good agreement is found between the experimental results of Al-Shayea [44] and our numerical predictions.Figure 12
A variation of crack initiation angleθ with the crack angle β.


## 3.2.1. Comparison of Numerical Predictions of Crack Initiation Angles with Experimental Results

In this section, we compared the numerical predictions of crack initiation angles with experimental results to verify the formulation. Al-Shayea [44] conducted uniaxial pressure test on limestone rock discs 98 mm and 84 mm in diameter and 22 mm in thickness with 30 mm notch.  The crack orientation angle β between the crack plane and the tensile stress was varied. Figure 12 shows the variation of the crack initiation angle θ with the crack angle β determined numerically and experimentally. A good agreement is found between the experimental results of Al-Shayea [44] and our numerical predictions.Figure 12
A variation of crack initiation angleθ with the crack angle β.


## 4. Comparison of Numerical Predictions of Crack Propagation Paths with Experimental Results

To demonstrate the proposed SDBEM procedure when predicting crack propagation in the anisotropic Bimaterials 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 (as shown in Figure13). 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. The anisotropic elastic properties for rocks of types (1) and (2) are given in Table 8. The ratios of E1/E2 and E1/G12 are 1.635 and 4.301, respectively. Since the value of E1/E2 = 1.635, this rock of type (1) can be classified as a slightly anisotropic rock. Two specimens with the material inclination angle ψ = 45° and 0°, defined as the SSD-1 and SSD-2, have crack angles β = 45° and 90°, respectively. After Brazilian tests with cracked bimaterial specimens, the paths of crack propagation for SSD-1 and SSD-2 are shown in Figures 14(a) and 15(a), 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 SDBEM 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 14(b) and 15(b) are the comparisons examples; it can be concluded that the proposed SDBEM is capable of predicting the crack propagation in anisotropic bimaterial rocks.Table 8
Elastic properties for the rocks of type (1) and (2).
Rock typesE1E2νν12GG12E1/E2E1/G12(1)28.04017.1500.1500.12012.1916.5201.6354.301(2)38.950—0.254—15.530———Geometry of a cracked straight through brazilian disc (CSTBD) specimen of anisotropic Bimaterial rock under diametrical loading.
(a)(b)Comparison between experimental observations and numerical predictions for specimen SSD-1: (a) photograph of specimen SSD-1 after failure (ψ = 45° and β = 45°) and (b) propagation of a crack at the center of a CSTBD specimen with ψ = 45° and β = 45°.
(a)(b)Comparison between experimental observations and numerical predictions for specimen SSD-2: (a) photograph of specimen SSD-2 after failure (ψ = 0° and β = 90°) and (b) propagation of an interfacial crack at the center of a CSTBD specimen with ψ = 0° and β = 90°.
(a)(b)


## 5. Conclusions

This paper shows that the mixed mode stress intensity factors of anisotropic Bimaterial rock under diametrical loading can be successfully determined by the SDBEM.A new SDBEM procedure based on the maximum tensile stress failure criterion was developed to predict the crack initiation direction and the crack propagation path in anisotropic Bimaterial rock discs under mixed mode loading. A good agreement was found between crack initiation angles and propagation paths predicted with the SDBEM and experimental observations reported by previous researchers on anisotropic materials. Numerical simulations of crack initiation and propagation in CSTBD specimens of  type (1)/(2) were also found to compare well with the experimental results. Since the present method is simple and can be used for curved cracks, it will be straightforward to extend the current SDBEM formulation to analyze fracture propagation in 2D anisotropic Bimaterials, which is currently under investigation by the authors.



---

*Source: 947818-2010-11-02.xml*