- Email: [email protected]
Dislocation-based semi-analytical method for calculating stress intensity factors of cracks: Two-dimensional cases
Engineering Fracture Mechanics 77 (2010) 3521–3531
Contents lists available at ScienceDirect
Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech
Dislocation-based semi-analytical method for calculating stress intensity factors of cracks: Two-dimensional cases Xi-Qiao Feng a,*, Yun-Fei Shi a, Xu-Yue Wang b, Bo Li a, Shou-Wen Yu a, Qiang Yang c a
AML, Department of Engineering Mechanics, Tsinghua University, Beijing 100084, China Department of Sciences, Harbin Institute of Technology, Shenzhen Graduate School, Shenzhen 518055, China c State Key Laboratory of Hydroscience and Engineering, Tsinghua University, Beijing 100084, China b
a r t i c l e
i n f o
Article history: Received 2 October 2009 Received in revised form 26 February 2010 Accepted 4 March 2010 Available online 11 March 2010 Keywords: Crack Stress intensity factor Dislocation Finite element method Singular integral equation Semi-analytical method
a b s t r a c t Determination of the stress intensity factors of cracks is a fundamental issue for assessing the performance safety and predicting the service lifetime of engineering structures. In the present paper, a dislocation-based semi-analytical method is presented by integrating the continuous dislocation model with the finite element method together. Using the superposition principle, a two-dimensional crack problem in a finite elastic body is reduced to the solution of a set of coupled singular integral equations and the calculation of the stress fields of a body which has the same shape as the original one but has no crack. It can easily solve crack problems of structures with arbitrary shape, and the calculated stress intensity factors show almost no dependence upon the finite element mesh. Some representative examples are given to illustrate the efficacy and accuracy of this novel numerical method. Only two-dimensional cases are addressed here, but this method can be extended to threedimensional problems. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction Despite the significant development of elastoplastic fracture mechanics methods in terms of such parameters as J-integral and crack opening displacement, the stress intensity factor (SIF) is still the most extensively utilized parameter for assessing the integrity and predicting the service life of engineering components and structures with cracks. Massive efforts have been directed towards developing analytical and numerical methods to calculate the SIFs under different conditions [1–9]. Analytical methods are available only for some simple geometry configurations and boundary conditions [1–3]. Therefore, various numerical approaches (e.g., finite element method [10–14], boundary element method [14–20], meshless method [21– 24], weight function method [25,26], and discrete element method [27–30]) have been established to calculate the SIFs of cracks. Among them, the finite element method has been used with great success for crack problems in a wide diversity of industrial fields. To determine the SIFs or other near-tip fields of a crack, very fine mesh should be generated in the vicinity of the crack front, making the numerical calculation time-consuming. An inappropriate choice of element types or sizes will lead to calculation results with considerable error and/or mesh dependence [10,11]. The boundary element method has been proved as another powerful technique to solve elastic solids with cracks. Different from the finite element method, the boundary element method discretizes only the body’s boundary [14–20]. Nevertheless, the reduced number of unknowns does not necessarily lead to improved efficiency because the boundary element method generally produces a full matrix, while the stiffness matrices in finite element simulations are usually sparse
* Corresponding author. Tel.: +86 10 62772934; fax: +86 10 62781824. E-mail address: [email protected] (X.-Q. Feng). 0013-7944/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2010.03.004
3522
X.-Q. Feng et al. / Engineering Fracture Mechanics 77 (2010) 3521–3531
[18]. The element-free or meshless method rose in 1990s and has been successfully applied in some fracture problems [21]. For example, the element-free Galerkin (EFG) method [22,23] is well suited to the modeling of two- and three-dimensional crack problems because arbitrary surfaces across which the displacement function is discontinuous can readily be incorporated in this technique. However, the EFG method is difficult in enforcing essential boundary conditions [24]. In addition, some semi-analytical methods have also been proposed for solving crack problems [31–35]. In the present paper, a semi-analytical method, also called dislocation-based finite element method (DFEM), is presented to solve the SIFs of cracks in a finite body of arbitrary shape. It combines the continuous dislocation model with finite element method. In Section 2, the two-dimensional cracked plate problem is formulated to show the basic idea of this method. In Section 3, several representative examples are given to illustrate its efficiency and accuracy.
2. Dislocation-based finite element method 2.1. Decomposition of a crack problem We consider a finite plate X containing a straight and thickness-through crack and subjected to the distributed tractions
rcrack and rexternal on the crack surface and its outer boundary C, respectively, as shown in Fig. 1a. The plate may have an arbitrary shape, and the crack may either intersect with the boundary C (e.g., a boundary or edge crack) or not (e.g., an interior crack). Assume that the material is isotropic and linear elastic. According to the superposition principle, the crack problem in Fig. 1a can be decomposed into two additive subproblems, A and B, as shown in Figs. 1b and c, respectively. In the former, the plate has no crack but is subjected to the distributed force rexternal on its boundary, whereas in the latter, the plate of the same shape is traction-free on its boundary but subjected to the traction rB on the crack surfaces. The superposition principle and the traction condition on the crack faces in Fig. 1a require that rB = rcrack rA, where rA denotes the stress distribution along the line at the original crack position in Fig. 1b. rA can be easily calculated either by an analytical or a numerical method (e.g., finite element method) because the subproblem A has no crack. Thus, the SIFs of the crack in the considered configuration in Fig. 1a are identical to those in the subproblem B. However, it is noticed that the subproblem B is still hard to solve because of the arbitrary shape and finite size of the plate. To manage these complexities, the superposition principle is utilized again to further decompose the subproblem B into C and D, neither of which contains a crack, as shown in Fig. 2. In the subproblem C, the plate is infinite in size and contains, instead of the crack, two families of continuously distributed dislocations, with Burgers vectors normal and parallel to the original crack faces, respectively. The dislocation density distribution functions along the crack length direction are yet to be determined. Let rdislocation and rC denote the dislocation-induced stresses along the original plate boundary C and along the original crack faces, respectively. In the subproblem D, the plate has the same shape and size as that in Fig. 1a and a distributed traction rimage is applied on its boundary in order to fulfill the boundary conditions. The superposition principle and the traction-free boundary condition on C in Fig. 1c require that the image force on the boundary C induced by the distributed dislocations should satisfy rimage = rdislocation. Thus, the problem in Fig. 1a becomes the solution of the
(a)
σ external
σ crack
(b)
(c)
σ external
σA
σ B = σ crack − σ A
Fig. 1. Decomposition of the crack problem: (a) the original crack configuration, (b) subproblem A, and (c) subproblem B.
3523
X.-Q. Feng et al. / Engineering Fracture Mechanics 77 (2010) 3521–3531
(a)
(b) C
D
σC
σD
σ dislocation σimage = − σdislocation Fig. 2. Decomposition of subproblem B into (a) subproblem C and (b) subproblem D.
dislocation density functions in the infinite body problem C coupled with the defect-free plate problem D, both of which are easier to be solved than the original one, as will be described below. Provided that the stress boundary condition on the crack surfaces in Fig. 1c can be satisfied, i.e., rC + rD = rB, then the stress field in B is exactly equal to the addition of those in C and D and the SIFs in B can be uniquely determined from the stress field in C. The dislocation density functions in the subproblem C will be solved approximately from the condition rC + rD = rB, as will be formulated in what follows. 2.2. Continuous dislocation model Refer to the Cartesian coordinate system (O–xy), as shown in Fig. 1a, where the origin O is located at the center of the crack and the x axis is along the crack direction. Without loss of any generality, the half-length of the crack is assumed to be unity, or in other words, all the length parameters are normalized by the crack half-length c. According to the continuous ; v gT can be equivalently replaced by two families of condislocation model, the crack opening displacement vector u ¼ fu tinuously distributed dislocations with the Burgers vector b = {bx, by}T, where the dislocations with the component bx and by can be considered as gliding and climbing edge dislocations, respectively [4,5,36]. The corresponding dislocation density functions are defined as [4]
@ðuþ u Þ @u ¼ ; @x @x @ v @ðv þ v Þ ; ¼ gðxÞ ¼ @x @x
f ðxÞ ¼
ð1Þ
where u+(t), u(t), v+(t), and v(t) denote the displacements of the two crack surfaces at position t along the x and y directions, respectively. Then the families of gliding and climbing dislocations with an infinitesimal segment dt have the intensities of f(t)dt and g(t)dt, respectively. The dislocation-induced stresses at a position (x, y = 0) on the crack faces can be expressed as
ryx ðxÞ ¼ ryy ðxÞ ¼
Z
1
1 Z 1 1
rðIÞ yx ðt; xÞgðtÞdt þ rðIÞ yy ðt; xÞgðtÞdt þ
Z
1
1 Z 1 1
rðIIÞ yx ðt; xÞf ðtÞdt; rðIIÞ yy ðt; xÞf ðtÞdt;
ð2Þ
where the superscripts (I) and (II) stand for the quantities induced by the first (climbing) and second (sliding) dislocations, ðcÞ respectively, rab ðt; xÞ (c = I, II) denote the stresses at the point (x, y = 0) induced by a unit-intensity dislocation of the c-th family located at position (t, y = 0). For the subproblem B, the stress boundary condition on the crack surfaces is then written as:
rBxy ðxÞ ¼ rByy ðxÞ ¼
Z
1
1 Z 1 1
rðIÞ yx ðt; xÞgðtÞdt þ rðIÞ yy ðt; xÞgðtÞdt þ
Z
1
1 Z 1 1
rðIIÞ yx ðt; xÞf ðtÞdt; rðIIÞ yy ðt; xÞf ðtÞdt:
In addition, the single-value conditions of displacements should also be satisfied by the dislocation functions, i.e.,
ð3Þ
3524
X.-Q. Feng et al. / Engineering Fracture Mechanics 77 (2010) 3521–3531
ðxÞjx¼þ1 u x¼1 ¼ ðxÞjx¼þ1 x¼1
v
Z
¼
1
1 1
Z
1
Z 1 ðtÞ @u dt ¼ f ðtÞdt ¼ 0; @t 1 Z 1 @ v ðtÞ gðtÞdt ¼ 0; dt ¼ @t 1
ð4Þ
indicating the vanishing of the opening displacements at the left (x = 1) and right (x = 1) tips of the crack. ðcÞ As aforementioned, however, no analytical solution of the functions rab ðt; xÞ in Eq. (2) is available for a body with an arbitrary shape. Therefore, the subproblem B has been decomposed into C and D in Fig. 2. Then the dislocation-induced stresses in Eq. (2) can be additively decomposed as
rðacbÞ ðt; xÞ ¼ r^ ðacbÞ ðt; xÞ þ r~ ðacbÞ ðt; xÞ; ða; b ¼ 1; 2; c ¼ I; IIÞ;
ð5Þ
ðcÞ
^ ab ðt; xÞ denotes the stresses at point x induced by a unit-intensity dislocation at point t in the subproblem C, and where r r~ ðacbÞ ðt; xÞ the stresses induced by the image force on the boundary C in the subproblem D due to the presence of the same ^ ðacbÞ ðt; xÞ are given by [37] dislocation. From the well-known solution of an edge dislocation in an infinite and isotropic plate, r
^ ðIÞ r^ ðIÞ xx ðt; xÞ ¼ ryy ðt; xÞ ¼ p0
1 ; xt
r^ ðIÞ yx ðt; xÞ ¼ 0; ^ ðIIÞ r^ ðIIÞ xx ðt; xÞ ¼ ryy ðt; xÞ ¼ 0; r^ ðIIÞ yx ðt; xÞ ¼ p0
1 ; xt
ð6Þ
where p0 = G/[2p (1 m)], and G and m are the shear modulus and Poisson’s ratio of the material, respectively. Substituting Eqs. (5) and (6) into (3) leads to
rByx ðxÞ ¼ rByy ðxÞ ¼
Z
1
1 Z 1 1
r~ ðIÞ yx ðt; xÞgðtÞdt þ r~ ðIÞ yy ðt; xÞgðtÞdt þ
Z
1
1 Z 1 1
r~ ðIIÞ yx ðt; xÞf ðtÞdt þ p0 r~ ðIIÞ yy ðt; xÞf ðtÞdt p0
Z
1
1 Z 1 1
f ðtÞ dt; xt gðtÞ dt; xt
ð7Þ
which are singular integral equations with the Cauchy kernel (x t)1. ~ ðacbÞ ðt; xÞ in Eq. (5), induced by the image forces of the two families of dislocations with the densities g(x) and The stresses r f(x), can be easily solved from the defect-free subproblem D by using the finite element method, the boundary element method or other numerical methods. The image traction vector T(x) on a point x along the boundary C in D can be formulated as
T a ðxÞ ¼ rimage ab ðxÞnb ðxÞ ¼
Z
1
1
rðIÞ 0ab ðt; xÞgðtÞdt þ
Z
1
1
rðIIÞ 0ab ðt; xÞf ðtÞdt nb ðxÞ;
ð8Þ
ðcÞ
where r0ab ðt; xÞ stands for the image stresses at point x induced by a unit-intensity dislocation of the c-th family located at t, and n(x) is the outward unit vector normal to the boundary at x. Eq. (8) provides the boundary conditions on the crack face of the subproblem D. Thus, the two dislocation density functions g(x) and f(x) will be solved from Eq. (7) in conjunction with Eqs. (4) and (8). 2.3. Solution of the system of singular integral equations The numerical technique developed by Erdogan [4] is efficient for solving the system of simultaneous singular integral equations in Eq. (7). Denote [4,7–9]
pffiffiffiffiffiffiffiffiffiffiffiffiffi GðtÞ ¼ gðtÞ 1 t 2 ; pffiffiffiffiffiffiffiffiffiffiffiffiffi FðtÞ ¼ f ðtÞ 1 t2 :
ð9Þ
Then expressing the functions G(t) and F(t) in the series forms in terms of the Chebyshev polynomial of the first kind and truncating them to include only the first n terms, we obtain
GðtÞ ¼
n1 X
Ai T i ðtÞ;
i¼0
FðtÞ ¼
n1 X
Bi T i ðtÞ;
ð10Þ
i¼0
where Ti(t) is the Chebyshev polynomial of the first kind, and the coefficients Ai and Bi (i = 0, 1, 2, . . ., n 1) are constants as yet to be determined. Ti(t) is expressed as [1,4]
X.-Q. Feng et al. / Engineering Fracture Mechanics 77 (2010) 3521–3531
T i ðtÞ ¼ cos½i arccosðtÞ ¼
½i=2 X ð1Þj j¼0
ðiÞ! x2j ð1 x2 Þi2j : ð2jÞ!ði 2jÞ!
3525
ð11Þ
After discretization, Eqs. (7) and (4) can be rewritten as:
1 Fðt k Þ ¼ rByx ðxr Þ; xr t k ~ ðIIÞ Gðt k Þ þ r ðt ; x ÞFðt Þ ¼ rByy ðxr Þ; r k k yy
n pX ðIÞ ðIIÞ ~ yx r~ yx ðtk ; xr ÞGðt k Þ þ r ðt k ; xr Þ þ p0
n
k¼1
n pX
n
r~ ðIÞ yy ðt k ; xr Þ p0
k¼1
n X k¼1 n X
1 xr t k
Fðtk Þ ¼ 0; Gðt k Þ ¼ 0;
ð12Þ
k¼1
where the discretization points tk and xr are defined by
2k 1 p ; ðk ¼ 1; 2; . . . ; nÞ; 2n r xr ¼ cos p ; ðr ¼ 1; 2; . . . ; n 1Þ: n
tk ¼ cos
ð13Þ
The equation system in (12) includes 2n linear algebraic equations in terms of G(tk) and F(tk), and it can be solved easily via a numerical method. With increasing the term number n, the calculation results rapidly approach the exact solution. In our calculation, we truncate the series to include a sufficient number of terms in order to yield a solution of high accuracy. Thus, the stresses, strains, and displacements in the crack system can all be determined. 2.4. Calculation of SIFs Finally, the mode-I and mode-II SIFs at the two crack tips will be calculated according to their definitions:
K I ð1Þ ¼ K II ð1Þ ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2pjx 1jryy ðx; 0Þ; pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2pjx 1jryx ðx; 0Þ: lim þ
lim
x!1þ ;1
ð14Þ
x!1 ;1
From Eqs. (7)–(10), (and) (14) and using the following integral property of the Chebyshev polynomial of the first kind:
1
p
Z
1 1
ffii 1 T i ðtÞdt jxj jxj pffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi x x2 1 ; t x 1 t2 x x x2 1
ðjxj > 1;
i ¼ 0; 1; 2; . . .Þ;
ð15Þ
one can calculate the SIFs by
K I ð1Þ ¼ K 0 Gð1Þ ¼ K 0
n1 X ð1Þi Ai ; i¼0
n1 X K II ð1Þ ¼ K 0 Fð1Þ ¼ K 0 ð1Þi Bi ;
ð16Þ
i¼0
where
K 0 ¼ r0
pffiffiffiffi
p; r0 ¼ pp0 ¼
G : 2ð1 mÞ
ð17Þ
Evidently, the SIFs at the crack tips depend linearly on the coefficients Ai and Bi in the expansion series of the functions G(t) and F(t). In summary of the present dislocation-based semi-analytical or DFEM method, we use the superposition principle to decompose a two-dimensional crack problem in a finite elastic plate of arbitrary shape into two additive subproblems, namely, an infinite plate containing two families of continuous dislocations, and a defect-free body which has the same shape as the original one but contains no crack. Thus, the problem becomes the solution of a set of coupled singular integral equations and the calculation of the stress fields in a defect-free body. The Erdogan’s numerical technique [4] is further uti~ ðacbÞ ðt; xÞ induced by the image forces of a lized to solve the system of coupled singular integral equations, and the stress fields r specific distribution of climbing or gliding dislocations [e.g., each of the first n terms in Eq. (10)] can be calculated by the finite element method or other numerical methods. The advantages of this method include the avoidance of the singular stresses in the finite element calculation and the high accuracy of the obtained SIF results.
3526
X.-Q. Feng et al. / Engineering Fracture Mechanics 77 (2010) 3521–3531
What we have formulated above is only for a single straight crack in a finite plate of arbitrary shape, as shown in Fig. 1a. The formulation can be extended for many other configurations. When a plate contains m straight and thickness-through cracks, for example, we can easily write a system of 2m coupled singular integral equations similar to Eq. (7) and a system of 2m single-value conditions of displacements similar to Eq. (4). For the case of a smooth curved crack, the formulation is also similar but the discretization should be made along the actual crack shape and the corresponding image stresses should be modified by considering the rotation of the Burgers vectors of the two family of climbing and gliding edge dislocations. Besides the two-dimensional problems addressed above, another interesting case is an axisymmetric crack in a finite or semi-infinite body. For example, axisymmetric Hertzian contact cracks are often observed in the measurement of the fracture toughness of such materials as ceramics by using the indentation technique. In such a case, one can reformulate the above method by using the solutions of Somigliana ring dislocations [38,39] instead of those in Eq. (6).
3. Examples In this section, we will give several representative examples to illustrate the efficacy and accuracy of the above semi-analytical method. Since the tractions rcrack on the crack surface can be equivalently replaced by loads on the boundary C, we will address only traction-free cracks located in a finite plate with distributed tractions rexternal on its boundary.
3.1. A central crack in a rectangular plate We first consider a rectangular plate of width 2b containing a central crack of length 2c, as shown in Fig. 3a. A uniform tensile stress along the Y direction is applied on the boundary. The calculated SIFs with the increase in the number n of Chebyshev polynomial terms are shown in Fig. 3b and compared with those in the SIF handbook [3]. It is found that the present method yields a high accuracy of SIFs.
(a)
(b)
3.0 a / b = 1.0 DFEM, n = 0 DFEM, n = 1 DFEM, n = 2 DFEM, n = 3 DFEM, n = 4 DFEM, n = 5 DFEM, n = 6 Handbook, n = 0 Handbook, n = 1 Handbook, n = 2 Handbook, n = 3 Handbook, n = 4 Handbook, n = 5 Handbook, n = 6
( ) K I n / K0
2.5 2.0 1.5 1.0 0.5 0.0 0.0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
c/b Fig. 3. (a) A rectangular plate containing a central crack, and (b) the comparison of the calculated SIFs with those in the handbook [3].
X.-Q. Feng et al. / Engineering Fracture Mechanics 77 (2010) 3521–3531
3527
3.2. An inclined crack in a rectangular plate Now consider a rectangular plate of size 2a 2b, internally containing an inclined thickness-through crack with the length 2c and the angle h measured from the X1 axis, as shown in Fig. 4. A uniaxial tensile stress r0 along the Y1 direction is applied on the boundary. We can also conveniently calculate the SIFs via the method presented in Section 2. The calculapffiffiffiffiffiffi tion results are given in Fig. 5, where a/b = 1.0 and K 0 ¼ r0 pc. It is seen from Fig. 5a, where we take h = p/4, that with the
θ
Fig. 4. A rectangular plate containing an inclined crack.
KI / K0, KII / K0, KI / KII
(a) 3.5
a / b = 1.0
3.0
θ = 45
O
KI / K0
2.5
KII / K0 KI / KII
2.0 1.5 1.0 0.5 0.0
0.2
0.4
0.6
0.8
1.0
1.2
c/b
(b) 1.2
a / b = 1.0 c / b = 0.2 KI / K0
KI / K0, KII / K0
1.0 0.8
KII / K0
0.6 0.4 0.2 0.0 0
10
20
30
40
θ
50
60
70
80
o
Fig. 5. Variations of the SIFs of an inclined center-crack with respect to (a) the crack length and (b) the incline angle.
3528
X.-Q. Feng et al. / Engineering Fracture Mechanics 77 (2010) 3521–3531
c
increase of c/b in the range of 0 < c=b 6 0:2 the normalized SIFs increase slowly and KI/K0 KII/K0 0.5. This indicates that when the central crack is much shorter than the plate size, the boundary effect can be neglected and the SIFs of the crack can be approximated by those in an infinite plate. When a/b = 1.0 and c/b = 0.2, the dependences of KI/K0 and KII/K0 on the incline angle h are shown in Fig. 5b. The results of this example are also in a good agreement with those in the handbook [3].
θ
Fig. 6. An inclined crack at the boundary of a rectangular plate.
(a) 10
KI / K0, KII / K0, KI / KII
a1 / b = 0.5 8
a2 / b = 0.5 O
θ = 45
KI / K0
6
KII / K0 KI / KII
4 2 0 0.0
0.2
0.6
0.4
c/b
(b) KI / K0 , KII / K0 , KII / KI
1.5
1.0
c / b = 0.2 a1 / b = a2 / b = 0.5 KI / K0 KII / K0
0.5
KII / KI
0.0 0
10
20
30
40
θ
50
60
70
80
o
Fig. 7. Variations of the SIFs with the increase in (a) the crack length and (b) the incline angle.
X.-Q. Feng et al. / Engineering Fracture Mechanics 77 (2010) 3521–3531
3529
3.3. An inclined edge crack in a rectangular plate The third example addresses a rectangular plate containing an inclined edge crack with length c and angle h measured from the X1 axis, as shown in Fig. 6. It is subjected to uniform tension along the Y1 direction. The obtained results of the SIFs are plotted in Fig. 7 with respect to the crack length c and the incline angle h, respectively.
Fig. 8. A rectangular plate containing two collinear cracks.
(a) KI, 1L / K1, KI, 1R / K1, KI, 1R / KI, 1L
2.4
c1 / c2 = a / b = 1.0
2.2
e1 / b = e2 / b = 0.5
2.0
KI, 1L / K1 KI, 1R / K1
1.8
KI, 1R / KI, 1L
1.6 1.4 1.2 1.0 0.0
0.1
0.2
0.3
0.4
0.5
c1 / b
KI, 1L / K1, KI, 1R / K1, KI, 2L / K1, KI, 2R / K1
(b) 2.5
a / b = 1.0, c1 / b = 0.2 e1 / b = e2 / b = 0.5 KI, 1L / K1
2.0
KI, 1R / K1 KI, 2L / K1
1.5
KI, 2R / K1 1.0 0.5 0.0
0.0
0.5
1.0
1.5
2.0
c2 / c1 Fig. 9. Variations of the SIFs of two collinear cracks with respect to (a) the normalized crack length and (b) the crack length ratio of the two cracks.
3530
X.-Q. Feng et al. / Engineering Fracture Mechanics 77 (2010) 3521–3531
3.4. Two collinear cracks in a rectangular plate Finally, we give an example in which a rectangular plate of size 2a 2b contains two collinear and thickness-through cracks, as shown in Fig. 8. A uniaxial tensile stress r0 is assumed along the Y direction. The lengths of cracks are 2c1 and 2c2, and their eccentric distances are e1 and e2, respectively. For comparison, the following parameters in Ref. [3] are used in
our
ðnÞ K I;1R =K 0
¼
calculations: ðnÞ K I;2L =K 0
e1/e2 = a/b = 1.0
and
c1/b = c2/b = 0.2.
We
obtains
ðnÞ
ðnÞ
K I;1L =K 0 ¼ K I;2R =K 0 ¼ 1:096210
and
¼ 1:094521 , where the subscripts 1 and 2 stand for the left and right cracks, L and R indicate the left pffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi ðnÞ pc1 ¼ r0 pc2 . The handbook [3] gives K ðnÞ I;1L =K 0 ¼ K I;2R =K 0 ¼ 1:10
and right crack tips of a crack, respectively, and K 0 ¼ r0 ðnÞ
ðnÞ
and K I;1R =K 0 ¼ K I;2L =K 0 ¼ 1:09. Again, this example demonstrates the high accuracy of this numerical method. The variations of the normalized SIFs at the left and right tips of the right crack (E1F1) are plotted in Fig. 9a with respect to pffiffiffiffiffiffiffiffi c1/b, where K 1 ¼ r0 pc1 , c1/c2 = a/b = 1.0, and e1/b = e2/b = 0.5. Fig. 9b shows the influence of the crack length ratio c2/c1 on the SIFs of the two cracks, where a/b = 1.0, c1/b = 0.2, and e1/b = e2/b = 0.5. It is found that the SIFs at the four crack tips increase with increasing c2/c1. From this figure, one can analyze the effects of crack–crack interaction and crack–boundary interaction. In most cases, the longer crack will propagate first. Depending upon the geometric parameters, it may either propagate toward the free edge, evolving into an edge crack, or coalesce with the shorter crack. 4. Conclusions Based on the superposition principle and the analytical solutions of dislocations, we have presented a novel semi-analytical method to efficiently solve the SIFs of one or multiple cracks in a plate of arbitrary shape. By this means, the accuracy of the calculated SIFs does not depend on the mesh and can be totally controlled. Though the attention of this paper is focused on two-dimensional cases, this method can be extended to three-dimensional problems because the basic solutions of dislocations in the subproblem C and the finite element solution of the subproblem D have no limitation to two dimensions. Acknowledgements Supports from the National Natural Science Foundation of China (Grant Nos. 10732050, 10772058 and 10772093), Tsinghua University (2009THZ01001) and 973 Program of MOST (2010CB631005 and 2007CB936803) are gratefully acknowledged. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25]
Muskhelishvili NI. Singular integral equations. Groningen: Noordhoff; 1953. Tada H, Paris PC, Irwin GR. Stress analysis of cracks handbook. 3rd ed. St. Luis: Paris Production Incorporated; 2000. Aerospace Institute of China. Handbook of stress intensity factors. Beijing: Science Press; 1993. Erdogan F. Complex function technique. Continuum physics, vol. II. New York: Academic Press; 1975. Erdogan F. Stress distributions in bonded dissimilar materials with cracks. J Appl Mech 1965;32:403–10. Huang M, Long Y. Calculation of stress intensity factors of cracked Reissner plates by the sub-region mixed finite element method. Comput Struct 1988;30:837–40. Feng XQ, Xu M, Wang XY, Gu B. Fracture mechanics analysis of three-dimensional ion cut technology. J Mech Mater Struct 2007;2:1831–52. Gu B, Liu HY, Mai YW, Feng XQ, Yu SW. Fracture mechanics analysis on Smart-CutÒ technology. Part 1 – effects of stiffening wafer and defect interaction. Acta Mech Sinica 2009;25:73–81. Gu B, Yu SW, Feng XQ. Transient response of an interface crack between dissimilar piezoelectric layers under mechanical impacts. Int J Solids Struct 2002;39:1743–56. Banks-Sills L. Application of the finite element method to linear elastic fracture mechanics. Appl Mech Rev 1991;44:447–61. Camacho GT, Ortiz M. Computational modeling of impact damage in brittle materials. Int J Solids Struct 1996;33:2899–938. Sukumar N, Moës N, Moran B, Belytschko T. Extended finite element method for three-dimensional crack modeling. Int J Numer Methods Engng 2000;48:1549–70. Wyart E, Coulon D, Pardoen T, Remacle JF, Lani F. Application of the substructured finite element/extended finite element method (S-FE/XFE) to the analysis of cracks in aircraft thin walled structures. Engng Fract Mech 2009;76:44–58. Qin QH. The Trefftz finite and boundary element method. Southampton: WIT Press; 2000. Pan E. A general boundary element analysis of 2-D linear elastic fracture mechanics. Int J Fract 1997;88:41–59. Blandford GE, Ingraffea AR, Liggett JA. Two-dimensional stress intensity factor computations using the boundary element method. Int J Numer Methods Engng 1981;17:387–404. Portela A, Aliabadi MH, Rooke DP. The dual boundary element method: effective implementation for crack problems. Int J Numer Methods Engng 1992;33:1269–87. Garcia-Sanchu F, Zhang C, Saez A. 2-D transient dynamic analysis of cracked piezoelectric solids by a time-domain BEM. Computer Methods Appl Mech Engng 2008;197:3108–21. Qin QH, Lu M. BEM for crack-inclusion problems of plane thermopiezoelectric solids. Int J Numer Methods Engng 2000;48:1071–88. Qin QH, Mai YW. BEM for crack-hole problems in thermopiezoelectric materials. Engng Fract Mech 2002;69:577–88. Belytschko T, Krongauz Y, Organ D, Fleming M, Krysl P. Meshless methods: an overview and recent developments. Comput Methods Appl Mech Engng 1996;139:3–47. Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. Int J Numer Methods Engng 1994;37:229–56. Krysl P, Belytschko T. The element-free Galerkin method for dynamic propagation of arbitrary 3-D cracks. Int J Numer Methods Engng 1999;44:767–800. Belytschko T, Organ D, Krongauz Y. A coupled finite element–element-free Galerkin method. Comput Mech 1995;17:186–95. Buckner HF. A novel principle for the computation of the stress intensity factor. Z Angew Math Mech 1970;50:529–46.
X.-Q. Feng et al. / Engineering Fracture Mechanics 77 (2010) 3521–3531
3531
[26] Guagliano M, Vergani L. Model I stress intensity factors for curved cracks in gears by a weight functions method. Fatigue Fract Engng Mater Struct 2001;23:41–52. [27] Cundall PA, Hart RD. Numerical modeling of discontinua. Engng Comp 1992;9:101–13. [28] Cook BK, Jensen RP. Discrete element methods: numerical modeling of discontinua. ASCE 2002;22:22–3. [29] Abraham FF. Very large scale simulations of materials failure. Philos Trans Roy Soc London A 2002;360:367–82. [30] Abraham FF, Walkup R, Gao H, Duchaineau M, De la Rubia TD, Seager M. Simulating materials failure by using up to one billion atoms and the world’s fastest computer: brittle fracture. Proc Natl Acad Sci USA 2002;99:5777–82. [31] Pu SL, Hussain MA, Lorensen WE. Collapsed cubic isoparametric element as a singular element for crack problems. Int J Numer Methods Engng 1978;12:1727–42. [32] Carbonell R, Desroches J, Detournay E. A comparison between a semi-analytical and a numerical solution of a two-dimensional hydraulic fracture. Int J Solids Struct 1999;36:4869–88. [33] Gross D. Stress intensity factors of systems of cracks. Ing Arch 1982;51:301–10. [34] Feng XQ, Li JY, Yu SW. A simple method for calculating interaction of numerous microcracks and its applications. Int J Solids Struct 2003;40:447–64. [35] Ma L, Wang XY, Feng XQ, Yu SW. Numerical analysis of interaction and coalescence of numerous microcracks. Engng Fract Mech 2005;72:1841–65. [36] Eshelby JD. The force on an elastic singularity. Philos Trans Roy Soc London A 1951;244:87–112. [37] Hirth JP, Lothe J. Theory of dislocations. New York: Wiley; 1982. [38] Korsunsky A. The Somigliana ring dislocation revisited. 1. Papkovich potential solutions for dislocations in an infinite space. J Elast 1996;44:97–114. [39] Wang XY, Li LKY, Mai YW, Shen YG. Theoretical analysis of Hertzian contact fracture: ring crack. Engng Fract Mech 2008;75:4247–56.
Contents lists available at ScienceDirect
Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech
Dislocation-based semi-analytical method for calculating stress intensity factors of cracks: Two-dimensional cases Xi-Qiao Feng a,*, Yun-Fei Shi a, Xu-Yue Wang b, Bo Li a, Shou-Wen Yu a, Qiang Yang c a
AML, Department of Engineering Mechanics, Tsinghua University, Beijing 100084, China Department of Sciences, Harbin Institute of Technology, Shenzhen Graduate School, Shenzhen 518055, China c State Key Laboratory of Hydroscience and Engineering, Tsinghua University, Beijing 100084, China b
a r t i c l e
i n f o
Article history: Received 2 October 2009 Received in revised form 26 February 2010 Accepted 4 March 2010 Available online 11 March 2010 Keywords: Crack Stress intensity factor Dislocation Finite element method Singular integral equation Semi-analytical method
a b s t r a c t Determination of the stress intensity factors of cracks is a fundamental issue for assessing the performance safety and predicting the service lifetime of engineering structures. In the present paper, a dislocation-based semi-analytical method is presented by integrating the continuous dislocation model with the finite element method together. Using the superposition principle, a two-dimensional crack problem in a finite elastic body is reduced to the solution of a set of coupled singular integral equations and the calculation of the stress fields of a body which has the same shape as the original one but has no crack. It can easily solve crack problems of structures with arbitrary shape, and the calculated stress intensity factors show almost no dependence upon the finite element mesh. Some representative examples are given to illustrate the efficacy and accuracy of this novel numerical method. Only two-dimensional cases are addressed here, but this method can be extended to threedimensional problems. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction Despite the significant development of elastoplastic fracture mechanics methods in terms of such parameters as J-integral and crack opening displacement, the stress intensity factor (SIF) is still the most extensively utilized parameter for assessing the integrity and predicting the service life of engineering components and structures with cracks. Massive efforts have been directed towards developing analytical and numerical methods to calculate the SIFs under different conditions [1–9]. Analytical methods are available only for some simple geometry configurations and boundary conditions [1–3]. Therefore, various numerical approaches (e.g., finite element method [10–14], boundary element method [14–20], meshless method [21– 24], weight function method [25,26], and discrete element method [27–30]) have been established to calculate the SIFs of cracks. Among them, the finite element method has been used with great success for crack problems in a wide diversity of industrial fields. To determine the SIFs or other near-tip fields of a crack, very fine mesh should be generated in the vicinity of the crack front, making the numerical calculation time-consuming. An inappropriate choice of element types or sizes will lead to calculation results with considerable error and/or mesh dependence [10,11]. The boundary element method has been proved as another powerful technique to solve elastic solids with cracks. Different from the finite element method, the boundary element method discretizes only the body’s boundary [14–20]. Nevertheless, the reduced number of unknowns does not necessarily lead to improved efficiency because the boundary element method generally produces a full matrix, while the stiffness matrices in finite element simulations are usually sparse
* Corresponding author. Tel.: +86 10 62772934; fax: +86 10 62781824. E-mail address: [email protected] (X.-Q. Feng). 0013-7944/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2010.03.004
3522
X.-Q. Feng et al. / Engineering Fracture Mechanics 77 (2010) 3521–3531
[18]. The element-free or meshless method rose in 1990s and has been successfully applied in some fracture problems [21]. For example, the element-free Galerkin (EFG) method [22,23] is well suited to the modeling of two- and three-dimensional crack problems because arbitrary surfaces across which the displacement function is discontinuous can readily be incorporated in this technique. However, the EFG method is difficult in enforcing essential boundary conditions [24]. In addition, some semi-analytical methods have also been proposed for solving crack problems [31–35]. In the present paper, a semi-analytical method, also called dislocation-based finite element method (DFEM), is presented to solve the SIFs of cracks in a finite body of arbitrary shape. It combines the continuous dislocation model with finite element method. In Section 2, the two-dimensional cracked plate problem is formulated to show the basic idea of this method. In Section 3, several representative examples are given to illustrate its efficiency and accuracy.
2. Dislocation-based finite element method 2.1. Decomposition of a crack problem We consider a finite plate X containing a straight and thickness-through crack and subjected to the distributed tractions
rcrack and rexternal on the crack surface and its outer boundary C, respectively, as shown in Fig. 1a. The plate may have an arbitrary shape, and the crack may either intersect with the boundary C (e.g., a boundary or edge crack) or not (e.g., an interior crack). Assume that the material is isotropic and linear elastic. According to the superposition principle, the crack problem in Fig. 1a can be decomposed into two additive subproblems, A and B, as shown in Figs. 1b and c, respectively. In the former, the plate has no crack but is subjected to the distributed force rexternal on its boundary, whereas in the latter, the plate of the same shape is traction-free on its boundary but subjected to the traction rB on the crack surfaces. The superposition principle and the traction condition on the crack faces in Fig. 1a require that rB = rcrack rA, where rA denotes the stress distribution along the line at the original crack position in Fig. 1b. rA can be easily calculated either by an analytical or a numerical method (e.g., finite element method) because the subproblem A has no crack. Thus, the SIFs of the crack in the considered configuration in Fig. 1a are identical to those in the subproblem B. However, it is noticed that the subproblem B is still hard to solve because of the arbitrary shape and finite size of the plate. To manage these complexities, the superposition principle is utilized again to further decompose the subproblem B into C and D, neither of which contains a crack, as shown in Fig. 2. In the subproblem C, the plate is infinite in size and contains, instead of the crack, two families of continuously distributed dislocations, with Burgers vectors normal and parallel to the original crack faces, respectively. The dislocation density distribution functions along the crack length direction are yet to be determined. Let rdislocation and rC denote the dislocation-induced stresses along the original plate boundary C and along the original crack faces, respectively. In the subproblem D, the plate has the same shape and size as that in Fig. 1a and a distributed traction rimage is applied on its boundary in order to fulfill the boundary conditions. The superposition principle and the traction-free boundary condition on C in Fig. 1c require that the image force on the boundary C induced by the distributed dislocations should satisfy rimage = rdislocation. Thus, the problem in Fig. 1a becomes the solution of the
(a)
σ external
σ crack
(b)
(c)
σ external
σA
σ B = σ crack − σ A
Fig. 1. Decomposition of the crack problem: (a) the original crack configuration, (b) subproblem A, and (c) subproblem B.
3523
X.-Q. Feng et al. / Engineering Fracture Mechanics 77 (2010) 3521–3531
(a)
(b) C
D
σC
σD
σ dislocation σimage = − σdislocation Fig. 2. Decomposition of subproblem B into (a) subproblem C and (b) subproblem D.
dislocation density functions in the infinite body problem C coupled with the defect-free plate problem D, both of which are easier to be solved than the original one, as will be described below. Provided that the stress boundary condition on the crack surfaces in Fig. 1c can be satisfied, i.e., rC + rD = rB, then the stress field in B is exactly equal to the addition of those in C and D and the SIFs in B can be uniquely determined from the stress field in C. The dislocation density functions in the subproblem C will be solved approximately from the condition rC + rD = rB, as will be formulated in what follows. 2.2. Continuous dislocation model Refer to the Cartesian coordinate system (O–xy), as shown in Fig. 1a, where the origin O is located at the center of the crack and the x axis is along the crack direction. Without loss of any generality, the half-length of the crack is assumed to be unity, or in other words, all the length parameters are normalized by the crack half-length c. According to the continuous ; v gT can be equivalently replaced by two families of condislocation model, the crack opening displacement vector u ¼ fu tinuously distributed dislocations with the Burgers vector b = {bx, by}T, where the dislocations with the component bx and by can be considered as gliding and climbing edge dislocations, respectively [4,5,36]. The corresponding dislocation density functions are defined as [4]
@ðuþ u Þ @u ¼ ; @x @x @ v @ðv þ v Þ ; ¼ gðxÞ ¼ @x @x
f ðxÞ ¼
ð1Þ
where u+(t), u(t), v+(t), and v(t) denote the displacements of the two crack surfaces at position t along the x and y directions, respectively. Then the families of gliding and climbing dislocations with an infinitesimal segment dt have the intensities of f(t)dt and g(t)dt, respectively. The dislocation-induced stresses at a position (x, y = 0) on the crack faces can be expressed as
ryx ðxÞ ¼ ryy ðxÞ ¼
Z
1
1 Z 1 1
rðIÞ yx ðt; xÞgðtÞdt þ rðIÞ yy ðt; xÞgðtÞdt þ
Z
1
1 Z 1 1
rðIIÞ yx ðt; xÞf ðtÞdt; rðIIÞ yy ðt; xÞf ðtÞdt;
ð2Þ
where the superscripts (I) and (II) stand for the quantities induced by the first (climbing) and second (sliding) dislocations, ðcÞ respectively, rab ðt; xÞ (c = I, II) denote the stresses at the point (x, y = 0) induced by a unit-intensity dislocation of the c-th family located at position (t, y = 0). For the subproblem B, the stress boundary condition on the crack surfaces is then written as:
rBxy ðxÞ ¼ rByy ðxÞ ¼
Z
1
1 Z 1 1
rðIÞ yx ðt; xÞgðtÞdt þ rðIÞ yy ðt; xÞgðtÞdt þ
Z
1
1 Z 1 1
rðIIÞ yx ðt; xÞf ðtÞdt; rðIIÞ yy ðt; xÞf ðtÞdt:
In addition, the single-value conditions of displacements should also be satisfied by the dislocation functions, i.e.,
ð3Þ
3524
X.-Q. Feng et al. / Engineering Fracture Mechanics 77 (2010) 3521–3531
ðxÞjx¼þ1 u x¼1 ¼ ðxÞjx¼þ1 x¼1
v
Z
¼
1
1 1
Z
1
Z 1 ðtÞ @u dt ¼ f ðtÞdt ¼ 0; @t 1 Z 1 @ v ðtÞ gðtÞdt ¼ 0; dt ¼ @t 1
ð4Þ
indicating the vanishing of the opening displacements at the left (x = 1) and right (x = 1) tips of the crack. ðcÞ As aforementioned, however, no analytical solution of the functions rab ðt; xÞ in Eq. (2) is available for a body with an arbitrary shape. Therefore, the subproblem B has been decomposed into C and D in Fig. 2. Then the dislocation-induced stresses in Eq. (2) can be additively decomposed as
rðacbÞ ðt; xÞ ¼ r^ ðacbÞ ðt; xÞ þ r~ ðacbÞ ðt; xÞ; ða; b ¼ 1; 2; c ¼ I; IIÞ;
ð5Þ
ðcÞ
^ ab ðt; xÞ denotes the stresses at point x induced by a unit-intensity dislocation at point t in the subproblem C, and where r r~ ðacbÞ ðt; xÞ the stresses induced by the image force on the boundary C in the subproblem D due to the presence of the same ^ ðacbÞ ðt; xÞ are given by [37] dislocation. From the well-known solution of an edge dislocation in an infinite and isotropic plate, r
^ ðIÞ r^ ðIÞ xx ðt; xÞ ¼ ryy ðt; xÞ ¼ p0
1 ; xt
r^ ðIÞ yx ðt; xÞ ¼ 0; ^ ðIIÞ r^ ðIIÞ xx ðt; xÞ ¼ ryy ðt; xÞ ¼ 0; r^ ðIIÞ yx ðt; xÞ ¼ p0
1 ; xt
ð6Þ
where p0 = G/[2p (1 m)], and G and m are the shear modulus and Poisson’s ratio of the material, respectively. Substituting Eqs. (5) and (6) into (3) leads to
rByx ðxÞ ¼ rByy ðxÞ ¼
Z
1
1 Z 1 1
r~ ðIÞ yx ðt; xÞgðtÞdt þ r~ ðIÞ yy ðt; xÞgðtÞdt þ
Z
1
1 Z 1 1
r~ ðIIÞ yx ðt; xÞf ðtÞdt þ p0 r~ ðIIÞ yy ðt; xÞf ðtÞdt p0
Z
1
1 Z 1 1
f ðtÞ dt; xt gðtÞ dt; xt
ð7Þ
which are singular integral equations with the Cauchy kernel (x t)1. ~ ðacbÞ ðt; xÞ in Eq. (5), induced by the image forces of the two families of dislocations with the densities g(x) and The stresses r f(x), can be easily solved from the defect-free subproblem D by using the finite element method, the boundary element method or other numerical methods. The image traction vector T(x) on a point x along the boundary C in D can be formulated as
T a ðxÞ ¼ rimage ab ðxÞnb ðxÞ ¼
Z
1
1
rðIÞ 0ab ðt; xÞgðtÞdt þ
Z
1
1
rðIIÞ 0ab ðt; xÞf ðtÞdt nb ðxÞ;
ð8Þ
ðcÞ
where r0ab ðt; xÞ stands for the image stresses at point x induced by a unit-intensity dislocation of the c-th family located at t, and n(x) is the outward unit vector normal to the boundary at x. Eq. (8) provides the boundary conditions on the crack face of the subproblem D. Thus, the two dislocation density functions g(x) and f(x) will be solved from Eq. (7) in conjunction with Eqs. (4) and (8). 2.3. Solution of the system of singular integral equations The numerical technique developed by Erdogan [4] is efficient for solving the system of simultaneous singular integral equations in Eq. (7). Denote [4,7–9]
pffiffiffiffiffiffiffiffiffiffiffiffiffi GðtÞ ¼ gðtÞ 1 t 2 ; pffiffiffiffiffiffiffiffiffiffiffiffiffi FðtÞ ¼ f ðtÞ 1 t2 :
ð9Þ
Then expressing the functions G(t) and F(t) in the series forms in terms of the Chebyshev polynomial of the first kind and truncating them to include only the first n terms, we obtain
GðtÞ ¼
n1 X
Ai T i ðtÞ;
i¼0
FðtÞ ¼
n1 X
Bi T i ðtÞ;
ð10Þ
i¼0
where Ti(t) is the Chebyshev polynomial of the first kind, and the coefficients Ai and Bi (i = 0, 1, 2, . . ., n 1) are constants as yet to be determined. Ti(t) is expressed as [1,4]
X.-Q. Feng et al. / Engineering Fracture Mechanics 77 (2010) 3521–3531
T i ðtÞ ¼ cos½i arccosðtÞ ¼
½i=2 X ð1Þj j¼0
ðiÞ! x2j ð1 x2 Þi2j : ð2jÞ!ði 2jÞ!
3525
ð11Þ
After discretization, Eqs. (7) and (4) can be rewritten as:
1 Fðt k Þ ¼ rByx ðxr Þ; xr t k ~ ðIIÞ Gðt k Þ þ r ðt ; x ÞFðt Þ ¼ rByy ðxr Þ; r k k yy
n pX ðIÞ ðIIÞ ~ yx r~ yx ðtk ; xr ÞGðt k Þ þ r ðt k ; xr Þ þ p0
n
k¼1
n pX
n
r~ ðIÞ yy ðt k ; xr Þ p0
k¼1
n X k¼1 n X
1 xr t k
Fðtk Þ ¼ 0; Gðt k Þ ¼ 0;
ð12Þ
k¼1
where the discretization points tk and xr are defined by
2k 1 p ; ðk ¼ 1; 2; . . . ; nÞ; 2n r xr ¼ cos p ; ðr ¼ 1; 2; . . . ; n 1Þ: n
tk ¼ cos
ð13Þ
The equation system in (12) includes 2n linear algebraic equations in terms of G(tk) and F(tk), and it can be solved easily via a numerical method. With increasing the term number n, the calculation results rapidly approach the exact solution. In our calculation, we truncate the series to include a sufficient number of terms in order to yield a solution of high accuracy. Thus, the stresses, strains, and displacements in the crack system can all be determined. 2.4. Calculation of SIFs Finally, the mode-I and mode-II SIFs at the two crack tips will be calculated according to their definitions:
K I ð1Þ ¼ K II ð1Þ ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2pjx 1jryy ðx; 0Þ; pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2pjx 1jryx ðx; 0Þ: lim þ
lim
x!1þ ;1
ð14Þ
x!1 ;1
From Eqs. (7)–(10), (and) (14) and using the following integral property of the Chebyshev polynomial of the first kind:
1
p
Z
1 1
ffii 1 T i ðtÞdt jxj jxj pffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi x x2 1 ; t x 1 t2 x x x2 1
ðjxj > 1;
i ¼ 0; 1; 2; . . .Þ;
ð15Þ
one can calculate the SIFs by
K I ð1Þ ¼ K 0 Gð1Þ ¼ K 0
n1 X ð1Þi Ai ; i¼0
n1 X K II ð1Þ ¼ K 0 Fð1Þ ¼ K 0 ð1Þi Bi ;
ð16Þ
i¼0
where
K 0 ¼ r0
pffiffiffiffi
p; r0 ¼ pp0 ¼
G : 2ð1 mÞ
ð17Þ
Evidently, the SIFs at the crack tips depend linearly on the coefficients Ai and Bi in the expansion series of the functions G(t) and F(t). In summary of the present dislocation-based semi-analytical or DFEM method, we use the superposition principle to decompose a two-dimensional crack problem in a finite elastic plate of arbitrary shape into two additive subproblems, namely, an infinite plate containing two families of continuous dislocations, and a defect-free body which has the same shape as the original one but contains no crack. Thus, the problem becomes the solution of a set of coupled singular integral equations and the calculation of the stress fields in a defect-free body. The Erdogan’s numerical technique [4] is further uti~ ðacbÞ ðt; xÞ induced by the image forces of a lized to solve the system of coupled singular integral equations, and the stress fields r specific distribution of climbing or gliding dislocations [e.g., each of the first n terms in Eq. (10)] can be calculated by the finite element method or other numerical methods. The advantages of this method include the avoidance of the singular stresses in the finite element calculation and the high accuracy of the obtained SIF results.
3526
X.-Q. Feng et al. / Engineering Fracture Mechanics 77 (2010) 3521–3531
What we have formulated above is only for a single straight crack in a finite plate of arbitrary shape, as shown in Fig. 1a. The formulation can be extended for many other configurations. When a plate contains m straight and thickness-through cracks, for example, we can easily write a system of 2m coupled singular integral equations similar to Eq. (7) and a system of 2m single-value conditions of displacements similar to Eq. (4). For the case of a smooth curved crack, the formulation is also similar but the discretization should be made along the actual crack shape and the corresponding image stresses should be modified by considering the rotation of the Burgers vectors of the two family of climbing and gliding edge dislocations. Besides the two-dimensional problems addressed above, another interesting case is an axisymmetric crack in a finite or semi-infinite body. For example, axisymmetric Hertzian contact cracks are often observed in the measurement of the fracture toughness of such materials as ceramics by using the indentation technique. In such a case, one can reformulate the above method by using the solutions of Somigliana ring dislocations [38,39] instead of those in Eq. (6).
3. Examples In this section, we will give several representative examples to illustrate the efficacy and accuracy of the above semi-analytical method. Since the tractions rcrack on the crack surface can be equivalently replaced by loads on the boundary C, we will address only traction-free cracks located in a finite plate with distributed tractions rexternal on its boundary.
3.1. A central crack in a rectangular plate We first consider a rectangular plate of width 2b containing a central crack of length 2c, as shown in Fig. 3a. A uniform tensile stress along the Y direction is applied on the boundary. The calculated SIFs with the increase in the number n of Chebyshev polynomial terms are shown in Fig. 3b and compared with those in the SIF handbook [3]. It is found that the present method yields a high accuracy of SIFs.
(a)
(b)
3.0 a / b = 1.0 DFEM, n = 0 DFEM, n = 1 DFEM, n = 2 DFEM, n = 3 DFEM, n = 4 DFEM, n = 5 DFEM, n = 6 Handbook, n = 0 Handbook, n = 1 Handbook, n = 2 Handbook, n = 3 Handbook, n = 4 Handbook, n = 5 Handbook, n = 6
( ) K I n / K0
2.5 2.0 1.5 1.0 0.5 0.0 0.0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
c/b Fig. 3. (a) A rectangular plate containing a central crack, and (b) the comparison of the calculated SIFs with those in the handbook [3].
X.-Q. Feng et al. / Engineering Fracture Mechanics 77 (2010) 3521–3531
3527
3.2. An inclined crack in a rectangular plate Now consider a rectangular plate of size 2a 2b, internally containing an inclined thickness-through crack with the length 2c and the angle h measured from the X1 axis, as shown in Fig. 4. A uniaxial tensile stress r0 along the Y1 direction is applied on the boundary. We can also conveniently calculate the SIFs via the method presented in Section 2. The calculapffiffiffiffiffiffi tion results are given in Fig. 5, where a/b = 1.0 and K 0 ¼ r0 pc. It is seen from Fig. 5a, where we take h = p/4, that with the
θ
Fig. 4. A rectangular plate containing an inclined crack.
KI / K0, KII / K0, KI / KII
(a) 3.5
a / b = 1.0
3.0
θ = 45
O
KI / K0
2.5
KII / K0 KI / KII
2.0 1.5 1.0 0.5 0.0
0.2
0.4
0.6
0.8
1.0
1.2
c/b
(b) 1.2
a / b = 1.0 c / b = 0.2 KI / K0
KI / K0, KII / K0
1.0 0.8
KII / K0
0.6 0.4 0.2 0.0 0
10
20
30
40
θ
50
60
70
80
o
Fig. 5. Variations of the SIFs of an inclined center-crack with respect to (a) the crack length and (b) the incline angle.
3528
X.-Q. Feng et al. / Engineering Fracture Mechanics 77 (2010) 3521–3531
c
increase of c/b in the range of 0 < c=b 6 0:2 the normalized SIFs increase slowly and KI/K0 KII/K0 0.5. This indicates that when the central crack is much shorter than the plate size, the boundary effect can be neglected and the SIFs of the crack can be approximated by those in an infinite plate. When a/b = 1.0 and c/b = 0.2, the dependences of KI/K0 and KII/K0 on the incline angle h are shown in Fig. 5b. The results of this example are also in a good agreement with those in the handbook [3].
θ
Fig. 6. An inclined crack at the boundary of a rectangular plate.
(a) 10
KI / K0, KII / K0, KI / KII
a1 / b = 0.5 8
a2 / b = 0.5 O
θ = 45
KI / K0
6
KII / K0 KI / KII
4 2 0 0.0
0.2
0.6
0.4
c/b
(b) KI / K0 , KII / K0 , KII / KI
1.5
1.0
c / b = 0.2 a1 / b = a2 / b = 0.5 KI / K0 KII / K0
0.5
KII / KI
0.0 0
10
20
30
40
θ
50
60
70
80
o
Fig. 7. Variations of the SIFs with the increase in (a) the crack length and (b) the incline angle.
X.-Q. Feng et al. / Engineering Fracture Mechanics 77 (2010) 3521–3531
3529
3.3. An inclined edge crack in a rectangular plate The third example addresses a rectangular plate containing an inclined edge crack with length c and angle h measured from the X1 axis, as shown in Fig. 6. It is subjected to uniform tension along the Y1 direction. The obtained results of the SIFs are plotted in Fig. 7 with respect to the crack length c and the incline angle h, respectively.
Fig. 8. A rectangular plate containing two collinear cracks.
(a) KI, 1L / K1, KI, 1R / K1, KI, 1R / KI, 1L
2.4
c1 / c2 = a / b = 1.0
2.2
e1 / b = e2 / b = 0.5
2.0
KI, 1L / K1 KI, 1R / K1
1.8
KI, 1R / KI, 1L
1.6 1.4 1.2 1.0 0.0
0.1
0.2
0.3
0.4
0.5
c1 / b
KI, 1L / K1, KI, 1R / K1, KI, 2L / K1, KI, 2R / K1
(b) 2.5
a / b = 1.0, c1 / b = 0.2 e1 / b = e2 / b = 0.5 KI, 1L / K1
2.0
KI, 1R / K1 KI, 2L / K1
1.5
KI, 2R / K1 1.0 0.5 0.0
0.0
0.5
1.0
1.5
2.0
c2 / c1 Fig. 9. Variations of the SIFs of two collinear cracks with respect to (a) the normalized crack length and (b) the crack length ratio of the two cracks.
3530
X.-Q. Feng et al. / Engineering Fracture Mechanics 77 (2010) 3521–3531
3.4. Two collinear cracks in a rectangular plate Finally, we give an example in which a rectangular plate of size 2a 2b contains two collinear and thickness-through cracks, as shown in Fig. 8. A uniaxial tensile stress r0 is assumed along the Y direction. The lengths of cracks are 2c1 and 2c2, and their eccentric distances are e1 and e2, respectively. For comparison, the following parameters in Ref. [3] are used in
our
ðnÞ K I;1R =K 0
¼
calculations: ðnÞ K I;2L =K 0
e1/e2 = a/b = 1.0
and
c1/b = c2/b = 0.2.
We
obtains
ðnÞ
ðnÞ
K I;1L =K 0 ¼ K I;2R =K 0 ¼ 1:096210
and
¼ 1:094521 , where the subscripts 1 and 2 stand for the left and right cracks, L and R indicate the left pffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi ðnÞ pc1 ¼ r0 pc2 . The handbook [3] gives K ðnÞ I;1L =K 0 ¼ K I;2R =K 0 ¼ 1:10
and right crack tips of a crack, respectively, and K 0 ¼ r0 ðnÞ
ðnÞ
and K I;1R =K 0 ¼ K I;2L =K 0 ¼ 1:09. Again, this example demonstrates the high accuracy of this numerical method. The variations of the normalized SIFs at the left and right tips of the right crack (E1F1) are plotted in Fig. 9a with respect to pffiffiffiffiffiffiffiffi c1/b, where K 1 ¼ r0 pc1 , c1/c2 = a/b = 1.0, and e1/b = e2/b = 0.5. Fig. 9b shows the influence of the crack length ratio c2/c1 on the SIFs of the two cracks, where a/b = 1.0, c1/b = 0.2, and e1/b = e2/b = 0.5. It is found that the SIFs at the four crack tips increase with increasing c2/c1. From this figure, one can analyze the effects of crack–crack interaction and crack–boundary interaction. In most cases, the longer crack will propagate first. Depending upon the geometric parameters, it may either propagate toward the free edge, evolving into an edge crack, or coalesce with the shorter crack. 4. Conclusions Based on the superposition principle and the analytical solutions of dislocations, we have presented a novel semi-analytical method to efficiently solve the SIFs of one or multiple cracks in a plate of arbitrary shape. By this means, the accuracy of the calculated SIFs does not depend on the mesh and can be totally controlled. Though the attention of this paper is focused on two-dimensional cases, this method can be extended to three-dimensional problems because the basic solutions of dislocations in the subproblem C and the finite element solution of the subproblem D have no limitation to two dimensions. Acknowledgements Supports from the National Natural Science Foundation of China (Grant Nos. 10732050, 10772058 and 10772093), Tsinghua University (2009THZ01001) and 973 Program of MOST (2010CB631005 and 2007CB936803) are gratefully acknowledged. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25]
Muskhelishvili NI. Singular integral equations. Groningen: Noordhoff; 1953. Tada H, Paris PC, Irwin GR. Stress analysis of cracks handbook. 3rd ed. St. Luis: Paris Production Incorporated; 2000. Aerospace Institute of China. Handbook of stress intensity factors. Beijing: Science Press; 1993. Erdogan F. Complex function technique. Continuum physics, vol. II. New York: Academic Press; 1975. Erdogan F. Stress distributions in bonded dissimilar materials with cracks. J Appl Mech 1965;32:403–10. Huang M, Long Y. Calculation of stress intensity factors of cracked Reissner plates by the sub-region mixed finite element method. Comput Struct 1988;30:837–40. Feng XQ, Xu M, Wang XY, Gu B. Fracture mechanics analysis of three-dimensional ion cut technology. J Mech Mater Struct 2007;2:1831–52. Gu B, Liu HY, Mai YW, Feng XQ, Yu SW. Fracture mechanics analysis on Smart-CutÒ technology. Part 1 – effects of stiffening wafer and defect interaction. Acta Mech Sinica 2009;25:73–81. Gu B, Yu SW, Feng XQ. Transient response of an interface crack between dissimilar piezoelectric layers under mechanical impacts. Int J Solids Struct 2002;39:1743–56. Banks-Sills L. Application of the finite element method to linear elastic fracture mechanics. Appl Mech Rev 1991;44:447–61. Camacho GT, Ortiz M. Computational modeling of impact damage in brittle materials. Int J Solids Struct 1996;33:2899–938. Sukumar N, Moës N, Moran B, Belytschko T. Extended finite element method for three-dimensional crack modeling. Int J Numer Methods Engng 2000;48:1549–70. Wyart E, Coulon D, Pardoen T, Remacle JF, Lani F. Application of the substructured finite element/extended finite element method (S-FE/XFE) to the analysis of cracks in aircraft thin walled structures. Engng Fract Mech 2009;76:44–58. Qin QH. The Trefftz finite and boundary element method. Southampton: WIT Press; 2000. Pan E. A general boundary element analysis of 2-D linear elastic fracture mechanics. Int J Fract 1997;88:41–59. Blandford GE, Ingraffea AR, Liggett JA. Two-dimensional stress intensity factor computations using the boundary element method. Int J Numer Methods Engng 1981;17:387–404. Portela A, Aliabadi MH, Rooke DP. The dual boundary element method: effective implementation for crack problems. Int J Numer Methods Engng 1992;33:1269–87. Garcia-Sanchu F, Zhang C, Saez A. 2-D transient dynamic analysis of cracked piezoelectric solids by a time-domain BEM. Computer Methods Appl Mech Engng 2008;197:3108–21. Qin QH, Lu M. BEM for crack-inclusion problems of plane thermopiezoelectric solids. Int J Numer Methods Engng 2000;48:1071–88. Qin QH, Mai YW. BEM for crack-hole problems in thermopiezoelectric materials. Engng Fract Mech 2002;69:577–88. Belytschko T, Krongauz Y, Organ D, Fleming M, Krysl P. Meshless methods: an overview and recent developments. Comput Methods Appl Mech Engng 1996;139:3–47. Belytschko T, Lu YY, Gu L. Element-free Galerkin methods. Int J Numer Methods Engng 1994;37:229–56. Krysl P, Belytschko T. The element-free Galerkin method for dynamic propagation of arbitrary 3-D cracks. Int J Numer Methods Engng 1999;44:767–800. Belytschko T, Organ D, Krongauz Y. A coupled finite element–element-free Galerkin method. Comput Mech 1995;17:186–95. Buckner HF. A novel principle for the computation of the stress intensity factor. Z Angew Math Mech 1970;50:529–46.
X.-Q. Feng et al. / Engineering Fracture Mechanics 77 (2010) 3521–3531
3531
[26] Guagliano M, Vergani L. Model I stress intensity factors for curved cracks in gears by a weight functions method. Fatigue Fract Engng Mater Struct 2001;23:41–52. [27] Cundall PA, Hart RD. Numerical modeling of discontinua. Engng Comp 1992;9:101–13. [28] Cook BK, Jensen RP. Discrete element methods: numerical modeling of discontinua. ASCE 2002;22:22–3. [29] Abraham FF. Very large scale simulations of materials failure. Philos Trans Roy Soc London A 2002;360:367–82. [30] Abraham FF, Walkup R, Gao H, Duchaineau M, De la Rubia TD, Seager M. Simulating materials failure by using up to one billion atoms and the world’s fastest computer: brittle fracture. Proc Natl Acad Sci USA 2002;99:5777–82. [31] Pu SL, Hussain MA, Lorensen WE. Collapsed cubic isoparametric element as a singular element for crack problems. Int J Numer Methods Engng 1978;12:1727–42. [32] Carbonell R, Desroches J, Detournay E. A comparison between a semi-analytical and a numerical solution of a two-dimensional hydraulic fracture. Int J Solids Struct 1999;36:4869–88. [33] Gross D. Stress intensity factors of systems of cracks. Ing Arch 1982;51:301–10. [34] Feng XQ, Li JY, Yu SW. A simple method for calculating interaction of numerous microcracks and its applications. Int J Solids Struct 2003;40:447–64. [35] Ma L, Wang XY, Feng XQ, Yu SW. Numerical analysis of interaction and coalescence of numerous microcracks. Engng Fract Mech 2005;72:1841–65. [36] Eshelby JD. The force on an elastic singularity. Philos Trans Roy Soc London A 1951;244:87–112. [37] Hirth JP, Lothe J. Theory of dislocations. New York: Wiley; 1982. [38] Korsunsky A. The Somigliana ring dislocation revisited. 1. Papkovich potential solutions for dislocations in an infinite space. J Elast 1996;44:97–114. [39] Wang XY, Li LKY, Mai YW, Shen YG. Theoretical analysis of Hertzian contact fracture: ring crack. Engng Fract Mech 2008;75:4247–56.