- Email: [email protected]
cracks and their interactions
Theoretical and Applied Fracture Mechanics 38 (2002) 81–101 www.elsevier.com/locate/tafmec
Modeling of voids/cracks and their interactions C.H. Yang, A.K. Soh
*
Department of Mechanical Engineering, The University of Hong Kong, Pokfulam Road, Hong Kong, China
Abstract Numerical modeling of a two-dimensional elastic body containing multiple voids/cracks to study the interaction between these defects can be significantly simplified by developing special finite elements, each containing an internal circular/elliptic hole or a slit crack. These finite elements are developed using complex potentials and the conformal mapping technique. The elements developed can be divided into two categories, namely, the semi-analytic-type and hybrid-type elements. The latter element type is an improved version of the former due to the implementation of displacement continuity along the inter-element boundary. All the proposed elements can be easily combined with the conventional displacement elements, such as isoparametric elements, to analyze the above-mentioned problems without using complicated finite element meshes. Numerical examples have been employed to illustrate the modeling of voids/ cracks and their interactions. The results obtained using the semi-analytic-type elements are in good agreement with the theoretical results, and the corresponding results obtained using the hybrid-type elements show an improvement of the agreement with the theoretical results. However, the former element type is much easier to construct. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Circular/elliptical hole; Crack; Semi-analytic element; Hybrid element; Complex potentials; Conformal mapping
1. Introduction For brittle/quasi-brittle materials, the existence of multiple holes, or voids, and cracks or other forms significantly affects their mechanical behavior due to high stress concentrations. Thus, the ability to perform accurate analysis of a general two-dimensional body containing the above-mentioned defects is very important. To-date, many researchers have studied the local behavior of such defects using theoretical approaches [1–18]. However, all the theoretical solutions obtained are for
*
Corresponding author. Fax: +852-2858-5415. E-mail address: [email protected] (A.K. Soh).
relatively simple cases, such as a crack among several regularly spaced voids or a void among several regularly spaced cracks embedded in an infinite plate subjected to remote loading. These theoretical approaches cannot be applied to solve the problems of practical structures in which random distributions of many voids and cracks exist. Meanwhile, the finite element method has also been widely employed to solve such problems on practical structures. However, the conventional elements cannot ensure high accuracy in the vicinity of holes or cracks even if a very large number of elements are used. This is because the assumed displacement functions do not satisfy the prevailing conditions at the free boundary. Proposed in [19,20] are some modified methods to solve the
0167-8442/02/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 1 6 7 - 8 4 4 2 ( 0 2 ) 0 0 0 8 3 - 6
82
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
said problem. But the main weakness of these methods is that many elements are required to model the vicinity of the holes or cracks. Another approach to solve such problem is to construct special elements containing internal holes or cracks. Tong et al. [21] applied the hybrid-element formulation and complex variable technique to construct a special super-element for the analysis of the elastic stress intensity factors for the case of a plate containing an edge crack. Constructed in [22] are some Treffz-type finite elements with internal circular/elliptic holes or cracks based on two different variational formulae and applied Trefftz method for the use of trial functions, by employing the complex potentials and the conformal mapping technique. Two-dimensional semi-analytic finite elements [23] are constructed for a circular hole using complex potentials to define stress and strain functions inside the element directly. In the formulation of these elements, the boundary conditions at the free boundary of the hole are implemented explicitly. The main shortcoming of the formulation in [23] is the possible discontinuity along the inter-element boundaries. Such problems can be resolved by employing the hybrid displacement method based on the modified minimum potential energy principle [24]. It is worth noting that the usage of the above-mentioned elements leads to the use of only one element for the vicinity of hole/crack, and high accuracy is achieved in the stress concentrated region. In the present paper, the above-mentioned special elements containing an internal circular/elliptical hole or an internal crack will be developed using complex potentials and the conformal mapping technique. The above-mentioned hybrid displacement method will then be employed to solve the problem of discontinuity along inter-element boundaries. The versatility and accuracy of the special finite elements developed will be demonstrated by analyzing some two-dimensional elastic bodies containing multiple voids and cracks.
2. Basic theory and formulas For linear elastic analysis of a three-dimensional solid, the governing equations are as follows:
(a) Equilibrium conditions rij;j þ F i ¼ 0
ð1Þ
in V
(b) Stress–strain relations rij ¼ Dijkl ekl
or
eij ¼ Cijkl rkl
in V
ð2Þ
(c) Strain–displacement relations eij ¼ 12ðui;j þ uj;i Þ
in V
ð3Þ
(d) Displacement boundary conditions ui ¼ ui
on Su
ð4Þ
(e) Traction boundary conditions Ti ¼ T i
on Sr
ð5Þ
In the above equations, rij , eij , ui are the stress, strain and displacement components, respectively; F i , T i , ui are the prescribed body force and boundary traction and displacement components, respectively; Dijkl , Cijkl are the elastic and compliance constants, respectively; V is the volume of the solid; and Sr , Su are the portion of the boundary on which surface tractions and boundary displacements are prescribed, respectively. To consider a two-dimensional plate with holes or cracks in plane elasticity, the fields of stress and displacement are best expressed in terms of the complex potentials proposed by Muskhelishvili [25] as follows: n o rx þ ry ¼ 2 w0 ðzÞ þ w0 ðzÞ ð6Þ n o rx ry þ 2isxy ¼ 2 zw00 ðzÞ þ v0 ðzÞ
ð7Þ
2lðux þ iuy Þ ¼ jwðzÞ zw0 ðzÞ vðzÞ
ð8Þ
h iB Fx þ iFy ¼ i wðzÞ þ zw0 ðzÞ þ vðzÞ A h i 0 ¼ i wðzB Þ þ zB w ðzB Þ þ vðzB Þ CA
ð9Þ
where z ¼ x þ iy; j ¼ 3 4m for plane strain and j ¼ ð3 mÞ=ð1 þ mÞ for plane stress; A, B are the starting and end points of integration respectively; CA is a complex constant at the datum point A; (0 ) and (00 ) denote the first and second orders of complex differentiation with respect to z, respectively. As the treatment of boundary conditions in
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
the original domain usually fails, it is convenient and advantageous to make use of conformal mapping, for which the mapping function is given as follows: m z ¼ f ð1Þ ¼ c 1 þ ð10Þ 1 where 1 ¼ n þ ig and its inverse function is 1 ¼ f 1 ðzÞ. By means of conformal mapping, the original domains with holes or cracks are mapped onto the alternative boundaries as unit circles. Note that the advantage of such mapping is that an important boundary of an element can be mapped to a boundary that is geometrically simpler, such as a unit circle. In the case of an elliptical hole with major axis a and minor axis b, the arbitrary constants of Eq. (10) are given by c ¼ ða þ bÞ=2 and m ¼ ða bÞ=ða þ bÞ. Obviously, the circle and the crack can be considered as special cases by letting a ¼ b ¼ r for the circle with radius r, and b ¼ 0 for the crack with length 2a. The inverse mapping function for Eq. (10) is pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
1 1¼ z z2 4c2 m 2c
ð11Þ
To map the exterior of an ellipse onto the exterior of the unit circle, the sign of the root in Eq. (11) is chosen such that j1j P 1:0. In using conformal mapping, the complex functions are not assumed in the original domain but in the transformed element domain. Therefore, the stress and displacement fields are expressed in terms of f ð1Þ and the new functions wð1Þ and vð1Þ, chosen in the 1-plane, as follows: ( rx þ ry ¼ 2
0
0
w ð1Þ w ð1Þ þ f 0 ð1Þ f 0 ð1Þ
Fx þ iFy ¼ i wð1B0 Þ þ f ð1B0 Þ
ð12Þ
rx ry þ 2isxy ¼ ( ! ) w00 ð1Þ f 00 ð1Þ v0 ð1Þ 0 2 f ð 1Þ w ð1Þ 03 þ 0 f 02 ð1Þ f ð1Þ f ð1Þ
f 0 ð1Þ
f 0 ð1B0 Þ
# þ vð1B0 Þ CA0 ð15Þ
where A0 , B0 correspond to the original integration points A and B; and CA0 is also a complex constant; (0 ) and (00 ) denote the first and second orders of differentiation with respect to 1, respectively. Due to the simple boundary of the unit circle, the polar coordinates (q; h) are introduced in the transformed domain, and Eqs. (12)–(14) can be rewritten as follows: ) ( w0 ð1Þ w0 ð1Þ rq þ rh ¼ 2 ð16Þ þ f 0 ð1Þ f 0 ð1Þ rh rq þ 2isqh ( ¼ 2e2ia
w00 ð1Þ f 00 ð1Þ 0 f ð1Þ ð1Þ w f 02 ð1Þ f 03 ð1Þ
!
v0 ð1Þ þ 0 f ð1Þ
vð1Þ
ð13Þ
ð14Þ
)
ð17Þ 2lðuq þ iuh Þ ¼ e
ia
jwð1Þ f ð1Þ
w0 ð1Þ f 0 ð1Þ
! vð1Þ ð18Þ
where a is the angle between the axes x and 1, eia ¼
1 f 0 ð1Þ q jf 0 ð1Þj
and e2ia ¼
12 f 0 ð1Þ q2 f 0 ð1Þ
For the free boundary of a hole of unit circle in a plate, boundary conditions to be satisfied at 1 ¼ 10 ¼ eih or j1j ¼ 1:0 are ð19Þ
and Fx ¼ Fy ¼ 0
w0 ð1Þ
w0 ð1B0 Þ
rq ¼ sqh ¼ 0 in the 1-plane
)
2lðux þ iuy Þ ¼ jwð1Þ f ð1Þ
"
83
in the z-plane
From Eqs. (16), (17) and (19) ) ( ( w0 ð10 Þ w0 ð10 Þ w00 ð10 Þ 2ia f ð10 Þ e þ 0 f ð10 Þ f 0 ð10 Þ f 02 ð10 Þ ! ) f 00 ð10 Þ v0 ð10 Þ 0 w ð10 Þ 03 þ 0 ¼0 f ð10 Þ f ð10 Þ
ð20Þ
ð21Þ
84
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
Substituting Eq. (20) into Eq. (15) and letting CA0 ¼ 0, it is found that wð10 Þ þ f ð10 Þ
c v0 þ c 0 1 þ
w0 ð10 Þ
þ vð10 Þ ¼ 0
f 0 ð10 Þ
wð1Þ ¼ w0 þ
k
ak 1 þ
k¼1
lU X
vð1Þ ¼ v0 þ
cp 1p þ
p¼1
bl 1
ð23Þ
dq 1q
ð24Þ
q¼1
or in accordance with the expressions provided in [3]: w0z ð1Þ ¼ a0 þ
kU X
ak 1 k þ
k¼1
( 0
v ð1Þ ¼ c c0 þ
lU X
bl 1l
ð25Þ
l¼1
pU X
p
cp 1 þ
p¼1
qU X
) dq 1
q
ð26Þ
q¼1
where w0z ð1Þ is the first order complex derivative of wð1Þ with respect to z; kU , lU , pU , and qU are the upper limits of the series, which are related to the total number of nodes in the special element; w0 , v0 , ak , bl , cp , and dq are complex coefficients of the v0 , ak ¼ a^k þ i~ ak , forms w0 ¼ w^0 þ iw~0 , v0 ¼ v^0 þ i~ bl ¼ b^l þ ib~l , cp ¼ c^p þ i~ cp , and dq ¼ d^q þ id~q . Note that the complex constants w0 and v0 in Eqs. (23) and (24) are also used as the integral constants for Eqs. (25) and (26), respectively. Thus, from Eqs. (25) and (26), we obtain wð1Þ ¼ ( c w0 þ
a0 1 þ
kU X k¼1
m
a0 11 þ
lU X 1 1 ak 1kþ1 b 1ðl1Þ kþ1 l 1 l l¼2
kU X k¼2
q
)
The following conditions are to be introduced in Eqs. (27) and (28) to ensure that ux and uy are single-valued displacements: b1 ma1 ¼ 0 ð29Þ d1 ¼ 0 In addition, from Eqs. (22)–(24), we obtain
l
l¼1 qU X
U X 1 1 cp 1pþ1 dq 1ðq1Þ pþ1 q 1 q¼2
ð28Þ
w0 mb2 þ 2 a2 þ v0 m c2 ¼ 0 or v ¼ w þ mb 2a2 þ mc 0
pU X
pU X p¼1
ð22Þ
Thus, the boundary conditions are expressed in the 1-plane. The holomorphic functions wð1Þ and vð1Þ can be approximated, using Laurent series of 1 in the region where j1j P 1:0, as follows: kU X
vð1Þ ¼ (
!
U X 1 1 ak 1k1 bl 1ðlþ1Þ k1 l þ 1 l¼1
l
!)
ð30Þ
2
and by using Eqs. (25)–(28), we obtain w0 þ mb1 þ a1 þ v0 ¼ 0 or v0 ¼ w0 ð1 þ m2 Þa1
ð31Þ
From Eqs. (30) and (31), the relationship between the two integral constants v0 and w0 is established. This leads to a reduction of the number of unknown parameters and, thus, enables the acquisition of higher order terms in those complex series. The above-mentioned complex potentials w and v in the form of Laurent series of 1 will be used to deduce the detailed expressions for the displacement distribution matrix [U] and the stress distribution matrix [P]. The only difference from the special semi-analytic-type element and the special hybrid-type element is that all the coefficients related to the rigid body movement will be set to zero in prior in the latter.
3. Stiffness matrix of a semi-analytic-type element with a central circular/elliptic hole or an internal crack Fig. 1(a)–(c) show the typical eight-node rectangular elements containing a central circular and an elliptic hole and an internal crack, respectively. e The nodal displacement vector fqg is defined as fqge ¼ ½ ðux Þ1
ð27Þ
2
0
ðuy Þ1
ðux Þ2
ðuy Þ2
ðux Þn
ðuy Þn T ð32Þ
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
85
Fig. 1. Eight-node rectangular membrane element with a circular/elliptic hole or an internal crack.
where ðux Þi and ðuy Þi are the x and y displacement components of node i ði ¼ 1; . . . ; nÞ, respectively; and n is the total number of the nodes of the element. In order to express wð1Þ and vð1Þ in terms of e fqg , the upper limits in Eqs. (22) and (23) are chosen such that the total number of independent terms in wð1Þ and vð1Þ is approximately equal to that of the independent displacement components. Substituting Eqs. (23) and (24) or Eqs. (25)–(28) into Eq. (19) or (21), respectively, two sets of equations relating the complex coefficients are obtained. By back-substituting these equations into Eqs. (23) and (24) or Eqs. (25)–(28) and through the use of Eq. (14), the strain field within the element is determined. The displacement vector fqg at any point in the element can be expressed in terms of the above-mentioned independent coefficients as follows: ux fug ¼ ¼ ½Ufbg ð33Þ uy where fbg is a matrix consisting of the unknown independent coefficients in Eqs. (23) and (24) or Eqs. (25)–(28), ½U is the displacement distribution matrix given by ½U ¼ fUx ðq; hÞ Uy ðq; hÞgT It is worth noting that the complex potentials in Eqs. (23) and (24) provide very high accuracy in
dealing with circular holes, but they create complications in handling non-circular holes. However, Eqs. (25)–(28) provide relatively simple expressions for all the components in both fbg and ½U, and maintain good accuracy in all cases. Substituting the nodal displacements at all the eight nodal points of the element into Eq. (33), we obtain e
e
fqg ¼ ½U fbg
ð34Þ
where ½Ue ¼ f ðUx Þ1
ðUy Þ1
ðUx Þ2
ðUy Þ2
...
ðUx Þn
ðUy Þn gT
Further, 1
fbg ¼ ½Ue fqg
e
ð35Þ
Substituting Eq. (35) into Eq. (33) fug ¼
ux uy
1
e
¼ ½U½Ue fqg ¼ ½Nfqg 1
e
ð36Þ
where ½N 22n ¼ ½U½Ue is termed the elementshape function in the conventional finite element method. For the standard finite element approach, the strain vector feg at any point in the element can be obtained as follows:
86
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
8 9 2o e > = 6 ox < x> 0 ¼6 feg ¼ ey 4 > ; : > o exy oy 2 oU 3
0
3
4. Stiffness matrix of a hybrid-type element with a central circular/elliptic hole or an internal crack
o 7 7 oy 5fqg o ox
x
6 6 ¼ 2l6 4
ox oUy oy oUx oy
oUy ox
þ
7 7 e 1 e e 7½U fqg ¼ ½Bfqg 5
ð37Þ
One important version of the modified minimum potential energy principle for linear elastic analysis of a solid body, developed in [24], is given as, Pep
where ½B is the element strain–displacement matrix. Since the sub-matrices Ux ðq; hÞ and Uy ðq; hÞ are in terms of the coordinates ðq; hÞ, there results " oU # " oU # " oUy # " oUy #
¼
oq
ox oUx oy
¼ ½J
oUx oh
where " ox 2 6 ¼4
ox oh
oUy oh
¼ ½J
ox oUy oy
ð38Þ
oy oh
cm 1 qm2 cos h
cm 1 þ qm2 q sin h
From Eq. (38), we have " oUx # " oUx # ¼ ½J
3
cm 1 þ qm2 sin h 7
5 m cm 1 q2 q cos h
oq
2
and
oUx oh
oUy 4 ox oUy oy
3
2
5 ¼ ½J1 4
oUy oq oUy oh
3 5
ð39Þ where ½J1 ¼
"
oy oq
ð41aÞ
Pep ¼
Z ve
Z
Z 1 T T T e Ee F u dv T ~u ds 2 ser TT ð~u uÞ ds
ð41bÞ
ove
1
T~i ðui u~i Þ ds
or in matrix form,
is the Jacobian matrix related to the mapping function. Moreover, except at the point (0,0), ! " m2 m 2 jJj ¼ cm q 1 þ 4 2 2 cos 2h > 0 q q
ox oUx oy
Z
Z 1 Dijkl eij ekl F i ui dv T i u~i ds 2 ser
ove
#
oy oq
oq
½J ¼
oq
and
ve
x
x
Z
1 jJj ox oh
oy oq ox oh
#
" oq ¼
ox oq oy
oh ox
#
oh oy
The element stiffness matrix is given by Z e T ½Ks ¼ ½B ½D½BjJjt dA
ð40Þ
Ae
where ½D is the element elasticity matrix, t is the plate thickness and Ae is the area of an element.
where ve is the element volume; ser is the surface traction boundary; and ove is the boundary of the element volume. For displacement-type finite elements, the displacements are normally required to be continuous not only inside the element but also along the inter-element boundary to ensure the stability and precision of calculations. The displacement continuity condition along the inter-element boundary SC is uðCÞ ¼ ~uðCÞ
on SC
ð42Þ
The assumed linear or quadratic displacement function along the inter-element boundary can be expressed as h
i u~ðCÞ x ~uðCÞ ¼ ¼ ½Lfqge u~ðCÞ y
ð43Þ
where [L] is the interpolation function matrix only defined on SC [21,26]. In principle, each quantity in Eqs. (41a) and (41b) should be assumed individually but in fact the appropriate constructions of r, u and T can make Eqs. (1)–(5) automatically satisfied in advance. Further by considering the divergence the-
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
87
orem, the simplified form of the functional (41a) and (41b) used in the proposed element with internal hole/crack is Z Z $ ðCÞ %T $ ðCÞ %T ðCÞ 1 e ~ Pp ¼ T u ds T ð44Þ u ds 2 SC SC
5. Calculation of stresses and stress intensity factors
From Eqs. (12) and (13),
where frg ¼ f rx ry sxy g and ½S is the element stress–displacement matrix given by ½S ¼ ½D½B. In fact, there is a direct method to calculate the stresses once the independent parameter vector fbg for the complex potentials w and v have been determined, see also Eq. (45). With reference to Fig. 1, by assuming that b ¼ 0, the elliptic hole degenerates to a slit crack with length 2a. In this case, the mode I and II stress intensity factors (SIFs), KI and KII , respectively, at the crack tip need to be determined. The formulae for calculating the complex stress intensity factor K at the crack tip has been given in [27]. In general, the stress field around the crack tip for the mixedmode I and II fracture based on the coordinate systems (x0 ; y 0 ) and (r0 ; h0 ), refer to Fig. 2, is given as follows:
frg ¼ ½Pfbg
frg ¼ ½Dfeg ¼ ½Sfqg
where ½P is the stress distribution matrix and T ½P ¼ f Px ðq; hÞ Py ðq; hÞ Pxy ðq; hÞ g . The traction boundary conditions (5) on SC can be written as ( ) & ' TxðCÞ ðCÞ ½T ¼ ð46Þ ¼ RðCÞ fbg on SC ðCÞ Ty ^ðCÞ is a 2 3 unit matrix nðCÞ P, n where RðCÞ ¼ ^ outward normal to the common boundary SC . Alternatively, " # " # ðCÞ PðCÞ RðCÞ ðCÞ x nx þ Pxy ny x ðq; hÞ ½R ¼ ¼ ðCÞ RðCÞ PðCÞ y ðq; hÞ y ny þ Pxy nx By substituting the above relations, the hybrid functional becomes
KI h0 rx0 ¼ pffiffiffiffiffiffiffiffiffi cos 2 2pr0
KI h0 ry 0 ¼ pffiffiffiffiffiffiffiffiffi cos 2 2pr0 0
SC
Furthermore, its stationary value with respect to b is given by ð48Þ
Since the element strain energy Vs:e: is Vs:e: ¼ 12bT Hb ¼ 12qT GT H1 Gq
e
T
1
sx0 y 0
h0 3h0 2 þ cos cos 2 2 ! h0 3h0 1 þ sin sin 2 2 0
KII h h 3h þ pffiffiffiffiffiffiffiffiffi sin cos cos 2 2 2 2pr0 0 0 0 KI h h 3h ¼ pffiffiffiffiffiffiffiffiffi cos sin cos 2 2 2 2pr0 KII h0 þ pffiffiffiffiffiffiffiffiffi cos 2 2pr0
! !
ð52Þ
0
h0 3h0 1 sin sin 2 2
!
Moreover, ð49Þ
Therefore, the element stiffness matrix ½Keh can be expressed as ½Kh ¼ ½G ½H ½G
ð51Þ
h0 3h0 1 sin sin 2 2
KII h0 pffiffiffiffiffiffiffiffiffi sin 2 2pr0
ð47Þ
where Z Z 1 T H¼ R U ds ðRT U þ UT RÞ ds 2 SC SC Z G¼ RT L ds
oPev ¼ 0 ) b ¼ H1 Gq ob
e
T
ð45Þ
Pep ¼ 12bT Hb bT Gq
The stress vector at a point in the element can be determined from
ð50Þ
KI h0 KII h0 rx0 þ ry 0 ¼ 2 pffiffiffiffiffiffiffiffiffi cos pffiffiffiffiffiffiffiffiffi sin 2 2 2pr0 2pr0
! ð53Þ
Eq. (53) can be rewritten in complex form as follows:
88
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
Thus, if the complex coefficients of the complex potential wð1Þ can be determined, the SIFs can be obtained from Eq. (56). Some theoretical solutions for different crack problems are provided below for easy reference. For the mixed-mode crack problem shown in Fig. 3(a), the complex potential wð1Þ can also be derived theoretically using complex variables and expressed as follows: wð1Þ ¼ Fig. 2. The local coordinate system at the crack tip.
( ! rx þ ry ¼ 2Re K
1 2pðz z0 Þ
"1=2 )
¼ 4Re½w0 ðzÞ ð54Þ
where K ¼ KI iKII and z0 ¼ a. Note that the points z ¼ a in the complex plane are at 1 ¼ 1 in the mapped plane. Eq. (54) gives ! "1=2 1 0 K ¼ 2w ðzÞ ð55Þ 2pðz z0 Þ Let z ¼ z0 or 1 ! 1, pffiffiffiffiffiffi 1=2 K ¼ 2 2p lim ðz z0 Þ w0 ðzÞ ¼ 2 z!z0
rffiffiffi p 0 w ð1Þ1¼1 a ð56Þ
r0 að1 e2ib Þ 41
ð57Þ
where b is the angle between the centerline of the crack and the axis of the applied load. Substituting Eq. (57) into Eq. (56), the result is % % 1 pffiffiffiffiffiffi$ 1 pffiffiffiffiffiffi$ K ¼ r0 pa 1 e2ib 12 1¼1 ¼ r0 pa 1 e2ib 2 2 ð58Þ Therefore, pffiffiffiffiffiffi KI ¼ r0 pa sin2 b and
pffiffiffiffiffiffi KII ¼ r0 pa sin b cos b ð59Þ
A more general expression for the stress intensity factor of mode I crack in a finite elastic body is [28] pffiffiffiffiffiffi KI ¼ r0 paf ðaW Þ ð60Þ where r0 is the applied stress in the x or y direction and f ðaW Þ is a configuration factor which accounts for proximity effects of boundary surfaces.
Fig. 3. A square plate with an inclined crack or a single edge crack subjected to tension in one direction.
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
(a)
89
(b)
Fig. 4. A square plate with single crack or two collinear cracks subjected to tension in one direction.
For the case of a wide plate containing a long edge crack, as shown in Fig. 3(b), the opening mode stress intensity factor KI is given by 2 W W f ðaW Þ ¼ 1:12 0:231 þ 10:55 a a 3 4 W W 21:72 þ 30:39 ð61Þ a a where W is the width of the plate. The above factor is highly accurate (0.5% error) for W =a < 0:6. For the case of a wide plate containing a crack or two collinear cracks, as depicted in Fig. 4(a) and (b), respectively, the opening mode stress intensity factor KI is given by 1=2 W pa tan f ðaW Þ ¼ ð62Þ pa W where W is either the width of the plate or the distance between the centerlines of two collinear cracks.
micro-defects, the stress field in the vicinity of a macro-crack is greatly disturbed such that the original mode I crack problem becomes a mixedmode crack problem. For mixed-mode crack problems, the initial crack growth direction is a very important parameter that has to be determined. Proposed in [29–32] is a way to solve mixedmode fracture problems using the concept of strain energy density, which can be used to predict the crack growth direction. The strain energy density dW =dV has a singularity of inverse r near the crack tip given by $ %1 dW 2 ¼ a11 KI2 þ 2a12 KI KII þ a22 KII2 þ a33 KIII dV r ð63Þ It is apparent that the coefficient of the 1/r term in the above equation is defined only if r 6¼ 0. Hence, the element volume is always kept at a finite distance away from the crack border. This coefficient is termed as the strain energy density factor, i.e., 2 S ¼ a11 KI2 þ 2a12 KI KII þ a22 KII2 þ a33 KIII
6. Crack growth criterion The possible change of crack growth direction due to the interaction of different defects in an elastic body is of great interest to many researchers. It is worth noting that due to the influences of
ð64Þ
where KI , KII and KIII are the mode I, II and III stress intensity factors, respectively, which vary from point to point along the crack border. The coefficients aij , which vary with the polar angle d measured from the crack tip, as shown in Fig. 2, are given by
90
a11 ¼
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
1 ½ð3 4m cos dÞð1 þ cos dÞ 16l
ð65Þ
7. Numerical examples 7.1. Circular holes
a12 ¼ a22 ¼
a33 ¼
1 2 sin d½ cos d ð1 2mÞ 16l 1 ½4ð1 mÞð1 sin dÞ 16l þ ð1 þ cos dÞð3 cos d 1Þ 1 16l
ð66Þ
ð67Þ ð68Þ
where l and m are the shear modulus of elasticity and Poisson’s ratio, respectively. For two-dimensional problems where the crack extends in the xyplane, the stress intensity factors do not vary along the crack front, and the strain energy density factor S is only the function of the angle d. It is assumed that the crack will begin to extend in a direction, for which the strain energy density factor possesses a relative minimum value, i.e., 2 oS oS ¼ 0; >0 ð69Þ od d¼d0 od2 d¼d0 Based on oS=od ¼ 0 at d0 , the direction of initial crack growth d0 can be determined. Note that the above total strain energy is evaluated along a constant circle with radius r around the crack tip based on the theory of elasticity.
(a)
Fig. 5(a) and (b) show a square plate with a central circular hole subjected to two different loadings. The model devised for the finite element analyses using the semi-analytic-type element is shown in Fig. 6(a). Note that the proposed element was employed to model the region surrounding the circular hole, and 24 eight-node isoparametric elements were used to model the remaining region of the plate. Tables 1–3 present the normalized hoop stresses, rh =r0 , along the free boundary of the hole when the plate was subjected to (1) uniform tensile stresses in the x direction; (2) uniform tensile stresses in both the x and y directions; and (3) uniform tensile and compressive stresses in the x and y directions, respectively. Ten analyses were performed, by varying the ratio of the radius of the hole to half the side length of the proposed element, r=A, for each of the load case studied. It is obvious that the numerical results obtained using the proposed element are in good agreement with the corresponding theoretical results for all the load cases studied. Note that the maximum discrepancy is less than 5% for all the three load cases when r=A 6 0:7. Moreover, the maximum discrepancy is still less than 8% for all the three load cases when 0:7 6 r=A 6 0:8.
(b)
Fig. 5. A square plate with a central circular hole subjected to prescribed loads.
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
91
Fig. 6. Mesh of a square plate with a small central hole/slit crack devised using a special element (occupying the central hatch area) and some isoparametric elements.
Table 1 Normalized stress rh =r0 at the hole in a square plate subjected to uniform tensile stress in the x direction by using the semi-analytictype special element r=A
0.01 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Exact
h (°) 0
15
30
45
60
75
90
1.0244 (2.4%) 1.0287 (2.9%) 1.0291 (1.9%) 1.0194 (1.9%) 1.0089 (0.9%) 0.9873 (1.3%) 0.9692 (3.0%) 0.9575 (4.3%) 0.9237 (7.0%) 0.8748 (12.5%) 1.0
0.7566 (3.3%) 0.7602 (3.8%) 0.7602 (3.8%) 0.7521 (2.7%) 0.7428 (1.5%) 0.7249 (1.0%) 0.7120 (2.7%) 0.6991 (4.5%) 0.6931 (5.3%) 0.6843 (6.5%) 0.7321
0.0249 (–) 0.0206 (–) 0.0257 (–) 0.0216 (–) 0.0151 (–) 0.0635 (–) 0.0047 (–) 0.0091 (–) 0.0366 (–) 0.1034 (–) 0
0.9747 (2.5%) 0.9756 (3.4%) 0.9778 (2.2%) 0.9767 (3.3%) 0.9802 (2.0%) 0.9794 (2.1%) 0.9724 (2.7%) 0.9599 (4.0%) 0.9227 (7.8%) 0.8327 (16.7%) 1.0
1.9743 (1.3%) 1.9777 (1.1%) 1.9813 (0.9%) 1.9754 (1.7%) 1.9772 (1.1%) 1.9699 (1.5%) 1.9619 (2.0%) 1.9594 (2.0%) 1.9539 (2.3%) 1.9331 (3.3%) 2.0
2.7060 (1.0%) 2.7113 (0.6%) 2.7160 (0.6%) 2.7067 (0.9%) 2.7078 (0.9%) 2.6977 (1.3%) 2.6938 (1.4%) 2.7097 (0.8%) 2.7535 (0.8%) 2.8417 (4.0%) 2.7321
2.9738 (0.8%) 2.9799 (0.7%) 2.9849 (0.5%) 2.9744 (0.9%) 2.9752 (0.8%) 2.9642 (1.2%) 2.9624 (1.3%) 2.9870 (0.4%) 3.054 (1.8%) 3.1941 (6.5%) 3.0
7.2. Cracks Fig. 3 shows a square plate with a crack of inclined angle b ¼ 45° subjected to uniform tensile stress in the y direction. The model devised for the finite element analysis by using the semi-analytic-
type special element is shown in Fig. 7. Note that the proposed element was employed to model the region surrounding the crack, and eight-node isoparametric elements were used to model the remaining region of the plate. Eleven analyses were performed, by varying the ratio of the crack length
92
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
Table 2 Normalized stress rh =r0 at the hole in a square plate subjected to uniform tensile stresses in both the x and y directions by using the semi-analytic-type special element r=A
0.01 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Exact
h (°) 0
15
30
45
60
75
90
1.9495 (2.5%) 1.9512 (2.4%) 1.9560 (2.2%) 1.9555 (2.2%) 1.9673 (1.1%) 1.9784 (1.1%) 1.9954 (0.2%) 2.0425 (2.1%) 2.1340 (6.7%) 2.3240 (16.2%) 2.0
1.9495 (2.5%) 1.9512 (2.4%) 1.9558 (2.2%) 1.9548 (2.3%) 1.9654 (1.7%) 1.9732 (1.3%) 1.9824 (0.9%) 2.0114 (0.6%) 2.0614 (3.1%) 2.1590 (7.5%) 2.0
1.9495 (2.5%) 1.9511 (2.4%) 1.9555 (2.2%) 1.9536 (2.3%) 1.9616 (1.9%) 1.9629 (1.9%) 1.9562 (2.2%) 1.9488 (2.6%) 1.9156 (4.2%) 1.8281 (8.9%) 2.0
1.9495 (2.5%) 1.9510 (2.5%) 1.9553 (2.2%) 1.9530 (2.4%) 1.9598 (2.0%) 1.9577 (2.1%) 1.9431 (2.8%) 1.9175 (4.1%) 1.8427 (7.9%) 1.6625 (16.9%) 2.0
1.9495 (2.5%) 1.9511 (2.4%) 1.9555 (2.2%) 1.9536 (2.3%) 1.9616 (1.9%) 1.9629 (1.9%) 1.9562 (2.2%) 1.9488 (2.6%) 1.9156 (4.2%) 1.8281 (8.9%) 2.0
1.9495 (2.5%) 1.9512 (2.4%) 1.9558 (2.2%) 1.9548 (2.3%) 1.9654 (1.7%) 1.9732 (1.3%) 1.9824 (0.9%) 2.0114 (0.6%) 2.0614 (3.1%) 2.1590 (7.5%) 2.0
1.9495 (2.5%) 1.9512 (2.4%) 1.9560 (2.2%) 1.9555 (2.4%) 1.9673 (1.1%) 1.9784 (1.1%) 1.9954 (0.2%) 2.0425 (2.1%) 2.1340 (6.7%) 2.3240 (16.2%) 2.0
Table 3 Normalized stress rh =r0 at the hole in a square plate subjected to uniform tensile and compressive stresses in the x and y directions by using the semi-analytic-type special element, respectively r=A
0.01 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Exact
h (°) 0
15
30
45
60
75
90
4.0003 (0.0%) 4.0081 (0.2%) 4.0123 (0.3%) 3.9905 (0.2%) 3.9785 (0.5%) 3.9456 (1.4%) 3.9290 (1.8%) 3.9210 (2.0%) 3.9607 (1.0%) 4.047 (1.2%) 4.0
3.4643 (0.0%) 3.4711 (0.2%) 3.4749 (0.3%) 3.4562 (0.2%) 3.4461 (0.5%) 3.4162 (1.4%) 3.3972 (1.9%) 3.3978 (1.9%) 3.4327 (0.9%) 3.5081 (1.3%) 3.4641
2.0001 (0.0%) 2.0040 (0.2%) 2.0062 (0.3%) 1.9955 (0.2%) 1.9898 (0.5%) 1.9727 (1.4%) 1.9619 (1.9%) 1.9624 (1.9%) 1.9828 (0.9%) 2.0264 (1.3%) 2.0
0 (–) 0 (–) 0 (–) 0 (–) 0 (–) 0 (–) 0 (–) 0 (–) 0 (–) 0 (–) 0.0
2.0001 (0.0%) 2.0040 (0.2%) 2.0062 (0.3%) 1.9955 (0.2%) 1.9898 (0.5%) 1.9727 (1.4%) 1.9619 (1.9%) 1.9624 (1.9%) 1.9828 (0.9%) 2.0264 (1.3%) 2.0
3.4643 (0.0%) 3.4711 (0.2%) 3.4749 (0.3%) 3.4562 (0.2%) 3.4461 (0.5%) 3.4162 (1.4%) 3.3972 (1.9%) 3.3978 (1.9%) 3.4327 (0.9%) 3.5081 (1.3%) 3.4641
4.0003 (0.0%) 4.0081 (0.2%) 4.0123 (0.3%) 3.9905 (0.2%) 3.9785 (0.5%) 3.9456 (1.4%) 3.9290 (1.8%) 3.9210 (2.0%) 3.9607 (1.0%) 4.047 (1.2%) 4.0
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
93
Fig. 7. Mesh of a square plate with an inclined crack (b ¼ 45°) devised using a semi-analytic-type element and some isoparametric elements.
Fig. 8. Mesh of a square plate with two collinear cracks devised using two semi-analytic-type elements and some isoparametric elements.
Table 4 Stress intensity factor KI and KII at the crack tip (z ¼ a) of the internal crack with length 2a and inclined angle b ¼ 45° in a square plate subjected to uniform tensile stress in the y direction by using the semi-analytic-type special element
Fig. 4(a) and (b) show a square plate containing a single crack and two collinear cracks, respectively. The model devised for the finite element analyses by using the semi-analytic-type special element is shown in Fig. 8. Note that for both the single and double cracks, the proposed element was again employed to model only the region
a/A
0.01 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Numerical solutions
Theoretical solutions
Error
KI
KII
KI
KII
KI
KII
0.0909 0.2035 0.2881 0.4084 0.5029 0.5842 0.6603 0.7343 0.8094 0.8879 0.9699
0.0907 0.2034 0.2880 0.4083 0.5028 0.5843 0.6603 0.7341 0.8092 0.8876 0.9701
0.0886 0.1982 0.2803 0.3964 0.4854 0.5605 0.6266 0.6865 0.7415 0.7927 0.8408
0.0886 0.1982 0.2803 0.3964 0.4854 0.5605 0.6266 0.6865 0.7415 0.7927 0.8408
2.60% 2.67% 2.78% 3.03% 3.61% 4.23% 5.38% 6.96% 9.16% 12.01% 15.35%
2.37% 2.62% 2.75% 3.00% 3.58% 4.25% 5.38% 6.93% 9.13% 11.97% 15.38%
to the side length of the proposed element, a=A, for the load case studied. Table 4 presents the stress intensity factor KI and KII at the crack tip (z ¼ a) of the inclined crack. It is obvious that the maximum discrepancy between the proposed and theoretical results is less than 5.4% for a=A 6 0:5 However, the said discrepancy increases rapidly with the increase of a=A.
Table 5 Stress intensity factor KI and KII at the crack tip (z ¼ a) of the internal crack with length 2a in a square plate subjected to uniform tensile stress in the y direction by using the semi-analytic-type special element a/A
0.01 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Numerical solutions
Exact solutions
Error
KI
KII
KI
KII
KI
KII
0.1820 0.4071 0.5761 0.8167 1.0047 1.1684 1.3204 1.4685 1.6188 1.7757 1.9397
1.295e5 4.258e4 2.385e3 1.285e2 3.258e2 5.993e2 9.203e2 1.250e1 1.572e1 1.860e1 2.102e1
0.1772 0.3963 0.5605 0.7927 0.9709 1.1212 1.2536 1.3734 1.4836 1.5862 1.6827
0 0 0 0 0 0 0 0 0 0 0
2.70% 2.73% 2.78% 3.03% 3.48% 4.21% 5.33% 6.92% 9.11% 11.95% 15.27%
– – – – – – – – – – –
94
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
Table 6 Stress intensity factor KI and KII at the crack tip (z ¼ a) of the two collinear cracks, each of length 2a, in a square plate (W =A ¼ 4) subject to uniform tensile stress in the y direction by using the semi-analytic-type special element a/A
Numerical solutions Crack 1
0.01 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Error Exact solutions
Crack 1
KI
KII
Crack 2 KI
KII
KI
KII
KI
KII
KI
KII
0.1818 0.4066 0.5754 0.8155 1.0028 1.1657 1.3166 1.4701 1.6102 1.7605 1.9128
1.255e5 4.771e4 2.745e3 1.481e2 3.759e2 6.916e2 1.063e1 1.258e1 1.839e1 2.189e1 2.494e1
0.1818 0.4066 0.5754 0.8155 1.0028 1.1657 1.3166 1.4701 1.6102 1.7605 1.9128
1.326e5 4.753e4 2.749e3 1.481e2 3.759e2 6.920e2 1.064e1 1.258e1 1.840e1 2.190e1 2.496e1
0.1772 0.3964 0.5606 0.7935 0.9731 1.1256 1.2615 1.3859 1.5021 1.6123 1.7180
0 0 0 0 0 0 0 0 0 0 0
2.60% 2.57% 2.64% 2.77% 3.05% 3.56% 4.37% 6.08% 7.20% 9.19% 11.3%
– – – – – – – – – – –
2.60% 2.57% 2.64% 2.77% 3.05% 3.56% 4.37% 6.08% 7.20% 9.19% 11.3%
– – – – – – – – – – –
Table 7 Normalized stresses rh =r0 along the free boundary of circular hole in a square plate subjected to uniform tensile stress in x direction with different type special elements r=A
h (°) 0
0.01 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Exact
Table 8 Stress Intensity factor KI and KII at the crack tip (z ¼ a) of the internal crack with length 2a in a square plate subjected to uniform tensile stress in y direction with different type special elements a=A
90
Semi-analytic type
Hybrid type
Semi-analytic type
Hybrid type
1.0244 (2.4%) 1.0287 (2.9%) 1.0291 (1.9%) 1.0194 (1.9%) 1.0089 (0.9%) 0.9873 (1.3%) 0.9692 (3.0%) 0.9575 (4.3%) 0.9237 (7.0%) 0.8748 (12.5%) 1.0000
1.0000 (0%) 1.0000 (0%) 0.9995 (0.05%) 0.9992 (0.08%) 0.9987 (0.13%) 0.9985 (0.15%) 0.9972 (0.28%) 0.9955 (0.45%) 0.9934 (0.66%) 0.9842 (1.6%) 1.0000
2.9738 (0.8%) 2.9799 (0.7%) 2.9849 (0.5%) 2.9744 (0.9%) 2.9752 (0.8%) 2.9642 (1.2%) 2.9624 (1.3%) 2.9870 (0.4%) 3.054 (1.8%) 3.1941 (6.5%) 3.0000
3.0000 (0%) 3.0000 (0%) 3.0003 (0.01%) 3.0013 (0.03%) 3.0030 (0.10%) 3.0052 (0.14%) 3.0060 (0.20%) 3.0105 (0.35%) 3.0150 (0.5%) 3.0270 (0.9%) 3.0000
surrounding the crack, and eight-node isoparametric elements were used to model the remaining
Crack 2
0.01 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Semi-analytic type
Hybrid type
KI
KII
KI
KII
KI
KII
0.1820 (2.70%) 0.5761 (2.78%) 0.8169 (3.05%) 1.0051 (3.53%) 1.1691 (4.29%) 1.3215 (5.44%) 1.4701 (7.08%) 1.6209 (9.31%) 1.7780 (12.16%) 1.9419 (15.49%)
1.295e5
0.1772 (0%) 0.5606 (0.01%) 0.7929 (0.02%) 0.9716 (0.08%) 1.1222 (0.12%) 1.2556 (0.2%) 1.3777 (0.35%) 1.4948 (0.8%) 1.6010 (1.0%) 1.7067 (1.5%)
0
0.1772
0
0
0.5605
0
0
0.7927
0
0
0.9708
0
0
1.1210
0
0
1.2533
0
0
1.3729
0
0
1.4829
0
0
1.5853
0
0
1.6815
0
2.385e3 1.285e2 3.258e2 5.993e2 9.203e2 1.257e1 1.583e1 1.875e1 2.121e1
Exact solutions
region of the plate. Tables 5 and 6 present the stress intensity factor KI and KII at the crack tip for the case of single and double cracks, respectively.
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
95
Fig. 9. Comparison of the analytical and numerical results of the normalized KI for the case of a crack aligned along the centerline of a hole.
It is obvious that the proposed and theoretical results are in good agreement. 7.3. Comparison of semi-analytic and hybrid-element types From Tables 7 and 8, it is obvious that the results obtained using the hybrid-type elements are in significantly better agreement with the theoretical solutions than those obtained using the semianalytic-type elements. This is due to the fact that
the conditions of displacement continuity along the inter-element boundaries have been implemented in the former element types. Note that the hybrid-type elements are to be employed together with four-node isoparametric elements in devising a finite element model. This is to ensure the same displacement distributions at the interface between the special and the traditional isoparametric elements. As for the semianalytic-type elements, these elements can be used together with the conventional elements of the
96
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
Fig. 10. Comparison of the analytical and numerical results of the normalized KI for the case of a crack situated at a distance t from two holes.
same number of nodes and provide reasonably good accuracy. Fig. 6(a) and (b) show the finite element meshes devised using 1 semi-analytic-type element together with 24 eight-node isoparametric elements and 1 hybrid-type element together with 32 four-node isoparametric elements, respectively.
7.4. Interactions of voids and cracks Figs. 9–11 provide some examples for the modeling of the interaction between two circu-
lar holes and a crack using the high precision hybrid-type elements. Fig. 9 presents the analytical and numerical results of the normalized KI , i.e., pffiffiffi KI =ðr0 aÞ, of a crack, aligned along the centerline of a hole at a distance t apart, embedded in an infinite plate subjected to uniaxial tension r0 . The results obtained using the proposed hybrid-type element are in excellent agreement with the analytic results obtained in [15]. Fig. 10 shows an infinite plate containing a crack and two holes subjected to uniaxial loading r0 . The results obtained using the proposed hybrid-type element are again in excellent agreement with the analytic re-
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
97
Fig. 11. Comparison of the analytical and numerical results of the normalized KI for the case of a crack separating two holes.
sults obtained by Hu et al. [15]. This excellent agreement is also clearly shown in the case of an infinite plate, which contains a crack separating two holes, subject to uniaxial loading r0 as shown in Fig. 11. Fig. 12(a)–(c) show an infinite plate containing a semi-infinite macro-crack and a micro-defect, i.e., circular hole, elliptic hole or crack, respectively, subjected to uniaxial loading r0 in the y direction. The finite element mesh devised for analysis is shown in Fig. 12(d). Note that the hy-
brid element containing an edge crack developed in [21] was employed to model the region surrounding the tip of the main crack. Whereas, the proposed hybrid-type special elements with internal defects were employed to model the region surrounding the micro-defects. With reference to Figs. 13–17, it is obvious that the results obtained using the proposed hybrid-type elements are again in excellent agreement with the corresponding analytical results obtained in [16,17]. Obviously, the existence of micro-defects has significant effects
98
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
Fig. 12. A large plate with a long edge macro-crack and a micro-defect (in the shaded zone) and the FE model devised.
Fig. 13. Comparison of the analytical and numerical results of the normalized KI for the case of a macro-crack and a micro-defect.
on the macro-crack. As a result, the original mode I crack problem became a mixed-mode crack problem and the initial crack growth direction of the macro-crack was also changed.
8. Conclusions Eight-node rectangular elements, each with a central circular or elliptical hole or an internal
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
99
Fig. 14. Comparison of the analytical and numerical results of the normalized KII for the case of a macro-crack and a micro-elliptichole.
Fig. 15. Comparison of the analytical and numerical results of the normalized KII for the case of a macro-crack and a micro-crack.
crack, have been developed for modeling multiple voids, cracks to study the interaction between these defects. The good accuracy of the semi-analytictype elements are achieved due to the exact defi-
nition of the stress and strain fields within the elements and the implementation of the free boundary conditions along the free boundary of the hole or crack. However, these semi-analytic-type elements
100
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
Fig. 16. Comparison of the analytical and numerical results of the initial crack growth direction d0 of the macro-crack for the case of a macro-crack and a micro-elliptic-hole.
Fig. 17. Comparison of the analytical and numerical results of the initial crack growth direction d0 of the macro-crack for the case of a macro-crack and a micro-crack.
have been improved by implementing displacement continuity at the inter-element boundary to establish the hybrid-type elements. The proposed ele-
ments can be easily combined with the conventional elements, such as isoparametric elements, without any modifications, to obtain good results
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
using simple finite element meshes. It is interesting to note that both the semi-analytic-type elements and hybrid-type elements are particularly useful for modeling the interaction between many voids and cracks. The former element type is less accurate compared with the latter but much easier to construct and use. Acknowledgements Supports from the Research Grants Council of the Hong Kong Special Administrative Region, China (Project No. HKU 7081/00E) is acknowledged. The authors would like to thank Professor Daining Fang of Tsinghua University for the useful discussions. References [1] M. Isida, D.H. Chen, H. Nisitani, Plane problems of an arbitrary array of crack emanating from the edge of an elliptical hole, Eng. Fract. Mech. 21 (1985) 983–995. [2] M. Isida, H. Igawa, Analysis of a zig-zag array of circular holes in an infinite solid under uniaxial tension, Int. J. Solids Struct. 27 (1991) 227–258. [3] M. Isida, S. Nemat-Nasser, A unified analysis of various problems relating to circular holes with edge cracks, Eng. Fract. Mech. 27 (1987) 849–864. [4] G.C. Sih, Handbook of stress intensity factor, Lehigh University, Bethlehem, PA, 1973. [5] C. Atkinson, The interaction between a crack and an inclusion, Int. J. Eng. Sci. 10 (1972) 127–136. [6] F. Edorgan, G.D. Gupta, M. Ratawani, Interaction between a circular inclusion and an arbitrary oriented crack, J. Appl. Mech. 41 (1974) 1007–1013. [7] Y.Z. Chen, General case of multiple crack problems in an infinite plate, Eng. Fract. Mech. 20 (1984) 591–597. [8] Y.Z. Chen, Solutions of multiple crack problems of a circular plate or an infinite plate containing a circular hole by using Fredholm integral equation approach, Int. J. Fract. 25 (1984) 155–168. [9] Y.Z. Chen, A special boundary-element formulation for multiple-circular-hole problems in an infinite plate, Comput. Meth. Appl. Eng. Mech. 50 (1985) 263–273. [10] H. Tada, P. Paris, G. Irwin, The Stress Analysis of Cracks Handbook, Paris Production Inc, St. Louis, MO, 1985. [11] M. Isida, Analysis of stress intensity factors for plates containing random array of cracks, Bull. Jpn. Soc. Mech. Eng. 13 (1970) 635–642. [12] M. Isida, On the determination of stress intensity factors for some common structural problems, Eng. Fract. Mech. 2 (1970) 61–79.
101
[13] M. Isida, Methods of laurent series expansion for internal crack problems, in: G.C. Sih (Ed.), Methods of Analysis and Solutions of Crack Problems, Noordhoff, Leyden, 1973, pp. 56–130. [14] G.J. Rodin, Y.-L. Hwang, On the problem of linear elasticity for an infinite region containing an infinite number of non-intersecting spherical inhomogeneities, Int. J. Solid Struct. 27 (1991) 145–159. [15] K.X. Hu, A. Chandra, Y. Huang, Multiple void-crack interaction, Int. J. Solids Struct. 30 (11) (1993) 1473–1489. [16] S.X. Gong, H. Horii, General solution to the problem of microcracks near the tip of a main crack, J. Mech. Phys. Solids 37 (1) (1989) 27–46. [17] S.X. Gong, S.A. Meguid, Microdefect interacting with a main crack: a general treatment, Int. J. Mech. Sci. 34 (12) (1992) 933–945. [18] S.A. Meguid, P.E. Gaultier, S.X. Gong, A comparison between analytical and finite element analysis of main crack-microcrack interaction, Eng. Fract. Mech. 38 (6) (1991) 451–465. [19] A.K. Soh, An improved method for determining free boundary stresses, J. Strain Anal. 27 (2) (1992) 93–99. [20] A.K. Soh, Development of special finite elements for free boundaries, Int. J. Solids Struct. 36 (6) (1999) 899–917. [21] P. Tong, T.H.H. Pian, S.J. Lasry, A hybrid-element approach to crack problems in plane elasticity, Int. J. Numer. Meth. Eng. 7 (1973) 297–308. [22] R. Piltner, Special finite elements with holes and internal cracks, Int. J. Numer. Meth. Eng. 21 (1985) 1471–1485. [23] A.K. Soh, Z.F. Long, A high precision element with a central circular hole, Int. J. Solids Struct. 36 (35) (1999) 5485–5497. [24] P. Tong, New displacement hybrid finite element models for solid continua, Int. J. Numer. Meth. Eng. 2 (1970) 73–83. [25] N.I. Muskhelishvili, Some Basic Problems of the Mathematical Theory of Elasticity, fourth ed., Noordhoff International Publishing, Leyden, Netherlands, 1954. [26] T.H.H. Pian, Derivation of element stiffness matrices by assumed stress distribution, AIAA J. 2 (1964) 1333–1336. [27] G.C. Sih, P.C. Paris, F. Erdogan, Crack-tip stress intensity factors for plane extension and plate bending problems, J. Appl. Mech. 29 (1962) 306–312. [28] S.A. Meguid, in: Engineering Fracture Mechanics, Elsevier Science Publisher LTD, Great Britain, 1989, pp. 134–136. [29] G.C. Sih, A special theory of crack propagation, in: G.C. Sih (Ed.), Mechanics of fracture 1: Methods of analysis and solutions of crack problems, Noordhoff Internatioinal Publishing, Leyden, Netherlands, 1973, pp. XXI–XLV. [30] G.C. Sih, Some basic problems in fracture mechanics and new concepts, Eng. Fract. Mech. 5 (1973) 365–377. [31] G.C. Sih, Strain-energy-density factor applied to mixed mode crack problem, Int. J. Fract. 10 (3) (1975) 305–321. [32] G.C. Sih, A three-dimensional strain density factor theory of crack propagation, in: M.K. Kassir, G.C. Sih (Eds.), Mechanics of fracture 2: Three dimensional crack problems, Noordhoff International Publishing, Leyden, Netherlands, 1975, pp. XV–LIII.
Modeling of voids/cracks and their interactions C.H. Yang, A.K. Soh
*
Department of Mechanical Engineering, The University of Hong Kong, Pokfulam Road, Hong Kong, China
Abstract Numerical modeling of a two-dimensional elastic body containing multiple voids/cracks to study the interaction between these defects can be significantly simplified by developing special finite elements, each containing an internal circular/elliptic hole or a slit crack. These finite elements are developed using complex potentials and the conformal mapping technique. The elements developed can be divided into two categories, namely, the semi-analytic-type and hybrid-type elements. The latter element type is an improved version of the former due to the implementation of displacement continuity along the inter-element boundary. All the proposed elements can be easily combined with the conventional displacement elements, such as isoparametric elements, to analyze the above-mentioned problems without using complicated finite element meshes. Numerical examples have been employed to illustrate the modeling of voids/ cracks and their interactions. The results obtained using the semi-analytic-type elements are in good agreement with the theoretical results, and the corresponding results obtained using the hybrid-type elements show an improvement of the agreement with the theoretical results. However, the former element type is much easier to construct. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Circular/elliptical hole; Crack; Semi-analytic element; Hybrid element; Complex potentials; Conformal mapping
1. Introduction For brittle/quasi-brittle materials, the existence of multiple holes, or voids, and cracks or other forms significantly affects their mechanical behavior due to high stress concentrations. Thus, the ability to perform accurate analysis of a general two-dimensional body containing the above-mentioned defects is very important. To-date, many researchers have studied the local behavior of such defects using theoretical approaches [1–18]. However, all the theoretical solutions obtained are for
*
Corresponding author. Fax: +852-2858-5415. E-mail address: [email protected] (A.K. Soh).
relatively simple cases, such as a crack among several regularly spaced voids or a void among several regularly spaced cracks embedded in an infinite plate subjected to remote loading. These theoretical approaches cannot be applied to solve the problems of practical structures in which random distributions of many voids and cracks exist. Meanwhile, the finite element method has also been widely employed to solve such problems on practical structures. However, the conventional elements cannot ensure high accuracy in the vicinity of holes or cracks even if a very large number of elements are used. This is because the assumed displacement functions do not satisfy the prevailing conditions at the free boundary. Proposed in [19,20] are some modified methods to solve the
0167-8442/02/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 1 6 7 - 8 4 4 2 ( 0 2 ) 0 0 0 8 3 - 6
82
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
said problem. But the main weakness of these methods is that many elements are required to model the vicinity of the holes or cracks. Another approach to solve such problem is to construct special elements containing internal holes or cracks. Tong et al. [21] applied the hybrid-element formulation and complex variable technique to construct a special super-element for the analysis of the elastic stress intensity factors for the case of a plate containing an edge crack. Constructed in [22] are some Treffz-type finite elements with internal circular/elliptic holes or cracks based on two different variational formulae and applied Trefftz method for the use of trial functions, by employing the complex potentials and the conformal mapping technique. Two-dimensional semi-analytic finite elements [23] are constructed for a circular hole using complex potentials to define stress and strain functions inside the element directly. In the formulation of these elements, the boundary conditions at the free boundary of the hole are implemented explicitly. The main shortcoming of the formulation in [23] is the possible discontinuity along the inter-element boundaries. Such problems can be resolved by employing the hybrid displacement method based on the modified minimum potential energy principle [24]. It is worth noting that the usage of the above-mentioned elements leads to the use of only one element for the vicinity of hole/crack, and high accuracy is achieved in the stress concentrated region. In the present paper, the above-mentioned special elements containing an internal circular/elliptical hole or an internal crack will be developed using complex potentials and the conformal mapping technique. The above-mentioned hybrid displacement method will then be employed to solve the problem of discontinuity along inter-element boundaries. The versatility and accuracy of the special finite elements developed will be demonstrated by analyzing some two-dimensional elastic bodies containing multiple voids and cracks.
2. Basic theory and formulas For linear elastic analysis of a three-dimensional solid, the governing equations are as follows:
(a) Equilibrium conditions rij;j þ F i ¼ 0
ð1Þ
in V
(b) Stress–strain relations rij ¼ Dijkl ekl
or
eij ¼ Cijkl rkl
in V
ð2Þ
(c) Strain–displacement relations eij ¼ 12ðui;j þ uj;i Þ
in V
ð3Þ
(d) Displacement boundary conditions ui ¼ ui
on Su
ð4Þ
(e) Traction boundary conditions Ti ¼ T i
on Sr
ð5Þ
In the above equations, rij , eij , ui are the stress, strain and displacement components, respectively; F i , T i , ui are the prescribed body force and boundary traction and displacement components, respectively; Dijkl , Cijkl are the elastic and compliance constants, respectively; V is the volume of the solid; and Sr , Su are the portion of the boundary on which surface tractions and boundary displacements are prescribed, respectively. To consider a two-dimensional plate with holes or cracks in plane elasticity, the fields of stress and displacement are best expressed in terms of the complex potentials proposed by Muskhelishvili [25] as follows: n o rx þ ry ¼ 2 w0 ðzÞ þ w0 ðzÞ ð6Þ n o rx ry þ 2isxy ¼ 2 zw00 ðzÞ þ v0 ðzÞ
ð7Þ
2lðux þ iuy Þ ¼ jwðzÞ zw0 ðzÞ vðzÞ
ð8Þ
h iB Fx þ iFy ¼ i wðzÞ þ zw0 ðzÞ þ vðzÞ A h i 0 ¼ i wðzB Þ þ zB w ðzB Þ þ vðzB Þ CA
ð9Þ
where z ¼ x þ iy; j ¼ 3 4m for plane strain and j ¼ ð3 mÞ=ð1 þ mÞ for plane stress; A, B are the starting and end points of integration respectively; CA is a complex constant at the datum point A; (0 ) and (00 ) denote the first and second orders of complex differentiation with respect to z, respectively. As the treatment of boundary conditions in
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
the original domain usually fails, it is convenient and advantageous to make use of conformal mapping, for which the mapping function is given as follows: m z ¼ f ð1Þ ¼ c 1 þ ð10Þ 1 where 1 ¼ n þ ig and its inverse function is 1 ¼ f 1 ðzÞ. By means of conformal mapping, the original domains with holes or cracks are mapped onto the alternative boundaries as unit circles. Note that the advantage of such mapping is that an important boundary of an element can be mapped to a boundary that is geometrically simpler, such as a unit circle. In the case of an elliptical hole with major axis a and minor axis b, the arbitrary constants of Eq. (10) are given by c ¼ ða þ bÞ=2 and m ¼ ða bÞ=ða þ bÞ. Obviously, the circle and the crack can be considered as special cases by letting a ¼ b ¼ r for the circle with radius r, and b ¼ 0 for the crack with length 2a. The inverse mapping function for Eq. (10) is pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
1 1¼ z z2 4c2 m 2c
ð11Þ
To map the exterior of an ellipse onto the exterior of the unit circle, the sign of the root in Eq. (11) is chosen such that j1j P 1:0. In using conformal mapping, the complex functions are not assumed in the original domain but in the transformed element domain. Therefore, the stress and displacement fields are expressed in terms of f ð1Þ and the new functions wð1Þ and vð1Þ, chosen in the 1-plane, as follows: ( rx þ ry ¼ 2
0
0
w ð1Þ w ð1Þ þ f 0 ð1Þ f 0 ð1Þ
Fx þ iFy ¼ i wð1B0 Þ þ f ð1B0 Þ
ð12Þ
rx ry þ 2isxy ¼ ( ! ) w00 ð1Þ f 00 ð1Þ v0 ð1Þ 0 2 f ð 1Þ w ð1Þ 03 þ 0 f 02 ð1Þ f ð1Þ f ð1Þ
f 0 ð1Þ
f 0 ð1B0 Þ
# þ vð1B0 Þ CA0 ð15Þ
where A0 , B0 correspond to the original integration points A and B; and CA0 is also a complex constant; (0 ) and (00 ) denote the first and second orders of differentiation with respect to 1, respectively. Due to the simple boundary of the unit circle, the polar coordinates (q; h) are introduced in the transformed domain, and Eqs. (12)–(14) can be rewritten as follows: ) ( w0 ð1Þ w0 ð1Þ rq þ rh ¼ 2 ð16Þ þ f 0 ð1Þ f 0 ð1Þ rh rq þ 2isqh ( ¼ 2e2ia
w00 ð1Þ f 00 ð1Þ 0 f ð1Þ ð1Þ w f 02 ð1Þ f 03 ð1Þ
!
v0 ð1Þ þ 0 f ð1Þ
vð1Þ
ð13Þ
ð14Þ
)
ð17Þ 2lðuq þ iuh Þ ¼ e
ia
jwð1Þ f ð1Þ
w0 ð1Þ f 0 ð1Þ
! vð1Þ ð18Þ
where a is the angle between the axes x and 1, eia ¼
1 f 0 ð1Þ q jf 0 ð1Þj
and e2ia ¼
12 f 0 ð1Þ q2 f 0 ð1Þ
For the free boundary of a hole of unit circle in a plate, boundary conditions to be satisfied at 1 ¼ 10 ¼ eih or j1j ¼ 1:0 are ð19Þ
and Fx ¼ Fy ¼ 0
w0 ð1Þ
w0 ð1B0 Þ
rq ¼ sqh ¼ 0 in the 1-plane
)
2lðux þ iuy Þ ¼ jwð1Þ f ð1Þ
"
83
in the z-plane
From Eqs. (16), (17) and (19) ) ( ( w0 ð10 Þ w0 ð10 Þ w00 ð10 Þ 2ia f ð10 Þ e þ 0 f ð10 Þ f 0 ð10 Þ f 02 ð10 Þ ! ) f 00 ð10 Þ v0 ð10 Þ 0 w ð10 Þ 03 þ 0 ¼0 f ð10 Þ f ð10 Þ
ð20Þ
ð21Þ
84
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
Substituting Eq. (20) into Eq. (15) and letting CA0 ¼ 0, it is found that wð10 Þ þ f ð10 Þ
c v0 þ c 0 1 þ
w0 ð10 Þ
þ vð10 Þ ¼ 0
f 0 ð10 Þ
wð1Þ ¼ w0 þ
k
ak 1 þ
k¼1
lU X
vð1Þ ¼ v0 þ
cp 1p þ
p¼1
bl 1
ð23Þ
dq 1q
ð24Þ
q¼1
or in accordance with the expressions provided in [3]: w0z ð1Þ ¼ a0 þ
kU X
ak 1 k þ
k¼1
( 0
v ð1Þ ¼ c c0 þ
lU X
bl 1l
ð25Þ
l¼1
pU X
p
cp 1 þ
p¼1
qU X
) dq 1
q
ð26Þ
q¼1
where w0z ð1Þ is the first order complex derivative of wð1Þ with respect to z; kU , lU , pU , and qU are the upper limits of the series, which are related to the total number of nodes in the special element; w0 , v0 , ak , bl , cp , and dq are complex coefficients of the v0 , ak ¼ a^k þ i~ ak , forms w0 ¼ w^0 þ iw~0 , v0 ¼ v^0 þ i~ bl ¼ b^l þ ib~l , cp ¼ c^p þ i~ cp , and dq ¼ d^q þ id~q . Note that the complex constants w0 and v0 in Eqs. (23) and (24) are also used as the integral constants for Eqs. (25) and (26), respectively. Thus, from Eqs. (25) and (26), we obtain wð1Þ ¼ ( c w0 þ
a0 1 þ
kU X k¼1
m
a0 11 þ
lU X 1 1 ak 1kþ1 b 1ðl1Þ kþ1 l 1 l l¼2
kU X k¼2
q
)
The following conditions are to be introduced in Eqs. (27) and (28) to ensure that ux and uy are single-valued displacements: b1 ma1 ¼ 0 ð29Þ d1 ¼ 0 In addition, from Eqs. (22)–(24), we obtain
l
l¼1 qU X
U X 1 1 cp 1pþ1 dq 1ðq1Þ pþ1 q 1 q¼2
ð28Þ
w0 mb2 þ 2 a2 þ v0 m c2 ¼ 0 or v ¼ w þ mb 2a2 þ mc 0
pU X
pU X p¼1
ð22Þ
Thus, the boundary conditions are expressed in the 1-plane. The holomorphic functions wð1Þ and vð1Þ can be approximated, using Laurent series of 1 in the region where j1j P 1:0, as follows: kU X
vð1Þ ¼ (
!
U X 1 1 ak 1k1 bl 1ðlþ1Þ k1 l þ 1 l¼1
l
!)
ð30Þ
2
and by using Eqs. (25)–(28), we obtain w0 þ mb1 þ a1 þ v0 ¼ 0 or v0 ¼ w0 ð1 þ m2 Þa1
ð31Þ
From Eqs. (30) and (31), the relationship between the two integral constants v0 and w0 is established. This leads to a reduction of the number of unknown parameters and, thus, enables the acquisition of higher order terms in those complex series. The above-mentioned complex potentials w and v in the form of Laurent series of 1 will be used to deduce the detailed expressions for the displacement distribution matrix [U] and the stress distribution matrix [P]. The only difference from the special semi-analytic-type element and the special hybrid-type element is that all the coefficients related to the rigid body movement will be set to zero in prior in the latter.
3. Stiffness matrix of a semi-analytic-type element with a central circular/elliptic hole or an internal crack Fig. 1(a)–(c) show the typical eight-node rectangular elements containing a central circular and an elliptic hole and an internal crack, respectively. e The nodal displacement vector fqg is defined as fqge ¼ ½ ðux Þ1
ð27Þ
2
0
ðuy Þ1
ðux Þ2
ðuy Þ2
ðux Þn
ðuy Þn T ð32Þ
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
85
Fig. 1. Eight-node rectangular membrane element with a circular/elliptic hole or an internal crack.
where ðux Þi and ðuy Þi are the x and y displacement components of node i ði ¼ 1; . . . ; nÞ, respectively; and n is the total number of the nodes of the element. In order to express wð1Þ and vð1Þ in terms of e fqg , the upper limits in Eqs. (22) and (23) are chosen such that the total number of independent terms in wð1Þ and vð1Þ is approximately equal to that of the independent displacement components. Substituting Eqs. (23) and (24) or Eqs. (25)–(28) into Eq. (19) or (21), respectively, two sets of equations relating the complex coefficients are obtained. By back-substituting these equations into Eqs. (23) and (24) or Eqs. (25)–(28) and through the use of Eq. (14), the strain field within the element is determined. The displacement vector fqg at any point in the element can be expressed in terms of the above-mentioned independent coefficients as follows: ux fug ¼ ¼ ½Ufbg ð33Þ uy where fbg is a matrix consisting of the unknown independent coefficients in Eqs. (23) and (24) or Eqs. (25)–(28), ½U is the displacement distribution matrix given by ½U ¼ fUx ðq; hÞ Uy ðq; hÞgT It is worth noting that the complex potentials in Eqs. (23) and (24) provide very high accuracy in
dealing with circular holes, but they create complications in handling non-circular holes. However, Eqs. (25)–(28) provide relatively simple expressions for all the components in both fbg and ½U, and maintain good accuracy in all cases. Substituting the nodal displacements at all the eight nodal points of the element into Eq. (33), we obtain e
e
fqg ¼ ½U fbg
ð34Þ
where ½Ue ¼ f ðUx Þ1
ðUy Þ1
ðUx Þ2
ðUy Þ2
...
ðUx Þn
ðUy Þn gT
Further, 1
fbg ¼ ½Ue fqg
e
ð35Þ
Substituting Eq. (35) into Eq. (33) fug ¼
ux uy
1
e
¼ ½U½Ue fqg ¼ ½Nfqg 1
e
ð36Þ
where ½N 22n ¼ ½U½Ue is termed the elementshape function in the conventional finite element method. For the standard finite element approach, the strain vector feg at any point in the element can be obtained as follows:
86
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
8 9 2o e > = 6 ox < x> 0 ¼6 feg ¼ ey 4 > ; : > o exy oy 2 oU 3
0
3
4. Stiffness matrix of a hybrid-type element with a central circular/elliptic hole or an internal crack
o 7 7 oy 5fqg o ox
x
6 6 ¼ 2l6 4
ox oUy oy oUx oy
oUy ox
þ
7 7 e 1 e e 7½U fqg ¼ ½Bfqg 5
ð37Þ
One important version of the modified minimum potential energy principle for linear elastic analysis of a solid body, developed in [24], is given as, Pep
where ½B is the element strain–displacement matrix. Since the sub-matrices Ux ðq; hÞ and Uy ðq; hÞ are in terms of the coordinates ðq; hÞ, there results " oU # " oU # " oUy # " oUy #
¼
oq
ox oUx oy
¼ ½J
oUx oh
where " ox 2 6 ¼4
ox oh
oUy oh
¼ ½J
ox oUy oy
ð38Þ
oy oh
cm 1 qm2 cos h
cm 1 þ qm2 q sin h
From Eq. (38), we have " oUx # " oUx # ¼ ½J
3
cm 1 þ qm2 sin h 7
5 m cm 1 q2 q cos h
oq
2
and
oUx oh
oUy 4 ox oUy oy
3
2
5 ¼ ½J1 4
oUy oq oUy oh
3 5
ð39Þ where ½J1 ¼
"
oy oq
ð41aÞ
Pep ¼
Z ve
Z
Z 1 T T T e Ee F u dv T ~u ds 2 ser TT ð~u uÞ ds
ð41bÞ
ove
1
T~i ðui u~i Þ ds
or in matrix form,
is the Jacobian matrix related to the mapping function. Moreover, except at the point (0,0), ! " m2 m 2 jJj ¼ cm q 1 þ 4 2 2 cos 2h > 0 q q
ox oUx oy
Z
Z 1 Dijkl eij ekl F i ui dv T i u~i ds 2 ser
ove
#
oy oq
oq
½J ¼
oq
and
ve
x
x
Z
1 jJj ox oh
oy oq ox oh
#
" oq ¼
ox oq oy
oh ox
#
oh oy
The element stiffness matrix is given by Z e T ½Ks ¼ ½B ½D½BjJjt dA
ð40Þ
Ae
where ½D is the element elasticity matrix, t is the plate thickness and Ae is the area of an element.
where ve is the element volume; ser is the surface traction boundary; and ove is the boundary of the element volume. For displacement-type finite elements, the displacements are normally required to be continuous not only inside the element but also along the inter-element boundary to ensure the stability and precision of calculations. The displacement continuity condition along the inter-element boundary SC is uðCÞ ¼ ~uðCÞ
on SC
ð42Þ
The assumed linear or quadratic displacement function along the inter-element boundary can be expressed as h
i u~ðCÞ x ~uðCÞ ¼ ¼ ½Lfqge u~ðCÞ y
ð43Þ
where [L] is the interpolation function matrix only defined on SC [21,26]. In principle, each quantity in Eqs. (41a) and (41b) should be assumed individually but in fact the appropriate constructions of r, u and T can make Eqs. (1)–(5) automatically satisfied in advance. Further by considering the divergence the-
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
87
orem, the simplified form of the functional (41a) and (41b) used in the proposed element with internal hole/crack is Z Z $ ðCÞ %T $ ðCÞ %T ðCÞ 1 e ~ Pp ¼ T u ds T ð44Þ u ds 2 SC SC
5. Calculation of stresses and stress intensity factors
From Eqs. (12) and (13),
where frg ¼ f rx ry sxy g and ½S is the element stress–displacement matrix given by ½S ¼ ½D½B. In fact, there is a direct method to calculate the stresses once the independent parameter vector fbg for the complex potentials w and v have been determined, see also Eq. (45). With reference to Fig. 1, by assuming that b ¼ 0, the elliptic hole degenerates to a slit crack with length 2a. In this case, the mode I and II stress intensity factors (SIFs), KI and KII , respectively, at the crack tip need to be determined. The formulae for calculating the complex stress intensity factor K at the crack tip has been given in [27]. In general, the stress field around the crack tip for the mixedmode I and II fracture based on the coordinate systems (x0 ; y 0 ) and (r0 ; h0 ), refer to Fig. 2, is given as follows:
frg ¼ ½Pfbg
frg ¼ ½Dfeg ¼ ½Sfqg
where ½P is the stress distribution matrix and T ½P ¼ f Px ðq; hÞ Py ðq; hÞ Pxy ðq; hÞ g . The traction boundary conditions (5) on SC can be written as ( ) & ' TxðCÞ ðCÞ ½T ¼ ð46Þ ¼ RðCÞ fbg on SC ðCÞ Ty ^ðCÞ is a 2 3 unit matrix nðCÞ P, n where RðCÞ ¼ ^ outward normal to the common boundary SC . Alternatively, " # " # ðCÞ PðCÞ RðCÞ ðCÞ x nx þ Pxy ny x ðq; hÞ ½R ¼ ¼ ðCÞ RðCÞ PðCÞ y ðq; hÞ y ny þ Pxy nx By substituting the above relations, the hybrid functional becomes
KI h0 rx0 ¼ pffiffiffiffiffiffiffiffiffi cos 2 2pr0
KI h0 ry 0 ¼ pffiffiffiffiffiffiffiffiffi cos 2 2pr0 0
SC
Furthermore, its stationary value with respect to b is given by ð48Þ
Since the element strain energy Vs:e: is Vs:e: ¼ 12bT Hb ¼ 12qT GT H1 Gq
e
T
1
sx0 y 0
h0 3h0 2 þ cos cos 2 2 ! h0 3h0 1 þ sin sin 2 2 0
KII h h 3h þ pffiffiffiffiffiffiffiffiffi sin cos cos 2 2 2 2pr0 0 0 0 KI h h 3h ¼ pffiffiffiffiffiffiffiffiffi cos sin cos 2 2 2 2pr0 KII h0 þ pffiffiffiffiffiffiffiffiffi cos 2 2pr0
! !
ð52Þ
0
h0 3h0 1 sin sin 2 2
!
Moreover, ð49Þ
Therefore, the element stiffness matrix ½Keh can be expressed as ½Kh ¼ ½G ½H ½G
ð51Þ
h0 3h0 1 sin sin 2 2
KII h0 pffiffiffiffiffiffiffiffiffi sin 2 2pr0
ð47Þ
where Z Z 1 T H¼ R U ds ðRT U þ UT RÞ ds 2 SC SC Z G¼ RT L ds
oPev ¼ 0 ) b ¼ H1 Gq ob
e
T
ð45Þ
Pep ¼ 12bT Hb bT Gq
The stress vector at a point in the element can be determined from
ð50Þ
KI h0 KII h0 rx0 þ ry 0 ¼ 2 pffiffiffiffiffiffiffiffiffi cos pffiffiffiffiffiffiffiffiffi sin 2 2 2pr0 2pr0
! ð53Þ
Eq. (53) can be rewritten in complex form as follows:
88
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
Thus, if the complex coefficients of the complex potential wð1Þ can be determined, the SIFs can be obtained from Eq. (56). Some theoretical solutions for different crack problems are provided below for easy reference. For the mixed-mode crack problem shown in Fig. 3(a), the complex potential wð1Þ can also be derived theoretically using complex variables and expressed as follows: wð1Þ ¼ Fig. 2. The local coordinate system at the crack tip.
( ! rx þ ry ¼ 2Re K
1 2pðz z0 Þ
"1=2 )
¼ 4Re½w0 ðzÞ ð54Þ
where K ¼ KI iKII and z0 ¼ a. Note that the points z ¼ a in the complex plane are at 1 ¼ 1 in the mapped plane. Eq. (54) gives ! "1=2 1 0 K ¼ 2w ðzÞ ð55Þ 2pðz z0 Þ Let z ¼ z0 or 1 ! 1, pffiffiffiffiffiffi 1=2 K ¼ 2 2p lim ðz z0 Þ w0 ðzÞ ¼ 2 z!z0
rffiffiffi p 0 w ð1Þ1¼1 a ð56Þ
r0 að1 e2ib Þ 41
ð57Þ
where b is the angle between the centerline of the crack and the axis of the applied load. Substituting Eq. (57) into Eq. (56), the result is % % 1 pffiffiffiffiffiffi$ 1 pffiffiffiffiffiffi$ K ¼ r0 pa 1 e2ib 12 1¼1 ¼ r0 pa 1 e2ib 2 2 ð58Þ Therefore, pffiffiffiffiffiffi KI ¼ r0 pa sin2 b and
pffiffiffiffiffiffi KII ¼ r0 pa sin b cos b ð59Þ
A more general expression for the stress intensity factor of mode I crack in a finite elastic body is [28] pffiffiffiffiffiffi KI ¼ r0 paf ðaW Þ ð60Þ where r0 is the applied stress in the x or y direction and f ðaW Þ is a configuration factor which accounts for proximity effects of boundary surfaces.
Fig. 3. A square plate with an inclined crack or a single edge crack subjected to tension in one direction.
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
(a)
89
(b)
Fig. 4. A square plate with single crack or two collinear cracks subjected to tension in one direction.
For the case of a wide plate containing a long edge crack, as shown in Fig. 3(b), the opening mode stress intensity factor KI is given by 2 W W f ðaW Þ ¼ 1:12 0:231 þ 10:55 a a 3 4 W W 21:72 þ 30:39 ð61Þ a a where W is the width of the plate. The above factor is highly accurate (0.5% error) for W =a < 0:6. For the case of a wide plate containing a crack or two collinear cracks, as depicted in Fig. 4(a) and (b), respectively, the opening mode stress intensity factor KI is given by 1=2 W pa tan f ðaW Þ ¼ ð62Þ pa W where W is either the width of the plate or the distance between the centerlines of two collinear cracks.
micro-defects, the stress field in the vicinity of a macro-crack is greatly disturbed such that the original mode I crack problem becomes a mixedmode crack problem. For mixed-mode crack problems, the initial crack growth direction is a very important parameter that has to be determined. Proposed in [29–32] is a way to solve mixedmode fracture problems using the concept of strain energy density, which can be used to predict the crack growth direction. The strain energy density dW =dV has a singularity of inverse r near the crack tip given by $ %1 dW 2 ¼ a11 KI2 þ 2a12 KI KII þ a22 KII2 þ a33 KIII dV r ð63Þ It is apparent that the coefficient of the 1/r term in the above equation is defined only if r 6¼ 0. Hence, the element volume is always kept at a finite distance away from the crack border. This coefficient is termed as the strain energy density factor, i.e., 2 S ¼ a11 KI2 þ 2a12 KI KII þ a22 KII2 þ a33 KIII
6. Crack growth criterion The possible change of crack growth direction due to the interaction of different defects in an elastic body is of great interest to many researchers. It is worth noting that due to the influences of
ð64Þ
where KI , KII and KIII are the mode I, II and III stress intensity factors, respectively, which vary from point to point along the crack border. The coefficients aij , which vary with the polar angle d measured from the crack tip, as shown in Fig. 2, are given by
90
a11 ¼
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
1 ½ð3 4m cos dÞð1 þ cos dÞ 16l
ð65Þ
7. Numerical examples 7.1. Circular holes
a12 ¼ a22 ¼
a33 ¼
1 2 sin d½ cos d ð1 2mÞ 16l 1 ½4ð1 mÞð1 sin dÞ 16l þ ð1 þ cos dÞð3 cos d 1Þ 1 16l
ð66Þ
ð67Þ ð68Þ
where l and m are the shear modulus of elasticity and Poisson’s ratio, respectively. For two-dimensional problems where the crack extends in the xyplane, the stress intensity factors do not vary along the crack front, and the strain energy density factor S is only the function of the angle d. It is assumed that the crack will begin to extend in a direction, for which the strain energy density factor possesses a relative minimum value, i.e., 2 oS oS ¼ 0; >0 ð69Þ od d¼d0 od2 d¼d0 Based on oS=od ¼ 0 at d0 , the direction of initial crack growth d0 can be determined. Note that the above total strain energy is evaluated along a constant circle with radius r around the crack tip based on the theory of elasticity.
(a)
Fig. 5(a) and (b) show a square plate with a central circular hole subjected to two different loadings. The model devised for the finite element analyses using the semi-analytic-type element is shown in Fig. 6(a). Note that the proposed element was employed to model the region surrounding the circular hole, and 24 eight-node isoparametric elements were used to model the remaining region of the plate. Tables 1–3 present the normalized hoop stresses, rh =r0 , along the free boundary of the hole when the plate was subjected to (1) uniform tensile stresses in the x direction; (2) uniform tensile stresses in both the x and y directions; and (3) uniform tensile and compressive stresses in the x and y directions, respectively. Ten analyses were performed, by varying the ratio of the radius of the hole to half the side length of the proposed element, r=A, for each of the load case studied. It is obvious that the numerical results obtained using the proposed element are in good agreement with the corresponding theoretical results for all the load cases studied. Note that the maximum discrepancy is less than 5% for all the three load cases when r=A 6 0:7. Moreover, the maximum discrepancy is still less than 8% for all the three load cases when 0:7 6 r=A 6 0:8.
(b)
Fig. 5. A square plate with a central circular hole subjected to prescribed loads.
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
91
Fig. 6. Mesh of a square plate with a small central hole/slit crack devised using a special element (occupying the central hatch area) and some isoparametric elements.
Table 1 Normalized stress rh =r0 at the hole in a square plate subjected to uniform tensile stress in the x direction by using the semi-analytictype special element r=A
0.01 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Exact
h (°) 0
15
30
45
60
75
90
1.0244 (2.4%) 1.0287 (2.9%) 1.0291 (1.9%) 1.0194 (1.9%) 1.0089 (0.9%) 0.9873 (1.3%) 0.9692 (3.0%) 0.9575 (4.3%) 0.9237 (7.0%) 0.8748 (12.5%) 1.0
0.7566 (3.3%) 0.7602 (3.8%) 0.7602 (3.8%) 0.7521 (2.7%) 0.7428 (1.5%) 0.7249 (1.0%) 0.7120 (2.7%) 0.6991 (4.5%) 0.6931 (5.3%) 0.6843 (6.5%) 0.7321
0.0249 (–) 0.0206 (–) 0.0257 (–) 0.0216 (–) 0.0151 (–) 0.0635 (–) 0.0047 (–) 0.0091 (–) 0.0366 (–) 0.1034 (–) 0
0.9747 (2.5%) 0.9756 (3.4%) 0.9778 (2.2%) 0.9767 (3.3%) 0.9802 (2.0%) 0.9794 (2.1%) 0.9724 (2.7%) 0.9599 (4.0%) 0.9227 (7.8%) 0.8327 (16.7%) 1.0
1.9743 (1.3%) 1.9777 (1.1%) 1.9813 (0.9%) 1.9754 (1.7%) 1.9772 (1.1%) 1.9699 (1.5%) 1.9619 (2.0%) 1.9594 (2.0%) 1.9539 (2.3%) 1.9331 (3.3%) 2.0
2.7060 (1.0%) 2.7113 (0.6%) 2.7160 (0.6%) 2.7067 (0.9%) 2.7078 (0.9%) 2.6977 (1.3%) 2.6938 (1.4%) 2.7097 (0.8%) 2.7535 (0.8%) 2.8417 (4.0%) 2.7321
2.9738 (0.8%) 2.9799 (0.7%) 2.9849 (0.5%) 2.9744 (0.9%) 2.9752 (0.8%) 2.9642 (1.2%) 2.9624 (1.3%) 2.9870 (0.4%) 3.054 (1.8%) 3.1941 (6.5%) 3.0
7.2. Cracks Fig. 3 shows a square plate with a crack of inclined angle b ¼ 45° subjected to uniform tensile stress in the y direction. The model devised for the finite element analysis by using the semi-analytic-
type special element is shown in Fig. 7. Note that the proposed element was employed to model the region surrounding the crack, and eight-node isoparametric elements were used to model the remaining region of the plate. Eleven analyses were performed, by varying the ratio of the crack length
92
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
Table 2 Normalized stress rh =r0 at the hole in a square plate subjected to uniform tensile stresses in both the x and y directions by using the semi-analytic-type special element r=A
0.01 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Exact
h (°) 0
15
30
45
60
75
90
1.9495 (2.5%) 1.9512 (2.4%) 1.9560 (2.2%) 1.9555 (2.2%) 1.9673 (1.1%) 1.9784 (1.1%) 1.9954 (0.2%) 2.0425 (2.1%) 2.1340 (6.7%) 2.3240 (16.2%) 2.0
1.9495 (2.5%) 1.9512 (2.4%) 1.9558 (2.2%) 1.9548 (2.3%) 1.9654 (1.7%) 1.9732 (1.3%) 1.9824 (0.9%) 2.0114 (0.6%) 2.0614 (3.1%) 2.1590 (7.5%) 2.0
1.9495 (2.5%) 1.9511 (2.4%) 1.9555 (2.2%) 1.9536 (2.3%) 1.9616 (1.9%) 1.9629 (1.9%) 1.9562 (2.2%) 1.9488 (2.6%) 1.9156 (4.2%) 1.8281 (8.9%) 2.0
1.9495 (2.5%) 1.9510 (2.5%) 1.9553 (2.2%) 1.9530 (2.4%) 1.9598 (2.0%) 1.9577 (2.1%) 1.9431 (2.8%) 1.9175 (4.1%) 1.8427 (7.9%) 1.6625 (16.9%) 2.0
1.9495 (2.5%) 1.9511 (2.4%) 1.9555 (2.2%) 1.9536 (2.3%) 1.9616 (1.9%) 1.9629 (1.9%) 1.9562 (2.2%) 1.9488 (2.6%) 1.9156 (4.2%) 1.8281 (8.9%) 2.0
1.9495 (2.5%) 1.9512 (2.4%) 1.9558 (2.2%) 1.9548 (2.3%) 1.9654 (1.7%) 1.9732 (1.3%) 1.9824 (0.9%) 2.0114 (0.6%) 2.0614 (3.1%) 2.1590 (7.5%) 2.0
1.9495 (2.5%) 1.9512 (2.4%) 1.9560 (2.2%) 1.9555 (2.4%) 1.9673 (1.1%) 1.9784 (1.1%) 1.9954 (0.2%) 2.0425 (2.1%) 2.1340 (6.7%) 2.3240 (16.2%) 2.0
Table 3 Normalized stress rh =r0 at the hole in a square plate subjected to uniform tensile and compressive stresses in the x and y directions by using the semi-analytic-type special element, respectively r=A
0.01 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Exact
h (°) 0
15
30
45
60
75
90
4.0003 (0.0%) 4.0081 (0.2%) 4.0123 (0.3%) 3.9905 (0.2%) 3.9785 (0.5%) 3.9456 (1.4%) 3.9290 (1.8%) 3.9210 (2.0%) 3.9607 (1.0%) 4.047 (1.2%) 4.0
3.4643 (0.0%) 3.4711 (0.2%) 3.4749 (0.3%) 3.4562 (0.2%) 3.4461 (0.5%) 3.4162 (1.4%) 3.3972 (1.9%) 3.3978 (1.9%) 3.4327 (0.9%) 3.5081 (1.3%) 3.4641
2.0001 (0.0%) 2.0040 (0.2%) 2.0062 (0.3%) 1.9955 (0.2%) 1.9898 (0.5%) 1.9727 (1.4%) 1.9619 (1.9%) 1.9624 (1.9%) 1.9828 (0.9%) 2.0264 (1.3%) 2.0
0 (–) 0 (–) 0 (–) 0 (–) 0 (–) 0 (–) 0 (–) 0 (–) 0 (–) 0 (–) 0.0
2.0001 (0.0%) 2.0040 (0.2%) 2.0062 (0.3%) 1.9955 (0.2%) 1.9898 (0.5%) 1.9727 (1.4%) 1.9619 (1.9%) 1.9624 (1.9%) 1.9828 (0.9%) 2.0264 (1.3%) 2.0
3.4643 (0.0%) 3.4711 (0.2%) 3.4749 (0.3%) 3.4562 (0.2%) 3.4461 (0.5%) 3.4162 (1.4%) 3.3972 (1.9%) 3.3978 (1.9%) 3.4327 (0.9%) 3.5081 (1.3%) 3.4641
4.0003 (0.0%) 4.0081 (0.2%) 4.0123 (0.3%) 3.9905 (0.2%) 3.9785 (0.5%) 3.9456 (1.4%) 3.9290 (1.8%) 3.9210 (2.0%) 3.9607 (1.0%) 4.047 (1.2%) 4.0
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
93
Fig. 7. Mesh of a square plate with an inclined crack (b ¼ 45°) devised using a semi-analytic-type element and some isoparametric elements.
Fig. 8. Mesh of a square plate with two collinear cracks devised using two semi-analytic-type elements and some isoparametric elements.
Table 4 Stress intensity factor KI and KII at the crack tip (z ¼ a) of the internal crack with length 2a and inclined angle b ¼ 45° in a square plate subjected to uniform tensile stress in the y direction by using the semi-analytic-type special element
Fig. 4(a) and (b) show a square plate containing a single crack and two collinear cracks, respectively. The model devised for the finite element analyses by using the semi-analytic-type special element is shown in Fig. 8. Note that for both the single and double cracks, the proposed element was again employed to model only the region
a/A
0.01 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Numerical solutions
Theoretical solutions
Error
KI
KII
KI
KII
KI
KII
0.0909 0.2035 0.2881 0.4084 0.5029 0.5842 0.6603 0.7343 0.8094 0.8879 0.9699
0.0907 0.2034 0.2880 0.4083 0.5028 0.5843 0.6603 0.7341 0.8092 0.8876 0.9701
0.0886 0.1982 0.2803 0.3964 0.4854 0.5605 0.6266 0.6865 0.7415 0.7927 0.8408
0.0886 0.1982 0.2803 0.3964 0.4854 0.5605 0.6266 0.6865 0.7415 0.7927 0.8408
2.60% 2.67% 2.78% 3.03% 3.61% 4.23% 5.38% 6.96% 9.16% 12.01% 15.35%
2.37% 2.62% 2.75% 3.00% 3.58% 4.25% 5.38% 6.93% 9.13% 11.97% 15.38%
to the side length of the proposed element, a=A, for the load case studied. Table 4 presents the stress intensity factor KI and KII at the crack tip (z ¼ a) of the inclined crack. It is obvious that the maximum discrepancy between the proposed and theoretical results is less than 5.4% for a=A 6 0:5 However, the said discrepancy increases rapidly with the increase of a=A.
Table 5 Stress intensity factor KI and KII at the crack tip (z ¼ a) of the internal crack with length 2a in a square plate subjected to uniform tensile stress in the y direction by using the semi-analytic-type special element a/A
0.01 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Numerical solutions
Exact solutions
Error
KI
KII
KI
KII
KI
KII
0.1820 0.4071 0.5761 0.8167 1.0047 1.1684 1.3204 1.4685 1.6188 1.7757 1.9397
1.295e5 4.258e4 2.385e3 1.285e2 3.258e2 5.993e2 9.203e2 1.250e1 1.572e1 1.860e1 2.102e1
0.1772 0.3963 0.5605 0.7927 0.9709 1.1212 1.2536 1.3734 1.4836 1.5862 1.6827
0 0 0 0 0 0 0 0 0 0 0
2.70% 2.73% 2.78% 3.03% 3.48% 4.21% 5.33% 6.92% 9.11% 11.95% 15.27%
– – – – – – – – – – –
94
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
Table 6 Stress intensity factor KI and KII at the crack tip (z ¼ a) of the two collinear cracks, each of length 2a, in a square plate (W =A ¼ 4) subject to uniform tensile stress in the y direction by using the semi-analytic-type special element a/A
Numerical solutions Crack 1
0.01 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Error Exact solutions
Crack 1
KI
KII
Crack 2 KI
KII
KI
KII
KI
KII
KI
KII
0.1818 0.4066 0.5754 0.8155 1.0028 1.1657 1.3166 1.4701 1.6102 1.7605 1.9128
1.255e5 4.771e4 2.745e3 1.481e2 3.759e2 6.916e2 1.063e1 1.258e1 1.839e1 2.189e1 2.494e1
0.1818 0.4066 0.5754 0.8155 1.0028 1.1657 1.3166 1.4701 1.6102 1.7605 1.9128
1.326e5 4.753e4 2.749e3 1.481e2 3.759e2 6.920e2 1.064e1 1.258e1 1.840e1 2.190e1 2.496e1
0.1772 0.3964 0.5606 0.7935 0.9731 1.1256 1.2615 1.3859 1.5021 1.6123 1.7180
0 0 0 0 0 0 0 0 0 0 0
2.60% 2.57% 2.64% 2.77% 3.05% 3.56% 4.37% 6.08% 7.20% 9.19% 11.3%
– – – – – – – – – – –
2.60% 2.57% 2.64% 2.77% 3.05% 3.56% 4.37% 6.08% 7.20% 9.19% 11.3%
– – – – – – – – – – –
Table 7 Normalized stresses rh =r0 along the free boundary of circular hole in a square plate subjected to uniform tensile stress in x direction with different type special elements r=A
h (°) 0
0.01 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Exact
Table 8 Stress Intensity factor KI and KII at the crack tip (z ¼ a) of the internal crack with length 2a in a square plate subjected to uniform tensile stress in y direction with different type special elements a=A
90
Semi-analytic type
Hybrid type
Semi-analytic type
Hybrid type
1.0244 (2.4%) 1.0287 (2.9%) 1.0291 (1.9%) 1.0194 (1.9%) 1.0089 (0.9%) 0.9873 (1.3%) 0.9692 (3.0%) 0.9575 (4.3%) 0.9237 (7.0%) 0.8748 (12.5%) 1.0000
1.0000 (0%) 1.0000 (0%) 0.9995 (0.05%) 0.9992 (0.08%) 0.9987 (0.13%) 0.9985 (0.15%) 0.9972 (0.28%) 0.9955 (0.45%) 0.9934 (0.66%) 0.9842 (1.6%) 1.0000
2.9738 (0.8%) 2.9799 (0.7%) 2.9849 (0.5%) 2.9744 (0.9%) 2.9752 (0.8%) 2.9642 (1.2%) 2.9624 (1.3%) 2.9870 (0.4%) 3.054 (1.8%) 3.1941 (6.5%) 3.0000
3.0000 (0%) 3.0000 (0%) 3.0003 (0.01%) 3.0013 (0.03%) 3.0030 (0.10%) 3.0052 (0.14%) 3.0060 (0.20%) 3.0105 (0.35%) 3.0150 (0.5%) 3.0270 (0.9%) 3.0000
surrounding the crack, and eight-node isoparametric elements were used to model the remaining
Crack 2
0.01 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Semi-analytic type
Hybrid type
KI
KII
KI
KII
KI
KII
0.1820 (2.70%) 0.5761 (2.78%) 0.8169 (3.05%) 1.0051 (3.53%) 1.1691 (4.29%) 1.3215 (5.44%) 1.4701 (7.08%) 1.6209 (9.31%) 1.7780 (12.16%) 1.9419 (15.49%)
1.295e5
0.1772 (0%) 0.5606 (0.01%) 0.7929 (0.02%) 0.9716 (0.08%) 1.1222 (0.12%) 1.2556 (0.2%) 1.3777 (0.35%) 1.4948 (0.8%) 1.6010 (1.0%) 1.7067 (1.5%)
0
0.1772
0
0
0.5605
0
0
0.7927
0
0
0.9708
0
0
1.1210
0
0
1.2533
0
0
1.3729
0
0
1.4829
0
0
1.5853
0
0
1.6815
0
2.385e3 1.285e2 3.258e2 5.993e2 9.203e2 1.257e1 1.583e1 1.875e1 2.121e1
Exact solutions
region of the plate. Tables 5 and 6 present the stress intensity factor KI and KII at the crack tip for the case of single and double cracks, respectively.
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
95
Fig. 9. Comparison of the analytical and numerical results of the normalized KI for the case of a crack aligned along the centerline of a hole.
It is obvious that the proposed and theoretical results are in good agreement. 7.3. Comparison of semi-analytic and hybrid-element types From Tables 7 and 8, it is obvious that the results obtained using the hybrid-type elements are in significantly better agreement with the theoretical solutions than those obtained using the semianalytic-type elements. This is due to the fact that
the conditions of displacement continuity along the inter-element boundaries have been implemented in the former element types. Note that the hybrid-type elements are to be employed together with four-node isoparametric elements in devising a finite element model. This is to ensure the same displacement distributions at the interface between the special and the traditional isoparametric elements. As for the semianalytic-type elements, these elements can be used together with the conventional elements of the
96
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
Fig. 10. Comparison of the analytical and numerical results of the normalized KI for the case of a crack situated at a distance t from two holes.
same number of nodes and provide reasonably good accuracy. Fig. 6(a) and (b) show the finite element meshes devised using 1 semi-analytic-type element together with 24 eight-node isoparametric elements and 1 hybrid-type element together with 32 four-node isoparametric elements, respectively.
7.4. Interactions of voids and cracks Figs. 9–11 provide some examples for the modeling of the interaction between two circu-
lar holes and a crack using the high precision hybrid-type elements. Fig. 9 presents the analytical and numerical results of the normalized KI , i.e., pffiffiffi KI =ðr0 aÞ, of a crack, aligned along the centerline of a hole at a distance t apart, embedded in an infinite plate subjected to uniaxial tension r0 . The results obtained using the proposed hybrid-type element are in excellent agreement with the analytic results obtained in [15]. Fig. 10 shows an infinite plate containing a crack and two holes subjected to uniaxial loading r0 . The results obtained using the proposed hybrid-type element are again in excellent agreement with the analytic re-
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
97
Fig. 11. Comparison of the analytical and numerical results of the normalized KI for the case of a crack separating two holes.
sults obtained by Hu et al. [15]. This excellent agreement is also clearly shown in the case of an infinite plate, which contains a crack separating two holes, subject to uniaxial loading r0 as shown in Fig. 11. Fig. 12(a)–(c) show an infinite plate containing a semi-infinite macro-crack and a micro-defect, i.e., circular hole, elliptic hole or crack, respectively, subjected to uniaxial loading r0 in the y direction. The finite element mesh devised for analysis is shown in Fig. 12(d). Note that the hy-
brid element containing an edge crack developed in [21] was employed to model the region surrounding the tip of the main crack. Whereas, the proposed hybrid-type special elements with internal defects were employed to model the region surrounding the micro-defects. With reference to Figs. 13–17, it is obvious that the results obtained using the proposed hybrid-type elements are again in excellent agreement with the corresponding analytical results obtained in [16,17]. Obviously, the existence of micro-defects has significant effects
98
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
Fig. 12. A large plate with a long edge macro-crack and a micro-defect (in the shaded zone) and the FE model devised.
Fig. 13. Comparison of the analytical and numerical results of the normalized KI for the case of a macro-crack and a micro-defect.
on the macro-crack. As a result, the original mode I crack problem became a mixed-mode crack problem and the initial crack growth direction of the macro-crack was also changed.
8. Conclusions Eight-node rectangular elements, each with a central circular or elliptical hole or an internal
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
99
Fig. 14. Comparison of the analytical and numerical results of the normalized KII for the case of a macro-crack and a micro-elliptichole.
Fig. 15. Comparison of the analytical and numerical results of the normalized KII for the case of a macro-crack and a micro-crack.
crack, have been developed for modeling multiple voids, cracks to study the interaction between these defects. The good accuracy of the semi-analytictype elements are achieved due to the exact defi-
nition of the stress and strain fields within the elements and the implementation of the free boundary conditions along the free boundary of the hole or crack. However, these semi-analytic-type elements
100
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
Fig. 16. Comparison of the analytical and numerical results of the initial crack growth direction d0 of the macro-crack for the case of a macro-crack and a micro-elliptic-hole.
Fig. 17. Comparison of the analytical and numerical results of the initial crack growth direction d0 of the macro-crack for the case of a macro-crack and a micro-crack.
have been improved by implementing displacement continuity at the inter-element boundary to establish the hybrid-type elements. The proposed ele-
ments can be easily combined with the conventional elements, such as isoparametric elements, without any modifications, to obtain good results
C.H. Yang, A.K. Soh / Theoretical and Applied Fracture Mechanics 38 (2002) 81–101
using simple finite element meshes. It is interesting to note that both the semi-analytic-type elements and hybrid-type elements are particularly useful for modeling the interaction between many voids and cracks. The former element type is less accurate compared with the latter but much easier to construct and use. Acknowledgements Supports from the Research Grants Council of the Hong Kong Special Administrative Region, China (Project No. HKU 7081/00E) is acknowledged. The authors would like to thank Professor Daining Fang of Tsinghua University for the useful discussions. References [1] M. Isida, D.H. Chen, H. Nisitani, Plane problems of an arbitrary array of crack emanating from the edge of an elliptical hole, Eng. Fract. Mech. 21 (1985) 983–995. [2] M. Isida, H. Igawa, Analysis of a zig-zag array of circular holes in an infinite solid under uniaxial tension, Int. J. Solids Struct. 27 (1991) 227–258. [3] M. Isida, S. Nemat-Nasser, A unified analysis of various problems relating to circular holes with edge cracks, Eng. Fract. Mech. 27 (1987) 849–864. [4] G.C. Sih, Handbook of stress intensity factor, Lehigh University, Bethlehem, PA, 1973. [5] C. Atkinson, The interaction between a crack and an inclusion, Int. J. Eng. Sci. 10 (1972) 127–136. [6] F. Edorgan, G.D. Gupta, M. Ratawani, Interaction between a circular inclusion and an arbitrary oriented crack, J. Appl. Mech. 41 (1974) 1007–1013. [7] Y.Z. Chen, General case of multiple crack problems in an infinite plate, Eng. Fract. Mech. 20 (1984) 591–597. [8] Y.Z. Chen, Solutions of multiple crack problems of a circular plate or an infinite plate containing a circular hole by using Fredholm integral equation approach, Int. J. Fract. 25 (1984) 155–168. [9] Y.Z. Chen, A special boundary-element formulation for multiple-circular-hole problems in an infinite plate, Comput. Meth. Appl. Eng. Mech. 50 (1985) 263–273. [10] H. Tada, P. Paris, G. Irwin, The Stress Analysis of Cracks Handbook, Paris Production Inc, St. Louis, MO, 1985. [11] M. Isida, Analysis of stress intensity factors for plates containing random array of cracks, Bull. Jpn. Soc. Mech. Eng. 13 (1970) 635–642. [12] M. Isida, On the determination of stress intensity factors for some common structural problems, Eng. Fract. Mech. 2 (1970) 61–79.
101
[13] M. Isida, Methods of laurent series expansion for internal crack problems, in: G.C. Sih (Ed.), Methods of Analysis and Solutions of Crack Problems, Noordhoff, Leyden, 1973, pp. 56–130. [14] G.J. Rodin, Y.-L. Hwang, On the problem of linear elasticity for an infinite region containing an infinite number of non-intersecting spherical inhomogeneities, Int. J. Solid Struct. 27 (1991) 145–159. [15] K.X. Hu, A. Chandra, Y. Huang, Multiple void-crack interaction, Int. J. Solids Struct. 30 (11) (1993) 1473–1489. [16] S.X. Gong, H. Horii, General solution to the problem of microcracks near the tip of a main crack, J. Mech. Phys. Solids 37 (1) (1989) 27–46. [17] S.X. Gong, S.A. Meguid, Microdefect interacting with a main crack: a general treatment, Int. J. Mech. Sci. 34 (12) (1992) 933–945. [18] S.A. Meguid, P.E. Gaultier, S.X. Gong, A comparison between analytical and finite element analysis of main crack-microcrack interaction, Eng. Fract. Mech. 38 (6) (1991) 451–465. [19] A.K. Soh, An improved method for determining free boundary stresses, J. Strain Anal. 27 (2) (1992) 93–99. [20] A.K. Soh, Development of special finite elements for free boundaries, Int. J. Solids Struct. 36 (6) (1999) 899–917. [21] P. Tong, T.H.H. Pian, S.J. Lasry, A hybrid-element approach to crack problems in plane elasticity, Int. J. Numer. Meth. Eng. 7 (1973) 297–308. [22] R. Piltner, Special finite elements with holes and internal cracks, Int. J. Numer. Meth. Eng. 21 (1985) 1471–1485. [23] A.K. Soh, Z.F. Long, A high precision element with a central circular hole, Int. J. Solids Struct. 36 (35) (1999) 5485–5497. [24] P. Tong, New displacement hybrid finite element models for solid continua, Int. J. Numer. Meth. Eng. 2 (1970) 73–83. [25] N.I. Muskhelishvili, Some Basic Problems of the Mathematical Theory of Elasticity, fourth ed., Noordhoff International Publishing, Leyden, Netherlands, 1954. [26] T.H.H. Pian, Derivation of element stiffness matrices by assumed stress distribution, AIAA J. 2 (1964) 1333–1336. [27] G.C. Sih, P.C. Paris, F. Erdogan, Crack-tip stress intensity factors for plane extension and plate bending problems, J. Appl. Mech. 29 (1962) 306–312. [28] S.A. Meguid, in: Engineering Fracture Mechanics, Elsevier Science Publisher LTD, Great Britain, 1989, pp. 134–136. [29] G.C. Sih, A special theory of crack propagation, in: G.C. Sih (Ed.), Mechanics of fracture 1: Methods of analysis and solutions of crack problems, Noordhoff Internatioinal Publishing, Leyden, Netherlands, 1973, pp. XXI–XLV. [30] G.C. Sih, Some basic problems in fracture mechanics and new concepts, Eng. Fract. Mech. 5 (1973) 365–377. [31] G.C. Sih, Strain-energy-density factor applied to mixed mode crack problem, Int. J. Fract. 10 (3) (1975) 305–321. [32] G.C. Sih, A three-dimensional strain density factor theory of crack propagation, in: M.K. Kassir, G.C. Sih (Eds.), Mechanics of fracture 2: Three dimensional crack problems, Noordhoff International Publishing, Leyden, Netherlands, 1975, pp. XV–LIII.