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.
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 , 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–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, of the complex variables 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 functionscan in general express (2) and (4) as follows [19–21]: where ; denotes the real part of a complex variable or function; a prime denotes the derivative; the complex number ; 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 and lower half-planes be occupied by materials and , respectively (as shown in Figure 1).
Considering a concentrated force acting at the source point in material , we express the complex vector function as [20] where the vector function is where the argument has the generic form .
In (7), is a singular solution corresponding to a point force acting at the point on an anisotropic infinite plane with the elastic properties of material . This singular solution can be expressed as [20, 21] where , is the magnitude of the point force in the -direction, and where , the overbar indicates the complex conjugate, and superscript −1 indicates the matrix inverse. Two unknown vector functions, and , are solved in (7). The former is analytic in material , and the latter is analytic in material . 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 and are used to denote that the corresponding matrix or vector is for material and material , respectively.
Similarly, for a point force in material #1 , the complex vector functions can be expressed as where vector function is the infinite-plane solution given in (9), but with the elastic properties of material .
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 and are used to denote the corresponding quantities for materials and , respectively [15].(i)For source point and field point on material , where the matrix is defined in (10) with the anisotropic elastic properties of material , and (ii)For source point and field point on material and field point on material , with (iii)For source point and field point on material , where matrix is defined in (10) with the anisotropic elastic properties of material , and (iv)For source point and field point on materials and field point on material , 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].
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 point on the uncracked boundary, the displacement integral equation applied to the outer boundary results in the following form ( only, as shown in Figure 1): where 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 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 and are equal to 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 being a smooth point on the crack) applied to one side of the cracked surface is ( only) where is the unit outward normal at cracked surface and 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].
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 . 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 is the SIF vector and is a matrix with elements related to the anisotropic properties in material , 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 (as shown in Figure 2(e)): where subscript denotes the RCD component and superscript denotes the RCDs at nodes (as shown in Figure 2(b)), respectively. The shape functions are those introduced by [25]
(a) Discontinuous quadratic element of Type I
(b) Crack surface quadratic element of Type II
(c) Discontinuous quadratic element of Type III
(d) Discontinuous quadratic element of Type IV
(e) Crack tip quadratic element of Type V
In this case, the relation of the RCDs at a distance behind the crack tip and the SIFs can be found as [26–28] where
Substituting the RCDs into (26) and (28), we may obtain a set of algebraic equations with which the SIFs and 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 and . 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 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 , 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.
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
A horizontal crack under uniform pressure is shown in Figure 3. The crack, whose length is , is located distance from the interface. The Poisson’s ratios for materials and are , and the ratio of the shear modulus 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.
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 Figure 4. For the comparison, the crack length is , , 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 and 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.
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, and are the normalized SIFs; for , the failure mechanism is tensile fracture, and for , 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, , and shear fracture, , 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 and , 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 , (ii) a horizontal crack within material , 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 , was located at distances (material ), −1 (material ), and (interfacial crack) from the interface, respectively. Material was transversely isotropic with the plane of transverse isotropy inclined at angle with respect to the -axis. Material had five independent elastic constants , and 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 , , and . In this problem, the anisotropic direction within material varies from to .
| (a) Horizontal crack within material |
| (b) Horizontal crack within material |
(c) Interfacial crack within bi-material
For all cases, three sets of dimensionless elastic constants are considered. They are defined in Table 4. Here, , which is the ratio of the Young’s modulus of material to the Young’s modulus of material , 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 (, , and ) are shown as solid lines in the figures for comparison.
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 , respectively. The same consider an infinite plate of bi-materials, having the elastic constants ratio , and for material and , , and for material . The Poisson ratios are , material orientation , and the Young modulus ratio . One considers an interfacial crack of length and a horizontal crack situated near the interfacial crack with ah. The ratio is equal to 1. The longitudinal distance , 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 and of interfacial crack and inclined crack according to the inclined angle , respectively.
(a)
(b)
4.3. Discussion
An analysis of Figures 7 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 is symmetric with respect to . 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 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 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 for bi-materials with a high degree of anisotropy .(v) Figure 7(a) 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., or ) for case 1 (crack within material ). Figure 9(b) shows that the minimum values of the SIFs (mode II) occur when the anisotropic direction angle is about for case 3 (interfacial crack). Analysis of Figures 10 and 11 reveals that several major trends can be underlined as follows.(vi) Figure 10(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 . 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 . 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 , the SIFs start decreasing and stabilizing when approaching .(vii) The same phenomenon was observed in Figure 10(b). It is shown that the increase in the inclined angle causes the reduction in SIFs (mode II) until reaching a value at , and starts an immediate sharp increasing. The maximum of SIFs (mode II) is reached for . 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 , the SIFs start a slight increasing and stabilizing when approaching .(viii) Figure 11 illustrates, respectively, the variations of the normalized SIFs of tip C and D according to inclined angle . For Figure 11(a), between and , one observes that mode I of the normalized SIFs decreases with . The minimal value is reached at . From , the two normalized SIFs increase until , and the normalized SIFs at the tip D reaches a maximum value. In the interval of at , 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 and , 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 and and decrease between and . In the interval of at the normalized SIFs of tip C are slightly higher than those of tip D, and the inverse phenomenon is noted for angles ranging between , and . These results permit to note on one hand that the orientation of the inclined crack has influence on the variations of and a weak influence on the variations of , 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 . The presence of the bi-material interface gives different of zero whatever the value of and at .(x) The ratio of the normalized SIFs is concerned; Figure 12 shows that the fracture of shear mode occurs to the inclined crack angle is between and (tip C). For tip D, the shear mode occurs to the is between and , and between and , respectively. And the opening mode occurs to the is between and . 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.
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
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.
Acknowledgments
The authors would like to thank the National Science Council of the Republic of China, Taiwan, for financially supporting this research under Contract no. NSC-100-2221-E-006-220-MY3. They also would like to thank two anonymous reviewers for their very constructive comments that greatly improved the paper.