Abstract
This paper presents a single-domain boundary element method (SDBEM) for linear elastic fracture mechanics analysis in the 2D anisotropic bimaterial. In this formulation, the displacement integral equation is collocated on the uncracked boundary only, and the traction integral equation is collocated on one side of the crack surface only. The complete fundamental solution (Green's function) for anisotropic bi-materials was also derived and implemented into the boundary integral formulation so the discretization along the interface can be avoided except for the interfacial crack part. A special crack-tip element was introduced to capture exactly the crack-tip behavior. A computer program with the FORTRAN code has been developed to effectively calculate the stress intensity factors, crack initiation angle, and propagation path of an anisotropic bi-material. This SDBEM program has been verified having a good accuracy with the previous researches. In addition, a rock of type (1)/(2) disk specimen with a central crack was made to conduct the Brazilian test under diametrical loading. The result shows that the numerical analysis can predict relatively well the direction of crack initiation and the path of crack propagation.
1. Introduction
In rock masses, the interbed construction suffers from cracking which is caused by various factors. Of greater concern are those cracks that develop as a result of initiation and propagation path, leading to significant change in the failure resistance of the structure. Many relative researches and discussions are started like wildfire and never stopped.
Because of the discontinuities of weak interbed on laccoliths (rock mass), there are many damages that occurred which are due to earth stress efforts or under geotechnical engineering, such as slope slip, tunnel collapse, deep excavation collapse, and crack openings which are due to deep well drilling. Crack is one of the fracture models to cause those damages because the crack openings and propagation on the field affect the rock mass structure stability. This paper will discuss the behavior of crack propagation on the basis of the theory of fracture mechanics. We defined the rock sample on interbed as “Bimaterial rock”.
Theoretically, interfacial crack problems in isotropic Bimaterials were studied [1–3] where the authors showed that the stresses possess the singularity of . Rice [4] re-examined the elastic fracture mechanics concepts for the isotropic interfacial crack and introduced an intrinsic material length scale so that the definition of the stress intensity factor (SIFs) possesses the same physical significance as those for the homogeneous cracks. Clements [5] and Willis [6] studied interfacial crack problems in dissimilar anisotropic materials. They showed that the oscillatory behavior of the stresses and the phenomenon of interpenetrating of the crack faces were also present near the crack tip for anisotropic interface cracks. Recent studies on interfacial cracks in anisotropic materials have been conducted by many authors [7–21], and different definitions for the stress intensity factor exist. By introducing a characteristic length, however, the definition given by Gao et al. [15], and Wu [12, 13] is consistent with Rice’s general definition [4] and appears to be more explicit than other definitions.
The study of fracture initiation and propagation in anisotropic rocks is subjected to Brazilian loads. A numerical procedure based on the SDBEM and maximum tensile stress criterion has been developed to predict the angle of crack initiation and the path of crack propagation in anisotropic rocks. Crack propagation in an anisotropic homogeneous rock disc under mixed mode I-II loading is simulated by an incremental crack extension with a piecewise linear discretization. A computer program, which can automatically generate a new mesh, has been developed to simulate the crack propagation process. Some experimental observations of crack initiation angles and crack propagation were obtained by conducting diametrical loading of initially cracked discs of a gypsum/cement. It was found that the numerical analysis could predict relatively well the direction of crack initiation and the crack propagation path.
2. Theoretical Background
2.1. Green’s Function in Anisotropic Bimaterials
With three complex analytical functions , one can, in general, express displacements, stresses, and tractions as follows [11, 22, 23]: where , denotes the real part of a complex variable or function, a prime denotes the derivative, the three complex numbers and the elements of the complex matrices B and A are functions of the elastic properties [11, 22, 23].
Assume that the medium is composed of two joined dissimilar anisotropic and elastic half-planes. We let the interface be along the x-axis and the upper and lower half-planes occupied by materials () and (), respectively (Figure 1).
For concentrated force acting at the point () in material () (), we express the complex vector function as [11] where the vector function with the argument having the generic form .
In (2.2), is a singular solution corresponding to a point force acting at the point () in an anisotropic infinite plane with the elastic properties of material (). This singular solution can be expressed as [11, 23] where , is the point force vector, and H is given by
There are two unknown vector functions to be solved in (2.2), that is, and . While the former is analytic in the upper (material ()) half-plane, the latter is analytic in the lower (material ()) half-plane. These expressions can be found by requiring continuity of the resultant traction and displacement across the interface, along with the standard analytic continuation arguments. Following this approach and after some complex algebraic manipulation, the complex vector functions in materials () and () are obtained as
In (2.6), the special subscripts () and () are used exclusively to denote that the corresponding matrix or vector is in material () () and material () (), respectively.
Similarly, for a point force in material (2.1) (), these complex functions can be found as where the vector functions are the infinite plane solution given in (2.4) but with the elastic properties of material ().
With the complex function given in (2.6) and (2.7), Green’s functions of the displacement and traction can be obtained by substituting these complex functions into (2.1). These Green’s functions have four different forms depending on the relative location of the field and source points.
It is noteworthy that these Green’s function can be used to solve both plane stress and plane strain problems in anisotropic Bimaterials. Although the isotropic solution cannot be analytically reduced from these Green’s functions one can numerically approximate it by selecting a very weak anisotropic (or nearly isotropic) medium [24, 25].
2.2. BEM Formulation for 2D Cracked Anisotropic Bimaterials
In this section, we present an SDBEM formulation in which neither the artificial boundary nor the discretization along the un-cracked interface is necessary. This SDBEM formulation was widely used recently by Chen et al. [26], for homogeneous materials and is now extended to anisotropic Bimaterials.
The displacement integral equation applied to the outer boundary results in the following form ( only, Figure 2): where = 1,2, and are Green’s tractions and displacements, and are the boundary displacement and tractions, respectively, are quantities that depend on the geometry of the boundary and are equal to for a smooth boundary, and and are the field points and the source points on the boundary of the domain. has the same outward normal as . Here, the subscripts and denote the outer boundary and the crack surface, respectively.
The traction integral equation (for being a smooth point on the crack) applied to one side of the crack surfaces is ( only) where is the outward normal at the crack surface and is the fourth-order stiffness tensor.
Equations (2.8) and (2.9) form a pair of boundary integral equations [24, 27, 28] and can be used for the calculation of SIFs in anisotropic Bimaterials. The main feature of the BEM formulation is that it is a single-domain formulation with displacement integral (2.8) being collocated on the un-cracked boundary only and traction integral (2.9) on one side of the crack surface only. For problems without cracks, one needs (2.8) only, with the integral on the crack surface being discarded. Equation (2.8) then reduces to the well-known displacement integral on the un-cracked boundary being discarded. We emphasize here that since bimaterial Green’s functions are included in (2.8), discretization along the interface can be avoided, with the exception of the interfacial crack part which will be treated by the traction integral equation presented by (2.9).
It is well known that a cracked domain poses certain difficulties for BEM modeling (Cruse, 1988 [29]). Previously, fracture mechanics problems in isotropic or anisotropic bimaterials were mostly handled by the multidomain method in which each side of the crack surface is put into different domains and artificial boundaries are introduced to connect the crack surface to the un-cracked boundary. For the bimaterial case, discretization along the interface is also required if one uses the Kelvin-type (infinite domain) Green’s functions.
2.3. Crack-Tip Modeling
In fracture mechanics analysis, especially in the calculation of the SIFs, one needs to know the asymptotic behavior of the displacements and stresses near the crack-tip. In our BEM analysis of the SIFs, we propose to use the extrapolation method of the crack tip displacements. We therefore need to know the exact asymptotic behavior of the relative crack displacement behind the cracktip. This asymptotic expression has different forms depending on the location of the cracktip. In this paper, two cases will be discussed, that is, a crack-tip within the homogeneous material and an interfacial crack-tip. Inclined cracks terminating at the interface will be discussed in a future paper.
2.3.1. A Crack Tip within a Homogeneous Material
The mixed mode stress intensity factors (SIFs) for anisotropic media can be determined by using the extrapolation method of the relative crack displacement (RCD), combined with a set of the shape functions. The RCD is defined as [24]
For this case, the relation of the RCDs at a distance behind the crack tip and the SIFs can be found as [24, 25, 30] where
Substituting the RCDs into (2.10) and (2.11), we obtain a set of algebraic equations in which the SIFs and can be solved.
2.3.2. An Interfacial Crack Tip
For this case, the relative crack displacements at a distance behind the interfacial crack tip can be expressed, in terms of the three SIFs, as [15] where , and are the relative parameters in material () and material (). Utilizing (2.5), we defined the matrix of material as where and are two real matrices, and then utilizing these two matrices, we define matrix as and the characteristic relative to material Then, we used the characteristic obtain to define oscillation index as where is a identity matrix.
Utilizing the relationship between the characteristic and oscillation index , we define constant as where is the characteristic distance along material interface to crack tip.
Comparing (2.10) with (2.13), we noticed that while the relative crack displacement behaves as a square root for a crack tip within a homogeneous medium, for an interfacial crack tip, its behavior is , a square-root feature multiplied by weak oscillatory behaviors.
Equation (2.13) can be recast into the following form, which is more convenient for the current numerical applications: where is the characteristic length and is a matrix function with its expression given by Again, in order to capture the square-root and the weak oscillatory behavior, we construct a crack-tip element with tip at in terms of which the relative crack displacement is expressed as
In meshing the 2D anisotropic bimaterial problem (as shown in Figure 3), we assume that the interface is along the X-axis and the upper () and lower () half-planes are occupied by materials () and (), respectively. The corner of outer boundary is processed by the discontinuous elements Type I and Type III; the continuous element Type V is to deal with all outer smooth boundary; internal crack surface is processed by crack surface elements Type II; and crack tip element Type VI is to process crack tip problem. In order to avoid the oscillatory behavior of the interface, we mesh an anisotropic problem (as shown in Figure 3()) of bimaterial in which neither the interfacial elements nor the discretization along the un-cracked interface is necessary, with the exception of the interfacial crack part which will be treated by traction integral Equation (2.9) and (2.10) of RCD.
2.4. Crack Initiation and Fracture Propagation
In fracture mechanics, there are three criteria commonly used to predict the crack initiation angle: the maximum tensile stress criterion, or -criterion [31], the maximum energy release rate criterion, or G-criterion [32], and the minimum strain energy density criterion, or S-criterion [33]. Among them, the -criterion has been found to predict well the directions of crack initiation compared to the experimental results for polymethylmethacrylate [34, 35] and brittle clay [36]. Because of its simplicity, the -criterion seems to be the most popular criterion in mixed mode I-II fracture studies [37]. Therefore, the -criterion was used in this paper to determine the crack initiation angle for anisotropic plates.
For anisotropic materials, the general form of the elastic stress field near the crack tip in the local Cartesian coordinates of Figure 4 can be expressed in terms of the two stress intensity factors and a follows [30]:
Using coordinate transformation, the stress fields near the crack tip in the polar coordinates () of Figure 4 are
If the maximum -criterion is used, the angle of crack initiation, , must satisfy
A numerical procedure was applied to find the angle when is a maximum for known values of the material elastic constants, the anisotropic orientation angle , and the crack geometry.
In this paper, the process of crack propagation in an anisotropic homogeneous plate under mixed mode I-II loading is simulated by incremental crack extension with a piecewise linear discretization. For each incremental analysis, crack extension is conveniently modeled by a new boundary element. A computer program has been developed to automatically generate new data required for analyzing sequentially the changing boundary configuration. Based on the calculation of the SIFs and crack initiation angle for each increment, the procedure of crack propagation can be simulated. The steps in the crack propagation process are summarized as follows (Figure 5):(1)compute the SIFs using the proposed BEM;(2)determine the angle of crack initiation based on the maximum tensile stress criterion;(3)extend the crack by a linear element (of length selected by the user) along the direction determined in step ;(4)automatically generate the new BEM mesh;(5)repeat all the above steps until the crack is near the outer boundary.
3. Numerical Examples
In this section, a lot of numerical examples are presented to verify the formulation and to show the efficiency and versatility of the present SDBEM for problems related to fracture in anisotropic Bimaterial.
3.1. Stress Intensity Factors (SIFs)
3.1.1. Horizontal Crack in Material ()
A horizontal crack under a uniform pressure is shown in Figure 6. The crack has a length , and is located at a distance to the interface. The Poisson ratios for both materials () and () the same, that is, = = 0.3, while the ratio of the shear module ratio varies. A plane stress condition is assumed. In order to calculate the SIFs at crack tip () or (), 20 quadratic elements were used to discretize the crack surface. The results are given in Table 1 for various values of the shear module ratio. They are compared to the results given by Isida and Noguchi [38], using a body force integral equation method and those by Yuuki and Cho [39], using a multidomain BEM formulation. As can be observed from this table, the results compare quite well.
3.1.2. Vertical Crack in Material ()
Consider a vertical crack in material () subjected to far-field horizontal stresses as shown in Figure 7. The crack has a length and is located at a distance to the interface. The Poisson ratios for both materials () and () are the same, that is, = 0.35, = 0.3, and the ratio of the shear module is the same, that is, = 23.077, while the ratio of the crack length located at a distance () varies. A plane stress condition is assumed. In order to calculate the SIFs at crack tip () or (), 20 quadratic elements were used to discretize the crack surface. The results are given in Table 2 for various values of the crack length located at a distance. They are compared to the results given by Isida and Noguchi [38] using a body force integral equation method and those by Cook and Erdogan [41] using a Wiener-Hopf technique and an asymptotic analysis. As can be observed from this table, the results compare quite well.
3.1.3. Vertical Crack Intersecting an Interface
Consider a vertical crack intersecting an interface and subjected to far-field horizontal stresses as shown in Figure 8. The horizontal far field stresses applied in materials () and () are, respectively, and (=). The Poisson ratio and are assumed to be equal to 0.3 and the shear modulus ratio is assumed to vary. The distance of the crack tip () and () to the interface are the same, that is, == , the half-length of the crack. Again, a plane stress condition is assumed and 20 quadratic elements were used to discretize the crack surface. The SIFs at the crack tips () and () are listed in Table 3 for several values of the shear modulus ratio, and are compare to those given by Isida and Noguchi [38]. Again, the results between the two numerical analyses compare quite well.
3.1.4. An Interface Kinked Crack in Infinite Bimaterials
Consider an interface kinked () crack in infinite Bimaterials subjected to far-field tensile stresses as shown in Figure 9. The far-field tensile stresses applied in materials () and () are, respectively, and (=). The kink crack length is ; the main crack length is . Kink crack’s tip is (). The Poisson ratios and are assumed to be equal to 0.3. The Young module in material () is 1 GPa, and according to the shear module ratio with material () and material () ( = 0.25), we can get the Young module of material (). Then according to equation , we can get the relation between kink length and stress intensity factors. After numerical analysis, we compared with the results of Isida and Noguchi [38], using a body integral equation force method, as shown in Table 4, the results the two numerical analyses and found that they compare quite well.
3.1.5. Interfacial Horizontal Crack in Infinite Bimaterials
An interfacial crack along the x-axis of length is shown in Figure 10. The crack surface is under a uniform pressure and the materials can be either isotropic or anisotropic. Twenty quadratic elements were used to discretize the crack, and the characteristic length is assumed as .
The SIFs at the crack tip of an interfacial crack were also calculated for the anisotropic Bimaterial case. The anisotropic elastic properties in material () were assumed to be those of glass/epoxy with = 48.26 GPa, = 17.24 GPa, = 6.89 GPa, and = 0.29. For material (), a graphite/epoxy with = 114.8 GPa, = 11.7 GPa, = 9.66 GPa, and = 0.21 was selected [42]. The material axis in material () and material () makes angles and , respectively, with respect to the horizontal direction (Figure 10). While the material axis in material () was assumed to be along the horizontal direction (i.e., = 0), the -axis in material () makes an angle with respect to the horizontal direction. The interfacial SIFs at crack tip () obtained by the present method are listed in Table 5 and compared to the exact solutions proposed by Wu [12]. A very good agreement is found between the numerical analysis and the exact solution.
3.1.6. Interfacial Horizontal Crack in Finite Bimaterials
The example is included as a comparison with the literature in order to demonstrate the accuracy of SDBEM approach for an interfacial crack in anisotropic bimaterial plate. The geometry is that of rectangular plate and is shown in Figure 11. For the comparison, crack length is taken as = 2, = , and = 0.4, and static tensile loading is applied on the upper and the lower boundary of the plate. Plane stress condition is assumed. The anisotropic elastic properties for materials () and () are given in Table 6. The normalized complex stress intensity factors at crack tip () or () are listed in Table 7 together with those from the work of Cho et al. [43], who used a multidomain BEM formulation and the results from Wünsche et al. [21] for a finite body. The outer boundary and interfacial crack surface were discretized with 40 continuous and 20 discontinuous quadratic elements, respectively. It is obvious from Table 7 that these are very close to those obtained by the other researchers.
3.2. Crack Initiation Angles
3.2.1. Comparison of Numerical Predictions of Crack Initiation Angles with Experimental Results
In this section, we compared the numerical predictions of crack initiation angles with experimental results to verify the formulation. Al-Shayea [44] conducted uniaxial pressure test on limestone rock discs 98 mm and 84 mm in diameter and 22 mm in thickness with 30 mm notch. The crack orientation angle between the crack plane and the tensile stress was varied. Figure 12 shows the variation of the crack initiation angle with the crack angle determined numerically and experimentally. A good agreement is found between the experimental results of Al-Shayea [44] and our numerical predictions.
4. Comparison of Numerical Predictions of Crack Propagation Paths with Experimental Results
To demonstrate the proposed SDBEM procedure when predicting crack propagation in the anisotropic Bimaterials under mixed mode I–II loading, the propagation path in a CSTBD specimen is numerically predicted and compared with the actual laboratory observations. In these experiments, a crack initially inclined with respect to the applied stress is allowed to grow under concentrated diametrical loading (as shown in Figure 13). The Brazilian tests on CSTBD specimens with a diameter of 7.4 cm, a thickness of 1 cm, and a crack length of 2.2 cm are conducted to observe the actual propagation paths and are compared with the numerical predictions. The anisotropic elastic properties for rocks of types () and () are given in Table 8. The ratios of and are 1.635 and 4.301, respectively. Since the value of = 1.635, this rock of type () can be classified as a slightly anisotropic rock. Two specimens with the material inclination angle = and , defined as the SSD-1 and SSD-2, have crack angles = and , respectively. After Brazilian tests with cracked bimaterial specimens, the paths of crack propagation for SSD-1 and SSD-2 are shown in Figures 14(a) and 15(a), respectively. It can be observed that the crack propagates nearly perpendicular to the crack surface in the initial stage and then rapidly approaches toward the loading point. The proposed SDBEM procedure is also used to simulate crack propagation in the CSTBD specimens. The outer boundary and crack surface are discretized with 28 continuous and 20 discontinuous quadratic elements, respectively. Figures 14(b) and 15(b) are the comparisons examples; it can be concluded that the proposed SDBEM is capable of predicting the crack propagation in anisotropic bimaterial rocks.
(a)
(b)
(a)
(b)
(a)
(b)
5. Conclusions
This paper shows that the mixed mode stress intensity factors of anisotropic Bimaterial rock under diametrical loading can be successfully determined by the SDBEM.
A new SDBEM procedure based on the maximum tensile stress failure criterion was developed to predict the crack initiation direction and the crack propagation path in anisotropic Bimaterial rock discs under mixed mode loading. A good agreement was found between crack initiation angles and propagation paths predicted with the SDBEM and experimental observations reported by previous researchers on anisotropic materials. Numerical simulations of crack initiation and propagation in CSTBD specimens of type ()/() were also found to compare well with the experimental results. Since the present method is simple and can be used for curved cracks, it will be straightforward to extend the current SDBEM formulation to analyze fracture propagation in 2D anisotropic Bimaterials, which is currently under investigation by the authors.