- Email: [email protected]
Numerical methods for the determination of multiple stress singularities and related stress intensity coefficients
Engineering Fracture Mechanics 63 (1999) 775±790
www.elsevier.com/locate/engfracmech
Numerical methods for the determination of multiple stress singularities and related stress intensity coecients Jin-Quan Xu a,*, Yi-Hua Liu b, Xiao-Gui Wang a a
b
Department of Mechanics, Zhejiang University, Hangzhou 310027, People's Republic of China Department of Applied Mathematics and Mechanics, Hefei University of Technology, Hefei 230009, People's Republic of China Received 3 December 1998; accepted 3 March 1999
Abstract This paper proposed numerical methods to determine the multiple stress singularities (including the oscillatory stress singularities) and the related stress intensity coecients, by the use of common numerical solutions (stresses or displacements) obtained by an ordinary numerical tool such as ®nite element method (FEM) or boundary element method (BEM). To verify the eciency of the present methods, two models of bonded dissimilar materials under the plane strain state are analyzed by BEM, and the orders of the stress singularities and the related stress intensity coecients are examined numerically. The results show that all the orders of the stress singularities at an interface edge can be determined precisely by the present method, and the related stress intensity coecients can be determined by the extrapolation method with a very good linearity. It is found that the methods presented in this paper are very simple and ecient. Moreover, they can be easily extended to any singular problem. # 1999 Elsevier Science Ltd. All rights reserved. Keywords: Multiple stress singularities; Stress intensity coecients; Common numerical solutions; Extrapolation method; Bonded dissimilar materials
1. Introduction It is well known that there may appear multiple singularities or oscillatory singularities at a stress singular point such as an interface edge, an interface corner, the tip of an interface crack * Corresponding author. Tel.: +86-571-7951769; fax: +86-571-7951464. E-mail address: [email protected] (J.-Q. Xu) 0013-7944/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 3 - 7 9 4 4 ( 9 9 ) 0 0 0 4 4 - 2
776
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
(or a crack terminating at the interface) in bonded dissimilar materials, as well as the tip of a V-notch in a homogenous material (e.g., Refs. [1±5]). When the stress singularity is not single (i.e., it is multiple), the orders of the stress singularities can be only determined generally by analytic methods. However, there are many cases to which the theoretical results are not available, or they are too dicult to analyze theoretically, so that numerical methods to determine the multiple stress singularities and the related stress intensity coecients are strongly expected in applied mechanics and engineering. Recently, some numerical methods [6± 9] have been proposed. For example, Somaratna and Ting [7] proposed a particularly constructed ®nite element method in which the orders of the multiple or oscillatory stress singularities can be determined by the eigenequation in the form of homogeneous linear algebraic equations. On the other hand, various numerical methods [10±16] as well as photoelastic [17] and strain gage [18] methods have been reported to determine the stress intensity factors of an interface crack. However, for the interface edge, it is not easy to determine stress intensity coecients related to the multiple singularities. Munz and Yang [19] presented a numerical method to calculate stress intensity coecients of the interface edge based on the known analytic angular functions. Leslie et al. [20] gave out an approach based on the in¯uence function. For a V-notch problem, Zhao and Hahn [21] proposed a method to determine the stress intensity factors from the results of a mixed mode crack. Chen [22] proposed a method based on the body force method. Some additional related investigations are also available in the literature. In general, the extrapolation procedure which was developed by Chan et al.[23] and extended later by Carpenter [24] could be considered as a useful tool to determine the stress intensity coecients for multiple singularities problems if some modi®cations are introduced. The purpose of this study is to set up the numerical methods by which the orders of the multiple stress singularities (or oscillatory singularities) and the related stress intensity coecients can be determined from the common numerical solutions (stresses or displacements) obtained by ordinary numerical tools such as ®nite element method (FEM) and boundary element method (BEM). The method would be applicable even if there is no information about the theoretical solutions. Two bi-material models under the plane strain state are analyzed by BEM to illustrate the application of the proposed methods. The numerical results show that all the orders of the multiple stress singularities (including oscillatory singularities) can be determined precisely by the present numerical method, and the related stress intensity coecients can also be determined by the extrapolation method proposed in this paper with a very good linearity. Moreover, a test problem is also analyzed to investigate the accuracy of the present extrapolation method. The results show that the present algorithm is ecient and accurate enough too.
2. Numerical analytical methods 2.1. Determination of multiple stress singularities For simplicity, the problem with two stress singularities is considered here. Problems with more than two singularities can be considered in the same way. By the way, the order of the
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
777
single stress singularity can be determined directly by the log±log (double logarithmic) diagram of the stress and distance from the singular point. In the plane problem, the singular stress ®eld near the singular point with two stress singularities can be described by (in the polar coordinates (r,y)) sij K1 rÿl1 fij1 y K2 rÿl2 fij2 y
1
where lk 0
2
where Hk Kk Lÿlk (k 1,2). In order to determine l1 and l2 , we select arbitrarily two stress components s1 and s2 , then in the selected direction y y0 , we have sl H1l r^ÿl1 H2l r^ÿl2 ,
l 1,2,
3
where Hkl Hk flk y0 ,
k,l 1,2:
4
After some algebraic computations of Eq. (3) we have s1 ÿ A3ÿk s2 Bk r^ÿlk , k 1,2,
5
where A1 H11 =H12 ,
A2 H21 =H22 , B1 H12 A1 ÿ A2 , B2 H22 A2 ÿ A1 ,
6
are coecients to be determined. From Eq. (5), it can be found that the double logarithmic distributions log s1 ÿ A3ÿk s2 0log r^ (k 1,2) should be linear in the dominant region of the singular stress ®eld (1), and the slopes are related to the orders of the stress singularities l1 and l2 , respectively. However, to draw the double logarithmic diagram, coecients A1 and A2 should be predetermined. In order to determine A1 and A2 numerically, we select three points distributed as an equally proportional series in the selected direction y y0 , i.e., rn rnÿ1 r1 ; n 1, 2, 3 (r is a proportional factor and r1 is the distance between the ®rst point and the singular point). The (l 1,2; n 1,2,3), which can be corresponding stress components are expressed by s n l obtained by the ordinary numerical procedure such as FEM or BEM. Substituting these stress k yields components into Eq. (5) and eliminating Bk r^ÿl n aA2k bAk c 0, k 1,2, where
7
778
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
2 a s 22 ÿs 21 s 23 , b s 11 s 23 s 13 s 21 ÿ 2s 12 s 22 , c
s 12
2
ÿs 11 s 13 ,
8
are coecients which can be determined by the numerical results. The roots of Eq. (7) may have the following three cases: 1. b2 ÿ 4ac 0 In this case A1 A2 ÿb= 2a. This leads to l1 l2 , so that it is the case with the single stress singularity. 2. b2 ÿ 4ac > 0 In this case, A1 and A2 are all real, i.e., k p b ÿ ÿ 1 b2 ÿ 4ac , k 1,2: 9 Ak ÿ 2a This leads to two dierent real singularities. 3. b2 ÿ 4ac<0 In this case A1 and A2 are a couple of conjugate complexes, i.e., p b ÿ i 4ac ÿ b2 A1 A 2 AR iAI ÿ , 10 2a p where i ÿ1. This leads to a couple of the complex conjugate orders of the stress singularities which mean oscillatory stress singularities. It should be noted that in this case lk , Bk , Kk , Hk , Hkl , flk y (k,l 1,2) are also complex conjugate, i.e., G1 G 2 GR iGI ,
11
in which Gk lk ,Bk ,Kk ,Hk ,Hkl ,flk y. Substituting Eq. (11) into Eq. (5) yields seics BeicB r^ÿ lR ilI ,
12
where
q s s1 ÿ AR s2 2 AI s2 2 ,
q B B 2R B 2I ,
cs tanÿ1 AI s2 = s1 ÿ AR s2 , cB tanÿ1 BI =BR :
13
From Eq. (12), it can be found that the double logarithmic distribution log s0log r^ and the single logarithmic distribution cs log e0log r^ should be linear in the dominant region of the
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
779
oscillatory singular stress ®eld, and their slopes are related to the real part lR and the imaginary part lI of the order l1 of the oscillatory stress singularity, respectively. To obtain accurate values of coecients A1 and A2 , a few sets with three points distributed as the equally proportional series should be considered. Due to the computational error, A1 and A2 obtained by each set may be slightly dierent. However, because A1 and A2 should be constants in the selected direction, they can be determined by the least square method from these slight dierent values. 2.2. Determination of stress intensity coecients 2.2.1. The case with two real singularities Assume that the stress intensity coecients K1 and K2 can be de®ned by the stresses sy and try in the direction y y0 (e.g., the interface), i.e., sy0 sy jyy0 K1 rÿl1 K2 rÿl2 fy2 y0 , try0 try jyy0 K1 rÿl1 fry1 y0 K2 rÿl2 :
14
This de®nition of the stress intensity coecients means the angular functions fy1 y and fry2 y are normalized in the direction y y0 , i.e., fy1 y0 fry2 y0 1. Taking s1 sy0 , s2 try0 , and according to Eq. (4), we have H11 H1 , H22 H2 . Therefore, from Eqs. (5) and (6), the stress intensity coecients can be determined by the following extrapolation method ÿ K1 Alimrl1 sy0 ÿ A2 try0 , r40
ÿ K2 Alimrl2 try0 ÿ Aÿ1 1 sy0 , r40
15
where A A1 = A1 ÿ A2 , the coecients A1 and A2 have been given by Eq. (9). The linearity of this extrapolation method would be investigated numerically in Section 3. 2.2.2. The case with a couple of complex conjugate singularities In this case, substituting Eq. (11) into (2) yields sy rÿlR k1 Fy1 y cos lI ln r^ k2 Fy2 y sin lI ln r^ , try rÿlR k1 Fry1 y sin lI ln r^ k2 Fry2 y cos lI ln r^ ,
16
where the stress intensity coecients k1 and k2 as well as the angular functions Fy1 y, Fy2 y, Fry1 y, and Fry2 y are k1 HR LlR , k2 HI LlR ,
17
780
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
Fy1 y 2 fyR y ÿ HfyI y ,
Fy2 y 2 fyR y H ÿ1 fyI y ,
Fry1 y 2 HfryR y fryI l ,
Fry2 y 2 H ÿ1 fryR y ÿ fryI y ,
18
in which H HI =HR . We also assume that k1 and k2 can be de®ned by the stresses sy and try in the direction y y0 , i.e., sy0 sy jyy0 rÿlR k1 cos lI ln r^ k2 Fy2 y0 sin lI ln r^ , try0 try jyy0 rÿlR k1 Fry1 y0 sin lI ln r^ k2 cos lI ln r^ :
19
This de®nition means the angular functions Fy1 y and Fry2 y are normalized in the direction y y0 , i.e., Fy1 y0 Fry2 y0 1. Taking s1 sy0 and s2 try0 , then from Eqs. (4), (6) and (11) we get BR A HR AI ÿ 2H1I AR ÿ2AI H2I , BI A HR AR 2H1I AI AI HI ,
20
where A AI = A2R A2I . After the substitution of Eqs. (11) and (20) into (5) and some algebraic computations, we can ®nally obtain the following extrapolation method: n lR ÿ k1 Aÿ1 sy0 ÿ AR try0 AR sin lI ln r^ AI cos lI ln r^ AI try0 AR cos lI lim r I r40
ln r^ ÿ AI sin lI ln r^
o ,
lR ÿ k2 Aÿ1 sy0 ÿ AR try0 sin lI ln r^ AI try0 cos lI ln r^ , I lim r r40
21
where the coecients AR and AI have been given by Eq. (10). The linearity of this extrapolation method would be discussed numerically in Section 3 too. For the interface crack, because Fry1 y0 ÿFy2 y0 1 (assume that the interface is in the direction y y0 ) and lR 1=2, from Eqs. (4), (6) and (11), we have AR 0 and AI ÿ1, so that Eq. (21) can be simpli®ed as p k1 lim r sy0 cos lI ln r^ try0 sin lI ln r^ , r40
p k2 lim r try0 cos lI ln r^ ÿ sy0 sin lI ln r^ : r40
22
It is obvious that the dierence between the de®nition of Eq. (22) and Rice's [25] de®nition p about the stress intensity factors of an interface crack is only a proportional factor 2p.
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
781
3. Numerical results To verify the eciency of the above numerical methods, the orders of the stress singularities and the related stress intensity coecients of two models of the interface edge under plane strain state are calculated numerically by using the solutions of the BEM. The program BEM2D [26] used here, is especially ecient for the analysis of bonded dissimilar materials, which was developed by the author. A test problem is also analyzed to check the accuracy of the present extrapolation method for stress intensity coecients. 3.1. Two real singularities The analytical model is shown in Fig. 1. The material constants are E1 108 GPa, n1 0:33; E2 304 GPa, n2 0:27. The boundary element mesh near the interface edge A is shown in Fig. 2. The nodes near the interface edge A are arranged in the equally proportional series distribution (the proportional factor r 2) and the minimum element length is about 2:0 10ÿ7 W. The double logarithmic distributions of the dimensionless stress components sy and txy on the interface near the interface edge A are shown in Fig. 3. From Fig. 3 it can be found that the double logarithmic distributions are not parallel to each other. If one takes the slopes as the stress singular order, dierent values can be obtained by adopting dierent stress components. Moreover, though the distribution seems to be a line, the slope determined by dierent regions shows that it is a curve in fact. Therefore, it is the case of multiple singularities. According to Bogy [1], the theoretical solutions of the singularity orders at the interface edge A are l1 0:491479 and l2 0:218652. Following the method proposed in Section 2, a group of the coecients A1 and A2 for the dierent distance r1 can be calculated by Eq. (9), as shown in Fig. 4. By the use of the least square method, the values of these
Fig. 1. Calculated model 1.
782
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
Fig. 2. Mesh division near the interface edge A.
coecients are determined as A1 1:398343 and A2 ÿ2:385579. Fig. 5 shows the double logarithmic distributions log sy ÿ A3ÿk txy =p0log r=W k 1,2. It can be found that each distribution behaves a good linearity, and the slopes which correspond to the orders of the stress singularities are l1 0:490515 and l2 0:222642. It is obvious that the above numerical results agree with the theoretical solutions very well. The relative errors are less than 0.2 and 1.8%, respectively. It is noted that the error of the weaker stress singularity is relatively large, but it still has good accuracy. According to Eq. (16), the extrapolating diagram of the stress intensity coecients is shown in Fig. 6. From the ®gure, it is found that the extrapolation method proposed in this paper has a good linearity. The dimensionless stress intensity coecients are extrapolated as K10 K1 = pW l1 0:360847 and K20 K2 = pW l2 0:167114. 3.2. A couple of conjugate complex singularities The analytical model is shown in Fig. 7. The material constants are E1 1000 GPa, n1 0:3; E2 1 GPa, n2 0:2. The boundary element mesh near the interface edge B is shown
Fig. 3. log sy =p and log txy =p vs. log r=W .
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
783
Fig. 4. Dispersity of coecients Ak .
in Fig. 8. Near the interface edge B, the nodes are also arranged in the equally proportional series (the proportional factor r 2) and the minimum element length is about 0:2 10ÿ7 W. The double logarithmic distributions of the dimensionless stress components sy and txy on the interface near the interface edge B are shown in Fig. 9. From the ®gure, it can be found that the double logarithmic distributions are not linear. According to Bogy [1], the theoretical solutions of the singularity orders at the interface edge B are a couple of conjugate complex, and its real part is lR 0:305859 and the imaginary part lI 0:103203. Following the method proposed in Section 2, we have b2 ÿ 4ac<0, so that the numerical results also show the oscillatory stress singularities. According to Eq. (10), a group of coecients AR and AI can be obtained, as shown in Fig. 10. By the least square method, we get AR 1:494505 and AI 0:877353. Figs. 11 and 12 show the distributions of log s=p0log r^ and cs log e0log r^ (L W), respectively. Their slopes which correspond to the real and imaginary parts of the oscillatory stress singularities are lR 0:306515 and lI 0:104132. It is found that these values
Fig. 5. log sy ÿ A3ÿk txy =p (k 1,2) vs. log r=W .
784
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
Fig. 6. Extrapolation of dimensionless stress intensity coecients Kk0 Kk = pW lk .
agree with the theoretical solutions very well, the relative errors are within 0.2 and 0.9%, respectively. According to Eq. (21), the extrapolating diagram to determine the stress intensity coecients can be drawn in Fig. 13. It can be also found that this extrapolation method has a very good linearity. The dimensionless stress intensity coecients are calculated as k10 k1 = pW lR 1:403879 and k20 k2 = pW lR ÿ0:227451.
Fig. 7. Calculated model 2.
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
785
Fig. 8. Mesh division near the interface edge B.
3.3. Test problem of the extrapolation method To investigate the accuracy and eciency of the present extrapolation method for the stress intensity coecients, the con®guration shown in Fig. 14 is analyzed by BEM too. According to Refs. [27,28], the singularity orders at point A are 1 ÿ l1 0:54484 and 1 ÿ l2 0:90853. The stress intensity coecient is K1 1:011682a1 1:011682 2580 2610 psi inl1 . By the algorithm described in Section 2, we obtain A1 0:544168, A2 ÿ0:157363, and l1 0:452873, l2 0:093041. The extrapolation diagram for K1 is shown in Fig. 15. From the ®gure, one obtains K1 2720 psi inl1 . The dierence between Carpenter's and this paper's results is only 4.2%. This fact means that the extrapolation proposed by this paper is ecient and accurate enough. By the way, it should be noted that the de®nition of KII in reference [28] is dierent from that of Eq. (14) or Eq. (15) in this paper, so that the comparison about the value of KII cannot be carried out.
Fig. 9. log sy =p and log txy =p vs. log r=W .
786
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
Fig. 10. Dispersity of coecients AR and AI .
4. Discussions 4.1. The case with more than two singularities For the problem with more than two singularities, the displacement and singular stress ®elds near the singular point can be expressed as ui
N X k1
Kk r1ÿlk gik y ,
sij
N X Kk rÿlk fijk y ,
23
k1
where N is the number of the singularity orders. Apparently, by the similar procedures in Section 2, a series of equations corresponding to Eq. (5) can be obtained as follows
Fig. 11. log s=p vs. log r=W .
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
Fig. 12. cs log e vs. log r=W .
Fig. 13. Extrapolation of complex dimensionless stress intensity coecients k0 k= pW lR k10 ik20 .
Fig. 14. Test problem model.
787
788
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
Fig. 15. Extrapolation of dimensionless stress intensity coecient K1.
s1 ÿ AN1ÿk1 s2 ÿ ÿ AN1ÿkNÿ1 sN Bk r^ÿlk ,
k 1,2, . . . ,N,
24
where sl (l 1,2, . . . ,N) are the selected stress components. If the number of the stress components is less than the number of the singularities, ui =r can be used as sl . The coecients Akm can be also determined by ordinary numerical results by selecting enough points distributed as the equally proportional series. The details of this deduction are omitted because it is only a simple algebraic computation. 4.2. General procedure In the engineering problem, the number of singularities is generally unknown. For this case, the following procedures are proposed: 1. Drawing the double logarithmic distributions of the stress components sij , if these distributions are all linear and parallel to each other, it means only single singularity exists and their slopes are the order of the stress singularity. 2. If the above distributions are not linear or not parallel to each other, the problem has multiple singularities. Then, it can be assumed that there exist 2,3, . . . ,N singularities step by step until all the double logarithmic distributions (according to Eq. (5) or (12) or (24)) behave linear. From the slopes, all the orders of singularity can be determined.
4.3. Note points It should be noted that in order to determine precisely the orders of the multiple stress singularities, a quite higher accuracy of the numerical results near the singular point is required. Moreover, it would be best if the nodes near the singular point could be arranged in the equally proportional series distribution. By the way, if a FEM code is used, the stress components on the interface of bonded dissimilar materials may be not accurate enough to
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
789
calculate the orders of the stress singularities for the sake of the stress discontinuity on the interface, so that the direction in which the analysis is carried out should be selected in the inside of the composed materials. The present methods are also suitable for anisotropic problems. For three-dimensional problems, if the points with an equally proportional series relation are selected in a radial direction y y0 and j j0 (in the global coordinates (r,y,j)), the above procedures are also applicable obviously.
5. Conclusions In this paper, the numerical methods are presented to determine the orders of the multiple or oscillatory stress singularities as well as the related stress intensity coecients from the common numerical solutions (stresses or displacements) obtained by an ordinary numerical code such as the FEM or BEM. Two models of the interface edge have been analyzed to verify the eectiveness. The results show that the multiple or oscillatory stress singularities can be determined precisely by the present method by using common numerical solutions. The extrapolation method to determine the stress intensity coecients presented in this paper has a very good linearity. The numerical results show that the methods are very convenient and ecient for the engineering problems. Moreover, these methods can be used to anisotropic problems and three-dimensional problems.
Acknowledgements The present work was supported by the National Natural Science Foundation of China under Grant No. 19502011.
References [1] Bogy DB. Two edge-bonded elastic wedges of dierent materials and wedge angles under surface tractions. Journal of Applied Mechanics 1971;38:377±86. [2] Bogy DB, Wang KC. Stress singularities at interface corners in bonded dissimilar isotropic elastic materials. International Journal of Solids and Structures 1971;7:993±1005. [3] Williams ML. The stresses around a fault or crack in dissimilar media. Bulletin of the Seismological Society of America 1959;49:199±204. [4] Bogy DB. On the plane elastostatic problem of a loaded crack terminating at a material interface. Journal of Applied Mechanics 1971;38:911±8. [5] Chen DH, Nisitani H. Singular stress ®elds near the tip of a V-notch in a semi-in®nite plate. Transactions of the Japan Society of Mechanical Engineers (in Japanese) 1991;A57:1406±11. [6] Yeh JR, Tadjbakhsh IG. Stress singularity in composite laminates by ®nite element method. Journal of Composite Materials 1986;20:347±64. [7] Somaratna N, Ting TCT. Three-dimensional stress singularities in anisotropic materials and composites. International Journal of Engineering Science 1986;24:1115±34.
790
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
[8] Gu L, Belytschko T. A numerical study of stress singularities in a two-material wedge. International Journal of Solids and Structures 1994;31:865±89. [9] Pageau SS, Biggers Jr. SB. A ®nite element approach to three-dimensional singular stress states in anisotropic multi-material wedges and junctions. International Journal of Solids and Structures 1996;33:33±47. [10] Yuuki R, Xu JQ. Stress intensity factors for the interface crack between dissimilar orthotropic materials. Transactions of the Japan Society of Mechanical Engineers (in Japanese) 1991;A57:1542±9. [11] Tan CL, Gao YL, Afagh FF. Boundary element analysis of interface cracks between dissimilar anisotropic materials. International Journal of Solids and Structures 1992;29:3201±20. [12] Leslie BS. Weight functions for interface cracks. International Journal of Fracture 1993;60:89±95. [13] Wu YL. A new method for evaluation of stress intensities for interface cracks. Engineering Fracture Mechanics 1994;48:755±61. [14] Charalambides PG, Zhang WB. An energy method for calculating the stress intensities in orthotropic bimaterial fracture. International Journal of Fracture 1996;76:97±120. [15] Yang XX, Kuang ZB. Contour integral method for stress intensity factors of interface crack. International Journal of Fracture 1996;78:299±313. [16] Dong YX, Wang ZM, Wang B. On the computation of stress intensity factors for interfacial cracks using quarter-point boundary elements. Engineering Fracture Mechanics 1997;57:335±42. [17] Lu H, Chiang FP. Photoelastic determination of stress intensity factor of an interface crack in a bi-material. Journal of Applied Mechanics 1993;60:93±100. [18] Ricci V, Shukla A, Singh RP. Evaluation of fracture mechanics parameters in bimaterial systems using strain gages. Engineering Fracture Mechanics 1997;58:273±83. [19] Munz D, Yang YY. Stresses near the edge of bonded dissimilar materials described by two stress intensity factors. International Journal of Fracture 1993;60:169±77. [20] Leslie BS, Yang YY, Munz D. An in¯uence function for stress intensity factors of bimaterial notched bodies. International Journal of Fracture 1997;85:333±50. [21] Zhao Z, Hahn HG. Determining the SIF of a V-notch from the results of a mixed-mode crack. Engineering Fracture Mechanics 1992;43:511±8. [22] Chen DH. Stress intensity factors for V-notched strip under tension or in-plane bending. International Journal of Fracture 1995;70:81±97. [23] Chan SK, Tuba IS, Wilson WK. On the ®nite element method in linear fracture mechanics. Engng Fracture Mech 1970;2:1±17. [24] William C Carpenter. Extrapolation techniques for determining stress intensity factors. Engng Fracture Mech 1983;2:325±32. [25] Rice JR. Elastic fracture mechanics concepts for interfacial cracks. Journal of Applied Mechanics 1988;55:98± 103. [26] Yuuki R, Xu JQ. Development of the boundary element method program (BEM2D/EWS) and analyses of dissimilar materials and interface cracks. Seisan-Kenkyu (in Japanese) 1991;43:306±9. [27] William C Carpenter. A collocation procedure for determining fracture mechanics parameters at a corner. International Journal of Fracture 1984;24:255±66. [28] William C Carpenter. Mode I and mode II stress intensities for plates with cracks of ®nite opening. International Journal of Fracture 1984;26:201±14.
www.elsevier.com/locate/engfracmech
Numerical methods for the determination of multiple stress singularities and related stress intensity coecients Jin-Quan Xu a,*, Yi-Hua Liu b, Xiao-Gui Wang a a
b
Department of Mechanics, Zhejiang University, Hangzhou 310027, People's Republic of China Department of Applied Mathematics and Mechanics, Hefei University of Technology, Hefei 230009, People's Republic of China Received 3 December 1998; accepted 3 March 1999
Abstract This paper proposed numerical methods to determine the multiple stress singularities (including the oscillatory stress singularities) and the related stress intensity coecients, by the use of common numerical solutions (stresses or displacements) obtained by an ordinary numerical tool such as ®nite element method (FEM) or boundary element method (BEM). To verify the eciency of the present methods, two models of bonded dissimilar materials under the plane strain state are analyzed by BEM, and the orders of the stress singularities and the related stress intensity coecients are examined numerically. The results show that all the orders of the stress singularities at an interface edge can be determined precisely by the present method, and the related stress intensity coecients can be determined by the extrapolation method with a very good linearity. It is found that the methods presented in this paper are very simple and ecient. Moreover, they can be easily extended to any singular problem. # 1999 Elsevier Science Ltd. All rights reserved. Keywords: Multiple stress singularities; Stress intensity coecients; Common numerical solutions; Extrapolation method; Bonded dissimilar materials
1. Introduction It is well known that there may appear multiple singularities or oscillatory singularities at a stress singular point such as an interface edge, an interface corner, the tip of an interface crack * Corresponding author. Tel.: +86-571-7951769; fax: +86-571-7951464. E-mail address: [email protected] (J.-Q. Xu) 0013-7944/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 3 - 7 9 4 4 ( 9 9 ) 0 0 0 4 4 - 2
776
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
(or a crack terminating at the interface) in bonded dissimilar materials, as well as the tip of a V-notch in a homogenous material (e.g., Refs. [1±5]). When the stress singularity is not single (i.e., it is multiple), the orders of the stress singularities can be only determined generally by analytic methods. However, there are many cases to which the theoretical results are not available, or they are too dicult to analyze theoretically, so that numerical methods to determine the multiple stress singularities and the related stress intensity coecients are strongly expected in applied mechanics and engineering. Recently, some numerical methods [6± 9] have been proposed. For example, Somaratna and Ting [7] proposed a particularly constructed ®nite element method in which the orders of the multiple or oscillatory stress singularities can be determined by the eigenequation in the form of homogeneous linear algebraic equations. On the other hand, various numerical methods [10±16] as well as photoelastic [17] and strain gage [18] methods have been reported to determine the stress intensity factors of an interface crack. However, for the interface edge, it is not easy to determine stress intensity coecients related to the multiple singularities. Munz and Yang [19] presented a numerical method to calculate stress intensity coecients of the interface edge based on the known analytic angular functions. Leslie et al. [20] gave out an approach based on the in¯uence function. For a V-notch problem, Zhao and Hahn [21] proposed a method to determine the stress intensity factors from the results of a mixed mode crack. Chen [22] proposed a method based on the body force method. Some additional related investigations are also available in the literature. In general, the extrapolation procedure which was developed by Chan et al.[23] and extended later by Carpenter [24] could be considered as a useful tool to determine the stress intensity coecients for multiple singularities problems if some modi®cations are introduced. The purpose of this study is to set up the numerical methods by which the orders of the multiple stress singularities (or oscillatory singularities) and the related stress intensity coecients can be determined from the common numerical solutions (stresses or displacements) obtained by ordinary numerical tools such as ®nite element method (FEM) and boundary element method (BEM). The method would be applicable even if there is no information about the theoretical solutions. Two bi-material models under the plane strain state are analyzed by BEM to illustrate the application of the proposed methods. The numerical results show that all the orders of the multiple stress singularities (including oscillatory singularities) can be determined precisely by the present numerical method, and the related stress intensity coecients can also be determined by the extrapolation method proposed in this paper with a very good linearity. Moreover, a test problem is also analyzed to investigate the accuracy of the present extrapolation method. The results show that the present algorithm is ecient and accurate enough too.
2. Numerical analytical methods 2.1. Determination of multiple stress singularities For simplicity, the problem with two stress singularities is considered here. Problems with more than two singularities can be considered in the same way. By the way, the order of the
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
777
single stress singularity can be determined directly by the log±log (double logarithmic) diagram of the stress and distance from the singular point. In the plane problem, the singular stress ®eld near the singular point with two stress singularities can be described by (in the polar coordinates (r,y)) sij K1 rÿl1 fij1 y K2 rÿl2 fij2 y
1
where lk 0
2
where Hk Kk Lÿlk (k 1,2). In order to determine l1 and l2 , we select arbitrarily two stress components s1 and s2 , then in the selected direction y y0 , we have sl H1l r^ÿl1 H2l r^ÿl2 ,
l 1,2,
3
where Hkl Hk flk y0 ,
k,l 1,2:
4
After some algebraic computations of Eq. (3) we have s1 ÿ A3ÿk s2 Bk r^ÿlk , k 1,2,
5
where A1 H11 =H12 ,
A2 H21 =H22 , B1 H12 A1 ÿ A2 , B2 H22 A2 ÿ A1 ,
6
are coecients to be determined. From Eq. (5), it can be found that the double logarithmic distributions log s1 ÿ A3ÿk s2 0log r^ (k 1,2) should be linear in the dominant region of the singular stress ®eld (1), and the slopes are related to the orders of the stress singularities l1 and l2 , respectively. However, to draw the double logarithmic diagram, coecients A1 and A2 should be predetermined. In order to determine A1 and A2 numerically, we select three points distributed as an equally proportional series in the selected direction y y0 , i.e., rn rnÿ1 r1 ; n 1, 2, 3 (r is a proportional factor and r1 is the distance between the ®rst point and the singular point). The (l 1,2; n 1,2,3), which can be corresponding stress components are expressed by s n l obtained by the ordinary numerical procedure such as FEM or BEM. Substituting these stress k yields components into Eq. (5) and eliminating Bk r^ÿl n aA2k bAk c 0, k 1,2, where
7
778
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
2 a s 22 ÿs 21 s 23 , b s 11 s 23 s 13 s 21 ÿ 2s 12 s 22 , c
s 12
2
ÿs 11 s 13 ,
8
are coecients which can be determined by the numerical results. The roots of Eq. (7) may have the following three cases: 1. b2 ÿ 4ac 0 In this case A1 A2 ÿb= 2a. This leads to l1 l2 , so that it is the case with the single stress singularity. 2. b2 ÿ 4ac > 0 In this case, A1 and A2 are all real, i.e., k p b ÿ ÿ 1 b2 ÿ 4ac , k 1,2: 9 Ak ÿ 2a This leads to two dierent real singularities. 3. b2 ÿ 4ac<0 In this case A1 and A2 are a couple of conjugate complexes, i.e., p b ÿ i 4ac ÿ b2 A1 A 2 AR iAI ÿ , 10 2a p where i ÿ1. This leads to a couple of the complex conjugate orders of the stress singularities which mean oscillatory stress singularities. It should be noted that in this case lk , Bk , Kk , Hk , Hkl , flk y (k,l 1,2) are also complex conjugate, i.e., G1 G 2 GR iGI ,
11
in which Gk lk ,Bk ,Kk ,Hk ,Hkl ,flk y. Substituting Eq. (11) into Eq. (5) yields seics BeicB r^ÿ lR ilI ,
12
where
q s s1 ÿ AR s2 2 AI s2 2 ,
q B B 2R B 2I ,
cs tanÿ1 AI s2 = s1 ÿ AR s2 , cB tanÿ1 BI =BR :
13
From Eq. (12), it can be found that the double logarithmic distribution log s0log r^ and the single logarithmic distribution cs log e0log r^ should be linear in the dominant region of the
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
779
oscillatory singular stress ®eld, and their slopes are related to the real part lR and the imaginary part lI of the order l1 of the oscillatory stress singularity, respectively. To obtain accurate values of coecients A1 and A2 , a few sets with three points distributed as the equally proportional series should be considered. Due to the computational error, A1 and A2 obtained by each set may be slightly dierent. However, because A1 and A2 should be constants in the selected direction, they can be determined by the least square method from these slight dierent values. 2.2. Determination of stress intensity coecients 2.2.1. The case with two real singularities Assume that the stress intensity coecients K1 and K2 can be de®ned by the stresses sy and try in the direction y y0 (e.g., the interface), i.e., sy0 sy jyy0 K1 rÿl1 K2 rÿl2 fy2 y0 , try0 try jyy0 K1 rÿl1 fry1 y0 K2 rÿl2 :
14
This de®nition of the stress intensity coecients means the angular functions fy1 y and fry2 y are normalized in the direction y y0 , i.e., fy1 y0 fry2 y0 1. Taking s1 sy0 , s2 try0 , and according to Eq. (4), we have H11 H1 , H22 H2 . Therefore, from Eqs. (5) and (6), the stress intensity coecients can be determined by the following extrapolation method ÿ K1 Alimrl1 sy0 ÿ A2 try0 , r40
ÿ K2 Alimrl2 try0 ÿ Aÿ1 1 sy0 , r40
15
where A A1 = A1 ÿ A2 , the coecients A1 and A2 have been given by Eq. (9). The linearity of this extrapolation method would be investigated numerically in Section 3. 2.2.2. The case with a couple of complex conjugate singularities In this case, substituting Eq. (11) into (2) yields sy rÿlR k1 Fy1 y cos lI ln r^ k2 Fy2 y sin lI ln r^ , try rÿlR k1 Fry1 y sin lI ln r^ k2 Fry2 y cos lI ln r^ ,
16
where the stress intensity coecients k1 and k2 as well as the angular functions Fy1 y, Fy2 y, Fry1 y, and Fry2 y are k1 HR LlR , k2 HI LlR ,
17
780
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
Fy1 y 2 fyR y ÿ HfyI y ,
Fy2 y 2 fyR y H ÿ1 fyI y ,
Fry1 y 2 HfryR y fryI l ,
Fry2 y 2 H ÿ1 fryR y ÿ fryI y ,
18
in which H HI =HR . We also assume that k1 and k2 can be de®ned by the stresses sy and try in the direction y y0 , i.e., sy0 sy jyy0 rÿlR k1 cos lI ln r^ k2 Fy2 y0 sin lI ln r^ , try0 try jyy0 rÿlR k1 Fry1 y0 sin lI ln r^ k2 cos lI ln r^ :
19
This de®nition means the angular functions Fy1 y and Fry2 y are normalized in the direction y y0 , i.e., Fy1 y0 Fry2 y0 1. Taking s1 sy0 and s2 try0 , then from Eqs. (4), (6) and (11) we get BR A HR AI ÿ 2H1I AR ÿ2AI H2I , BI A HR AR 2H1I AI AI HI ,
20
where A AI = A2R A2I . After the substitution of Eqs. (11) and (20) into (5) and some algebraic computations, we can ®nally obtain the following extrapolation method: n lR ÿ k1 Aÿ1 sy0 ÿ AR try0 AR sin lI ln r^ AI cos lI ln r^ AI try0 AR cos lI lim r I r40
ln r^ ÿ AI sin lI ln r^
o ,
lR ÿ k2 Aÿ1 sy0 ÿ AR try0 sin lI ln r^ AI try0 cos lI ln r^ , I lim r r40
21
where the coecients AR and AI have been given by Eq. (10). The linearity of this extrapolation method would be discussed numerically in Section 3 too. For the interface crack, because Fry1 y0 ÿFy2 y0 1 (assume that the interface is in the direction y y0 ) and lR 1=2, from Eqs. (4), (6) and (11), we have AR 0 and AI ÿ1, so that Eq. (21) can be simpli®ed as p k1 lim r sy0 cos lI ln r^ try0 sin lI ln r^ , r40
p k2 lim r try0 cos lI ln r^ ÿ sy0 sin lI ln r^ : r40
22
It is obvious that the dierence between the de®nition of Eq. (22) and Rice's [25] de®nition p about the stress intensity factors of an interface crack is only a proportional factor 2p.
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
781
3. Numerical results To verify the eciency of the above numerical methods, the orders of the stress singularities and the related stress intensity coecients of two models of the interface edge under plane strain state are calculated numerically by using the solutions of the BEM. The program BEM2D [26] used here, is especially ecient for the analysis of bonded dissimilar materials, which was developed by the author. A test problem is also analyzed to check the accuracy of the present extrapolation method for stress intensity coecients. 3.1. Two real singularities The analytical model is shown in Fig. 1. The material constants are E1 108 GPa, n1 0:33; E2 304 GPa, n2 0:27. The boundary element mesh near the interface edge A is shown in Fig. 2. The nodes near the interface edge A are arranged in the equally proportional series distribution (the proportional factor r 2) and the minimum element length is about 2:0 10ÿ7 W. The double logarithmic distributions of the dimensionless stress components sy and txy on the interface near the interface edge A are shown in Fig. 3. From Fig. 3 it can be found that the double logarithmic distributions are not parallel to each other. If one takes the slopes as the stress singular order, dierent values can be obtained by adopting dierent stress components. Moreover, though the distribution seems to be a line, the slope determined by dierent regions shows that it is a curve in fact. Therefore, it is the case of multiple singularities. According to Bogy [1], the theoretical solutions of the singularity orders at the interface edge A are l1 0:491479 and l2 0:218652. Following the method proposed in Section 2, a group of the coecients A1 and A2 for the dierent distance r1 can be calculated by Eq. (9), as shown in Fig. 4. By the use of the least square method, the values of these
Fig. 1. Calculated model 1.
782
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
Fig. 2. Mesh division near the interface edge A.
coecients are determined as A1 1:398343 and A2 ÿ2:385579. Fig. 5 shows the double logarithmic distributions log sy ÿ A3ÿk txy =p0log r=W k 1,2. It can be found that each distribution behaves a good linearity, and the slopes which correspond to the orders of the stress singularities are l1 0:490515 and l2 0:222642. It is obvious that the above numerical results agree with the theoretical solutions very well. The relative errors are less than 0.2 and 1.8%, respectively. It is noted that the error of the weaker stress singularity is relatively large, but it still has good accuracy. According to Eq. (16), the extrapolating diagram of the stress intensity coecients is shown in Fig. 6. From the ®gure, it is found that the extrapolation method proposed in this paper has a good linearity. The dimensionless stress intensity coecients are extrapolated as K10 K1 = pW l1 0:360847 and K20 K2 = pW l2 0:167114. 3.2. A couple of conjugate complex singularities The analytical model is shown in Fig. 7. The material constants are E1 1000 GPa, n1 0:3; E2 1 GPa, n2 0:2. The boundary element mesh near the interface edge B is shown
Fig. 3. log sy =p and log txy =p vs. log r=W .
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
783
Fig. 4. Dispersity of coecients Ak .
in Fig. 8. Near the interface edge B, the nodes are also arranged in the equally proportional series (the proportional factor r 2) and the minimum element length is about 0:2 10ÿ7 W. The double logarithmic distributions of the dimensionless stress components sy and txy on the interface near the interface edge B are shown in Fig. 9. From the ®gure, it can be found that the double logarithmic distributions are not linear. According to Bogy [1], the theoretical solutions of the singularity orders at the interface edge B are a couple of conjugate complex, and its real part is lR 0:305859 and the imaginary part lI 0:103203. Following the method proposed in Section 2, we have b2 ÿ 4ac<0, so that the numerical results also show the oscillatory stress singularities. According to Eq. (10), a group of coecients AR and AI can be obtained, as shown in Fig. 10. By the least square method, we get AR 1:494505 and AI 0:877353. Figs. 11 and 12 show the distributions of log s=p0log r^ and cs log e0log r^ (L W), respectively. Their slopes which correspond to the real and imaginary parts of the oscillatory stress singularities are lR 0:306515 and lI 0:104132. It is found that these values
Fig. 5. log sy ÿ A3ÿk txy =p (k 1,2) vs. log r=W .
784
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
Fig. 6. Extrapolation of dimensionless stress intensity coecients Kk0 Kk = pW lk .
agree with the theoretical solutions very well, the relative errors are within 0.2 and 0.9%, respectively. According to Eq. (21), the extrapolating diagram to determine the stress intensity coecients can be drawn in Fig. 13. It can be also found that this extrapolation method has a very good linearity. The dimensionless stress intensity coecients are calculated as k10 k1 = pW lR 1:403879 and k20 k2 = pW lR ÿ0:227451.
Fig. 7. Calculated model 2.
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
785
Fig. 8. Mesh division near the interface edge B.
3.3. Test problem of the extrapolation method To investigate the accuracy and eciency of the present extrapolation method for the stress intensity coecients, the con®guration shown in Fig. 14 is analyzed by BEM too. According to Refs. [27,28], the singularity orders at point A are 1 ÿ l1 0:54484 and 1 ÿ l2 0:90853. The stress intensity coecient is K1 1:011682a1 1:011682 2580 2610 psi inl1 . By the algorithm described in Section 2, we obtain A1 0:544168, A2 ÿ0:157363, and l1 0:452873, l2 0:093041. The extrapolation diagram for K1 is shown in Fig. 15. From the ®gure, one obtains K1 2720 psi inl1 . The dierence between Carpenter's and this paper's results is only 4.2%. This fact means that the extrapolation proposed by this paper is ecient and accurate enough. By the way, it should be noted that the de®nition of KII in reference [28] is dierent from that of Eq. (14) or Eq. (15) in this paper, so that the comparison about the value of KII cannot be carried out.
Fig. 9. log sy =p and log txy =p vs. log r=W .
786
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
Fig. 10. Dispersity of coecients AR and AI .
4. Discussions 4.1. The case with more than two singularities For the problem with more than two singularities, the displacement and singular stress ®elds near the singular point can be expressed as ui
N X k1
Kk r1ÿlk gik y ,
sij
N X Kk rÿlk fijk y ,
23
k1
where N is the number of the singularity orders. Apparently, by the similar procedures in Section 2, a series of equations corresponding to Eq. (5) can be obtained as follows
Fig. 11. log s=p vs. log r=W .
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
Fig. 12. cs log e vs. log r=W .
Fig. 13. Extrapolation of complex dimensionless stress intensity coecients k0 k= pW lR k10 ik20 .
Fig. 14. Test problem model.
787
788
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
Fig. 15. Extrapolation of dimensionless stress intensity coecient K1.
s1 ÿ AN1ÿk1 s2 ÿ ÿ AN1ÿkNÿ1 sN Bk r^ÿlk ,
k 1,2, . . . ,N,
24
where sl (l 1,2, . . . ,N) are the selected stress components. If the number of the stress components is less than the number of the singularities, ui =r can be used as sl . The coecients Akm can be also determined by ordinary numerical results by selecting enough points distributed as the equally proportional series. The details of this deduction are omitted because it is only a simple algebraic computation. 4.2. General procedure In the engineering problem, the number of singularities is generally unknown. For this case, the following procedures are proposed: 1. Drawing the double logarithmic distributions of the stress components sij , if these distributions are all linear and parallel to each other, it means only single singularity exists and their slopes are the order of the stress singularity. 2. If the above distributions are not linear or not parallel to each other, the problem has multiple singularities. Then, it can be assumed that there exist 2,3, . . . ,N singularities step by step until all the double logarithmic distributions (according to Eq. (5) or (12) or (24)) behave linear. From the slopes, all the orders of singularity can be determined.
4.3. Note points It should be noted that in order to determine precisely the orders of the multiple stress singularities, a quite higher accuracy of the numerical results near the singular point is required. Moreover, it would be best if the nodes near the singular point could be arranged in the equally proportional series distribution. By the way, if a FEM code is used, the stress components on the interface of bonded dissimilar materials may be not accurate enough to
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
789
calculate the orders of the stress singularities for the sake of the stress discontinuity on the interface, so that the direction in which the analysis is carried out should be selected in the inside of the composed materials. The present methods are also suitable for anisotropic problems. For three-dimensional problems, if the points with an equally proportional series relation are selected in a radial direction y y0 and j j0 (in the global coordinates (r,y,j)), the above procedures are also applicable obviously.
5. Conclusions In this paper, the numerical methods are presented to determine the orders of the multiple or oscillatory stress singularities as well as the related stress intensity coecients from the common numerical solutions (stresses or displacements) obtained by an ordinary numerical code such as the FEM or BEM. Two models of the interface edge have been analyzed to verify the eectiveness. The results show that the multiple or oscillatory stress singularities can be determined precisely by the present method by using common numerical solutions. The extrapolation method to determine the stress intensity coecients presented in this paper has a very good linearity. The numerical results show that the methods are very convenient and ecient for the engineering problems. Moreover, these methods can be used to anisotropic problems and three-dimensional problems.
Acknowledgements The present work was supported by the National Natural Science Foundation of China under Grant No. 19502011.
References [1] Bogy DB. Two edge-bonded elastic wedges of dierent materials and wedge angles under surface tractions. Journal of Applied Mechanics 1971;38:377±86. [2] Bogy DB, Wang KC. Stress singularities at interface corners in bonded dissimilar isotropic elastic materials. International Journal of Solids and Structures 1971;7:993±1005. [3] Williams ML. The stresses around a fault or crack in dissimilar media. Bulletin of the Seismological Society of America 1959;49:199±204. [4] Bogy DB. On the plane elastostatic problem of a loaded crack terminating at a material interface. Journal of Applied Mechanics 1971;38:911±8. [5] Chen DH, Nisitani H. Singular stress ®elds near the tip of a V-notch in a semi-in®nite plate. Transactions of the Japan Society of Mechanical Engineers (in Japanese) 1991;A57:1406±11. [6] Yeh JR, Tadjbakhsh IG. Stress singularity in composite laminates by ®nite element method. Journal of Composite Materials 1986;20:347±64. [7] Somaratna N, Ting TCT. Three-dimensional stress singularities in anisotropic materials and composites. International Journal of Engineering Science 1986;24:1115±34.
790
J.-Q. Xu et al. / Engineering Fracture Mechanics 63 (1999) 775±790
[8] Gu L, Belytschko T. A numerical study of stress singularities in a two-material wedge. International Journal of Solids and Structures 1994;31:865±89. [9] Pageau SS, Biggers Jr. SB. A ®nite element approach to three-dimensional singular stress states in anisotropic multi-material wedges and junctions. International Journal of Solids and Structures 1996;33:33±47. [10] Yuuki R, Xu JQ. Stress intensity factors for the interface crack between dissimilar orthotropic materials. Transactions of the Japan Society of Mechanical Engineers (in Japanese) 1991;A57:1542±9. [11] Tan CL, Gao YL, Afagh FF. Boundary element analysis of interface cracks between dissimilar anisotropic materials. International Journal of Solids and Structures 1992;29:3201±20. [12] Leslie BS. Weight functions for interface cracks. International Journal of Fracture 1993;60:89±95. [13] Wu YL. A new method for evaluation of stress intensities for interface cracks. Engineering Fracture Mechanics 1994;48:755±61. [14] Charalambides PG, Zhang WB. An energy method for calculating the stress intensities in orthotropic bimaterial fracture. International Journal of Fracture 1996;76:97±120. [15] Yang XX, Kuang ZB. Contour integral method for stress intensity factors of interface crack. International Journal of Fracture 1996;78:299±313. [16] Dong YX, Wang ZM, Wang B. On the computation of stress intensity factors for interfacial cracks using quarter-point boundary elements. Engineering Fracture Mechanics 1997;57:335±42. [17] Lu H, Chiang FP. Photoelastic determination of stress intensity factor of an interface crack in a bi-material. Journal of Applied Mechanics 1993;60:93±100. [18] Ricci V, Shukla A, Singh RP. Evaluation of fracture mechanics parameters in bimaterial systems using strain gages. Engineering Fracture Mechanics 1997;58:273±83. [19] Munz D, Yang YY. Stresses near the edge of bonded dissimilar materials described by two stress intensity factors. International Journal of Fracture 1993;60:169±77. [20] Leslie BS, Yang YY, Munz D. An in¯uence function for stress intensity factors of bimaterial notched bodies. International Journal of Fracture 1997;85:333±50. [21] Zhao Z, Hahn HG. Determining the SIF of a V-notch from the results of a mixed-mode crack. Engineering Fracture Mechanics 1992;43:511±8. [22] Chen DH. Stress intensity factors for V-notched strip under tension or in-plane bending. International Journal of Fracture 1995;70:81±97. [23] Chan SK, Tuba IS, Wilson WK. On the ®nite element method in linear fracture mechanics. Engng Fracture Mech 1970;2:1±17. [24] William C Carpenter. Extrapolation techniques for determining stress intensity factors. Engng Fracture Mech 1983;2:325±32. [25] Rice JR. Elastic fracture mechanics concepts for interfacial cracks. Journal of Applied Mechanics 1988;55:98± 103. [26] Yuuki R, Xu JQ. Development of the boundary element method program (BEM2D/EWS) and analyses of dissimilar materials and interface cracks. Seisan-Kenkyu (in Japanese) 1991;43:306±9. [27] William C Carpenter. A collocation procedure for determining fracture mechanics parameters at a corner. International Journal of Fracture 1984;24:255±66. [28] William C Carpenter. Mode I and mode II stress intensities for plates with cracks of ®nite opening. International Journal of Fracture 1984;26:201±14.