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

1 Department of Resources Engineering, National Cheng Kung University, No. 1 University Road, Tainan City 70101, Taiwan 2Graduate Institute of Applied Geology, National Central University, No. 300 Jhongda Road, Jhongli City 32001, Taiwan 3 Graduate Institute of Geophysics, National Central University, No. 300 Jhongda Road, Jhongli City 32001, Taiwan 4Geotechnical Engineering Research Center, Sinotech Engineering Consultants, Inc., Taipei 110, Taiwan


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  −1/2± , where  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][8][9][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][12][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 singledomain 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][17][18].

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,   (  ) of the complex variables   =  +    ( = 1,2) where   is the root of the following characteristic equation [19]: where the coefficient   is the compliance component calculated using the - coordinate system.The detailed relationships of these components with material elasticity can be found in Chen et al. [14].If the root   of ( 1) is assumed to be distinct, the general expressions for the stress and displacements are where the elements of the complex matrices  are The traction components in the  and  directions are Here, the complex analytical functions   (  ) can in general express ( 2) and ( 4) as follows [19][20][21]: where   =  +   ; Re denotes the real part of a complex variable or function; a prime denotes the derivative; the complex number   ( = 1, 2); the elements of complex matrices  are defined in (3); and the elements of complex matrices  can be defined as Assume that the medium is composed of two joined dissimilar anisotropic and elastic half-planes.Let the interface be along the -axis and let the upper ( > 0) and lower ( < 0) half-planes be occupied by materials #1 and #2, respectively (as shown in Figure 1).
Considering a concentrated force acting at the source point ( 0 ,  0 ) in material #2( 0 < 0), we express the complex vector function as [20] where the vector function is where the argument has the generic form  =  + .In (7),  0 #2 is a singular solution corresponding to a point force acting at the point ( 0 ,  0 ) on an anisotropic infinite plane with the elastic properties of material #2.This singular solution can be expressed as [20,21] where   =  0 +    0 ,   is the magnitude of the point force in the -direction, and where  = √ −1, the overbar indicates the complex conjugate, and superscript −1 indicates the matrix inverse.Two unknown vector functions,  #1 () and  #2 (), 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 Therefore, the complex vector functions can be expressed as In (12), subscripts #1 and #2 are used to denote that the corresponding matrix or vector is for material #1 ( > 0) and material #2 ( < 0), respectively.
Similarly, for a point force in material #1 ( 0 > 0), the complex vector functions can be expressed as where vector function  0 #1 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,  *  is the Green's function for displacement and  *  is the Green's function for traction.Superscripts and subscripts #1 and #2 are used to denote the corresponding quantities for materials #1 ( > 0) and #2 ( < 0), respectively [15].
(i) For source point () and field point () on material #1 ( > 0), where the matrix  is defined in (10) with the anisotropic elastic properties of material #1, and (ii) For source point () and field point () on material #1 ( > 0) and field point () on material #2 ( < 0), with (iii) For source point () and field point () on material #2 ( < 0), where matrix  is defined in (10) with the anisotropic elastic properties of material #2, and (iv) For source point () and field point () on materials #2 ( < 0) and field point () on material #1 ( > 0), with 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].

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 singledomain 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 point  0 , on the uncracked boundary, the displacement integral equation applied to the outer boundary results in the following form ( 0 , ∈ Γ  only, as shown in Figure 1): where , ,  = 1, 2;  *  and  *  are the Green's tractions and displacements given in ( 14), ( 16), (18), and (20);   and   are the boundary displacements and tractions, respectively; and   and  0 , are the field points and the source points on boundary Γ of the domain, respectively.Γ  has the same outward normal as Γ + .  are coefficients that depend only on the local geometry of the uncracked boundary at point  0 , and are equal to   /2 for a smooth boundary.Here, subscripts  and  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 (for  0  being a smooth point on the crack) applied to one side of the cracked surface is ( where   is the unit outward normal at cracked surface  0 , and  lim  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 (  ) are determined using the following expression: 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].

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 cracktip 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  behind the crack tip can be expressed in terms of the three SIFs as [17] where  = [ II ,  I ,  III ]  is the SIF vector and  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)): where subscript  denotes the RCD component and superscript  (1, 2, 3) denotes the RCDs at nodes  = −2/3, 0, 2/3 (as shown in Figure 2(b)), respectively.The shape functions   are those introduced by [25] In this case, the relation of the RCDs at a distance  behind the crack tip and the SIFs can be found as [26][27][28] where Substituting the RCDs into ( 26) and ( 28), we may obtain a set of algebraic equations with which the SIFs  I and  II can be solved.
(ii) Interfacial Crack Tip.In this case, the relative crack displacements at a distance  behind the interfacial crack tip can be expressed in terms of the three SIFs as [29] where   ,   ,   , and  are the relative parameters in materials #1 and #2.Utilizing (10), we defined the matrix of a material as where  and  are two real matrices.These two matrices are used to define matrix  as The characteristic  relative to material is as follows: We use the characteristic  to define oscillation index  as where  is a 3 × 3 identity matrix.
The relationship between characteristic  and oscillation index  is used to define constant   as where  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  (1/2)+ , 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: where  is the characteristic length and  is a matrix function expressed as 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 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 Figure 2, during analysis.For these elements, both the geometry and the boundary quantities are approximated by intrinsic coordinate , and shape function   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.

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.  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.

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 bimaterial plate.The geometry is that of a rectangular plate, as shown in Figure 4.For the comparison, the crack length is 2, ℎ = 2, / = 0.4, and static tensile loading  is applied to the upper and the lower boundaries of the plate.
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.

Numerical Analysis
After its accuracy was checked, the BEM program was applied to more challenging problems involving an anisotropic bimaterial 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,  I and  II are the normalized SIFs; for  I > 0, the failure mechanism is tensile fracture, and for  II > 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,  I / II > 1, and shear fracture,  I / II < 1, resp.).

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  I and  II , are equal to with 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 farfield vertical tensile stresses, as shown in Figure 5.The crack, which had a length of 2, was located at distances / = 1 (material #1), −1 (material #2), and / = 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 -axis.Material #1 had five independent elastic constants ( #1 ,   #1 , ] #1 , ]  #1 , and  #1 ) in the local coordinate system;  # and   # are the Young's moduli in the plane of transverse isotropy and in a direction normal to it, respectively; ] # and ]  # 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   # is the shear modulus in the planes  normal to the plane of transverse isotropy.For material #2,  #2 /  #2 = 1, #2 /  #2 = 2.5, and ] #2 = ]  #2 = 0.25.In this problem, the anisotropic direction () within material #1 varies from 0 ∘ to 180 ∘ .For all cases, three sets of dimensionless elastic constants are considered.They are defined in Table 4. Here,  #1 / #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 ( #1 /  #1 = 1,  #1 /  #1 = 2.5, and ]  #1 = 0.25) are shown as solid lines in the figures for comparison.

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   / = 0.25, respectively.The same consider an infinite plate of bimaterials, having the elastic constants ratio  #1 /  #1 = 1,  #1 /  #1 = 0.4, and ]  #1 = 0.2499 for material #1 and  #2 /  #2 = 1/3,  #2 /  #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  #1 / #2 = 1/5.One considers an interfacial crack of length 2 and a horizontal crack situated near the interfacial crack with ah.The ratio / ℎ is equal to 1.The longitudinal distance   / = 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.(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.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 Figure 11(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  I and a weak influence on the variations of  II , 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  I different of zero whatever the value of  and  II = 0 at  = 80 ∘ .
(x) The ratio of the normalized SIFs ( I / II ) is concerned; Figure 12 shows that the fracture of shear mode ( I / II < 1) occurs to the inclined crack angle is between 30 ∘ and 150 ∘ (tip C).For tip D, the shear mode ( I / II < 1) occurs to the  is between 50 ∘ and 95 ∘ , and between 105 ∘ and 120 ∘ , respectively.And the opening mode ( I / II > 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.

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.

Figure 1 :
Figure 1: Definition of the coordinate systems in an anisotropic bimaterial.

Figure 2 :
Figure 2: Three-geometric-node quadratic elements used to approximate the uncracked boundary and the crack for two-dimensional problems.

Figure 3 :Figure 4 :
Figure 3: Horizontal crack under uniform pressure in material #1 of an infinite bi-material.

Figure 6 :
Figure 6: A crack situated near the interfacial crack within infinite bi-materials.

Figures 10 and 11
Figures 10 and 11 represent the variations of normalized SIFs ( I and  II ) of interfacial crack and inclined crack according to the inclined angle , respectively.

4. 3 .
Discussion.An analysis of Figures 7 to 9 reveals the following major trends.

Figure 10 :Figure 11 : 1 Figure 12 :
Figure 10: Variation of normalized SIFs of an interfacial crack with an inclined crack angle .
Material #1.A horizontal crack under uniform pressure  is shown in Figure3.The crack, whose length is 2, is located distance  from the interface.The Poisson's ratios for materials #1 and #2 are ] #1 = ] #2 = 0.3, and the ratio of the shear modulus  #2 / #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

Table 3 :
Comparison of the normalized complex SIFs for finite anisotropic problem (interfacial crack).