# Two-Dimensional Stress Intensity Factor Analysis of Cracks in Anisotropic Bimaterial

Authors: Chia-Huei Tu; Jia-Jyun Dong; Chao-Shi Chen; Chien-Chung Ke; Jyun-Yong Jhan; Hsien Jui Yu
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

This paper presents a 2D numerical technique based on the boundary element method (BEM) for the analysis of linear elastic fracture mechanics (LEFM) problems on stress intensity factors (SIFs) involving anisotropic bimaterials. The most outstanding feature of this analysis is that it is a singledomain method, yet it is very accurate, efficient, and versatile (i.e., the material properties of the medium can be anisotropic as well as isotropic). A computer program using the BEM formula translation (FORTRAN 90) code was developed to effectively calculate the stress intensity factors (SIFs) in an anisotropic bi-material. This BEM program has been verified and showed good accuracy compared with the previous studies. Numerical examples of stress intensity factor calculation for a straight crack with various locations in both finite and infinite bimaterials are presented. It was found that very accurate results can be obtained using the proposed method, even with relatively simple discretization. The results of the numerical analysis also show that material anisotropy can greatly affect the stress intensity factor.


---

## Body

## 1. Introduction

The problem of cracks between two dissimilar materials has been widely studied over the past several decades, stemming mainly from the desire to understand the failure modes of composites, including structures, rocks, and concrete. Williams [1] presented the first study of the plane problem of cracks between dissimilar isotropic materials. Williams showed that stresses possess the singularity of r-1/2±iε, where r is the radius distance from the crack tip and ε is a bi-material constant. England [2] investigated the problem of finite cracks between dissimilar isotropic materials. Rice and Sih [3] studied similar problems and derived the expressions of the stress fields near crack tips. Rice [4] reexamined the elastic fracture mechanics concepts of the isotropic interfacial crack and introduced an intrinsic material length scale so that the definition of the stress intensity factors (SIFs) possessed the same physical significance as those for homogeneous cracks. Clements [5] and Willis [6] extended the problem studied by England [2] for dissimilar anisotropic materials. They also showed the oscillatory behavior of the stresses and the phenomenon of interpenetration of the crack face near the crack tip in anisotropic interface cracks. Wu [7] extended the problem studied by Rice [4] for anisotropic bi-material cracks. Recently, many authors [7–10] have conducted studies on interfacial cracks in anisotropic materials and different definitions for the stress intensity factor existence.  However, the data on different crack locations, crack lengths, and degree of material anisotropy are scarce in the literature.Recently, several BEMs have been proposed for the study of cracked media [11–13]. It involves two sets of boundary integral equations and is, in general, superior to the aforementioned BEMs. Consequently, general mixed mode crack problems can be solved in a BEM formulation. The single-domain analysis can eliminate remeshing problems, which are typical of the FEM and the subregional BEM.In this paper, we develop a technique for single-domain BEM formulation in which neither the artificial boundary nor the discretization along the uncracked interface is necessary. Chen et al. [14] presented this single-domain BEM formulation for homogeneous materials. We combined it with Green’s functions of bi-materials [15] to extend it to anisotropic bi-materials. The BEM formulation is such that the displacement integral equation is applied only to the outer boundary, and the traction integral equation is applied only to one side of the cracked surface. A decoupling technique can be used to determine mixed mode SIFs based on the relative displacement at the crack tip. Numerical cases for mixed mode SIFs for anisotropic bi-materials with different crack locations and degrees of material anisotropy are presented. The numerical results obtained using the BEM formulation are verified by several previous studies [8, 16–18].


## 2. Methodology

### 2.1. Basic Equations for Anisotropic Elasticity

For 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 is the root of the following characteristic equation [19]:
(1)a11μ4-2a16μ3+(2a12+a66)μ2-2a26μ+a22=0,
where the coefficient aij is the compliance component calculated using the x-y coordinate system. The detailed relationships of these components with material elasticity can be found in Chen et al. [14]. If the root μj of (1) is assumed to be distinct, the general expressions for the stress and displacements are
(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 the elements of the complex matrices A are
(3)A1j=a11μj2+a22-a16μj,A2j=a12μj+a22μj-a26(j=1,2).
The traction components in the x and y directions are
(4)Tx=2Re[μ1ϕ1(z1)+μ2ϕ2(z2)],Ty=-2Re[ϕ1(z1)+ϕ2(z2)].
Here, the complex analytical functionsϕi(zi)can in general express (2) and (4) as follows [19–21]:
(5)ui=2Re[A1iϕ1(z1)+A2iϕ2(z2)],Ti=-2Re[B1iϕ1(z1)+B2iϕ2(z2)],σ2i=2Re[B1iϕ1′(z1)+B2iϕ2′(z2)],σ1i=-2Re[B1iϕ1′(z1)+B2iϕ2′(z2)],(i=1,2),
where zj=x+μjy; Re denotes the real part of a complex variable or function; a prime denotes the derivative; the complex number μj(j=1,2); the elements of complex matrices A are defined in (3); and the elements of complex matrices B can be defined as
(6)Bij=[-μ1-μ211].Assume that the medium is composed of two joined dissimilar anisotropic and elastic half-planes. Let the interface be along thex-axis and let the upper (y>0) and lower (y<0) half-planes be occupied by materials #1 and #2, respectively (as shown in Figure 1).Figure 1
Definition of the coordinate systems in an anisotropic bi-material.Considering a concentrated force acting at the source point(x0,y0) in material #2(y0<0), we express the complex vector function as [20]
(7)ϕ(z)={ϕ#1(z),z∈(1),ϕ#2(z)+ϕ#20(z),z∈(2),
where the vector function is
(8)ϕ(z)=[ϕ1(z),ϕ2(z)]T,
where the argument has the generic form z=x+μy.In (7), ϕ#20 is a singular solution corresponding to a point force acting at the point (x0,y0) on an anisotropic infinite plane with the elastic properties of material #2. This singular solution can be expressed as [20, 21]
(9)ϕk0(zk)=-12π[Hk1P1ln(zk-sk)+Hk2P2ln(zk-sk)],
where sk=x0+μky0, Pk is the magnitude of the point force in the k-direction, and
(10)H=A-1(Y-1+Y¯-1)-1,Y=iAB-1,
where i=-1, the overbar indicates the complex conjugate, and superscript −1 indicates the matrix inverse. Two unknown vector functions, ϕ#1(z) and ϕ#2(z), are solved in (7). The former is analytic in material #1, and the latter is analytic in material #2. These expressions can be found by requiring the continuity of the resultant traction and displacement across the interface in addition to the standard analytic continuation arguments. Substituting (9) for (7), we obtain
(11)ϕ#1(Z)=B#1-1(Y#1+Y¯#2)-1(Y¯#2+Y#2)B#2ϕ#20(z),z∈#1,ϕ#2(Z)=B#2-1(Y¯#1+Y#2)-1(Y¯#2-Y¯#1)B¯#2ϕ#20(z),z∈#2.
Therefore, the complex vector functions can be expressed as
(12)ϕ(z)={B#1-1(Y#1+Y¯#2)-1×(Y¯#2+Y#2)B#2ϕ#20(z),z∈#1,B#2-1(Y¯#1+Y#2)-1(Y¯#2-Y¯#1)×B¯#2ϕ#20(z)+ϕ#20(z),z∈#2.In (12), subscripts #1 and #2 are used to denote that the corresponding matrix or vector is for material #1(y>0) and material #2(y<0), respectively.Similarly, for a point force in material #1(y0>0), the complex vector functions can be expressed as
(13)ϕ(z)={B#1-1(Y¯#2+Y#1)-1×(Y¯#1-Y¯#2)B¯#1ϕ¯#20(z)+ϕ#10(z),z∈#1,B#2-1(Y#2+Y¯#1)-1×(Y¯#1+Y#1)B#1ϕ#10(z),z∈#2,
where vector function ϕ#10 is the infinite-plane solution given in (9), but with the elastic properties of material #1.The Green’s functions of the displacement and traction can be obtained by substituting the complex vector functions in (12) and (13) for (5). Here, Ukl* is the Green’s function for displacement and Tkl* is the Green’s function for traction. Superscripts and subscripts #1 and #2 are used to denote the corresponding quantities for materials #1(y>0) and #2(y<0), respectively [15].(i)
For source point(s) and field point (z) on material #1(y>0),
(14)Ukl*=-1πRe{∑j=12Alj#1[∑j=12Alj#1ln(zj#1-sj#1)Hjk#1+∑i=12Wji11ln(zj#1-si#1)H¯ik#1∑]∑j=12Alj#1},Tkl*=1πRe{+∑i=12Wji11μj#1nx-nyzj#1-si#1H¯ik#1∑]∑j=12Blj#1[∑μj#1nx-nyzj#1-sj#1Hjk#1+∑i=12Wji11μj#1nx-nyzj#1-si#1H¯ik#1∑]},
where the matrix H is defined in (10) with the anisotropic elastic properties of material #1, and
(15)W11=B#1-1(Y#1+Y¯#2)-1(Y¯#1-Y¯#2)B¯#1.(ii)
For source point(s) and field point (z) on material #1(y>0) and field point (z) on material #2(y<0),
(16)Ukl*=-1πRe{∑j=12Alj#2[∑i=12Wji12ln(zj#2-si#1)Hik#1]},Tkl*=1πRe{∑j=12Blj#2[∑i=12Wji12μj#2nx-nyzj#2-si#1H¯ik#1]},
with
(17)W12=B#2-1(Y#2+Y¯#1)-1(Y¯#1-Y#1)B#1.(iii)
For source point(s) and field point (z) on material #2(y<0),
(18)Ukl*=-1πRe{∑j=12Alj#2[∑ln(zj#2-sj#2)Hjk#2+∑i=12Wji22ln(zj#2-si#2)H¯ik#2∑]},Tkl*=1πRe{∑j=12Blj#2[∑μj#2nx-nyzj#2-sj#2Hjk#2+∑i=12Wji22μj#2nx-nyzj#2-si#2H¯ik#2∑]},
where matrix H is defined in (10) with the anisotropic elastic properties of material #2, and
(19)W22=B#2-1(Y#2+Y¯#1)-1(Y¯#2-Y¯#1)B¯#2.(iv)
For source point(s) and field point (z) on materials #2(y<0) and field point (z) on material #1(y>0),
(20)Ukl*=-1πRe{∑j=12Alj#1[∑i=12Wji21ln(zj#1-si#2)Hik#2]},Tkl*=1πRe{∑j=12Blj#1[∑i=12Wji21μj#1nx-nyzj#1-si#2H¯ik#2]},
with
(21)W21=B#1-1(Y#1+Y¯#2)-1(Y¯#2-Y#2)B#2.It should be noted that these Green’s functions can be used to solve both plane stress and plane strain problems for anisotropic bi-materials. 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 [22, 23].


### 2.2. Single-Domain Boundary Integral Equations

In this section, we present single-domain boundary integral equations in which neither the artificial boundary nor the discretization along the uncracked interface is necessary. This single-domain boundary integral equation was used recently by Chen et al. [14] for homogeneous materials, and it is now extended to anisotropic bi-materials.For a pointsk,B0 on the uncracked boundary, the displacement integral equation applied to the outer boundary results in the following form (sk,B0∈ΓB only, as shown in Figure 1):
(22)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 the Green’s tractions and displacements given in (14), (16), (18), and (20); uj and tj are the boundary displacements and tractions, respectively; and zk and sk,B0 are the field points and the source points on boundary Γ of the domain, respectively. ΓC has the same outward normal as ΓC+. bij are coefficients that depend only on the local geometry of the uncracked boundary at point sk,B0 and are equal to δij/2 for a smooth boundary. Here, subscripts B and C denote the outer boundary and the cracked surface, respectively. In deriving (22), we have assumed that the tractions on the two faces of a crack are equal and opposite. We emphasize here that since the bi-material Green’s functions are included in (22), discretization along the interface can be avoided, with the exception of the interfacial crack part, which is treated by the traction integral equation presented in the following.It should be noted that all the terms on the right-hand side of (22) have weak singularities and thus are integrable. Although the second term on the left-hand side of (22) has a strong singularity, it can be treated using the rigid-body motion method.The traction integral equation (forsk0 being a smooth point on the crack) applied to one side of the cracked surface is (sk,C0∈ΓC+ only)
(23)0.5tj(sk,C0)+nm(sk,C0)∫ΓΒClimkTij,k*(sk,C0,zk,B)×uj(zk,B)dΓ(zk,B)+nm(sk,C0)∫ΓCClimkTij,k*(sk,C0,zk,C)×[uj(zk,C+)-uj(zk,C-)]dΓ(zk,C)=nm(sk,C0)∫ΓBClimkUij,k*(sk,C0,zk,B)tj(zk,B)dΓ(zk,B),
where nm is the unit outward normal at cracked surface sk,C0 and Climk is the fourth-order stiffness tensor.Equations (22) and (23) form a pair of boundary integral equations [14, 21, 24] and can be used to calculate SIFs in anisotropic bi-materials. The main feature of the BEM formulation is that it is a single-domain formulation, with the displacement integral equation (22) being collocated only on the uncracked boundary and the traction integral equation (23) only on one side of the crack surface. For problems without cracks, one only needs (22), with the integral on the cracked surface being discarded. Equation (22) then reduces to the well-known displacement integral on the uncracked boundary.The internal stressesσ(zk) are determined using the following expression:
(24)σlm(zk)+∫ΓΒClimkTij,k*(sk,C0,zk,B)uj(zk,B)dΓ(zk,B)+∫ΓCClimkTij,k*(sk,C0,zk,C)×[uj(zk,C+)-uj(zk,C-)]dΓ(zk,C)=∫ΓBClimkUij,k*(sk,C0,zk,B)tj(zk,B)dΓ(zk,B).For given solutions (related to the body force of gravity, rotational forces, and the far-field stresses), the boundary integral equations (22) and (23) can be discretized and solved numerically for the unknown boundary displacements (or displacement discontinuities on the cracked surface) and tractions. In solving these equations, the hypersingular integral term in (23) can be handled using an accurate and efficient Gauss quadrature formula, which is similar to the traditional weighted Gauss quadrature but has a different weight [25].


### 2.3. Evaluation of Stress Intensity Factor

In fracture mechanics analysis, especially in the calculation of SIFs, one needs to know the asymptotic behavior of the displacements and stresses near the crack tip. In our BEM analysis of SIFs, we propose to use the extrapolation method of crack-tip displacement. We, therefore, need to know the exact asymptotic behavior of the relative crack displacement (RCD) behind the crack tip. The form of the asymptotic expression depends on the location of the crack tip. In this study, two cases are discussed: a crack tip in homogeneous material and an interfacial crack tip.(i)Crack Tip in Homogeneous Material. Assume that the crack tip is in material #1. The asymptotic behavior of the relative displacement at a distance r behind the crack tip can be expressed in terms of the three SIFs as [17]
(25)Δu(r)=22rπRe(Y#1)K,
where K=[KII,KI,KIII]T is the SIF vector and Y is a matrix with elements related to the anisotropic properties in material #1, as defined in (10).In order to calculate the square-root characteristic of the RCD near the crack tip, we constructed the following new crack tip element with the tip atξ=-1 (as shown in Figure 2(e)):
(26)Δui=∑k=13fkΔuik,
where subscript i denotes the RCD component and superscript k(1,2,3) denotes the RCDs at nodes ξ=-2/3,0,2/3 (as shown in Figure 2(b)), respectively. The shape functions fk are those introduced by [25]
(27)f1=338ξ+1[5-8(ξ+1)+3(ξ+1)2],f2=14ξ+1[-5+18(ξ+1)+9(ξ+1)2],f3=3385ξ+1[1-4(ξ+1)+3(ξ+1)2].Three-geometric-node quadratic elements used to approximate the uncracked boundary and the crack for two-dimensional problems.
(a)
Discontinuous quadratic element ofType I
(b)
Crack surface quadratic element ofType II
(c)
Discontinuous quadratic element ofType III
(d)
Discontinuous quadratic element ofType IV
(e)
Crack tip quadratic element ofType VIn this case, the relation of the RCDs at a distancer behind the crack tip and the SIFs can be found as [26–28]
(28)Δu1=22rπ(H11KI+H12KII),Δu2=22rπ(H21KI+H22KII),
where
(29)H11=Im(μ2P11-μ2P12μ1-μ2),H12=Im(P11-P12μ1-μ2),H21=Im(μ2P21-μ2P22μ1-μ2),H22=Im(P21-P22μ1-μ2).Substituting the RCDs into (26) and (28), we may obtain a set of algebraic equations with which the SIFs KI and KII can be solved.(ii)Interfacial Crack Tip. In this case, the relative crack displacements at a distance r behind the interfacial crack tip can be expressed in terms of the three SIFs as [29]
(30)Δu(r)=(∑j=13cjDQje-πδjr(1/2)+δj)K,
where cj, δj, Qj, and D are the relative parameters in materials #1 and #2. Utilizing (10), we defined the matrix of a material as
(31)Y#1+Y¯#2=D-iV,
where D and V are two real matrices. These two matrices are used to define matrix P as
(32)P=-D-1V.
The characteristic β relative to material is as follows:
(33)β=-12tr(P2).We use the characteristicβ to define oscillation index ε as
(34)ε=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.The relationship between characteristicβ and oscillation index ε is used to define constant cj as
(35)c1=22πβ2,c2=-e-πεdiε2π(1+2iε)β2cosh(πε),c3=-e-πεdiε2π(1-2iε)β2cosh(πε),
where d is the characteristic distance along the material interface to the crack tip.Comparing (25) with (30), we noticed that while the relative crack displacement behaves as a square root for a crack tip in a homogeneous medium, an interfacial crack-tip’s behavior is r(1/2)+iδ, a square root feature multiplied by a weak oscillation.Equation (30) can be written in the following form, which is more convenient for the current numerical applications:
(36)Δu(r)=2rπM(rd)K,
where d is the characteristic length and M is a matrix function expressed as
(37)M(x)=Dβ2{∑(P2+β2I)-(∑[cos(εlnx)+2εsin(εlnx)]P2+β[sin(εlnx)-2εcos(εlnx)]P∑)×((1+4ε2)cosh(πε))-1∑}.Again, in order to capture the square root and weak oscillatory behavior, we constructed a crack-tip element with the tip atξ = −1 (as shown in Figure 2(e)) in terms of which the relative crack displacement is expressed as
(38)Δu(r)=M(rd)[f1Δu11f2Δu12f3Δu13f1Δu21f2Δu22f3Δu23f1Δu31f2Δu31f3Δu33].


### 2.4. Construction of the Numerical Model

Quadratic elements are often used to model problems with curvilinear geometry. Three-geometric-node quadratic elements can be subdivided into five types, as shown in Figure2, during analysis. For these elements, both the geometry and the boundary quantities are approximated by intrinsic coordinate ξ, and shape function fk is obtained from (27). Generally, the three geometric nodes of an element are used as collocation points. Such an element is referred to as a continuous quadratic element of Type IV (as shown in Figure 2(d)). To model the corner points or the points where there is a change in the boundary conditions, one-sided discontinuous quadratic elements (commonly referred to as partially discontinuous elements) are used. These elements can be one of four types depending on which extreme collocation point has been shifted inside the element to model the geometric or physical singularity. If the third collocation node is shifted inside, then the resulting element is referred to as a partially discontinuous quadratic element of Type I (as shown in Figure 2(a)). If the first node is shifted inside, then the element is called a partially discontinuous quadratic element of Type III (as shown in Figure 2(c)). One may opt for moving both extreme nodes inside, which results in discontinuous elements of Type II and Type V (to model the cracked surface or the points where there is in the crack tip positions) as shown in Figures 2(b) and 2(e), respectively. In the graphical representation of these elements, we use an “×” for a geometric node, a “○” for a collocation node, and a “◊” for the crack tip position.


## 2.1. Basic Equations for Anisotropic Elasticity

For 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 is the root of the following characteristic equation [19]:
(1)a11μ4-2a16μ3+(2a12+a66)μ2-2a26μ+a22=0,
where the coefficient aij is the compliance component calculated using the x-y coordinate system. The detailed relationships of these components with material elasticity can be found in Chen et al. [14]. If the root μj of (1) is assumed to be distinct, the general expressions for the stress and displacements are
(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 the elements of the complex matrices A are
(3)A1j=a11μj2+a22-a16μj,A2j=a12μj+a22μj-a26(j=1,2).
The traction components in the x and y directions are
(4)Tx=2Re[μ1ϕ1(z1)+μ2ϕ2(z2)],Ty=-2Re[ϕ1(z1)+ϕ2(z2)].
Here, the complex analytical functionsϕi(zi)can in general express (2) and (4) as follows [19–21]:
(5)ui=2Re[A1iϕ1(z1)+A2iϕ2(z2)],Ti=-2Re[B1iϕ1(z1)+B2iϕ2(z2)],σ2i=2Re[B1iϕ1′(z1)+B2iϕ2′(z2)],σ1i=-2Re[B1iϕ1′(z1)+B2iϕ2′(z2)],(i=1,2),
where zj=x+μjy; Re denotes the real part of a complex variable or function; a prime denotes the derivative; the complex number μj(j=1,2); the elements of complex matrices A are defined in (3); and the elements of complex matrices B can be defined as
(6)Bij=[-μ1-μ211].Assume that the medium is composed of two joined dissimilar anisotropic and elastic half-planes. Let the interface be along thex-axis and let the upper (y>0) and lower (y<0) half-planes be occupied by materials #1 and #2, respectively (as shown in Figure 1).Figure 1
Definition of the coordinate systems in an anisotropic bi-material.Considering a concentrated force acting at the source point(x0,y0) in material #2(y0<0), we express the complex vector function as [20]
(7)ϕ(z)={ϕ#1(z),z∈(1),ϕ#2(z)+ϕ#20(z),z∈(2),
where the vector function is
(8)ϕ(z)=[ϕ1(z),ϕ2(z)]T,
where the argument has the generic form z=x+μy.In (7), ϕ#20 is a singular solution corresponding to a point force acting at the point (x0,y0) on an anisotropic infinite plane with the elastic properties of material #2. This singular solution can be expressed as [20, 21]
(9)ϕk0(zk)=-12π[Hk1P1ln(zk-sk)+Hk2P2ln(zk-sk)],
where sk=x0+μky0, Pk is the magnitude of the point force in the k-direction, and
(10)H=A-1(Y-1+Y¯-1)-1,Y=iAB-1,
where i=-1, the overbar indicates the complex conjugate, and superscript −1 indicates the matrix inverse. Two unknown vector functions, ϕ#1(z) and ϕ#2(z), are solved in (7). The former is analytic in material #1, and the latter is analytic in material #2. These expressions can be found by requiring the continuity of the resultant traction and displacement across the interface in addition to the standard analytic continuation arguments. Substituting (9) for (7), we obtain
(11)ϕ#1(Z)=B#1-1(Y#1+Y¯#2)-1(Y¯#2+Y#2)B#2ϕ#20(z),z∈#1,ϕ#2(Z)=B#2-1(Y¯#1+Y#2)-1(Y¯#2-Y¯#1)B¯#2ϕ#20(z),z∈#2.
Therefore, the complex vector functions can be expressed as
(12)ϕ(z)={B#1-1(Y#1+Y¯#2)-1×(Y¯#2+Y#2)B#2ϕ#20(z),z∈#1,B#2-1(Y¯#1+Y#2)-1(Y¯#2-Y¯#1)×B¯#2ϕ#20(z)+ϕ#20(z),z∈#2.In (12), subscripts #1 and #2 are used to denote that the corresponding matrix or vector is for material #1(y>0) and material #2(y<0), respectively.Similarly, for a point force in material #1(y0>0), the complex vector functions can be expressed as
(13)ϕ(z)={B#1-1(Y¯#2+Y#1)-1×(Y¯#1-Y¯#2)B¯#1ϕ¯#20(z)+ϕ#10(z),z∈#1,B#2-1(Y#2+Y¯#1)-1×(Y¯#1+Y#1)B#1ϕ#10(z),z∈#2,
where vector function ϕ#10 is the infinite-plane solution given in (9), but with the elastic properties of material #1.The Green’s functions of the displacement and traction can be obtained by substituting the complex vector functions in (12) and (13) for (5). Here, Ukl* is the Green’s function for displacement and Tkl* is the Green’s function for traction. Superscripts and subscripts #1 and #2 are used to denote the corresponding quantities for materials #1(y>0) and #2(y<0), respectively [15].(i)
For source point(s) and field point (z) on material #1(y>0),
(14)Ukl*=-1πRe{∑j=12Alj#1[∑j=12Alj#1ln(zj#1-sj#1)Hjk#1+∑i=12Wji11ln(zj#1-si#1)H¯ik#1∑]∑j=12Alj#1},Tkl*=1πRe{+∑i=12Wji11μj#1nx-nyzj#1-si#1H¯ik#1∑]∑j=12Blj#1[∑μj#1nx-nyzj#1-sj#1Hjk#1+∑i=12Wji11μj#1nx-nyzj#1-si#1H¯ik#1∑]},
where the matrix H is defined in (10) with the anisotropic elastic properties of material #1, and
(15)W11=B#1-1(Y#1+Y¯#2)-1(Y¯#1-Y¯#2)B¯#1.(ii)
For source point(s) and field point (z) on material #1(y>0) and field point (z) on material #2(y<0),
(16)Ukl*=-1πRe{∑j=12Alj#2[∑i=12Wji12ln(zj#2-si#1)Hik#1]},Tkl*=1πRe{∑j=12Blj#2[∑i=12Wji12μj#2nx-nyzj#2-si#1H¯ik#1]},
with
(17)W12=B#2-1(Y#2+Y¯#1)-1(Y¯#1-Y#1)B#1.(iii)
For source point(s) and field point (z) on material #2(y<0),
(18)Ukl*=-1πRe{∑j=12Alj#2[∑ln(zj#2-sj#2)Hjk#2+∑i=12Wji22ln(zj#2-si#2)H¯ik#2∑]},Tkl*=1πRe{∑j=12Blj#2[∑μj#2nx-nyzj#2-sj#2Hjk#2+∑i=12Wji22μj#2nx-nyzj#2-si#2H¯ik#2∑]},
where matrix H is defined in (10) with the anisotropic elastic properties of material #2, and
(19)W22=B#2-1(Y#2+Y¯#1)-1(Y¯#2-Y¯#1)B¯#2.(iv)
For source point(s) and field point (z) on materials #2(y<0) and field point (z) on material #1(y>0),
(20)Ukl*=-1πRe{∑j=12Alj#1[∑i=12Wji21ln(zj#1-si#2)Hik#2]},Tkl*=1πRe{∑j=12Blj#1[∑i=12Wji21μj#1nx-nyzj#1-si#2H¯ik#2]},
with
(21)W21=B#1-1(Y#1+Y¯#2)-1(Y¯#2-Y#2)B#2.It should be noted that these Green’s functions can be used to solve both plane stress and plane strain problems for anisotropic bi-materials. 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 [22, 23].


## 2.2. Single-Domain Boundary Integral Equations

In this section, we present single-domain boundary integral equations in which neither the artificial boundary nor the discretization along the uncracked interface is necessary. This single-domain boundary integral equation was used recently by Chen et al. [14] for homogeneous materials, and it is now extended to anisotropic bi-materials.For a pointsk,B0 on the uncracked boundary, the displacement integral equation applied to the outer boundary results in the following form (sk,B0∈ΓB only, as shown in Figure 1):
(22)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 the Green’s tractions and displacements given in (14), (16), (18), and (20); uj and tj are the boundary displacements and tractions, respectively; and zk and sk,B0 are the field points and the source points on boundary Γ of the domain, respectively. ΓC has the same outward normal as ΓC+. bij are coefficients that depend only on the local geometry of the uncracked boundary at point sk,B0 and are equal to δij/2 for a smooth boundary. Here, subscripts B and C denote the outer boundary and the cracked surface, respectively. In deriving (22), we have assumed that the tractions on the two faces of a crack are equal and opposite. We emphasize here that since the bi-material Green’s functions are included in (22), discretization along the interface can be avoided, with the exception of the interfacial crack part, which is treated by the traction integral equation presented in the following.It should be noted that all the terms on the right-hand side of (22) have weak singularities and thus are integrable. Although the second term on the left-hand side of (22) has a strong singularity, it can be treated using the rigid-body motion method.The traction integral equation (forsk0 being a smooth point on the crack) applied to one side of the cracked surface is (sk,C0∈ΓC+ only)
(23)0.5tj(sk,C0)+nm(sk,C0)∫ΓΒClimkTij,k*(sk,C0,zk,B)×uj(zk,B)dΓ(zk,B)+nm(sk,C0)∫ΓCClimkTij,k*(sk,C0,zk,C)×[uj(zk,C+)-uj(zk,C-)]dΓ(zk,C)=nm(sk,C0)∫ΓBClimkUij,k*(sk,C0,zk,B)tj(zk,B)dΓ(zk,B),
where nm is the unit outward normal at cracked surface sk,C0 and Climk is the fourth-order stiffness tensor.Equations (22) and (23) form a pair of boundary integral equations [14, 21, 24] and can be used to calculate SIFs in anisotropic bi-materials. The main feature of the BEM formulation is that it is a single-domain formulation, with the displacement integral equation (22) being collocated only on the uncracked boundary and the traction integral equation (23) only on one side of the crack surface. For problems without cracks, one only needs (22), with the integral on the cracked surface being discarded. Equation (22) then reduces to the well-known displacement integral on the uncracked boundary.The internal stressesσ(zk) are determined using the following expression:
(24)σlm(zk)+∫ΓΒClimkTij,k*(sk,C0,zk,B)uj(zk,B)dΓ(zk,B)+∫ΓCClimkTij,k*(sk,C0,zk,C)×[uj(zk,C+)-uj(zk,C-)]dΓ(zk,C)=∫ΓBClimkUij,k*(sk,C0,zk,B)tj(zk,B)dΓ(zk,B).For given solutions (related to the body force of gravity, rotational forces, and the far-field stresses), the boundary integral equations (22) and (23) can be discretized and solved numerically for the unknown boundary displacements (or displacement discontinuities on the cracked surface) and tractions. In solving these equations, the hypersingular integral term in (23) can be handled using an accurate and efficient Gauss quadrature formula, which is similar to the traditional weighted Gauss quadrature but has a different weight [25].


## 2.3. Evaluation of Stress Intensity Factor

In fracture mechanics analysis, especially in the calculation of SIFs, one needs to know the asymptotic behavior of the displacements and stresses near the crack tip. In our BEM analysis of SIFs, we propose to use the extrapolation method of crack-tip displacement. We, therefore, need to know the exact asymptotic behavior of the relative crack displacement (RCD) behind the crack tip. The form of the asymptotic expression depends on the location of the crack tip. In this study, two cases are discussed: a crack tip in homogeneous material and an interfacial crack tip.(i)Crack Tip in Homogeneous Material. Assume that the crack tip is in material #1. The asymptotic behavior of the relative displacement at a distance r behind the crack tip can be expressed in terms of the three SIFs as [17]
(25)Δu(r)=22rπRe(Y#1)K,
where K=[KII,KI,KIII]T is the SIF vector and Y is a matrix with elements related to the anisotropic properties in material #1, as defined in (10).In order to calculate the square-root characteristic of the RCD near the crack tip, we constructed the following new crack tip element with the tip atξ=-1 (as shown in Figure 2(e)):
(26)Δui=∑k=13fkΔuik,
where subscript i denotes the RCD component and superscript k(1,2,3) denotes the RCDs at nodes ξ=-2/3,0,2/3 (as shown in Figure 2(b)), respectively. The shape functions fk are those introduced by [25]
(27)f1=338ξ+1[5-8(ξ+1)+3(ξ+1)2],f2=14ξ+1[-5+18(ξ+1)+9(ξ+1)2],f3=3385ξ+1[1-4(ξ+1)+3(ξ+1)2].Three-geometric-node quadratic elements used to approximate the uncracked boundary and the crack for two-dimensional problems.
(a)
Discontinuous quadratic element ofType I
(b)
Crack surface quadratic element ofType II
(c)
Discontinuous quadratic element ofType III
(d)
Discontinuous quadratic element ofType IV
(e)
Crack tip quadratic element ofType VIn this case, the relation of the RCDs at a distancer behind the crack tip and the SIFs can be found as [26–28]
(28)Δu1=22rπ(H11KI+H12KII),Δu2=22rπ(H21KI+H22KII),
where
(29)H11=Im(μ2P11-μ2P12μ1-μ2),H12=Im(P11-P12μ1-μ2),H21=Im(μ2P21-μ2P22μ1-μ2),H22=Im(P21-P22μ1-μ2).Substituting the RCDs into (26) and (28), we may obtain a set of algebraic equations with which the SIFs KI and KII can be solved.(ii)Interfacial Crack Tip. In this case, the relative crack displacements at a distance r behind the interfacial crack tip can be expressed in terms of the three SIFs as [29]
(30)Δu(r)=(∑j=13cjDQje-πδjr(1/2)+δj)K,
where cj, δj, Qj, and D are the relative parameters in materials #1 and #2. Utilizing (10), we defined the matrix of a material as
(31)Y#1+Y¯#2=D-iV,
where D and V are two real matrices. These two matrices are used to define matrix P as
(32)P=-D-1V.
The characteristic β relative to material is as follows:
(33)β=-12tr(P2).We use the characteristicβ to define oscillation index ε as
(34)ε=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.The relationship between characteristicβ and oscillation index ε is used to define constant cj as
(35)c1=22πβ2,c2=-e-πεdiε2π(1+2iε)β2cosh(πε),c3=-e-πεdiε2π(1-2iε)β2cosh(πε),
where d is the characteristic distance along the material interface to the crack tip.Comparing (25) with (30), we noticed that while the relative crack displacement behaves as a square root for a crack tip in a homogeneous medium, an interfacial crack-tip’s behavior is r(1/2)+iδ, a square root feature multiplied by a weak oscillation.Equation (30) can be written in the following form, which is more convenient for the current numerical applications:
(36)Δu(r)=2rπM(rd)K,
where d is the characteristic length and M is a matrix function expressed as
(37)M(x)=Dβ2{∑(P2+β2I)-(∑[cos(εlnx)+2εsin(εlnx)]P2+β[sin(εlnx)-2εcos(εlnx)]P∑)×((1+4ε2)cosh(πε))-1∑}.Again, in order to capture the square root and weak oscillatory behavior, we constructed a crack-tip element with the tip atξ = −1 (as shown in Figure 2(e)) in terms of which the relative crack displacement is expressed as
(38)Δu(r)=M(rd)[f1Δu11f2Δu12f3Δu13f1Δu21f2Δu22f3Δu23f1Δu31f2Δu31f3Δu33].


## 2.4. Construction of the Numerical Model

Quadratic elements are often used to model problems with curvilinear geometry. Three-geometric-node quadratic elements can be subdivided into five types, as shown in Figure2, during analysis. For these elements, both the geometry and the boundary quantities are approximated by intrinsic coordinate ξ, and shape function fk is obtained from (27). Generally, the three geometric nodes of an element are used as collocation points. Such an element is referred to as a continuous quadratic element of Type IV (as shown in Figure 2(d)). To model the corner points or the points where there is a change in the boundary conditions, one-sided discontinuous quadratic elements (commonly referred to as partially discontinuous elements) are used. These elements can be one of four types depending on which extreme collocation point has been shifted inside the element to model the geometric or physical singularity. If the third collocation node is shifted inside, then the resulting element is referred to as a partially discontinuous quadratic element of Type I (as shown in Figure 2(a)). If the first node is shifted inside, then the element is called a partially discontinuous quadratic element of Type III (as shown in Figure 2(c)). One may opt for moving both extreme nodes inside, which results in discontinuous elements of Type II and Type V (to model the cracked surface or the points where there is in the crack tip positions) as shown in Figures 2(b) and 2(e), respectively. In the graphical representation of these elements, we use an “×” for a geometric node, a “○” for a collocation node, and a “◊” for the crack tip position.


## 3. Verification of the BEM Program

The Green’s functions and the particular crack-tip elements were incorporated into the boundary integral equations, and the results were programmed using FORTRAN code. In this section, the following numerical examples are presented to verify the formulation and to show the efficiency and versatility of the present BEM method for problems related to fractures in bi-material.


### 3.1. Horizontal Crack in Material#1

A horizontal crack under uniform pressureP is shown in Figure 3. The crack, whose length is 2a, is located distance d from the interface. The Poisson’s ratios for materials #1 and #2 are ν#1=ν#2=0.3, and the ratio of the shear modulus G#2/G#1 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 modulus ratio. They are compared to the results of Isida and Noguchi [16], who used a body force integral equation method, and those of Ryoji and Sang-Bong [17], who used a multidomain BEM formulation. It can be seen in this table that the results are in agreement.Table 1
Comparison of SIFs (horizontal crack).
G
(
2
)
/
G
(
1
)
d
/
2
a
Proposed approach
Isida and Noguchi [16]
Diff. (%)
Ryoji and Sang-Bong [17]
Diff. (%)
K
I
/
(
p
π
a
) of tip A (or B)
0.25
0.05
1.476
1.468
−0.57
1.468
−0.54
0.25
0.5
1.198
1.197
−0.09
1.197
−0.12
2.0
0.05
0.871
0.872
0.14
0.869
−0.17
2.0
0.5
0.936
0.935
−0.06
0.934
−0.16
K
II
/
(
p
π
a
) of tip A (or B)
0.25
0.05
0.285
0.286
0.35
0.292
2.50
0.25
0.5
0.071
0.071
0.70
0.072
1.67
2.0
0.05
−0.088
−0.087
−1.38
−0.085
−4.01
2.0
0.5
−0.023
−0.024
2.50
−0.023
−3.54Figure 3
Horizontal crack under uniform pressure in material#1 of an infinite bi-material.


### 3.2. Interfacial Crack in an Anisotropic Bi-Material Plate

The following example is used for comparison with results from the literature in order to demonstrate the accuracy of the BEM approach for an interfacial crack in an anisotropic bi-material plate. The geometry is that of a rectangular plate, as shown in Figure4. For the comparison, the crack length is 2a,h=2w,  a/w=0.4, and static tensile loading σ is applied to the upper and the lower boundaries of the plate.Figure 4
Interfacial crack within a finite rectangular plate of a bi-material.A plane stress condition is assumed. The anisotropic elastic properties for materials#1 and #2 are given in Table 2. The normalized complex stress intensity factors at crack tips A and B are listed in Table 3 together with those from Sang-Bong et al. [18], who used a multidomain BEM formulation, and the results from Wünsche et al. [8] for a finite body. The outer boundary and interfacial crack surface were discretized with 20 continuous and 20 discontinuous quadratic elements, respectively. Table 3 shows that the BEM approach produces results that agree well with those obtained by other researchers.Table 2
Elastic properties for materials #1 and #2.
Material
E (MPa)
E
′ (MPa)
ν
′
G
′
Material #1
100
50
0.3
10.009
Material #2
(i)
100
45
0.3
9.525
(ii)
100
40
0.3
9.010
(iii)
100
30
0.3
7.860
(iv)
100
10
0.3
4.630Table 3
Comparison of the normalized complex SIFs for finite anisotropic problem (interfacial crack).
Material #2E#2′/E#2
|
K
|
/
(
σ
π
a
) of tip A (or B)
Sang-Bong et al. [18]
Wünsche et al. [8]
Diff. (%)
Proposed approach
Diff. (%)
(i) 0.45
1.317
1.312
0.38
1.3132
0.29
(ii) 0.40
1.337
1.333
0.30
1.3351
0.14
(iii) 0.30
1.392
1.386
0.43
1.3922
−0.01
(iv) 0.10
1.697
1.689
0.47
1.6968
−0.06
*|K|=KI2+KII2.


## 3.1. Horizontal Crack in Material#1

A horizontal crack under uniform pressureP is shown in Figure 3. The crack, whose length is 2a, is located distance d from the interface. The Poisson’s ratios for materials #1 and #2 are ν#1=ν#2=0.3, and the ratio of the shear modulus G#2/G#1 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 modulus ratio. They are compared to the results of Isida and Noguchi [16], who used a body force integral equation method, and those of Ryoji and Sang-Bong [17], who used a multidomain BEM formulation. It can be seen in this table that the results are in agreement.Table 1
Comparison of SIFs (horizontal crack).
G
(
2
)
/
G
(
1
)
d
/
2
a
Proposed approach
Isida and Noguchi [16]
Diff. (%)
Ryoji and Sang-Bong [17]
Diff. (%)
K
I
/
(
p
π
a
) of tip A (or B)
0.25
0.05
1.476
1.468
−0.57
1.468
−0.54
0.25
0.5
1.198
1.197
−0.09
1.197
−0.12
2.0
0.05
0.871
0.872
0.14
0.869
−0.17
2.0
0.5
0.936
0.935
−0.06
0.934
−0.16
K
II
/
(
p
π
a
) of tip A (or B)
0.25
0.05
0.285
0.286
0.35
0.292
2.50
0.25
0.5
0.071
0.071
0.70
0.072
1.67
2.0
0.05
−0.088
−0.087
−1.38
−0.085
−4.01
2.0
0.5
−0.023
−0.024
2.50
−0.023
−3.54Figure 3
Horizontal crack under uniform pressure in material#1 of an infinite bi-material.


## 3.2. Interfacial Crack in an Anisotropic Bi-Material Plate

The following example is used for comparison with results from the literature in order to demonstrate the accuracy of the BEM approach for an interfacial crack in an anisotropic bi-material plate. The geometry is that of a rectangular plate, as shown in Figure4. For the comparison, the crack length is 2a,h=2w,  a/w=0.4, and static tensile loading σ is applied to the upper and the lower boundaries of the plate.Figure 4
Interfacial crack within a finite rectangular plate of a bi-material.A plane stress condition is assumed. The anisotropic elastic properties for materials#1 and #2 are given in Table 2. The normalized complex stress intensity factors at crack tips A and B are listed in Table 3 together with those from Sang-Bong et al. [18], who used a multidomain BEM formulation, and the results from Wünsche et al. [8] for a finite body. The outer boundary and interfacial crack surface were discretized with 20 continuous and 20 discontinuous quadratic elements, respectively. Table 3 shows that the BEM approach produces results that agree well with those obtained by other researchers.Table 2
Elastic properties for materials #1 and #2.
Material
E (MPa)
E
′ (MPa)
ν
′
G
′
Material #1
100
50
0.3
10.009
Material #2
(i)
100
45
0.3
9.525
(ii)
100
40
0.3
9.010
(iii)
100
30
0.3
7.860
(iv)
100
10
0.3
4.630Table 3
Comparison of the normalized complex SIFs for finite anisotropic problem (interfacial crack).
Material #2E#2′/E#2
|
K
|
/
(
σ
π
a
) of tip A (or B)
Sang-Bong et al. [18]
Wünsche et al. [8]
Diff. (%)
Proposed approach
Diff. (%)
(i) 0.45
1.317
1.312
0.38
1.3132
0.29
(ii) 0.40
1.337
1.333
0.30
1.3351
0.14
(iii) 0.30
1.392
1.386
0.43
1.3922
−0.01
(iv) 0.10
1.697
1.689
0.47
1.6968
−0.06
*|K|=KI2+KII2.


## 4. Numerical Analysis

After its accuracy was checked, the BEM program was applied to more challenging problems involving an anisotropic bi-material to determine the SIFs of crack tips. A generalized plane stress condition was assumed in all the problems. In the analysis, zero traction on the crack surface was assumed. Here,FI and FII are the normalized SIFs; for FI>0, the failure mechanism is tensile fracture, and for FII>0, the crack is subjected to an anticlockwise shear force. For a 2D problem under mixed mode loading, the failure mechanism may be defined by the SIFs of crack tips (tensile fracture, FI/FII>1, and shear fracture, FI/FII<1, resp.).


### 4.1. Effect of the Degree of Anisotropy on the Stress Intensity Factor

In this chapter, we analyze the SIFs for an anisotropic bi-material with various anisotropic directions,ψ, and degrees of material anisotropy using the BEM program. The normalized SIFs, defined as FI and FII, are equal to
(39)FI=KIK0,FII=KIIK0,
with
(40)K0=σπa.In order to evaluate the influence of bi-material anisotropy on the SIFs, we considered the following three cases: (i) a horizontal crack within material#1, (ii) a horizontal crack within material #2, and (iii) a crack at the interface of the materials. The infinite bi-material was subjected to far-field vertical tensile stresses, as shown in Figure 5. The crack, which had a length of 2a, was located at distances d/a=1 (material #1), −1 (material #2), and d/a=0 (interfacial crack) from the interface, respectively. Material #1 was transversely isotropic with the plane of transverse isotropy inclined at angle ψ#1 with respect to the x-axis. Material #1 had five independent elastic constants (E#1,E#1′,ν#1,ν#1′, andG#1) in the local coordinate system; E#i and E#i′ are the Young’s moduli in the plane of transverse isotropy and in a direction normal to it, respectively; ν#i and ν#i′ are the Poisson’s ratios that characterize the lateral strain response in the plane of transverse isotropy to a stress acting parallel and normal to it, respectively; and G#i′ is the shear modulus in the planes normal to the plane of transverse isotropy. For material #2, E#2/E#2′=1,E#2/G#2′=2.5, and ν#2=ν#2′=0.25. In this problem, the anisotropic direction (ψ) within material #1 varies from 0° to 180°.Horizontal cracks in an infinite bi-material under far-field stresses.
(a)
Horizontal crack within material#1
(b)
Horizontal crack within material#2
(c)
Interfacial crack within bi-materialFor all cases, three sets of dimensionless elastic constants are considered. They are defined in Table4. Here, E#1/E#2, which is the ratio of the Young’s modulus of material #1 to the Young’s modulus of material #2, equals one. In order to calculate the SIFs at the crack tip, 20 quadratic elements were used to discretize the crack surface. The numerical results are plotted in Figures 7 to 9 for cases I, II, and III, respectively. The results for the isotropic case (E#1/E#1′=1, E#1/G#1′=2.5, and ν#1′=0.25) are shown as solid lines in the figures for comparison.Table 4
Sets of dimensionless elastic constants(ν#1=0.25).
E
#
1
/
E
#
1
′
E
#
1
/
G
#
1
′
ν
#
1
′
1/3, 1/2, 1, 2, 3
2.5
0.25
*Subscript #1 is applied in material #1.


### 4.2. An Inclined Crack Situated Near the Interfacial Crack

This case treats interaction between the interfacial crack and an inclined crack with an angleβ. The ratio dt/a=0.25, respectively. The same consider an infinite plate of bi-materials, having the elastic constants ratio E#1/E#1′=1,E#1/G#1′=0.4, and ν#1′=0.2499 for material #1 and E#2/E#2′=1/3, E#2/G#2′=0.4, and ν#2′=0.2499 for material #2. The Poisson ratios are ν#1=ν#2=0.25, material orientation ψ#1=ψ#2=0°, and the Young modulus ratio E#1/E#2=1/5. One considers an interfacial crack of length 2a and a horizontal crack situated near the interfacial crack with ah. The ratio a/ah is equal to 1. The longitudinal distance dl/a=1, and the plate is subjected to a far-field vertical tensile stresses as shown in Figure 6. For each crack, surface is meshed by 20 quadratic elements. The plane stress conditions were supposed. Calculations were carried out by our BEM code. Figures 10 and 11 represent the variations of normalized SIFs (FI and FII) of interfacial crack and inclined crack according to the inclined angle β, respectively.Figure 6
A crack situated near the interfacial crack within infinite bi-materials.Variation of normalized SIFs of the crack tip in material#1 with the anisotropic directions ψ#1 for case I (E#1/G#1′=2.5,v#1′=0.25,E#1/E#2=1).
(a)
(b)


### 4.3. Discussion

An analysis of Figures7 to 9 reveals the following major trends.(i)
For the SIFs in mode I, the orientation of the planes of material anisotropy with respect to the horizontal plane has a strong influence on the value of the SIFs for case I and case II. The influence is small for case III. This means that the effects of mode I on the crack in the homogeneous material are greater than those of the interface.(ii)
For all cases, the variation of the SIFs with the anisotropic directionψ#1 is symmetric with respect to ψ#1=90°. However, this type of symmetry was not observed for the SIFs (mode II) in case 1 and case 2 (Figures 7(b) and 8(b)).(iii)
For the interfacial crack (case III), the anisotropic directionψ#1 has a strong influence on the value of the SIFs for mode II (Figure 9(b)). The influence is small for mode I. It should be noted that the effects of E#1/E#1′ on the SIFs of mode II are greater than those of mode I.(iv)
There is a greater variation in the SIFs with the anisotropic directionψ#1 for bi-materials with a high degree of anisotropy (E#1/E#1′=3or1/3).(v)
Figure7(a)  (E#1/E#1′=1/3) indicates that the maximum values of the SIFs (mode I) occur when the far-field stresses are perpendicular to the plane of transverse isotropy (i.e., ψ#1=0° or 180°) for case 1 (crack within material #1). Figure 9(b)  (E#1/E#1′=1/3) shows that the minimum values of the SIFs (mode II) occur when the anisotropic direction angle is about 90° for case 3 (interfacial crack).Analysis of Figures10 and 11 reveals that several major trends can be underlined as follows.(vi)
Figure10(a) presents the SIFs (mode I) variation as an inclined angle β of the inclined crack for different values of the interfacial crack. It is shown that the increase in the inclined angle β causes a slight increase in SIFs (mode I) until reaching a value at β=90°. From these angles, the tip D meets at the interface of bi-materials and the normalized SIFs start an immediate sharp increasing. The rise continued until it reached β=105°. The maximum of SIFs (mode I) is reached for this last angle. In this position, the distance between tip D and the interfacial crack is minimal which increases the stress interaction between the two crack tips. For β>120°, the SIFs start decreasing and stabilizing when approaching 180°.(vii)
The same phenomenon was observed in Figure10(b). It is shown that the increase in the inclined angle β causes the reduction in SIFs (mode II) until reaching a value at β=100°, and starts an immediate sharp increasing. The maximum of SIFs (mode II) is reached for β=110°. In this position, the distance between tip D and the interfacial crack is minimal which increases the stress interaction between the two crack tips. For β>120°, the SIFs start a slight increasing and stabilizing when approaching 180°.(viii)
Figure11 illustrates, respectively, the variations of the normalized SIFs of tip C and D according to inclined angle β. For Figure 11(a), between 0° and 90°, one observes that mode I of the normalized SIFs decreases with β. The minimal value is reached at β=90°. From 90°, the two normalized SIFs increase until 138°, and the normalized SIFs at the tip D reaches a maximum value. In the interval of 90° at 138°, the normalized SIFs of tip D are higher than those of tip C, because the tip D approaches the tip of the interfacial crack. Between 138° and 180°, tip D moves away from the tip A of interfacial crack.(ix)
For Figure11(b), the normalized SIFs (mode II) of the two tips grow in interval of the angles ranging between 0° and 35° and decrease between 35° and 80°. In the interval of 0° at 32° the normalized SIFs of tip C are slightly higher than those of tip D, and the inverse phenomenon is noted for angles ranging between 32°, and 100°. These results permit to note on one hand that the orientation of the inclined crack has influence on the variations of FI and a weak influence on the variations of FII, and in another hand, when the angle of this orientation increases, the normalized SIFs of mode I increase and automatically generate the reduction in mode II SIFs. For homogeneous and isotropic materials, the two SIFs for inclined crack are nulls at β=90°. The presence of the bi-material interface gives FI different of zero whatever the value of β and FII=0 at β=80°.(x)
The ratio of the normalized SIFs(FI/FII) is concerned; Figure 12 shows that the fracture of shear mode (FI/FII<1) occurs to the inclined crack angle is between 30° and 150° (tip C). For tip D, the shear mode (FI/FII<1) occurs to the β is between 50° and 95°, and between 105° and 120°, respectively. And the opening mode (FI/FII>1) occurs to the β is between 95° and 105°. In this position, the distance between tip D and tip A of the interfacial crack is minimal which increases the stress interaction between the two crack tips.Variation of Normalized SIFs of crack tip in material#2 with the anisotropic directions ψ#1 for case II (E#1/G#1′=2.5,ν#1′=0.25,E#1/E#2=1).
(a)
(b)Variation of normalized SIFs of interfacial crack tip with the anisotropic directionsψ#1 for case III (E#1/G#1′=2.5,ν#1′=0.25,E#1/E#2=1).
(a)
(b)Variation of normalized SIFs of an interfacial crack with an inclined crack angleβ.
(a)
(b)Variation of normalized SIFs of an inclined crack with the inclined angleβ.
(a)
(b)Figure 12
Variation of the ratio ofFI/FII of an inclined crack with the inclined angle β.


## 4.1. Effect of the Degree of Anisotropy on the Stress Intensity Factor

In this chapter, we analyze the SIFs for an anisotropic bi-material with various anisotropic directions,ψ, and degrees of material anisotropy using the BEM program. The normalized SIFs, defined as FI and FII, are equal to
(39)FI=KIK0,FII=KIIK0,
with
(40)K0=σπa.In order to evaluate the influence of bi-material anisotropy on the SIFs, we considered the following three cases: (i) a horizontal crack within material#1, (ii) a horizontal crack within material #2, and (iii) a crack at the interface of the materials. The infinite bi-material was subjected to far-field vertical tensile stresses, as shown in Figure 5. The crack, which had a length of 2a, was located at distances d/a=1 (material #1), −1 (material #2), and d/a=0 (interfacial crack) from the interface, respectively. Material #1 was transversely isotropic with the plane of transverse isotropy inclined at angle ψ#1 with respect to the x-axis. Material #1 had five independent elastic constants (E#1,E#1′,ν#1,ν#1′, andG#1) in the local coordinate system; E#i and E#i′ are the Young’s moduli in the plane of transverse isotropy and in a direction normal to it, respectively; ν#i and ν#i′ are the Poisson’s ratios that characterize the lateral strain response in the plane of transverse isotropy to a stress acting parallel and normal to it, respectively; and G#i′ is the shear modulus in the planes normal to the plane of transverse isotropy. For material #2, E#2/E#2′=1,E#2/G#2′=2.5, and ν#2=ν#2′=0.25. In this problem, the anisotropic direction (ψ) within material #1 varies from 0° to 180°.Horizontal cracks in an infinite bi-material under far-field stresses.
(a)
Horizontal crack within material#1
(b)
Horizontal crack within material#2
(c)
Interfacial crack within bi-materialFor all cases, three sets of dimensionless elastic constants are considered. They are defined in Table4. Here, E#1/E#2, which is the ratio of the Young’s modulus of material #1 to the Young’s modulus of material #2, equals one. In order to calculate the SIFs at the crack tip, 20 quadratic elements were used to discretize the crack surface. The numerical results are plotted in Figures 7 to 9 for cases I, II, and III, respectively. The results for the isotropic case (E#1/E#1′=1, E#1/G#1′=2.5, and ν#1′=0.25) are shown as solid lines in the figures for comparison.Table 4
Sets of dimensionless elastic constants(ν#1=0.25).
E
#
1
/
E
#
1
′
E
#
1
/
G
#
1
′
ν
#
1
′
1/3, 1/2, 1, 2, 3
2.5
0.25
*Subscript #1 is applied in material #1.


## 4.2. An Inclined Crack Situated Near the Interfacial Crack

This case treats interaction between the interfacial crack and an inclined crack with an angleβ. The ratio dt/a=0.25, respectively. The same consider an infinite plate of bi-materials, having the elastic constants ratio E#1/E#1′=1,E#1/G#1′=0.4, and ν#1′=0.2499 for material #1 and E#2/E#2′=1/3, E#2/G#2′=0.4, and ν#2′=0.2499 for material #2. The Poisson ratios are ν#1=ν#2=0.25, material orientation ψ#1=ψ#2=0°, and the Young modulus ratio E#1/E#2=1/5. One considers an interfacial crack of length 2a and a horizontal crack situated near the interfacial crack with ah. The ratio a/ah is equal to 1. The longitudinal distance dl/a=1, and the plate is subjected to a far-field vertical tensile stresses as shown in Figure 6. For each crack, surface is meshed by 20 quadratic elements. The plane stress conditions were supposed. Calculations were carried out by our BEM code. Figures 10 and 11 represent the variations of normalized SIFs (FI and FII) of interfacial crack and inclined crack according to the inclined angle β, respectively.Figure 6
A crack situated near the interfacial crack within infinite bi-materials.Variation of normalized SIFs of the crack tip in material#1 with the anisotropic directions ψ#1 for case I (E#1/G#1′=2.5,v#1′=0.25,E#1/E#2=1).
(a)
(b)


## 4.3. Discussion

An analysis of Figures7 to 9 reveals the following major trends.(i)
For the SIFs in mode I, the orientation of the planes of material anisotropy with respect to the horizontal plane has a strong influence on the value of the SIFs for case I and case II. The influence is small for case III. This means that the effects of mode I on the crack in the homogeneous material are greater than those of the interface.(ii)
For all cases, the variation of the SIFs with the anisotropic directionψ#1 is symmetric with respect to ψ#1=90°. However, this type of symmetry was not observed for the SIFs (mode II) in case 1 and case 2 (Figures 7(b) and 8(b)).(iii)
For the interfacial crack (case III), the anisotropic directionψ#1 has a strong influence on the value of the SIFs for mode II (Figure 9(b)). The influence is small for mode I. It should be noted that the effects of E#1/E#1′ on the SIFs of mode II are greater than those of mode I.(iv)
There is a greater variation in the SIFs with the anisotropic directionψ#1 for bi-materials with a high degree of anisotropy (E#1/E#1′=3or1/3).(v)
Figure7(a)  (E#1/E#1′=1/3) indicates that the maximum values of the SIFs (mode I) occur when the far-field stresses are perpendicular to the plane of transverse isotropy (i.e., ψ#1=0° or 180°) for case 1 (crack within material #1). Figure 9(b)  (E#1/E#1′=1/3) shows that the minimum values of the SIFs (mode II) occur when the anisotropic direction angle is about 90° for case 3 (interfacial crack).Analysis of Figures10 and 11 reveals that several major trends can be underlined as follows.(vi)
Figure10(a) presents the SIFs (mode I) variation as an inclined angle β of the inclined crack for different values of the interfacial crack. It is shown that the increase in the inclined angle β causes a slight increase in SIFs (mode I) until reaching a value at β=90°. From these angles, the tip D meets at the interface of bi-materials and the normalized SIFs start an immediate sharp increasing. The rise continued until it reached β=105°. The maximum of SIFs (mode I) is reached for this last angle. In this position, the distance between tip D and the interfacial crack is minimal which increases the stress interaction between the two crack tips. For β>120°, the SIFs start decreasing and stabilizing when approaching 180°.(vii)
The same phenomenon was observed in Figure10(b). It is shown that the increase in the inclined angle β causes the reduction in SIFs (mode II) until reaching a value at β=100°, and starts an immediate sharp increasing. The maximum of SIFs (mode II) is reached for β=110°. In this position, the distance between tip D and the interfacial crack is minimal which increases the stress interaction between the two crack tips. For β>120°, the SIFs start a slight increasing and stabilizing when approaching 180°.(viii)
Figure11 illustrates, respectively, the variations of the normalized SIFs of tip C and D according to inclined angle β. For Figure 11(a), between 0° and 90°, one observes that mode I of the normalized SIFs decreases with β. The minimal value is reached at β=90°. From 90°, the two normalized SIFs increase until 138°, and the normalized SIFs at the tip D reaches a maximum value. In the interval of 90° at 138°, the normalized SIFs of tip D are higher than those of tip C, because the tip D approaches the tip of the interfacial crack. Between 138° and 180°, tip D moves away from the tip A of interfacial crack.(ix)
For Figure11(b), the normalized SIFs (mode II) of the two tips grow in interval of the angles ranging between 0° and 35° and decrease between 35° and 80°. In the interval of 0° at 32° the normalized SIFs of tip C are slightly higher than those of tip D, and the inverse phenomenon is noted for angles ranging between 32°, and 100°. These results permit to note on one hand that the orientation of the inclined crack has influence on the variations of FI and a weak influence on the variations of FII, and in another hand, when the angle of this orientation increases, the normalized SIFs of mode I increase and automatically generate the reduction in mode II SIFs. For homogeneous and isotropic materials, the two SIFs for inclined crack are nulls at β=90°. The presence of the bi-material interface gives FI different of zero whatever the value of β and FII=0 at β=80°.(x)
The ratio of the normalized SIFs(FI/FII) is concerned; Figure 12 shows that the fracture of shear mode (FI/FII<1) occurs to the inclined crack angle is between 30° and 150° (tip C). For tip D, the shear mode (FI/FII<1) occurs to the β is between 50° and 95°, and between 105° and 120°, respectively. And the opening mode (FI/FII>1) occurs to the β is between 95° and 105°. In this position, the distance between tip D and tip A of the interfacial crack is minimal which increases the stress interaction between the two crack tips.Variation of Normalized SIFs of crack tip in material#2 with the anisotropic directions ψ#1 for case II (E#1/G#1′=2.5,ν#1′=0.25,E#1/E#2=1).
(a)
(b)Variation of normalized SIFs of interfacial crack tip with the anisotropic directionsψ#1 for case III (E#1/G#1′=2.5,ν#1′=0.25,E#1/E#2=1).
(a)
(b)Variation of normalized SIFs of an interfacial crack with an inclined crack angleβ.
(a)
(b)Variation of normalized SIFs of an inclined crack with the inclined angleβ.
(a)
(b)Figure 12
Variation of the ratio ofFI/FII of an inclined crack with the inclined angle β.


## 5. Conclusions

In this study, we developed a single-domain BEM formulation in which neither the artificial boundary nor the discretization along the uncracked interface is necessary. Fedeliński [11] presented the single-domain BEM formulation for homogeneous materials. We combined Chen’s formulation with the Green’s functions of bi-materials derived by Ke et al. [12] to extend it to anisotropic bi-materials. The major achievements of the research work are summarized in the following. (i)
A decoupling technique was used to determine the SIFs of the mixed mode and the oscillation on the interfacial crack based on the relative displacements at the crack tip. Five types of three-node quadratic elements were utilized to approximate the crack tip and the outer boundary. A crack surface and an interfacial crack surface were evaluated using the asymptotical relation between the SIFs. Since the interfacial crack has an oscillation singular behavior, we used a special crack-tip element [26] to capture this behavior.(ii)
Calculation of the SIFs was conducted for several situations, including cracks along and away from the interface. Numerical results show that the proposed method is very accurate, even with relatively coarse mesh discretization. In addition, the use of proposed BEM program is also very fast. The calculation time for each sample was typically 10 s on a PC with an Intel Core i7-2600 CPU at 3.4 GHz and 4 GB of RAM.(iii)
The study of the fracture behavior of a crack is of practical importance due to the increasing application of anisotropic bi-material. Therefore, the fracture mechanisms for anisotropic bi-material need to be investigated.



---

*Source: 721656-2013-03-11.xml*