- Email: [email protected]
Determining the stress intensity factors KI, KII and the T-term via the conservation laws using the boundary element method
Engineering Fracture Mdwnicr
0013-7944(94)EOO88-X
Printed
49. No.
Vol.
Copyright
Pergamon
‘c
in Great
I, pp.49-60.
1994 Else&r Britain.
All
0013.7944/94
1994
Science Ltd rights
reserved
$7.00 + 0.00
DETERMINING THE STRESS INTENSITY FACTORS K,, K,, AND THE T-TERM VIA THE CONSERVATION LAWS USING THE BOUNDARY ELEMENT METHOD P. C. OLSEN Department
of Structural
Engineering, Lyngby,
Technical Denmark
University
of Denmark,
2800.
Abstract-A new method for determining the stress intensity factors K, and K,, is presented. The method is based on the conservation laws, from which the property that the J,, J, and L integrals are zero for any closed path of integration is utilized. The integrals above are evaluated numerically for a circular path surrounding the crack tip and along the crack faces, ending, however, a distance A from the tip. At the distance A, all fields are replaced by the singular first-order expansion series, including the T-term for r~,, Based on the J,, J, and L integrals, numerically evaluated as described above, equations are derived, from which both the stress intensity factors K, and K,,, as well as the T-term are determined. Numerical tests based on the boundary element method are presented, in which the accuracy of the method is demonstrated. The total energy release rate. represented by the J, integral, is typically determined with an accuracy of within l-2%. The distribution of the total energy release rate between mode I and mode II is, however, less accurately determined. The T-term is determined with a rather large uncertainty.
1. INTRODUCTION MANY possibilities exist for evaluating the total energy release rate due to a change in the crack surface geometry. The difficult part is deciding which portion of the energy release rate has to be subscribed to which mode. In mixed mode problems, it is essential to distinguish between the different modes. Evaluating the J integral due to Rice [l], one can easily determine the total energy release rate. Budiansky and Rice [2] later realised that the J integral is actually the first component of a vector and that all components can be interpreted as energy release rates, when a cavity or a flat crack is given a translation. Also in [2], the L and M integrals were interpreted as energy release rates for rotation and expansion, respectively, of cavities and flat cracks. Numerically it is rather attractive to determine the energy release rate via the J integral because it can be evaluated in a stress field, far away from the crack tip and thus avoiding the singular stresses. It is not necessary, numerically, to model the singular behaviour near the crack tip correctly. With a sufficiently fine mesh, standard shape functions can be applied. However, in using the second component of the J integral in order to determine both K, and K,, , terms involving near field stresses must be included, whereby the attractiveness of the procedure becomes less pronounced. Eischen [4] applied the J integrals in order to determine both K, and K,, , as well as the constant T-term of the first-order stress expansion. The finite element method was applied and both standard elements, as well as singular crack tip elements were used. Both K,, and the T-term were reported to be rather sensitive to the types of element used at the crack tip. In this paper the determination of both K, and K,,, as well as the T-term in plane isotropic elasticity by means of evaluating the J integrals and the L integral is demonstrated. The integrals have to be evaluated in near field stresses. However, the integrands are replaced by singular terms from the first-order expansion closest to the crack tip, where the numerical values are unreliable. The boundary element method is applied.
2. THE
CONSERVATION
LAWS
The conservation laws are the direct consequence of the basic equations of the theory of elasticity. Knowles and Sternberg [3] have shown that, when these equations are all satisfied, the 49
50
following
P. C. OLSEN
integrals:
(1)
(2)
(3)
(4)
vanish for arbitrary closed integration paths S surrounding a homogeneous singly connected isotropic elastic domain. In eqs (l-4) W is the strain energy density, n is the unit outward normal vector of S, u is the deformation vector, L is the first-order strain measure, p is the traction vector at S corresponding to the unit outward normal vector and finally, x and y are coordinates of points at S. The subscripts refer to components in a global Cartesian reference frame. The conservation laws are now applied for a path enclosing a crack tip as shown in Fig. 1. This path is composed of three parts as shown: firstly, an arbitrarily shaped curve starting at the crack face at a given distance from the crack tip, surrounding the crack tip and ending at the opposite crack face. In practice this curve is often chosen as a circle; secondly, a circle with centre at the crack tip and a small radius r surrounding the crack tip; and thirdly, lines coinciding with both crack faces, connecting the first two curves. The complete path is thus closed, resulting in the following equations being valid: Jc”’ + J.y
+ J;”
= 0
J,v
+
JTb ?
=
0
L Face +
LTip
=
0,
+
L Far +
where the superscripts this paper.
JFe ?
refer to the parts shown in Fig. 1. The M integral
Fig. 1. Particular
integration
path for calculating
the conservation
(3
(6)
(7) is not further
integrals.
applied
in
51
Determining K,, K,, and the T-term using BEM
3. J AND L AT THE Displacements and expansion series, i.e.
stresses
close
CRACK
to the crack
TIP
tip are assumed
-sin(8/2)(2
to follow
the first-order
+ cos(8/2)cos(38/2))
T
sin(e/2)cos(e/2)cos(38/2)
+
i c0s(e/2)(i - sin(B/2)sin(38/2))
0 ,
(9)
(i 0
j
in which G is the shear modulus. ICis equal to (3 - 4v) in plane strain and (3 - v)/(l + v) for plane stress (see [6]). K, and K,, are the stress intensity factors. Regarding the inclusion of the T-term, see for example Eftis [5]. Introducing these expressions into the J and L integrals, performing the integration for the circle surrounding the crack tip and evaluating the results in the limit of vanishing radius, the following results may be derived (see [7]):
(10)
4. J AND L AT THE The J and L integrals
CRACK
FACES
at the crack faces are given by: JFaCe .r
Jc”“=
0
‘(W+ s
LFace
=
=
-
(13) Wm)dr
(14)
0
A (W’
-
W-jr
&,
(15)
s 0
where W+ and Wm are the strain energy at the upper and lower faces, respectively. Assuming displacements and stresses to follow the first-order expansion series eqs (8) (9) at the crack faces close to the crack tip, the stress o,, along the crack faces is given by: T - 2 J(::i),
8 = z
(Upper)
(T,,.=
(16) 8 = -7t (Lower)
resulting
in a difference
in strain
w+_w-=-&
energy
between
the two crack faces:
=
4
K,, T E’J(2w)
’
(17)
P. C. OLSEN
52
E’ equals E, the modulus of elasticity in plane stress; in plane strain E’ = E/( 1 - v’) where 1’is the Poisson ratio. JX and L can now be evaluated explicitly, by introducing eq. (17) into eqs (14), (I 5): JF””
=
_
8
K1t
(18)
T&
E'JO~)
(19)
5. DETERMINING
K,,K,,AND THE T-TERM
Both the stress intensity factors K, and K,,, as well as the T-term can now be determined by means of the J and L integrals evaluated numerically in stress fields far from the crack tip. It is only a matter of simple manipulation to arrive at the following results:
(20) K,, = J(E’J, T =
- K;)
(21)
3E'JW W,fi
(22)
In eqs (20) and (21) the superscript “Far” of J,. J,. and L has been suppressed. Equations (20)-(22) are the basic equations to be thoroughly investigated next. In the case of K,, becoming zero, i.e. for pure mode I, the method above is not applicable. In this case the strain energy difference between the opposite crack faces is zero, which means that .lVand L are path independent in the sense of J, path independence. However, the author cannot convince himself as to whether the T-term must be zero in pure mode I or not.
6. THE
BOUNDARY
ELEMENT
METHOD
The numerical analysis is done by means of the boundary element method. There is no need to reiterate this method here and the reader is referred to the literature, see for example Brebbia [S]. When modelling flat cracks, the method considered most suitable is the multi-domain boundary element method (Blandford [9]), where the domain in front of each crack is divided up into subdomains, in order to prevent the singularities in the system matrix otherwise arising. When determining the J and L integrals, one needs to determine very precisely displacements and stresses close to the boundary, due to the subdivision. The boundary element method is not famous for being very accurate in determining stresses near the boundary. On the contrary, it is almost a standard rule not to determine stresses at internal points which are closer to the boundary than the length of an element. However, it is possible to achieve very accurate results at points even extremely close to the boundary, provided that integration is performed analytically. This demands that either the boundary is not curved or that the curvature is small enough so that the boundary can be approximated with straight line segments. The program being used supports two classes of elements, Lagrange elements and cubic spline elements. The cubic spline elements have been implemented in order to achieve continuity in the normal stresses, acting tangential to the boundary. This stress component is derived from the first derivatives of the boundary displacements, and thus differentiability is needed for the boundary displacements. The cubic spline elements satisfy this requirement internally. The same effect is achieved by using high-order Lagrange elements. However, high-order Lagrange polynomiums are very poor interpolants. The nodes in the cubic spline element may be positioned with variable where si is the spacing between spacing. A simple power law has been implemented, i.e. s,+ , = fsi, two consecutive nodes. Normally quadratic Lagrange elements are used throughout, except for at the crack faces and the elements adjacent to the crack tips, where the cubic spline elements are used.
Determining
K,, K,, and the T-term
53
using BEM
The partial cross derivatives au,/+ and &,./ax appearing in the conservation integrals are determined as finite differences, i.e. by determining the displacements at closely spaced points, calculating the differences and dividing by the distance. One can achieve very accurate results, which are easy to check, as the sum of the derivatives must equal twice the shear strain. However, it is advisable to perform the calculation by determining the kernels in the boundary element method corresponding to the cross derivatives (see Aliabadi and Rooke [6]), analytically, just as the displacements and stresses are determined analytically. Boundary stress discontinuities are modelled by means of the double node concept, i.e. two nodes are placed at boundary points, where stress discontinuities may be expected, which is the case at a crack tip. For double nodes, where unknown tractions appear at both nodes, discontinuous elements are used. 7. INTRODUCTORY
NUMERICAL
TEST
The example of a slanted crack subjected to a uniform pressure is used for the introductory numerical calculations. The geometry of the specimen is shown in Fig. 2. The dots indicate the nodal points of the boundary element discretization. In the enlargement, the nodal points close to the crack tip are shown. The nodal points along the crack face and at the internal boundaries next to the crack have been positioned corresponding to a factor f = 1.25 in the power law of nodal spacing. Cubic spline elements, each containing 11 nodes, are used at the crack face and at the internal boundaries. Quadratic Lagrange elements are used at the external boundaries. The length of half the crack is a = 0.1 m whereas the angle of the crack with the vertical is a = 60. For the sheet of infinite size, analytical solutions exist: K, = oj(na)
sin’s,
K,, = oJ(na)
cos a sin x
and
T = g cos 2~,
(23)
where CJ is the undisturbed far field stress (Eftis [5]). Firstly it has been tested as to what extent one can expect the conservation laws to be satisfied when integrating along a closed path not enclosing any defects, in which case the integrals must turn out to be zero. A circular path was chosen with a radius of r = 0.05 m and with the centre placed at two different positions Pl and P2 as shown in Fig. 2. The first position is characterized as having almost uniform stresses around the entire periphery of the circle, whereas the second position is so close to one of the crack tips that large variations in the stresses are expected. In Fig. 3 the variations of the integrands of J and L along the periphery are shown. In Fig. 3a, it is seen that for position Pl the variation is almost harmonic, indicating an almost constant stress field, for which the variation is solely due to the projection of the normal vector. For an integration path in a non-uniform stress field, like the one at position P2, the integrands vary irregularly as shown in Fig. 3b.
Mesh
Fig. 2. The example
of a slanted
crack
at
the
used for testing
crack
the algorithms
face
P. c‘. OLSEN
54
A)
Position
B)
Pl
-3.00 [ -180.00
Position
P2
I -90.00
90.00
Fig. 3. The variation
0.00
180.00
Angle
Angle of the integrands
of the conservation
integrals
at the positions
Pl and P2
In order to perform the integration, the curves above are interpolated by means of natural cubic splines with zero curvature at the interval ends. The spline is determined by means of 61 regularly spaced points. In Table 1 the values of the J and L integrals are shown. It is seen that the accuracy with which the J and L integrals are determined using the boundary element method is astonishingly high. It would be interesting to see if the finite element method can achieve the same degree of accuracy. However, to the author’s knowledge similar investigations for the finite element method have not been published. Having seen that the boundary element method is an excellent tool for determining the J and L integrals, the integration path is now positioned in order to determine the stress intensity factors, i.e. the centre of the circle is positioned at the crack tip, the periphery starting at one crack face and ending at the opposite crack face. In Fig. 4 the variations of the integrands of J,, Jr and L are plotted for a circle with r = 0.08 m at this position. In Table 2 the calculated values of J,, J, and L are shown for different radii of the integration circle. J, is, as it should be, practically independent of the choice of radius, with a deviation on the fourth decimal only. J? is path dependent, as has been discussed by others (Hermann and Hermann[7]). L is likewise path dependent, whereby path dependence means that J? and L depend on the choice of radius. To check path independency in the case of the integration path surrounding the crack tip, an integration path is now chosen, composed of (1) a circle with the centre positioned at the crack tip, the periphery starting at one crack surface and ending at the opposite crack surface and (2) two lines coinciding with the crack surfaces, starting at the points where the circle intersects the crack surfaces and ending a distance A from the crack tip (see Fig. 5). For a given distance A, J, and L must be independent of the choice of radius. In Table 3 the .l,, integral is tabulated for different radii and distances to the crack tip. It is seen that for a given distance A, TVis practically independent of the choice of radius. Similarly, the L integral is tabulated in Table 4. It is seen that path independence cannot be achieved. However, results improve the further away from the crack tip the integration path ends. It is also remarkable that the order of magnitude of the L integral Table 2. The J,, J, and L integrals for different radii in the case of a circle of integration surrounding the crack tip
Table 1. The J and L integrals determined for a closed path of integration positioned in a regular stress field (PI) and very near to a crack tip in a highly irregular stress field (P2) Position J., 4 L Pl P2
6.603 x 1O-6 2.807 x 10-j
-4.911 -8.629
x lo-* x IO-’
2.080 x IO-’ 2.319 x IO-’
Radius
J.
J.
0.04 0.06 0.08 0.10 0.12 0.14 0.16
0.24630 0.24622 0.24643 0.24637 0.24637 0.24641 0.24648
0.13826 0.12700 0.12120 0.11926 0.12123 0.12723 0.13763
L -0.0007071 -0.0012116 -0.0015833 -0.0017723 -0.0015023 -0.0006614 0.0009726
Determining K,, K,, and the T-term using BEM
55
Angle Fig. 4. The variations of the integrands of the conservation integrals for a circle positioned at the crack tip.
in this example is more than one hundred times smaller than that of the J integrals, although all integrands have the same order of magnitude. This of course demands a higher degree of accuracy.
8. EVALUATING
K,, K,, AND THE
Z--TERM
Remembering, that eqs (20)-(22) rely on a first-order stress expansion closest to the crack tip corresponding to a distance A, it is obvious that the smaller the distance A is, the better the approximation is and the more accurate the stress intensity factors and the T-term can be determined. However, if one comes too close to the crack tip, the numerical values become unreliable. This is especially true for the L integral, as can be seen in Table 4. In Fig. 6 the normal stresses acting tangential to the crack surface are plotted. It is seen that the numerical representation becomes unreliable at a distance of approximately 0.04m from the crack tips. This figure explains the very inaccurate results for the L integral for distances A less than approximately 0.04 m. The idea is therefore to choose a distance A as small as possible and still have reliable numerical data. In Table 5 eqs (20)-(22) are evaluated for A > 0.04 m, for which the L integral has been determined with reasonable accuracy. Also in this table, the effect of choosing different radii can be seen.
Fields represented expansion terms tance A next to
by singular over a disthe crack tip
Fig. 5. Integration path around the crack tip composed of a circle and lines along the crack surfaces.
56
P. C. OLSEN Table 3. Path mdependence
J, J, for drfferent radn and distances
of
to the crack
tip
Radius
A = 0.01
A = 0.02
A = 0.03
A = 0.04
A = 0.05
A = 0.06
A = 0.07
A = 0.08
0.04
0.1794x
0.06
0.17986
0.08 0.10 0.12 0.14 0.16
0. I7992 0.17992 0. I7997 0.17997 0.17923
0.15515 0. I5554 0. I5559 0. I5560 0. IS564 0.1 SS64 0.15490
0.14585 0.14623 0. I4628 0. I4629 0.14633 0.14634 0.14560
0. I381-6 0. 13X64 0. 13870 0.13x7 I 0.13875 0.13875 0.13x01
U.13178 0. I3 I84 0.1318s 0.13189 0.13189 0.131 I5
0. I2700 0.12706 0. I2706 0.12710 0.12711 0.12637
0.12364 0.12365 0.12369 0.12369 0.12295
0. I z IX 0.12121 0.12125 0.12125 0.12051
Table 4. Path independence
of L. L for different
radii and distances
to the crack
tip. Values of L are multiplied
_
by 100.
Radius
A = 0.01
A = 0.02
A = 0.03
A = 0.04
A = 0.05
A = 0.06
A = 0.07
A = 0.08
0.04 0.06 0.08 0.10 0.12 0.14 0.16
0.01393 0.02030 0.02354 0.02142 0.02654 0.03198 0.02724
-0.02157 -0.01520 -0.01196 -0.01408 -0.00896 -0.00353 - 0.00826
-0.0440x -0.03771 -0.03447 -0.03660 -0.03147 -0.02604 - 0.03078
-0.07071 - 0.06434 -0.061 IO -0.06322 -0.05810 -0.05267 - 0.05740
-0.09500 -0.09176 -0.09389 -0.08876 -0.08333 -0.08807
-0.12116 -0.11792 -0.12005 -0.11492 -0.10949 -0.11422
-0.14033 -0.14215 -0.13703 -0.1316 -0.13633
-0.15823 -0.16036 -0.15523 -0.14980 - 0.15454
The relative variations of the stress intensity factors determined in Table 5 are 1.5% for K, and 7.5% for K,, , respectively. The T-term is less accurately determined, with a relative variation of 13.5%. The analytical solution for the sheet of infinite size is [5] K, = 0.4204
K,, = 0.2427
T = -0.5.
(24)
The results obtained numerically deviate from these analytical values and can be explained partly by the fact that the sheet has a finite size, and partly by the fact that the support is completely fixed, preventing any movement. It is also observed that there is a difference in the total energy
Lower 4.00
-4.00
Upper
,
I
I I 0.04
I 0.00
I 0.08
Distance Fig. 6. The tangential Table 5. The stress intensity
normal
factors
along
K,, K,, and the T-term
A = 0.05 T
from
stresses
4
A = 0.06 ~ .~~~~ K,, T
I I 0.16
I 0.12
crack
I 0.20
tip
the upper and lower crack for different
surfaces
radii and distances
to the crack
A = 0.07
Radius
K,
41
4
41
0.06
0.4496
0.2100
-0.3805
0.4504
0.2082
-0.3721
-
-
0.08 0.10 0.12 0.14 0.16
0.4511 0.4502 0.4522 0.4542 0.4530
0.2071 0.2090 0.2047 0.2002 0.2031
-0.3725 -0.3777 -0.3645 -0.3499 -0.3646
0.4517 0.4509 0.4525 0.4542 0.4534
0.2059 0.2074 0.2039 0.2002 0.2024
-0.3663 -0.3701 -0.3604 -0.3498 -0.3610
0.4532 0.4525 0.4539 0.4553 0.4546
0.2026 0.2039 0.2010 0.1978 0.1995
tip
A = 0.08 T -0.3508 -0.3538 -0.3461 -0.3377 -0.3469
4
4,
-
-
0.4551 0.4545 0.4556 0.4568 0.4563
0.1984 0.1995 0.1969 0.1942 0.1955
T -0.3314 -0.3339 -0.3274 -0.3205 -0.3283
57
Determining K,, K,, and the T-term using BEM
Boundary
element
mesh
Normal
Displacements
Tangential
tractions
Fig. 7. Comparison with Eischen [4]. Boundary element mesh, displacements
normal
stresses
and stresses.
release rate G. The analytical values of eq. (24) result in G = 0.2356, whereas the numerical evaluation of J, results in a value of G = 0.2464, a result which is considered to be very accurate. 9. ADDITIONAL
RESULTS
Very few articles are concerned with the determination of the T-term. However, Eischen [4] developed a method to determine both K,,K,, ,as well as the T-term and closed his paper with a numerical example. This example deals with a rectangular sheet subject to uniform tension and penetrated by a crack, with a crack mouth starting at the left-hand side of the vertical boundary and extending upwards at an angle of 45”. The size of the sheet is 1 x 2 and the length of the crack is a = 0.5656. In Fig. 7 the boundary element representation of the sheet is shown. Also in this figure the distorted shape of the geometry, the traction normal to the boundary and the tangential normal stresses are shown. The circular path of integration is chosen with r = 0.2828, for which the variations of J,, JI and L are shown in Fig. 8(a). Integrating these curves, one gets J, = 4.323, JJ = - 3.922 and L = -0.06561. In Fig. 8(b) the tangential normal stresses along the crack faces are plotted. These stresses seem to be reliable at distances larger than approximately one-quarter of the crack length.
A)
Variation
of
integrands
B)
Stress
-5.04
‘oh
+ L o.io
’ 0.40
Distance
Angle
of the integrands
crack
Jx = 4.323 Jy =-3.775
Jx= 4.323 Jy=-3.922 L =-0.06561
Fig. 8. The variation
along
(A) and the stress distribution
I
I
I
0.30
0.40
0.40
from
along
crack
0. 0
tip
the crack
faces (B)
face
P. C. OLSEN Table 6. K,, K,, and T for varying values of A m the case of the example shown in Fig. 6. A 0.05657 0.08485 0.1131 0.1414 0.1697 0.1980 0.2263
K,
K,,
T
1.961 I .949 I.936 I.928
0.69 I2 0.7244 0.7589 0.7782 0.7948 0.8100 0.8240
-0.5890 ~ I .2444 ~ 0.9600 -0.8293 PO.7325 -0.6534 -0.5882
I.921 1.915 1.909
In Table 6 the variations of K,, K,, and T for varying distances A are shown. The most probable values are those obtained with the smallest A, for which the numerical stress representation is still reliable and the author concludes that K, = 1.936, K,, = 0.7589 and T = -0.96 are the most probable values. Eischen [4] obtained K, = 1.917, K,, = 0.8065 and T = 0.822, whereas Bowie [lo] obtained K, = 1.933 and K,, = 0.8265. The latter did not determine the T-term. It is very interesting that Eischen and the author agree upon the total energy release rate, both having a value of G = 4.325, whereas the energy release rate resulting from Bowie’s calculations is G = 4.420. On the other hand, Bowie and the author agree upon the mode I stress intensity factor, having almost identical results. Thus, if one accepts firstly that the total energy release rate has been correctly determined by Eischen and the author and secondly, that the identical values of K, obtained by Bowie and the author are correct, then K,, = 0.7589 as determined by the author must be the correct value for the mode II stress intensity factor. This article concludes with two examples concerned with the bending of beams. The beam in the first example is a four point bending problem and is a pure mode I fracture problem. Results for this beam have been presented by Blandford et al. [9]. The boundary element mesh is shown in Fig. 9. Also in the figure the distorted shape, the normal traction and the tangential normal stresses of the boundary element analysis are shown. The conservation integrals are determined for a circle surrounding the crack tip and having a radius r = 0.05. This results in J, = 0.3609 whereas JI and L are of an order of magnitude three smaller (J,. = -5.54 x 10. 5 and L = 4.74 x 10-j). The mode I SIF factor is thus determined as K, = 0.6007. Blandford et al. reports of a value K, = 3.818J(2.54/100) = 0.6085. The authors result is thus within an accuracy of 1.3%. A) Boundary
C) Normal
element
mesh.
traction.
B) Distorted
D) Tangential
Fig. 9. Mode I fracture
of beam.
shape
normal
stresses.
Determining
A) Boundary
C) Normal
element
K,, K,, and the T-term
mesh
B) Distorted
D) Tangential
traction.
Fig. 10. Mixed mode I and mode II fracture
Fig. 11. Shear stress distribution
along
vertical
59
using BEM
shape
normal
stresses.
of a beam
sections
of the beam.
Finally a beam with dimensions as before, however, with a notch placed at the quarter-point and with a point load at the centre is analysed. In Fig. 10 the boundary element mesh is shown. Also in this figure the distorted shape, the normal traction and the tangential normal stresses of the boundary element analysis are shown. The distribution of the shear stresses along vertical sections of the beam at different locations are shown in Fig. 11. The well known parabolic distribution is recognised at one end of the beam, whereas the shear stresses are wildly disturbed by the notch at the other end. It is also recognisable that the shear stresses behave singularly in front of the crack tip. The conservation integrals are determined for a circle surrounding the crack tip and having a radius r = 0.05. This results in J, = 0.2241, JY = -0.07755 and L = -0.0002011. Using a distance of A = 0.025, the correction amounts to JJ = -0.08071 and L = -0.0002917. Thus from eqns (20)-(22), KI = 0.4666, K,, = -0.07985 and T = -0.06455.
10. CONCLUSIONS A new method for deriving the stress intensity factors has been presented. The method is based on the conservation laws, from which the property that the J,, TVand L integrals are zero for any closed path of integration is utilized. The integrals above are evaluated numerically for a circular path surrounding the crack tip and along the crack faces from the point of intersection of the circle towards the crack tip, ending, however, a distance A from the tip. At the distance A, all fields are replaced by the singular first-order expansion series, including the T-term for c,,. Based on the J,, Jr and L integrals, numerically evaluated as described above, equations are derived, from which both the stress intensity factors K, and K,, , as well as the T-term are determined. Numerical tests, based on the boundary element method are presented, in which the accuracy of the method is demonstrated. The total energy release rate, represented by the J, integral is typically determined with an accuracy of within l-2%. The distribution of the total energy release rate between mode I and mode II is, however, less accurately determined. The T-term is determined with a rather large
60
P. C. OLSEN
uncertainty. This has to be seen as a result of the magnitude of the L integral, which is an order of two smaller than both J, and J, However, by carefully observing the quality of the numerical data, one may come to a conclusion with regard to the most probable values. From the numerical tests, it is concluded that the best results are achieved by using a radius of half the crack length for the circular integration path and by extending the integration along the crack faces to a distance A from the crack tip, corresponding to a quarter of the crack length.
REFERENCES [l] J. R. Rice, Mathematical analysis in the mechanics of fracture, in Fracture, An Advanced Treatise, Volume II. Mulhetical Fundamenfals (Edited by H. Liebowitz) Academic Press, New York (1968). [2] B. Budiansky and J. R. Rice, Conservation laws and energy-release rates. J. uppl. Mech. 40, 201-203 (1973). [3] J. K. Knowles and E. Sternberg, On a class of conservation laws in linearized and finite elastostatics. Arch. Rational Mech. Anal. 44, 187-211 (1972). [4] J. W. Eischen, An improved method for computing the J, integral. Engng Fracfure Mech. 26(S), 691-700 (1987). [S] J. Eftis, Load biaxiality and fracture: a two-sided history of complementing errors. Engng Fracture Mech. 26(4), 567-592 (I 987). [6] M. H. Aliabadi and D. P. Rooke, Numerical Fracture Mechanics. Computational Mechanics Publications/Kluwer Academic Publishers, Boston, MA (1991). [7] A. G. Hermann and G. Hermann, On energy-relase rates for a plane crack. J. uppI. Mech. 48, Q-528 (1981). [8] C. A. Brebbia. The Boundary Elemem Methodfor Engineers. PenTech Press (1978). [9] G. E. Blandford, A. E. IngratTea and J. A. Liggett, Two-dimensional stress intensity factor computations using the boundary element method. Int. J. numer. Methods Engng 17, 3877404 (1981). [IO] 0. L. Bowie, Solution of crack problems by mapping technique, in Mechanics of Fracture, Vol. 1, pp. I-55. Noordhoff, Leyden (1978). (Received 7 July 1993)
0013-7944(94)EOO88-X
Printed
49. No.
Vol.
Copyright
Pergamon
‘c
in Great
I, pp.49-60.
1994 Else&r Britain.
All
0013.7944/94
1994
Science Ltd rights
reserved
$7.00 + 0.00
DETERMINING THE STRESS INTENSITY FACTORS K,, K,, AND THE T-TERM VIA THE CONSERVATION LAWS USING THE BOUNDARY ELEMENT METHOD P. C. OLSEN Department
of Structural
Engineering, Lyngby,
Technical Denmark
University
of Denmark,
2800.
Abstract-A new method for determining the stress intensity factors K, and K,, is presented. The method is based on the conservation laws, from which the property that the J,, J, and L integrals are zero for any closed path of integration is utilized. The integrals above are evaluated numerically for a circular path surrounding the crack tip and along the crack faces, ending, however, a distance A from the tip. At the distance A, all fields are replaced by the singular first-order expansion series, including the T-term for r~,, Based on the J,, J, and L integrals, numerically evaluated as described above, equations are derived, from which both the stress intensity factors K, and K,,, as well as the T-term are determined. Numerical tests based on the boundary element method are presented, in which the accuracy of the method is demonstrated. The total energy release rate. represented by the J, integral, is typically determined with an accuracy of within l-2%. The distribution of the total energy release rate between mode I and mode II is, however, less accurately determined. The T-term is determined with a rather large uncertainty.
1. INTRODUCTION MANY possibilities exist for evaluating the total energy release rate due to a change in the crack surface geometry. The difficult part is deciding which portion of the energy release rate has to be subscribed to which mode. In mixed mode problems, it is essential to distinguish between the different modes. Evaluating the J integral due to Rice [l], one can easily determine the total energy release rate. Budiansky and Rice [2] later realised that the J integral is actually the first component of a vector and that all components can be interpreted as energy release rates, when a cavity or a flat crack is given a translation. Also in [2], the L and M integrals were interpreted as energy release rates for rotation and expansion, respectively, of cavities and flat cracks. Numerically it is rather attractive to determine the energy release rate via the J integral because it can be evaluated in a stress field, far away from the crack tip and thus avoiding the singular stresses. It is not necessary, numerically, to model the singular behaviour near the crack tip correctly. With a sufficiently fine mesh, standard shape functions can be applied. However, in using the second component of the J integral in order to determine both K, and K,, , terms involving near field stresses must be included, whereby the attractiveness of the procedure becomes less pronounced. Eischen [4] applied the J integrals in order to determine both K, and K,, , as well as the constant T-term of the first-order stress expansion. The finite element method was applied and both standard elements, as well as singular crack tip elements were used. Both K,, and the T-term were reported to be rather sensitive to the types of element used at the crack tip. In this paper the determination of both K, and K,,, as well as the T-term in plane isotropic elasticity by means of evaluating the J integrals and the L integral is demonstrated. The integrals have to be evaluated in near field stresses. However, the integrands are replaced by singular terms from the first-order expansion closest to the crack tip, where the numerical values are unreliable. The boundary element method is applied.
2. THE
CONSERVATION
LAWS
The conservation laws are the direct consequence of the basic equations of the theory of elasticity. Knowles and Sternberg [3] have shown that, when these equations are all satisfied, the 49
50
following
P. C. OLSEN
integrals:
(1)
(2)
(3)
(4)
vanish for arbitrary closed integration paths S surrounding a homogeneous singly connected isotropic elastic domain. In eqs (l-4) W is the strain energy density, n is the unit outward normal vector of S, u is the deformation vector, L is the first-order strain measure, p is the traction vector at S corresponding to the unit outward normal vector and finally, x and y are coordinates of points at S. The subscripts refer to components in a global Cartesian reference frame. The conservation laws are now applied for a path enclosing a crack tip as shown in Fig. 1. This path is composed of three parts as shown: firstly, an arbitrarily shaped curve starting at the crack face at a given distance from the crack tip, surrounding the crack tip and ending at the opposite crack face. In practice this curve is often chosen as a circle; secondly, a circle with centre at the crack tip and a small radius r surrounding the crack tip; and thirdly, lines coinciding with both crack faces, connecting the first two curves. The complete path is thus closed, resulting in the following equations being valid: Jc”’ + J.y
+ J;”
= 0
J,v
+
JTb ?
=
0
L Face +
LTip
=
0,
+
L Far +
where the superscripts this paper.
JFe ?
refer to the parts shown in Fig. 1. The M integral
Fig. 1. Particular
integration
path for calculating
the conservation
(3
(6)
(7) is not further
integrals.
applied
in
51
Determining K,, K,, and the T-term using BEM
3. J AND L AT THE Displacements and expansion series, i.e.
stresses
close
CRACK
to the crack
TIP
tip are assumed
-sin(8/2)(2
to follow
the first-order
+ cos(8/2)cos(38/2))
T
sin(e/2)cos(e/2)cos(38/2)
+
i c0s(e/2)(i - sin(B/2)sin(38/2))
0 ,
(9)
(i 0
j
in which G is the shear modulus. ICis equal to (3 - 4v) in plane strain and (3 - v)/(l + v) for plane stress (see [6]). K, and K,, are the stress intensity factors. Regarding the inclusion of the T-term, see for example Eftis [5]. Introducing these expressions into the J and L integrals, performing the integration for the circle surrounding the crack tip and evaluating the results in the limit of vanishing radius, the following results may be derived (see [7]):
(10)
4. J AND L AT THE The J and L integrals
CRACK
FACES
at the crack faces are given by: JFaCe .r
Jc”“=
0
‘(W+ s
LFace
=
=
-
(13) Wm)dr
(14)
0
A (W’
-
W-jr
&,
(15)
s 0
where W+ and Wm are the strain energy at the upper and lower faces, respectively. Assuming displacements and stresses to follow the first-order expansion series eqs (8) (9) at the crack faces close to the crack tip, the stress o,, along the crack faces is given by: T - 2 J(::i),
8 = z
(Upper)
(T,,.=
(16) 8 = -7t (Lower)
resulting
in a difference
in strain
w+_w-=-&
energy
between
the two crack faces:
=
4
K,, T E’J(2w)
’
(17)
P. C. OLSEN
52
E’ equals E, the modulus of elasticity in plane stress; in plane strain E’ = E/( 1 - v’) where 1’is the Poisson ratio. JX and L can now be evaluated explicitly, by introducing eq. (17) into eqs (14), (I 5): JF””
=
_
8
K1t
(18)
T&
E'JO~)
(19)
5. DETERMINING
K,,K,,AND THE T-TERM
Both the stress intensity factors K, and K,,, as well as the T-term can now be determined by means of the J and L integrals evaluated numerically in stress fields far from the crack tip. It is only a matter of simple manipulation to arrive at the following results:
(20) K,, = J(E’J, T =
- K;)
(21)
3E'JW W,fi
(22)
In eqs (20) and (21) the superscript “Far” of J,. J,. and L has been suppressed. Equations (20)-(22) are the basic equations to be thoroughly investigated next. In the case of K,, becoming zero, i.e. for pure mode I, the method above is not applicable. In this case the strain energy difference between the opposite crack faces is zero, which means that .lVand L are path independent in the sense of J, path independence. However, the author cannot convince himself as to whether the T-term must be zero in pure mode I or not.
6. THE
BOUNDARY
ELEMENT
METHOD
The numerical analysis is done by means of the boundary element method. There is no need to reiterate this method here and the reader is referred to the literature, see for example Brebbia [S]. When modelling flat cracks, the method considered most suitable is the multi-domain boundary element method (Blandford [9]), where the domain in front of each crack is divided up into subdomains, in order to prevent the singularities in the system matrix otherwise arising. When determining the J and L integrals, one needs to determine very precisely displacements and stresses close to the boundary, due to the subdivision. The boundary element method is not famous for being very accurate in determining stresses near the boundary. On the contrary, it is almost a standard rule not to determine stresses at internal points which are closer to the boundary than the length of an element. However, it is possible to achieve very accurate results at points even extremely close to the boundary, provided that integration is performed analytically. This demands that either the boundary is not curved or that the curvature is small enough so that the boundary can be approximated with straight line segments. The program being used supports two classes of elements, Lagrange elements and cubic spline elements. The cubic spline elements have been implemented in order to achieve continuity in the normal stresses, acting tangential to the boundary. This stress component is derived from the first derivatives of the boundary displacements, and thus differentiability is needed for the boundary displacements. The cubic spline elements satisfy this requirement internally. The same effect is achieved by using high-order Lagrange elements. However, high-order Lagrange polynomiums are very poor interpolants. The nodes in the cubic spline element may be positioned with variable where si is the spacing between spacing. A simple power law has been implemented, i.e. s,+ , = fsi, two consecutive nodes. Normally quadratic Lagrange elements are used throughout, except for at the crack faces and the elements adjacent to the crack tips, where the cubic spline elements are used.
Determining
K,, K,, and the T-term
53
using BEM
The partial cross derivatives au,/+ and &,./ax appearing in the conservation integrals are determined as finite differences, i.e. by determining the displacements at closely spaced points, calculating the differences and dividing by the distance. One can achieve very accurate results, which are easy to check, as the sum of the derivatives must equal twice the shear strain. However, it is advisable to perform the calculation by determining the kernels in the boundary element method corresponding to the cross derivatives (see Aliabadi and Rooke [6]), analytically, just as the displacements and stresses are determined analytically. Boundary stress discontinuities are modelled by means of the double node concept, i.e. two nodes are placed at boundary points, where stress discontinuities may be expected, which is the case at a crack tip. For double nodes, where unknown tractions appear at both nodes, discontinuous elements are used. 7. INTRODUCTORY
NUMERICAL
TEST
The example of a slanted crack subjected to a uniform pressure is used for the introductory numerical calculations. The geometry of the specimen is shown in Fig. 2. The dots indicate the nodal points of the boundary element discretization. In the enlargement, the nodal points close to the crack tip are shown. The nodal points along the crack face and at the internal boundaries next to the crack have been positioned corresponding to a factor f = 1.25 in the power law of nodal spacing. Cubic spline elements, each containing 11 nodes, are used at the crack face and at the internal boundaries. Quadratic Lagrange elements are used at the external boundaries. The length of half the crack is a = 0.1 m whereas the angle of the crack with the vertical is a = 60. For the sheet of infinite size, analytical solutions exist: K, = oj(na)
sin’s,
K,, = oJ(na)
cos a sin x
and
T = g cos 2~,
(23)
where CJ is the undisturbed far field stress (Eftis [5]). Firstly it has been tested as to what extent one can expect the conservation laws to be satisfied when integrating along a closed path not enclosing any defects, in which case the integrals must turn out to be zero. A circular path was chosen with a radius of r = 0.05 m and with the centre placed at two different positions Pl and P2 as shown in Fig. 2. The first position is characterized as having almost uniform stresses around the entire periphery of the circle, whereas the second position is so close to one of the crack tips that large variations in the stresses are expected. In Fig. 3 the variations of the integrands of J and L along the periphery are shown. In Fig. 3a, it is seen that for position Pl the variation is almost harmonic, indicating an almost constant stress field, for which the variation is solely due to the projection of the normal vector. For an integration path in a non-uniform stress field, like the one at position P2, the integrands vary irregularly as shown in Fig. 3b.
Mesh
Fig. 2. The example
of a slanted
crack
at
the
used for testing
crack
the algorithms
face
P. c‘. OLSEN
54
A)
Position
B)
Pl
-3.00 [ -180.00
Position
P2
I -90.00
90.00
Fig. 3. The variation
0.00
180.00
Angle
Angle of the integrands
of the conservation
integrals
at the positions
Pl and P2
In order to perform the integration, the curves above are interpolated by means of natural cubic splines with zero curvature at the interval ends. The spline is determined by means of 61 regularly spaced points. In Table 1 the values of the J and L integrals are shown. It is seen that the accuracy with which the J and L integrals are determined using the boundary element method is astonishingly high. It would be interesting to see if the finite element method can achieve the same degree of accuracy. However, to the author’s knowledge similar investigations for the finite element method have not been published. Having seen that the boundary element method is an excellent tool for determining the J and L integrals, the integration path is now positioned in order to determine the stress intensity factors, i.e. the centre of the circle is positioned at the crack tip, the periphery starting at one crack face and ending at the opposite crack face. In Fig. 4 the variations of the integrands of J,, Jr and L are plotted for a circle with r = 0.08 m at this position. In Table 2 the calculated values of J,, J, and L are shown for different radii of the integration circle. J, is, as it should be, practically independent of the choice of radius, with a deviation on the fourth decimal only. J? is path dependent, as has been discussed by others (Hermann and Hermann[7]). L is likewise path dependent, whereby path dependence means that J? and L depend on the choice of radius. To check path independency in the case of the integration path surrounding the crack tip, an integration path is now chosen, composed of (1) a circle with the centre positioned at the crack tip, the periphery starting at one crack surface and ending at the opposite crack surface and (2) two lines coinciding with the crack surfaces, starting at the points where the circle intersects the crack surfaces and ending a distance A from the crack tip (see Fig. 5). For a given distance A, J, and L must be independent of the choice of radius. In Table 3 the .l,, integral is tabulated for different radii and distances to the crack tip. It is seen that for a given distance A, TVis practically independent of the choice of radius. Similarly, the L integral is tabulated in Table 4. It is seen that path independence cannot be achieved. However, results improve the further away from the crack tip the integration path ends. It is also remarkable that the order of magnitude of the L integral Table 2. The J,, J, and L integrals for different radii in the case of a circle of integration surrounding the crack tip
Table 1. The J and L integrals determined for a closed path of integration positioned in a regular stress field (PI) and very near to a crack tip in a highly irregular stress field (P2) Position J., 4 L Pl P2
6.603 x 1O-6 2.807 x 10-j
-4.911 -8.629
x lo-* x IO-’
2.080 x IO-’ 2.319 x IO-’
Radius
J.
J.
0.04 0.06 0.08 0.10 0.12 0.14 0.16
0.24630 0.24622 0.24643 0.24637 0.24637 0.24641 0.24648
0.13826 0.12700 0.12120 0.11926 0.12123 0.12723 0.13763
L -0.0007071 -0.0012116 -0.0015833 -0.0017723 -0.0015023 -0.0006614 0.0009726
Determining K,, K,, and the T-term using BEM
55
Angle Fig. 4. The variations of the integrands of the conservation integrals for a circle positioned at the crack tip.
in this example is more than one hundred times smaller than that of the J integrals, although all integrands have the same order of magnitude. This of course demands a higher degree of accuracy.
8. EVALUATING
K,, K,, AND THE
Z--TERM
Remembering, that eqs (20)-(22) rely on a first-order stress expansion closest to the crack tip corresponding to a distance A, it is obvious that the smaller the distance A is, the better the approximation is and the more accurate the stress intensity factors and the T-term can be determined. However, if one comes too close to the crack tip, the numerical values become unreliable. This is especially true for the L integral, as can be seen in Table 4. In Fig. 6 the normal stresses acting tangential to the crack surface are plotted. It is seen that the numerical representation becomes unreliable at a distance of approximately 0.04m from the crack tips. This figure explains the very inaccurate results for the L integral for distances A less than approximately 0.04 m. The idea is therefore to choose a distance A as small as possible and still have reliable numerical data. In Table 5 eqs (20)-(22) are evaluated for A > 0.04 m, for which the L integral has been determined with reasonable accuracy. Also in this table, the effect of choosing different radii can be seen.
Fields represented expansion terms tance A next to
by singular over a disthe crack tip
Fig. 5. Integration path around the crack tip composed of a circle and lines along the crack surfaces.
56
P. C. OLSEN Table 3. Path mdependence
J, J, for drfferent radn and distances
of
to the crack
tip
Radius
A = 0.01
A = 0.02
A = 0.03
A = 0.04
A = 0.05
A = 0.06
A = 0.07
A = 0.08
0.04
0.1794x
0.06
0.17986
0.08 0.10 0.12 0.14 0.16
0. I7992 0.17992 0. I7997 0.17997 0.17923
0.15515 0. I5554 0. I5559 0. I5560 0. IS564 0.1 SS64 0.15490
0.14585 0.14623 0. I4628 0. I4629 0.14633 0.14634 0.14560
0. I381-6 0. 13X64 0. 13870 0.13x7 I 0.13875 0.13875 0.13x01
U.13178 0. I3 I84 0.1318s 0.13189 0.13189 0.131 I5
0. I2700 0.12706 0. I2706 0.12710 0.12711 0.12637
0.12364 0.12365 0.12369 0.12369 0.12295
0. I z IX 0.12121 0.12125 0.12125 0.12051
Table 4. Path independence
of L. L for different
radii and distances
to the crack
tip. Values of L are multiplied
_
by 100.
Radius
A = 0.01
A = 0.02
A = 0.03
A = 0.04
A = 0.05
A = 0.06
A = 0.07
A = 0.08
0.04 0.06 0.08 0.10 0.12 0.14 0.16
0.01393 0.02030 0.02354 0.02142 0.02654 0.03198 0.02724
-0.02157 -0.01520 -0.01196 -0.01408 -0.00896 -0.00353 - 0.00826
-0.0440x -0.03771 -0.03447 -0.03660 -0.03147 -0.02604 - 0.03078
-0.07071 - 0.06434 -0.061 IO -0.06322 -0.05810 -0.05267 - 0.05740
-0.09500 -0.09176 -0.09389 -0.08876 -0.08333 -0.08807
-0.12116 -0.11792 -0.12005 -0.11492 -0.10949 -0.11422
-0.14033 -0.14215 -0.13703 -0.1316 -0.13633
-0.15823 -0.16036 -0.15523 -0.14980 - 0.15454
The relative variations of the stress intensity factors determined in Table 5 are 1.5% for K, and 7.5% for K,, , respectively. The T-term is less accurately determined, with a relative variation of 13.5%. The analytical solution for the sheet of infinite size is [5] K, = 0.4204
K,, = 0.2427
T = -0.5.
(24)
The results obtained numerically deviate from these analytical values and can be explained partly by the fact that the sheet has a finite size, and partly by the fact that the support is completely fixed, preventing any movement. It is also observed that there is a difference in the total energy
Lower 4.00
-4.00
Upper
,
I
I I 0.04
I 0.00
I 0.08
Distance Fig. 6. The tangential Table 5. The stress intensity
normal
factors
along
K,, K,, and the T-term
A = 0.05 T
from
stresses
4
A = 0.06 ~ .~~~~ K,, T
I I 0.16
I 0.12
crack
I 0.20
tip
the upper and lower crack for different
surfaces
radii and distances
to the crack
A = 0.07
Radius
K,
41
4
41
0.06
0.4496
0.2100
-0.3805
0.4504
0.2082
-0.3721
-
-
0.08 0.10 0.12 0.14 0.16
0.4511 0.4502 0.4522 0.4542 0.4530
0.2071 0.2090 0.2047 0.2002 0.2031
-0.3725 -0.3777 -0.3645 -0.3499 -0.3646
0.4517 0.4509 0.4525 0.4542 0.4534
0.2059 0.2074 0.2039 0.2002 0.2024
-0.3663 -0.3701 -0.3604 -0.3498 -0.3610
0.4532 0.4525 0.4539 0.4553 0.4546
0.2026 0.2039 0.2010 0.1978 0.1995
tip
A = 0.08 T -0.3508 -0.3538 -0.3461 -0.3377 -0.3469
4
4,
-
-
0.4551 0.4545 0.4556 0.4568 0.4563
0.1984 0.1995 0.1969 0.1942 0.1955
T -0.3314 -0.3339 -0.3274 -0.3205 -0.3283
57
Determining K,, K,, and the T-term using BEM
Boundary
element
mesh
Normal
Displacements
Tangential
tractions
Fig. 7. Comparison with Eischen [4]. Boundary element mesh, displacements
normal
stresses
and stresses.
release rate G. The analytical values of eq. (24) result in G = 0.2356, whereas the numerical evaluation of J, results in a value of G = 0.2464, a result which is considered to be very accurate. 9. ADDITIONAL
RESULTS
Very few articles are concerned with the determination of the T-term. However, Eischen [4] developed a method to determine both K,,K,, ,as well as the T-term and closed his paper with a numerical example. This example deals with a rectangular sheet subject to uniform tension and penetrated by a crack, with a crack mouth starting at the left-hand side of the vertical boundary and extending upwards at an angle of 45”. The size of the sheet is 1 x 2 and the length of the crack is a = 0.5656. In Fig. 7 the boundary element representation of the sheet is shown. Also in this figure the distorted shape of the geometry, the traction normal to the boundary and the tangential normal stresses are shown. The circular path of integration is chosen with r = 0.2828, for which the variations of J,, JI and L are shown in Fig. 8(a). Integrating these curves, one gets J, = 4.323, JJ = - 3.922 and L = -0.06561. In Fig. 8(b) the tangential normal stresses along the crack faces are plotted. These stresses seem to be reliable at distances larger than approximately one-quarter of the crack length.
A)
Variation
of
integrands
B)
Stress
-5.04
‘oh
+ L o.io
’ 0.40
Distance
Angle
of the integrands
crack
Jx = 4.323 Jy =-3.775
Jx= 4.323 Jy=-3.922 L =-0.06561
Fig. 8. The variation
along
(A) and the stress distribution
I
I
I
0.30
0.40
0.40
from
along
crack
0. 0
tip
the crack
faces (B)
face
P. C. OLSEN Table 6. K,, K,, and T for varying values of A m the case of the example shown in Fig. 6. A 0.05657 0.08485 0.1131 0.1414 0.1697 0.1980 0.2263
K,
K,,
T
1.961 I .949 I.936 I.928
0.69 I2 0.7244 0.7589 0.7782 0.7948 0.8100 0.8240
-0.5890 ~ I .2444 ~ 0.9600 -0.8293 PO.7325 -0.6534 -0.5882
I.921 1.915 1.909
In Table 6 the variations of K,, K,, and T for varying distances A are shown. The most probable values are those obtained with the smallest A, for which the numerical stress representation is still reliable and the author concludes that K, = 1.936, K,, = 0.7589 and T = -0.96 are the most probable values. Eischen [4] obtained K, = 1.917, K,, = 0.8065 and T = 0.822, whereas Bowie [lo] obtained K, = 1.933 and K,, = 0.8265. The latter did not determine the T-term. It is very interesting that Eischen and the author agree upon the total energy release rate, both having a value of G = 4.325, whereas the energy release rate resulting from Bowie’s calculations is G = 4.420. On the other hand, Bowie and the author agree upon the mode I stress intensity factor, having almost identical results. Thus, if one accepts firstly that the total energy release rate has been correctly determined by Eischen and the author and secondly, that the identical values of K, obtained by Bowie and the author are correct, then K,, = 0.7589 as determined by the author must be the correct value for the mode II stress intensity factor. This article concludes with two examples concerned with the bending of beams. The beam in the first example is a four point bending problem and is a pure mode I fracture problem. Results for this beam have been presented by Blandford et al. [9]. The boundary element mesh is shown in Fig. 9. Also in the figure the distorted shape, the normal traction and the tangential normal stresses of the boundary element analysis are shown. The conservation integrals are determined for a circle surrounding the crack tip and having a radius r = 0.05. This results in J, = 0.3609 whereas JI and L are of an order of magnitude three smaller (J,. = -5.54 x 10. 5 and L = 4.74 x 10-j). The mode I SIF factor is thus determined as K, = 0.6007. Blandford et al. reports of a value K, = 3.818J(2.54/100) = 0.6085. The authors result is thus within an accuracy of 1.3%. A) Boundary
C) Normal
element
mesh.
traction.
B) Distorted
D) Tangential
Fig. 9. Mode I fracture
of beam.
shape
normal
stresses.
Determining
A) Boundary
C) Normal
element
K,, K,, and the T-term
mesh
B) Distorted
D) Tangential
traction.
Fig. 10. Mixed mode I and mode II fracture
Fig. 11. Shear stress distribution
along
vertical
59
using BEM
shape
normal
stresses.
of a beam
sections
of the beam.
Finally a beam with dimensions as before, however, with a notch placed at the quarter-point and with a point load at the centre is analysed. In Fig. 10 the boundary element mesh is shown. Also in this figure the distorted shape, the normal traction and the tangential normal stresses of the boundary element analysis are shown. The distribution of the shear stresses along vertical sections of the beam at different locations are shown in Fig. 11. The well known parabolic distribution is recognised at one end of the beam, whereas the shear stresses are wildly disturbed by the notch at the other end. It is also recognisable that the shear stresses behave singularly in front of the crack tip. The conservation integrals are determined for a circle surrounding the crack tip and having a radius r = 0.05. This results in J, = 0.2241, JY = -0.07755 and L = -0.0002011. Using a distance of A = 0.025, the correction amounts to JJ = -0.08071 and L = -0.0002917. Thus from eqns (20)-(22), KI = 0.4666, K,, = -0.07985 and T = -0.06455.
10. CONCLUSIONS A new method for deriving the stress intensity factors has been presented. The method is based on the conservation laws, from which the property that the J,, TVand L integrals are zero for any closed path of integration is utilized. The integrals above are evaluated numerically for a circular path surrounding the crack tip and along the crack faces from the point of intersection of the circle towards the crack tip, ending, however, a distance A from the tip. At the distance A, all fields are replaced by the singular first-order expansion series, including the T-term for c,,. Based on the J,, Jr and L integrals, numerically evaluated as described above, equations are derived, from which both the stress intensity factors K, and K,, , as well as the T-term are determined. Numerical tests, based on the boundary element method are presented, in which the accuracy of the method is demonstrated. The total energy release rate, represented by the J, integral is typically determined with an accuracy of within l-2%. The distribution of the total energy release rate between mode I and mode II is, however, less accurately determined. The T-term is determined with a rather large
60
P. C. OLSEN
uncertainty. This has to be seen as a result of the magnitude of the L integral, which is an order of two smaller than both J, and J, However, by carefully observing the quality of the numerical data, one may come to a conclusion with regard to the most probable values. From the numerical tests, it is concluded that the best results are achieved by using a radius of half the crack length for the circular integration path and by extending the integration along the crack faces to a distance A from the crack tip, corresponding to a quarter of the crack length.
REFERENCES [l] J. R. Rice, Mathematical analysis in the mechanics of fracture, in Fracture, An Advanced Treatise, Volume II. Mulhetical Fundamenfals (Edited by H. Liebowitz) Academic Press, New York (1968). [2] B. Budiansky and J. R. Rice, Conservation laws and energy-release rates. J. uppl. Mech. 40, 201-203 (1973). [3] J. K. Knowles and E. Sternberg, On a class of conservation laws in linearized and finite elastostatics. Arch. Rational Mech. Anal. 44, 187-211 (1972). [4] J. W. Eischen, An improved method for computing the J, integral. Engng Fracfure Mech. 26(S), 691-700 (1987). [S] J. Eftis, Load biaxiality and fracture: a two-sided history of complementing errors. Engng Fracture Mech. 26(4), 567-592 (I 987). [6] M. H. Aliabadi and D. P. Rooke, Numerical Fracture Mechanics. Computational Mechanics Publications/Kluwer Academic Publishers, Boston, MA (1991). [7] A. G. Hermann and G. Hermann, On energy-relase rates for a plane crack. J. uppI. Mech. 48, Q-528 (1981). [8] C. A. Brebbia. The Boundary Elemem Methodfor Engineers. PenTech Press (1978). [9] G. E. Blandford, A. E. IngratTea and J. A. Liggett, Two-dimensional stress intensity factor computations using the boundary element method. Int. J. numer. Methods Engng 17, 3877404 (1981). [IO] 0. L. Bowie, Solution of crack problems by mapping technique, in Mechanics of Fracture, Vol. 1, pp. I-55. Noordhoff, Leyden (1978). (Received 7 July 1993)