- Email: [email protected]
A path independent integral for the crack extension force of the circular arc crack
Engineering Fracture Mechanics 66 (2000) 423±439
www.elsevier.com/locate/engfracmech
A path independent integral for the crack extension force of the circular arc crack M. Lorentzon*, K. Eriksson Department of Solid Mechanics, Lulea University of Technology, SE-971 87 Luleas, Sweden Received 25 October 1999; received in revised form 31 January 2000; accepted 25 February 2000
Abstract A path independent integral expression for the crack extension force of a two-dimensional circular arc crack is presented. The integral expression, which consists of a contour and an area integral, is derived from the principle of virtual work. It is implemented into a FEM post-processing program and the crack extension force is calculated for a circular arc crack in a linear elastic material. Comparison with exact solutions by Cotterell and Rice [10] for the eective elastic stress intensity factor shows acceptable accuracy for the numerical procedure used. 7 2000 Elsevier Science Ltd. All rights reserved. Keywords: Path-independent integral; Principle of virtual work; Circular arc crack; Crack extension force; Stress intensity factor; FEM-analysis
1. Introduction
Two dierent methods, the J-integral and the virtual crack extension method are at present widely used for calculating the energy release rate or crack extension force for a body of a rate independent material with a stationary crack. The J-integral, which was introduced by Rice [1] in 1968, is the basis of non-linear fracture mechanics theory. The virtual crack extension method was independently formulated by Parks [2] and Hellen [3] for the ®nite element method and was the ®rst approach to evaluate J in three-dimensional problems. The domain integral method, which is very similar to the virtual crack extension method, was developed by Shih et al. [4,5] and can be applied to both quasistatic and dynamic problems with elastic, plastic, or visco-plastic material response. However, these * Corresponding author. Tel.: +46-920-91977; fax: +46-920-91047. E-mail address: [email protected] (M. Lorentzon). 0013-7944/00/$ - see front matter 7 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 3 - 7 9 4 4 ( 0 0 ) 0 0 0 3 0 - 8
424
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
methods can only be used to predict initiation of growth of a crack that is straight in the direction of crack extension and do not allow a crack of a general shape. Eshelby [6] discovered that the force on a point defect in an elastic body can be determined with a path independent integral. The path independence property follows from zero divergence of Eshelby's energy±momentum tensor for a homogeneous material in a region, free from singularities. Later, Rice discovered, independently, the path independent J-integral, whose integrand is a normal component of Eshelby's tensor. The particular property that the J-integral, is zero on stress-free crack borders allows application to crack problems. Since the introduction of the J-integral much interest has been focused on path independent integrals in fracture mechanics applications. Only a few attempts have been made to formulate conservation integrals for cracks of a general shape, e.g. [7±8], and no general approach to this problem has been presented in the literature so far. As the divergence of Eshelby's energy±momentum tensor vanishes as well in admissible curvilinear coordinate systems as in Cartesian coordinates, it is reasonable to assume that path independent integrals exist for a crack of a general shape. In this paper, we derive a path independent integral expression of the crack extension force for a curved crack directly from a general energy balance equation, namely, the principle of virtual work. In our opinion, this principle is perhaps the most straightforward to obtain a path independent integral for a curved crack and at the same time capable of yielding a crack extension force. The method presented is fairly general and can be applied to dierent crack geometries. In [9] the conical crack is treated and further applications are under way. In this work attention is focused on the two-dimensional circular arc crack. The derived integral expression is implemented into a separate FEM post-processing program in order to verify numerical path independence and accuracy to calculate the crack extension force. Some of the major steps of the implementation procedure are brie¯y described below. The results obtained are compared with an exact analytical solution of the eective stress intensity factor for a circular arc crack by Cotterell and Rice [10]. 2. A path independent integral expression Consider an arbitrary homogeneous plane body of area A enclosed by a contour G which is taken positive counterclockwise, Fig. 1. Let dA be a surface element and dG an element of arc. Superscripts denote contravariant and subscripts covariant tensor properties. We will employ the usual summation
Fig. 1. An arbitrary plane homogeneous body of area A enclosed by a contour G:
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
425
convention for repeated indices, unless stated otherwise. The body is free of body forces and in static equilibrium under the action of speci®ed surface tractions ti. Assume further a continuous, small deformation displacement ®eld which does not violate any of the constraints of the problem under consideration. From the principle of virtual work, which is valid for any constitutive law, we know that the virtual work done by the surface tractions ti is equal to the corresponding internal virtual work. We can write this principle ti dui dG dW dA 1 G
A
where dui is a virtual displacement and dW virtual work per unit volume. Thus the integral dI dW dA ÿ ti dui dG A
G
2
is zero for any regular region, i.e. a region free from singularities. The strategy is now to reformulate the contour integral in Eq. (1) by using ordinary continuum mechanics assumptions and tensor analysis. Some of the details are omitted for brevity. Let yi , i 1, 2 be general curvilinear coordinates in the plane. The virtual displacement dui is chosen as minus the total dierential of the actual displacement ui (i.e. the solution to the problem studied), dui ÿui, k dyk
3
where the subscript `comma' (,) denotes covariant dierentiation. The increment dyk can be arbitrarily chosen and allows modelling of crack extension in the direction of increasing yi and also, conveniently, variation of the position of the crack tip. Substitution of the virtual displacement (3) and Cauchy's formula ti tij nj on the left-hand side of the virtual work balance Eq. (1), results in ÿ tij ui, k ÿ dyk nj dG 4 G
where tij is the stress tensor and nj the outward positive unit normal of G: By applying the divergence theorem and the equations of equilibrium, tij,j 0, we convert the contour integral (4) into an area integral as follows 5 ÿ tij ui, kj dyk dA ÿ tij ui, k dyk,j dA A
A
Let us now assume that there exists a strain energy function W eij with the property tij
@W , @ eij
6
where eij ui; j uj; i =2 is the strain tensor. Because of symmetry of the stress tensor and that the second covariant derivative is commutative ui, kj ui, jk ), the ®rst area integral in Eq. (5) can be written as ÿ @W @W k eij, k ÿ dyk dA ÿ dy dA 7 k @e ij A A @y Finally, substituting Eq. (7) into (5) and (4) into the left-hand side of Eq. (1) yields
426
dI
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
@W
A @y
dyk dA k
G
ti ui, k dyk dG
A
tij ui, k dyk,j dA
8
which is a general path independent integral expression for a regular region in 2D that holds in curvilinear coordinates. However, so far we have not considered any particular crack geometry. In order to keep the mathematical treatment of the path independent integral expression (8) as simple as possible, we select a crack with a fairly simple geometry. In the next section, we derive a path independent integral for a two-dimensional circular arc crack in polar coordinates. 2.1. The circular arc crack Consider a circular arc crack in polar coordinates r, a at r rc and ja ÿ p=2jRC where rc and C are constants (Fig. 2). The curved crack borders are assumed traction-free. To model crack extension, the increment dr, da 0, rc dac =r
9
is used, where rc dac is the virtual displacement of the crack tip. The components of the metric tensor gij are grr 1,
gaa r 2
10a,b
and all other gij 0: Below, gij is the associated metric tensor and in orthogonal coordinates gij 1=gij : The only non-zero Christoel symbols are Graa ÿr,
1 Gara Gaar : r
11a,b
p p By using the metric scaling factors gi i or gi i , i not summed, we can relate tensor components to their corresponding physical components, which are here written as `non-italics' sij , ui , etc.). For each term of the path independent integral expression (8) we have @W k 1 @W rda drr da dy dA Wna rda dG 12a k A @y A r @a G
Fig. 2. A circular arc crack in polar coordinates r, a:
G
A
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
ti ui, k dyk dG
G
tia ui, r dyr, a dA
pp p p ti gi i ui, k gi i gk k dyk gk k dG
G
ti ui, a dar dG
1 pp pp pp tia gi i ga a ui, r gi i gr r dyr, a gr r ga a dA sia ui; r ÿ rda dA r A A
427
12b
12c
where a subscript semicolon (;) denotes the physical component of the covariant derivative. Finally, working out the covariant derivative in (8) using the Christoel symbol (11b) and the identities (12(a)± (c)), yield after simpli®cation, the path independent integral expression in physical components ÿ 1 13 dI ÿ Wna ÿ ti ui; a dG rda sia ui; r dA rda r G A Consider the energy variation dI 0 in association with the increment da, (9). By taking rc dac outside the integrals, we obtain ÿ 1 14 Wna ÿ ti ui; a dG ÿ sia ui; r dA rc dac 0, dI ÿ r G A where r > 0 for all points on G: As rc dac is non-zero, it follows that the expression within the square bracket vanishes for any regular region enclosed by the contour G: Path independence of the integral expression in (14) was veri®ed analytically in some well-known situations: (a) a concentrated force at a point on a straight boundary, Flamant's solution [11], (b) a force at a point of an in®nite plate, Timoshenko's solution [12], and (c) a circular hole in an in®nite plate under remote uniaxial stress, Kirch's solution [13]. However, for the Flamant and Timoshenko solutions the area integral vanishes identically. Only in the last case the stress, strain and displacement ®elds are complex enough to yield a non-trivial result for the area integral. The contour G has four parts, so that G G0 Ge Ga Gb , where G0 and Ge are the contours, remote and in®nitely close to the crack tip, respectively, Ga and Gb , the paths along the upper and lower crack borders, respectively (Fig. 3). In the absence of crack border tractions ti and the fact that na 0, the contour integral in (14) is zero on Ga and Gb : Thus for a closed contour G, not containing a crack tip, we have
Fig. 3. A closed contour G is formed of Ga , Gb , G0 and Ge : rc dac is the virtual displacement of the crack tip.
428
Ge
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
ÿ
Wna ÿ ti ui; a dG
G0
ÿ Wna ÿ ti ui; a dG ÿ
A0
sia ui; r
1 dA 0 r
When Ge is closed around the crack tip and shrunk towards zero, we can identify ÿ Wna ÿ ti ui; a dG F limGe 4 0 ÿGe
15
16
as the crack extension force per unit crack length. Note that Ge is taken positive clockwise. It follows from (2) that dI 0, is the variation of the total energy of the domain considered, and F thus a force associated to the virtual displacement of the crack tip rc dac : From (15) and (16) we obtain ÿ 1 sia ui; r dA 17 F Wna ÿ ti ui; a dG ÿ r G0 A0 As the contour integral is zero on the crack borders, G0 can be closed around the crack tip. Because of path independence, any arbitrary path G0 and enclosed area A0 will yield the same value of F. This means that the strength of the singularity can be evaluated on and inside a path at any distance from the crack tip. The area integral can be regarded as a correction term for crack curvature. When the crack radius rc tends to in®nity, the contribution from the area integral term vanishes, and (17) is identical to Rice's well-known J-integral in two dimensions [1].
3. FEM post-processing routine The ®nite element method is used to evaluate the crack extension force F for the circular arc crack (17). The contour and area integrals are evaluated numerically in a separately developed FEM postprocessing program. The main part of the program comprises two subroutines; one for the contour integral and the other for the area integral. All input data to the post-processing routine are obtained from the ®nite element analysis output database. The numerical technique follows essentially a standard procedure [14]. Gaussian 3 3 quadrature is used for element integration in both subroutines and is performed in the natural coordinate system x, Z associated with an element. The total value of F is obtained by summing the contributions from all samplings points of all elements on the integration path and, when appropriate, in the integration domain. 3.1. The contour integral The contour integral subroutine evaluates the contour integral in (17) on the integration path G0 (Fig. 3). The discretized form of the contour integral term is given by FGo
X elements in G
3 X ÿ k1
Wna ÿ tr ur; a ta ua; a S k wk
18
where S is the arc length scaling factor. The quantities within the square brackets are evaluated along an integration path with either x = constant, or Z = constant within an element and at the Gaussian sampling points xp ,Zp and wk is the corresponding weight. All output data from the ®nite element analysis are in Cartesian form. The strain energy density W is a scalar quantity, and therefore, invariant. The unit normal vector na , the traction vector ti and displacement components ui in polar coordinates
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
429
are obtained from the corresponding Cartesian components trough tensor transformation and reduction to physical components. The contribution from a corner element in a rectangular mesh requires scaling. Two typical cases for a rectangular element are shown in Fig. 4. The scaling factor sc for case (a) and (b) is de®ned as sc 1 Pg =2 and sc 1 ÿ Pg =2, respectively, where Pg is the natural coordinate of the corresponding Gaussian point. The element contribution is evaluated as described above by integrating across the entire element according to (18) and then scaled to account for reduced path length. 3.2. The area integral The area integral subroutine evaluates the area integral in (17) in the integration domain A0. The discretized form of the area integral is 9 X X 1 19 FA 0 sra ur; r saa ua; r J qk r k elements in A k1 0
where J is the Jacobian. The quantities within the square brackets are evaluated at the nine Gaussian sampling points of an element, and qk is the corresponding weight. The stress and displacement components are obtained from the corresponding Cartesian components as described above. The outermost elements in the integration domain requires scaling. Two cases appears for a rectangular element, Fig. 5. The scaling factor used for case (a) and (b) is de®ned as sc Ae ÿ Aex =Ae and sc Ae ÿ Aex1 Aex2 =Ae , respectively, where Ae is the total area of the element and Aex is the area of the element outside the integration path. An element contribution is obtained from (19) and then scaled prior to summation. 4. Veri®cation of the F-integral routine In order to verify numerical accuracy, the contour integral routine has been applied to a standard linear elastic fracture mechanics problem with a straight crack, the double edge notched strip in tension. A general purpose ®nite element program, ABAQUS 5.7 [15], is used for the stress and displacement
Fig. 4. Contour paths through corner elements. The Gaussian points are marked with dark quadratic symbols.
430
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
analysis. The results obtained are compared with the analytical solution p and also, with ABAQUS inherent domain integral method routine. To obtain the desired 1= r strain singularity in crack problems, all nodes in the crack tip node set are tied together using ABAQUS inherent multi-point constraints. On element edges radial to the crack tip, the mid-side nodes are moved to the 1/4 point position. A 3 3 Gaussian quadrature is used in all calculations. The static stress analysis is performed in one single step with automatic time incrementation. 4.1. Analysis of the straight crack The stress intensity factor KI for the double edge notched strip subjected to tension is according to the standard solution p 20 KI s0 paf a=W where s0 is remote tensile stress, a crack length, W width of strip and f a=W a geometry function [16]. The crack extension force JKI , which is related to KI through the Rice±Irwin relation JKI
KI2 1 ÿ v 2 , E
21
where E is Young's modulus and n Poisson's ratio, is used to obtain a reference value of the standard solution. The ®nite element mesh of one quarter of the double-edge-notched tensile (DENT) specimen is composed of 88 eight-node iso-parametric elements (Fig. 6(a) and (b)). The ratio of the radial length of the smallest element to the width of the entire model is less than 0.01. The total number of degrees of freedom is 618. The load is far-®eld tension applied to the boundaries of the model (Fig. 7a). Eight dierent integration contour/domains are chosen to compute the crack extension force numerically for the F-integral and the domain method J-integral. The result is shown in Table 1 and a comparison with the standard solution in Fig. 7. The F-integral and the J-integral values (domain method) are both slightly lower than the standard solution. A small oscillation is obtained for the second contour with the F-integral method. Since the ®nite element mesh is fairly coarse, this result is most probably due to the numerical solution being less sensitive to the singular behaviour which dominates in the crack tip region. It is well known that this eect becomes
Fig. 5. Contour paths through the outermost elements in a domain. Ae is the total area of the element and Aex1 Aex2 is the area of the element outside the integration path.
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
431
Fig. 6. Geometry of the problem: (a) double edge-cracked specimen (DENT) subjected to remote tension, a=W 0:5, h=W 10, (b) the crack tip region.
Fig. 7. Normalized stress intensity factor values KFI =KI and KD I =KI for the contours 1±8. (DENT specimen). Dark symbols are used for results obtained with the domain method and light symbols for the F-integral method.
432
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
Table 1 Energy release rate values: double-edge notched tensile specimen case; a=W 0:5, h=W 5, s0 100 MPa, E 210 GPa and p n 0:3: KI 2:0719e 08 (Pa m according to [16]a Contour
JD (N/m)
p KD I (Pa m)
DF (%)
F (N/m)
p KFI (Pa m)
DF (%)
1 2 3 4 5 6 7 8
1:8876e 05 1:8531e 05 1:8591e 05 1:8601e 05 1:8603e 05 1:8603e 05 1:8603e 05 1:8601e 05
2:0871e 08 2:0679e 08 2:0713e 08 2:0718e 08 2:0720e 08 2:0720e 08 2:0720e 08 2:0718e 08
0.7 0.2 0.0 0.0 0.0 0.0 0.0 0.0
1:8869e 05 1:7792e 05 1:8225e 05 1:8391e 05 1:8468e 05 1:8515e 05 1:8555e 05 1:8580e 05
2:0867e 08 2:0263e 08 2:0508e 08 2:0601e 08 2:0644e 08 2:0670e 08 2:0693e 08 2:0707e 08
0.7 2.2 1.0 0.6 0.4 0.2 0.1 0.1
Average 2±8b
1:8590e 05
2:0713e 08
0.0
1:8361e 05
2:0584e 08
0.7
a b
JD the domain integral routine, F the integral routine (17). Contour 1 is omitted from the average value calculations.
Fig. 8. A circular arc crack in a two-dimensional solid subjected to a biaxial stress ®eld. rc 0:12 m and b 458:
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
433
smaller as the mesh is re®ned. However, path independence is still well preserved with both methods. The maximum obtained dierence from the standard solution is typically less than 0.7% for both methods, and has been judged sucient, as far as the F-integral is concerned, for further studies. 4.2. Analysis of the circular arc crack The problem analysed is a circular arc crack in a remote biaxial stress ®eld. The borders of the curved crack are traction free. Due to the symmetry of the problem, only one-half of the body is considered (Fig. 8). The boundary conditions for the ®nite element model are: u0
for x 0,
ÿ 10rc RyR10rc ,
22a
syy s1
for 0Rx < 10rc , y 210rc ,
22b
sxx s1
for x 10rc ,
22c
and ÿ 10rc RyR10rc
where x and y are Cartesian coordinates and u displacement in the x-direction. The ®nite element model is composed of 3040 eight-node isoparametric elements and the total number of degrees of freedom is 9773. In this case eight-node elements allow modelling of curved crack borders. The radial length D of the elements adjacent to the crack tip is D rc =24: Seventeen dierent integration domains are chosen to compute the crack extension force F numerically (Fig. 9). The contours 1±7 are ``symmetrically'' and 8±11 asymmetrically located around the
Fig. 9. De®nition of the integration domains 1±17.
434
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
crack. The contours 12±17 are arranged in order to gradually increase the in¯uence of the crack curvature. The contour 17 encloses the entire half-crack. The computed results are compared with the exact analytical solution for the eective stress intensity factor Ke obtained from the solution for KI and KII by Cotterell and Rice [10]. For a circular arc crack under a uniform biaxial stress ®eld we have: # (" syy sxx syy ÿ sxx cos 2 b=2 1=2 2 2 ÿ sin b=2cos b=2 KI pa 2 2 1 sin 2 b=2 ) syy ÿ sxx 2 3 cos 3b=2 ÿ sxy sin 3b=2 sin b=2 2
and
(" 1=2
KII pa
syy sxx 2
23a
# syy ÿ sxx sin b=2 2 2 ÿ sin b=2cos b=2 2 1 sin 2 b=2
) syy ÿ sxx 2 sin 3b=2 sxy cos 3b=2 cos b=2sin b=2 2
23b
where the crack length a and the angle b are de®ned in Fig. 8. The eective stress intensity factor is de®ned as q 24 Keff KI2 KII2 The corresponding eective stress intensity factor KFeff calculated from the F-integral using the Rice± Irwin relationship is
Fig. 10. Comparison of the normalized stress intensity factor values KFeff =Keff and KD eff =Keff for the circular arc crack, contours 1±7. Dark symbols: the domain method, light symbols: the F-integral method.
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
435
Table 2 Eective stress intensity factor values: circular arc crack; rc 0:12 m, b 458, s1 1000 MPa, E 210 GPa and n 0:3: p Keff 4:5035e 08 Pa m according to [10]a Contour
JD (N/m)a
p KD e (Pa m)
DD (%)
F (N/m)
p KFe (Pa m)
DF (%)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
8:8763e 05 8:8716e 05 8:8662e 05 8:8600e 05 8:8531e 05 8:8455e 05 8:8373e 05 ± ± ± ± ± ± ± ± ± ±
4:5259e 08 4:5247e 08 4:5233e 08 4:5217e 08 4:5200e 08 4:5180e 08 4:5159e 08 ± ± ± ± ± ± ± ± ± ±
0.5 0.5 0.4 0.4 0.4 0.3 0.3 ± ± ± ± ± ± ± ± ± ±
8:8556e 05 8:7767e 05 8:7366e 05 8:7093e 05 8:6323e 05 8:5841e 05 8:5203e 05 8:7143e 05 8:6760e 05 8:7151e 05 8:7161e 05 9:1406e 05 9:3562e 05 9:4824e 05 9:5485e 05 9:6002e 05 9:6195e 05
4:5206e 08 4:5004e 08 4:4901e 08 4:4831e 08 4:4633e 08 4:4508e 08 4:4342e 08 4:4844e 08 4:4745e 08 4:4846e 08 4:4849e 08 4:5928e 08 4:6466e 08 4:6779e 08 4:6941e 08 4:7068e 08 4:7116e 08
0.4 0.1 0.3 0.4 0.9 1.2 1.5 0.4 0.6 0.4 0.4 2.0 3.2 3.8 4.2 4.5 4.6
Averageb
8:8586e 05
4:5206e 08
0.4
8:9705e 05
4:5488e 08
1.0
a
JD the domain integral routine, F the integral routine according to (17). Contour 1 is omitted from the average value calculations. Contour 8±11 are radially asymmetric integration path with respect to the crack tip. b
Fig. 11. Normalized stress intensity factor values KFeff =Keff obtained with the F-integral method for the circular arc crack, contours 8±11.
436
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
KFeff
r FE , 1 ÿ n 2
25
which is used for comparison with Keff (24). 4.3. Numerical results KFe for all integration domains are shown in Table 2 and compared with Ke in Figs. 10±12. It is noted that there is a weak decreasing trend in KFe, for the contours 1±7 and a opposite trend for the contours 12±17. The values for the contours 12±17 are somewhat higher than those for the contours 1± 11. However, no oscillation in the value of KFe is obtained for the selected ®nite element discretization, and path independence is well preserved for all contours. Comparison of the results obtained and the analytical solution (24) indicates that the average dierence is of the order 1% and the maximum dierence is not greater than 4.6%, Table 2. 5. Discussion In this work a path independent integral expression of the crack extension force for the circular arc crack, composed of a contour integral and an area integral, is derived from the principle of virtual work. The analysis is fairly general and can be applied to curved cracks of various shapes. The increment dyk , which governs the virtual displacement ®eld, can be arbitrarily chosen and allows modelling of crack extension. The area integral arises due to the fact that the covariant derivative dyi,j is linear in dyk , so that in the present case the component dyr,a is non-zero in polar coordinates. The area integral can be regarded as a correction term due to crack curvature. A similar term appears for a plane crack with a curved crack front, e.g. the penny-shaped crack [17]. The area integral term tends to zero as the radius of curvature of the crack tends to in®nity and the crack itself approaches a straight crack. In the limit Rice's J-integral for a plane crack is obtained as a special case of the F-integral. The principle of virtual work is valid for any constitutive law and applies to all materials within the limitations of small deformation. However, in this work we limit our attention to linear or non-linear
Fig. 12. Normalized stress intensity factor values KFeff =Keff obtained with the F-integral method for the circular arc crack, contours 12±17.
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
437
elastic materials. That means that there exists a strain energy function with the property tij @W=@ eij and a unique stress±strain relationship. Under the restriction that unloading may not occur in any part of a body, a non-linear elastic behaviour of a material can be used to describe elastic±plastic behaviour. For a vanishingly small contour Ge around the crack tip we obtained the expression ÿ Wny ÿ ti ui; a dG, 26 F ÿlimGe 4 0 ÿGe
which can be interpreted as the crack extension force or the energy release rate for a linear or non-linear elastic material. The integrand in (16) is thus of order 1/rÄ (where rÄ is the radial distance from the crack tip) which is a necessary condition for F to be locally path independent and ®nite in the limit. This condition was pointed out by Moran and Shih [5] who applied a general balance law to derive a number of contour integrals for static and dynamic crack problems. In crack problems it is proper in certain situations to distinguish between remote boundary conditions and those on the crack borders. Remote boundary conditions may to some extent be only approximately satis®ed and still yield a crack extension force of reasonable accuracy, for instance in the case when the resultant load on a segment of the boundary is correct but not the exact stress distribution. On the other hand, it is essential for accuracy to satisfy the boundary conditions on the crack borders, particularly in the region close to the crack tip. All this follows from St. Venant's principle of disturbances decaying with distance. Thus, it is essential for a crack of a general shape to ful®l the boundary conditions on curved crack borders. On stress±free crack borders the contour integral should be zero unconditionally. For current formulations of the contour integral, it is not possible to ful®l this boundary condition in Cartesian coordinates for a curved crack. Satisfaction of the crack border boundary conditions for a curved crack requires a coordinate system in which the crack can be aligned with a suitable coordinate axis, in analogy with the Cartesian case. In the present solution of the crack extension force for the circular arc crack, the boundary condition of stress-free crack borders is satis®ed through the use of polar coordinates. 5.1. Numerical results Two plane crack problems, the double-edge-notched tensile (DENT) specimen, and the circular arc crack in a biaxial stress ®eld, were selected to verify the numerical accuracy and path independence of the derived expression (17). For the DENT specimen the J-integral was also calculated with the domain integral method. The J- and F-integrals were evaluated within and on similar, but not exactly identical contours. The J-integral contour enclosed entire elements and followed element borders and the Fintegral contour, a path through the outermost Gaussian sampling points of the integration area border elements [14]. Comparison of J and F with a standard solution [16] for a plane crack in a DENT specimen shows that the stress intensity factor KI can be determined within 0.7% with both routines. This is encouraging since the domain integral method is known in general to yield better accuracy than contour integrals. The eective stress intensity factor KFe for the circular arc crack obtained with the present solution was determined within 4.6% and less (mean 1%) of Cotterell's and Rice's [10] analytical solution. This numerical accuracy is in agreement with results from similar investigations using various methods for plane cracks [14]. Also, the greatest deviations were found for integration paths with some part close to the mesh border, as observed by Bakker [18], who also points out that better accuracy can be obtained with integration paths taken along element boundaries on the ground that a single element is the smallest unit in which equilibrium and other equations are satis®ed. By using entire elements the scaling procedure inherent in the method used [14], is avoided.
438
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
The Cotterell±Rice analytical solution is exact for a linear elastic material and its boundary conditions are ideally simple and applied at in®nity. The F-integral allows treatment of ®nite body, circular arc crack problems with boundary conditions of practical interest, i.e. loading or displacement prescribed on ®nite boundaries. The domain integral method, although strictly applicable to a straight crack only, also oers these features for a linear elastic material, if a vanishingly small integration domain around the crack tip is used. The F-integral on the other hand, allows a remote integration contour enclosing a ®nite part of a curved crack. The F-integral can thus be used to investigate for what portion of a curved crack where the straight crack the domain integral method still yields a reasonable result. Both methods can be used with non-linear materials. The major advantage of the domain integral method is that accurate results are usually obtained even with fairly coarse meshes. This feature was also observed in the present investigation, Table 1. Although de®ned only for a crack that is straight in the direction of crack extension, path independence is fairly well preserved for the innermost contours of the curved crack, inside which the crack is almost straight, Table 2. For the seven innermost and comparable but not identical integration paths, the domain integral J and F dier on the mean 0.4 and 0.7%, respectively, from the Cotterell±Rice solution. For contours enclosing a more curved part of the crack, it was not possible to calculate a J with the present mesh and ®nite element code. It has been reported that calculation of the J-integral in plane problems is less accurate for paths traversing the ®rst rings of element around a crack tip. Such results are usually due to the approximate nature of the numerical solution and are believed to indicate that the mesh should be re®ned in the crack tip region. The values of the F-integral obtained in the present investigation for all integration paths are very close to the analytical solution (24). Thus, the results indicate that the selected mesh and element topology and size were adequate for the present purpose. For a non-linear elastic material (and a deformation theory of plasticity material) loading is most strongly non-proportional close to the crack tip. Path independence is lost for an integration contour/ domain very near the crack tip and requires therefore a remote integration path, that may include such a portion of the crack that the domain integral method no longer applies. A remote integration path is not a problem with the F-integral because it is designed for curved cracks.
6. Conclusion A path independent integral expression for the crack extension force of a curved crack has been derived from the principle of virtual work. The derived integral for the crack extension force can conveniently be implemented into a FEM-post processing routine for numerical calculation. Results from FEM analyses show that the eective stress intensity factor of the circular arc crack can be determined within 4.6% from the exact analytical solution.
References [1] Rice JR. A path-independent integral and approximate analysis of strain concentration by notches and cracks. ASME Journal of Applied mechanics 1968;35:379±86. [2] Parks DM. A stiness derivative ®nite element technique for determination of crack tip stress intensity factors. International Journal of Fracture 1974;10:487±502. [3] Hellen TK. On the method of virtual crack extension. International Journal for Numerical Methods in Engineering 1975;9:187±207.
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
439
[4] Shih CF, Moran B, Nakamura T. Energy release rate along a three-dimensional crack front in a thermally stressed body. International Journal of Fracture 1986;30:79±102. [5] Moran B, Shih CF. A general treatment of crack tip contour integrals. International Journal of Fracture 1987;35:295±310. [6] Eshelby JD. The continuum theory of lattice defects. In: Seitz F, Turnbull D, editors. Solid State Physics, vol. 3. New York: Academic Press, 1956. [7] Beom HG, Earmme YY, Choi SY. Energy release rates for an arbitrarily curved interface crack. Mechanics of Materials 1994;18:195±204. [8] Choi SY, Beom HG, Earmme YY. A new conservation integral for an arbitrarily kinked interfacial crack. International Journal of Fracture 1996;80:45±57. [9] Eriksson K. The crack extension force of a curved crack derived from the principle of virtual work. International Journal of Fracture, in press. [10] Cotterell B, Rice JR. Slightly curved or kinked cracks. International Journal of Fracture 1980;16:155±69. [11] Flamant. Compt. Rend. 1892;114:1465. [12] Timoshenko SP, Goodier JN. Theory of elasticity. Tokyo: McGraw-Hill Kogakusha, 1970. [13] Kirch G. Verein Deutscher Ingenieure 1898;42. [14] Owen DRJ, Fawkes AJ. Engineering fracture mechanics: numerical methods and applications. Swansea: Pineridge Press Limited, 1983. [15] ABAQUS Users's Manual, Version 5.7. Providence, RI: Hibbit, Karlsson and Sorensen, Inc., 1997. [16] Wu X, Carlsson AJ. Weight functions and stress intensity factor solutions. Oxford: Pergamon Press, 1991. [17] Eriksson K. On the point-wise J-value of axisymmetric plane cracks. International Journal of Fracture 1998;91:L31±L36. [18] Bakker A. On the numerical evaluation of J in three dimensions. In: Sih GC, Sommer E, Dahl W, editors. Proceedings IC Applications of Fracture Mechanics to Materials and Structures. Martinus Nijho Publishers: The Hauge, Boston, Lancaster, 1984.
www.elsevier.com/locate/engfracmech
A path independent integral for the crack extension force of the circular arc crack M. Lorentzon*, K. Eriksson Department of Solid Mechanics, Lulea University of Technology, SE-971 87 Luleas, Sweden Received 25 October 1999; received in revised form 31 January 2000; accepted 25 February 2000
Abstract A path independent integral expression for the crack extension force of a two-dimensional circular arc crack is presented. The integral expression, which consists of a contour and an area integral, is derived from the principle of virtual work. It is implemented into a FEM post-processing program and the crack extension force is calculated for a circular arc crack in a linear elastic material. Comparison with exact solutions by Cotterell and Rice [10] for the eective elastic stress intensity factor shows acceptable accuracy for the numerical procedure used. 7 2000 Elsevier Science Ltd. All rights reserved. Keywords: Path-independent integral; Principle of virtual work; Circular arc crack; Crack extension force; Stress intensity factor; FEM-analysis
1. Introduction
Two dierent methods, the J-integral and the virtual crack extension method are at present widely used for calculating the energy release rate or crack extension force for a body of a rate independent material with a stationary crack. The J-integral, which was introduced by Rice [1] in 1968, is the basis of non-linear fracture mechanics theory. The virtual crack extension method was independently formulated by Parks [2] and Hellen [3] for the ®nite element method and was the ®rst approach to evaluate J in three-dimensional problems. The domain integral method, which is very similar to the virtual crack extension method, was developed by Shih et al. [4,5] and can be applied to both quasistatic and dynamic problems with elastic, plastic, or visco-plastic material response. However, these * Corresponding author. Tel.: +46-920-91977; fax: +46-920-91047. E-mail address: [email protected] (M. Lorentzon). 0013-7944/00/$ - see front matter 7 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 3 - 7 9 4 4 ( 0 0 ) 0 0 0 3 0 - 8
424
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
methods can only be used to predict initiation of growth of a crack that is straight in the direction of crack extension and do not allow a crack of a general shape. Eshelby [6] discovered that the force on a point defect in an elastic body can be determined with a path independent integral. The path independence property follows from zero divergence of Eshelby's energy±momentum tensor for a homogeneous material in a region, free from singularities. Later, Rice discovered, independently, the path independent J-integral, whose integrand is a normal component of Eshelby's tensor. The particular property that the J-integral, is zero on stress-free crack borders allows application to crack problems. Since the introduction of the J-integral much interest has been focused on path independent integrals in fracture mechanics applications. Only a few attempts have been made to formulate conservation integrals for cracks of a general shape, e.g. [7±8], and no general approach to this problem has been presented in the literature so far. As the divergence of Eshelby's energy±momentum tensor vanishes as well in admissible curvilinear coordinate systems as in Cartesian coordinates, it is reasonable to assume that path independent integrals exist for a crack of a general shape. In this paper, we derive a path independent integral expression of the crack extension force for a curved crack directly from a general energy balance equation, namely, the principle of virtual work. In our opinion, this principle is perhaps the most straightforward to obtain a path independent integral for a curved crack and at the same time capable of yielding a crack extension force. The method presented is fairly general and can be applied to dierent crack geometries. In [9] the conical crack is treated and further applications are under way. In this work attention is focused on the two-dimensional circular arc crack. The derived integral expression is implemented into a separate FEM post-processing program in order to verify numerical path independence and accuracy to calculate the crack extension force. Some of the major steps of the implementation procedure are brie¯y described below. The results obtained are compared with an exact analytical solution of the eective stress intensity factor for a circular arc crack by Cotterell and Rice [10]. 2. A path independent integral expression Consider an arbitrary homogeneous plane body of area A enclosed by a contour G which is taken positive counterclockwise, Fig. 1. Let dA be a surface element and dG an element of arc. Superscripts denote contravariant and subscripts covariant tensor properties. We will employ the usual summation
Fig. 1. An arbitrary plane homogeneous body of area A enclosed by a contour G:
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
425
convention for repeated indices, unless stated otherwise. The body is free of body forces and in static equilibrium under the action of speci®ed surface tractions ti. Assume further a continuous, small deformation displacement ®eld which does not violate any of the constraints of the problem under consideration. From the principle of virtual work, which is valid for any constitutive law, we know that the virtual work done by the surface tractions ti is equal to the corresponding internal virtual work. We can write this principle ti dui dG dW dA 1 G
A
where dui is a virtual displacement and dW virtual work per unit volume. Thus the integral dI dW dA ÿ ti dui dG A
G
2
is zero for any regular region, i.e. a region free from singularities. The strategy is now to reformulate the contour integral in Eq. (1) by using ordinary continuum mechanics assumptions and tensor analysis. Some of the details are omitted for brevity. Let yi , i 1, 2 be general curvilinear coordinates in the plane. The virtual displacement dui is chosen as minus the total dierential of the actual displacement ui (i.e. the solution to the problem studied), dui ÿui, k dyk
3
where the subscript `comma' (,) denotes covariant dierentiation. The increment dyk can be arbitrarily chosen and allows modelling of crack extension in the direction of increasing yi and also, conveniently, variation of the position of the crack tip. Substitution of the virtual displacement (3) and Cauchy's formula ti tij nj on the left-hand side of the virtual work balance Eq. (1), results in ÿ tij ui, k ÿ dyk nj dG 4 G
where tij is the stress tensor and nj the outward positive unit normal of G: By applying the divergence theorem and the equations of equilibrium, tij,j 0, we convert the contour integral (4) into an area integral as follows 5 ÿ tij ui, kj dyk dA ÿ tij ui, k dyk,j dA A
A
Let us now assume that there exists a strain energy function W eij with the property tij
@W , @ eij
6
where eij ui; j uj; i =2 is the strain tensor. Because of symmetry of the stress tensor and that the second covariant derivative is commutative ui, kj ui, jk ), the ®rst area integral in Eq. (5) can be written as ÿ @W @W k eij, k ÿ dyk dA ÿ dy dA 7 k @e ij A A @y Finally, substituting Eq. (7) into (5) and (4) into the left-hand side of Eq. (1) yields
426
dI
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
@W
A @y
dyk dA k
G
ti ui, k dyk dG
A
tij ui, k dyk,j dA
8
which is a general path independent integral expression for a regular region in 2D that holds in curvilinear coordinates. However, so far we have not considered any particular crack geometry. In order to keep the mathematical treatment of the path independent integral expression (8) as simple as possible, we select a crack with a fairly simple geometry. In the next section, we derive a path independent integral for a two-dimensional circular arc crack in polar coordinates. 2.1. The circular arc crack Consider a circular arc crack in polar coordinates r, a at r rc and ja ÿ p=2jRC where rc and C are constants (Fig. 2). The curved crack borders are assumed traction-free. To model crack extension, the increment dr, da 0, rc dac =r
9
is used, where rc dac is the virtual displacement of the crack tip. The components of the metric tensor gij are grr 1,
gaa r 2
10a,b
and all other gij 0: Below, gij is the associated metric tensor and in orthogonal coordinates gij 1=gij : The only non-zero Christoel symbols are Graa ÿr,
1 Gara Gaar : r
11a,b
p p By using the metric scaling factors gi i or gi i , i not summed, we can relate tensor components to their corresponding physical components, which are here written as `non-italics' sij , ui , etc.). For each term of the path independent integral expression (8) we have @W k 1 @W rda drr da dy dA Wna rda dG 12a k A @y A r @a G
Fig. 2. A circular arc crack in polar coordinates r, a:
G
A
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
ti ui, k dyk dG
G
tia ui, r dyr, a dA
pp p p ti gi i ui, k gi i gk k dyk gk k dG
G
ti ui, a dar dG
1 pp pp pp tia gi i ga a ui, r gi i gr r dyr, a gr r ga a dA sia ui; r ÿ rda dA r A A
427
12b
12c
where a subscript semicolon (;) denotes the physical component of the covariant derivative. Finally, working out the covariant derivative in (8) using the Christoel symbol (11b) and the identities (12(a)± (c)), yield after simpli®cation, the path independent integral expression in physical components ÿ 1 13 dI ÿ Wna ÿ ti ui; a dG rda sia ui; r dA rda r G A Consider the energy variation dI 0 in association with the increment da, (9). By taking rc dac outside the integrals, we obtain ÿ 1 14 Wna ÿ ti ui; a dG ÿ sia ui; r dA rc dac 0, dI ÿ r G A where r > 0 for all points on G: As rc dac is non-zero, it follows that the expression within the square bracket vanishes for any regular region enclosed by the contour G: Path independence of the integral expression in (14) was veri®ed analytically in some well-known situations: (a) a concentrated force at a point on a straight boundary, Flamant's solution [11], (b) a force at a point of an in®nite plate, Timoshenko's solution [12], and (c) a circular hole in an in®nite plate under remote uniaxial stress, Kirch's solution [13]. However, for the Flamant and Timoshenko solutions the area integral vanishes identically. Only in the last case the stress, strain and displacement ®elds are complex enough to yield a non-trivial result for the area integral. The contour G has four parts, so that G G0 Ge Ga Gb , where G0 and Ge are the contours, remote and in®nitely close to the crack tip, respectively, Ga and Gb , the paths along the upper and lower crack borders, respectively (Fig. 3). In the absence of crack border tractions ti and the fact that na 0, the contour integral in (14) is zero on Ga and Gb : Thus for a closed contour G, not containing a crack tip, we have
Fig. 3. A closed contour G is formed of Ga , Gb , G0 and Ge : rc dac is the virtual displacement of the crack tip.
428
Ge
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
ÿ
Wna ÿ ti ui; a dG
G0
ÿ Wna ÿ ti ui; a dG ÿ
A0
sia ui; r
1 dA 0 r
When Ge is closed around the crack tip and shrunk towards zero, we can identify ÿ Wna ÿ ti ui; a dG F limGe 4 0 ÿGe
15
16
as the crack extension force per unit crack length. Note that Ge is taken positive clockwise. It follows from (2) that dI 0, is the variation of the total energy of the domain considered, and F thus a force associated to the virtual displacement of the crack tip rc dac : From (15) and (16) we obtain ÿ 1 sia ui; r dA 17 F Wna ÿ ti ui; a dG ÿ r G0 A0 As the contour integral is zero on the crack borders, G0 can be closed around the crack tip. Because of path independence, any arbitrary path G0 and enclosed area A0 will yield the same value of F. This means that the strength of the singularity can be evaluated on and inside a path at any distance from the crack tip. The area integral can be regarded as a correction term for crack curvature. When the crack radius rc tends to in®nity, the contribution from the area integral term vanishes, and (17) is identical to Rice's well-known J-integral in two dimensions [1].
3. FEM post-processing routine The ®nite element method is used to evaluate the crack extension force F for the circular arc crack (17). The contour and area integrals are evaluated numerically in a separately developed FEM postprocessing program. The main part of the program comprises two subroutines; one for the contour integral and the other for the area integral. All input data to the post-processing routine are obtained from the ®nite element analysis output database. The numerical technique follows essentially a standard procedure [14]. Gaussian 3 3 quadrature is used for element integration in both subroutines and is performed in the natural coordinate system x, Z associated with an element. The total value of F is obtained by summing the contributions from all samplings points of all elements on the integration path and, when appropriate, in the integration domain. 3.1. The contour integral The contour integral subroutine evaluates the contour integral in (17) on the integration path G0 (Fig. 3). The discretized form of the contour integral term is given by FGo
X elements in G
3 X ÿ k1
Wna ÿ tr ur; a ta ua; a S k wk
18
where S is the arc length scaling factor. The quantities within the square brackets are evaluated along an integration path with either x = constant, or Z = constant within an element and at the Gaussian sampling points xp ,Zp and wk is the corresponding weight. All output data from the ®nite element analysis are in Cartesian form. The strain energy density W is a scalar quantity, and therefore, invariant. The unit normal vector na , the traction vector ti and displacement components ui in polar coordinates
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
429
are obtained from the corresponding Cartesian components trough tensor transformation and reduction to physical components. The contribution from a corner element in a rectangular mesh requires scaling. Two typical cases for a rectangular element are shown in Fig. 4. The scaling factor sc for case (a) and (b) is de®ned as sc 1 Pg =2 and sc 1 ÿ Pg =2, respectively, where Pg is the natural coordinate of the corresponding Gaussian point. The element contribution is evaluated as described above by integrating across the entire element according to (18) and then scaled to account for reduced path length. 3.2. The area integral The area integral subroutine evaluates the area integral in (17) in the integration domain A0. The discretized form of the area integral is 9 X X 1 19 FA 0 sra ur; r saa ua; r J qk r k elements in A k1 0
where J is the Jacobian. The quantities within the square brackets are evaluated at the nine Gaussian sampling points of an element, and qk is the corresponding weight. The stress and displacement components are obtained from the corresponding Cartesian components as described above. The outermost elements in the integration domain requires scaling. Two cases appears for a rectangular element, Fig. 5. The scaling factor used for case (a) and (b) is de®ned as sc Ae ÿ Aex =Ae and sc Ae ÿ Aex1 Aex2 =Ae , respectively, where Ae is the total area of the element and Aex is the area of the element outside the integration path. An element contribution is obtained from (19) and then scaled prior to summation. 4. Veri®cation of the F-integral routine In order to verify numerical accuracy, the contour integral routine has been applied to a standard linear elastic fracture mechanics problem with a straight crack, the double edge notched strip in tension. A general purpose ®nite element program, ABAQUS 5.7 [15], is used for the stress and displacement
Fig. 4. Contour paths through corner elements. The Gaussian points are marked with dark quadratic symbols.
430
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
analysis. The results obtained are compared with the analytical solution p and also, with ABAQUS inherent domain integral method routine. To obtain the desired 1= r strain singularity in crack problems, all nodes in the crack tip node set are tied together using ABAQUS inherent multi-point constraints. On element edges radial to the crack tip, the mid-side nodes are moved to the 1/4 point position. A 3 3 Gaussian quadrature is used in all calculations. The static stress analysis is performed in one single step with automatic time incrementation. 4.1. Analysis of the straight crack The stress intensity factor KI for the double edge notched strip subjected to tension is according to the standard solution p 20 KI s0 paf a=W where s0 is remote tensile stress, a crack length, W width of strip and f a=W a geometry function [16]. The crack extension force JKI , which is related to KI through the Rice±Irwin relation JKI
KI2 1 ÿ v 2 , E
21
where E is Young's modulus and n Poisson's ratio, is used to obtain a reference value of the standard solution. The ®nite element mesh of one quarter of the double-edge-notched tensile (DENT) specimen is composed of 88 eight-node iso-parametric elements (Fig. 6(a) and (b)). The ratio of the radial length of the smallest element to the width of the entire model is less than 0.01. The total number of degrees of freedom is 618. The load is far-®eld tension applied to the boundaries of the model (Fig. 7a). Eight dierent integration contour/domains are chosen to compute the crack extension force numerically for the F-integral and the domain method J-integral. The result is shown in Table 1 and a comparison with the standard solution in Fig. 7. The F-integral and the J-integral values (domain method) are both slightly lower than the standard solution. A small oscillation is obtained for the second contour with the F-integral method. Since the ®nite element mesh is fairly coarse, this result is most probably due to the numerical solution being less sensitive to the singular behaviour which dominates in the crack tip region. It is well known that this eect becomes
Fig. 5. Contour paths through the outermost elements in a domain. Ae is the total area of the element and Aex1 Aex2 is the area of the element outside the integration path.
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
431
Fig. 6. Geometry of the problem: (a) double edge-cracked specimen (DENT) subjected to remote tension, a=W 0:5, h=W 10, (b) the crack tip region.
Fig. 7. Normalized stress intensity factor values KFI =KI and KD I =KI for the contours 1±8. (DENT specimen). Dark symbols are used for results obtained with the domain method and light symbols for the F-integral method.
432
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
Table 1 Energy release rate values: double-edge notched tensile specimen case; a=W 0:5, h=W 5, s0 100 MPa, E 210 GPa and p n 0:3: KI 2:0719e 08 (Pa m according to [16]a Contour
JD (N/m)
p KD I (Pa m)
DF (%)
F (N/m)
p KFI (Pa m)
DF (%)
1 2 3 4 5 6 7 8
1:8876e 05 1:8531e 05 1:8591e 05 1:8601e 05 1:8603e 05 1:8603e 05 1:8603e 05 1:8601e 05
2:0871e 08 2:0679e 08 2:0713e 08 2:0718e 08 2:0720e 08 2:0720e 08 2:0720e 08 2:0718e 08
0.7 0.2 0.0 0.0 0.0 0.0 0.0 0.0
1:8869e 05 1:7792e 05 1:8225e 05 1:8391e 05 1:8468e 05 1:8515e 05 1:8555e 05 1:8580e 05
2:0867e 08 2:0263e 08 2:0508e 08 2:0601e 08 2:0644e 08 2:0670e 08 2:0693e 08 2:0707e 08
0.7 2.2 1.0 0.6 0.4 0.2 0.1 0.1
Average 2±8b
1:8590e 05
2:0713e 08
0.0
1:8361e 05
2:0584e 08
0.7
a b
JD the domain integral routine, F the integral routine (17). Contour 1 is omitted from the average value calculations.
Fig. 8. A circular arc crack in a two-dimensional solid subjected to a biaxial stress ®eld. rc 0:12 m and b 458:
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
433
smaller as the mesh is re®ned. However, path independence is still well preserved with both methods. The maximum obtained dierence from the standard solution is typically less than 0.7% for both methods, and has been judged sucient, as far as the F-integral is concerned, for further studies. 4.2. Analysis of the circular arc crack The problem analysed is a circular arc crack in a remote biaxial stress ®eld. The borders of the curved crack are traction free. Due to the symmetry of the problem, only one-half of the body is considered (Fig. 8). The boundary conditions for the ®nite element model are: u0
for x 0,
ÿ 10rc RyR10rc ,
22a
syy s1
for 0Rx < 10rc , y 210rc ,
22b
sxx s1
for x 10rc ,
22c
and ÿ 10rc RyR10rc
where x and y are Cartesian coordinates and u displacement in the x-direction. The ®nite element model is composed of 3040 eight-node isoparametric elements and the total number of degrees of freedom is 9773. In this case eight-node elements allow modelling of curved crack borders. The radial length D of the elements adjacent to the crack tip is D rc =24: Seventeen dierent integration domains are chosen to compute the crack extension force F numerically (Fig. 9). The contours 1±7 are ``symmetrically'' and 8±11 asymmetrically located around the
Fig. 9. De®nition of the integration domains 1±17.
434
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
crack. The contours 12±17 are arranged in order to gradually increase the in¯uence of the crack curvature. The contour 17 encloses the entire half-crack. The computed results are compared with the exact analytical solution for the eective stress intensity factor Ke obtained from the solution for KI and KII by Cotterell and Rice [10]. For a circular arc crack under a uniform biaxial stress ®eld we have: # (" syy sxx syy ÿ sxx cos 2 b=2 1=2 2 2 ÿ sin b=2cos b=2 KI pa 2 2 1 sin 2 b=2 ) syy ÿ sxx 2 3 cos 3b=2 ÿ sxy sin 3b=2 sin b=2 2
and
(" 1=2
KII pa
syy sxx 2
23a
# syy ÿ sxx sin b=2 2 2 ÿ sin b=2cos b=2 2 1 sin 2 b=2
) syy ÿ sxx 2 sin 3b=2 sxy cos 3b=2 cos b=2sin b=2 2
23b
where the crack length a and the angle b are de®ned in Fig. 8. The eective stress intensity factor is de®ned as q 24 Keff KI2 KII2 The corresponding eective stress intensity factor KFeff calculated from the F-integral using the Rice± Irwin relationship is
Fig. 10. Comparison of the normalized stress intensity factor values KFeff =Keff and KD eff =Keff for the circular arc crack, contours 1±7. Dark symbols: the domain method, light symbols: the F-integral method.
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
435
Table 2 Eective stress intensity factor values: circular arc crack; rc 0:12 m, b 458, s1 1000 MPa, E 210 GPa and n 0:3: p Keff 4:5035e 08 Pa m according to [10]a Contour
JD (N/m)a
p KD e (Pa m)
DD (%)
F (N/m)
p KFe (Pa m)
DF (%)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
8:8763e 05 8:8716e 05 8:8662e 05 8:8600e 05 8:8531e 05 8:8455e 05 8:8373e 05 ± ± ± ± ± ± ± ± ± ±
4:5259e 08 4:5247e 08 4:5233e 08 4:5217e 08 4:5200e 08 4:5180e 08 4:5159e 08 ± ± ± ± ± ± ± ± ± ±
0.5 0.5 0.4 0.4 0.4 0.3 0.3 ± ± ± ± ± ± ± ± ± ±
8:8556e 05 8:7767e 05 8:7366e 05 8:7093e 05 8:6323e 05 8:5841e 05 8:5203e 05 8:7143e 05 8:6760e 05 8:7151e 05 8:7161e 05 9:1406e 05 9:3562e 05 9:4824e 05 9:5485e 05 9:6002e 05 9:6195e 05
4:5206e 08 4:5004e 08 4:4901e 08 4:4831e 08 4:4633e 08 4:4508e 08 4:4342e 08 4:4844e 08 4:4745e 08 4:4846e 08 4:4849e 08 4:5928e 08 4:6466e 08 4:6779e 08 4:6941e 08 4:7068e 08 4:7116e 08
0.4 0.1 0.3 0.4 0.9 1.2 1.5 0.4 0.6 0.4 0.4 2.0 3.2 3.8 4.2 4.5 4.6
Averageb
8:8586e 05
4:5206e 08
0.4
8:9705e 05
4:5488e 08
1.0
a
JD the domain integral routine, F the integral routine according to (17). Contour 1 is omitted from the average value calculations. Contour 8±11 are radially asymmetric integration path with respect to the crack tip. b
Fig. 11. Normalized stress intensity factor values KFeff =Keff obtained with the F-integral method for the circular arc crack, contours 8±11.
436
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
KFeff
r FE , 1 ÿ n 2
25
which is used for comparison with Keff (24). 4.3. Numerical results KFe for all integration domains are shown in Table 2 and compared with Ke in Figs. 10±12. It is noted that there is a weak decreasing trend in KFe, for the contours 1±7 and a opposite trend for the contours 12±17. The values for the contours 12±17 are somewhat higher than those for the contours 1± 11. However, no oscillation in the value of KFe is obtained for the selected ®nite element discretization, and path independence is well preserved for all contours. Comparison of the results obtained and the analytical solution (24) indicates that the average dierence is of the order 1% and the maximum dierence is not greater than 4.6%, Table 2. 5. Discussion In this work a path independent integral expression of the crack extension force for the circular arc crack, composed of a contour integral and an area integral, is derived from the principle of virtual work. The analysis is fairly general and can be applied to curved cracks of various shapes. The increment dyk , which governs the virtual displacement ®eld, can be arbitrarily chosen and allows modelling of crack extension. The area integral arises due to the fact that the covariant derivative dyi,j is linear in dyk , so that in the present case the component dyr,a is non-zero in polar coordinates. The area integral can be regarded as a correction term due to crack curvature. A similar term appears for a plane crack with a curved crack front, e.g. the penny-shaped crack [17]. The area integral term tends to zero as the radius of curvature of the crack tends to in®nity and the crack itself approaches a straight crack. In the limit Rice's J-integral for a plane crack is obtained as a special case of the F-integral. The principle of virtual work is valid for any constitutive law and applies to all materials within the limitations of small deformation. However, in this work we limit our attention to linear or non-linear
Fig. 12. Normalized stress intensity factor values KFeff =Keff obtained with the F-integral method for the circular arc crack, contours 12±17.
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
437
elastic materials. That means that there exists a strain energy function with the property tij @W=@ eij and a unique stress±strain relationship. Under the restriction that unloading may not occur in any part of a body, a non-linear elastic behaviour of a material can be used to describe elastic±plastic behaviour. For a vanishingly small contour Ge around the crack tip we obtained the expression ÿ Wny ÿ ti ui; a dG, 26 F ÿlimGe 4 0 ÿGe
which can be interpreted as the crack extension force or the energy release rate for a linear or non-linear elastic material. The integrand in (16) is thus of order 1/rÄ (where rÄ is the radial distance from the crack tip) which is a necessary condition for F to be locally path independent and ®nite in the limit. This condition was pointed out by Moran and Shih [5] who applied a general balance law to derive a number of contour integrals for static and dynamic crack problems. In crack problems it is proper in certain situations to distinguish between remote boundary conditions and those on the crack borders. Remote boundary conditions may to some extent be only approximately satis®ed and still yield a crack extension force of reasonable accuracy, for instance in the case when the resultant load on a segment of the boundary is correct but not the exact stress distribution. On the other hand, it is essential for accuracy to satisfy the boundary conditions on the crack borders, particularly in the region close to the crack tip. All this follows from St. Venant's principle of disturbances decaying with distance. Thus, it is essential for a crack of a general shape to ful®l the boundary conditions on curved crack borders. On stress±free crack borders the contour integral should be zero unconditionally. For current formulations of the contour integral, it is not possible to ful®l this boundary condition in Cartesian coordinates for a curved crack. Satisfaction of the crack border boundary conditions for a curved crack requires a coordinate system in which the crack can be aligned with a suitable coordinate axis, in analogy with the Cartesian case. In the present solution of the crack extension force for the circular arc crack, the boundary condition of stress-free crack borders is satis®ed through the use of polar coordinates. 5.1. Numerical results Two plane crack problems, the double-edge-notched tensile (DENT) specimen, and the circular arc crack in a biaxial stress ®eld, were selected to verify the numerical accuracy and path independence of the derived expression (17). For the DENT specimen the J-integral was also calculated with the domain integral method. The J- and F-integrals were evaluated within and on similar, but not exactly identical contours. The J-integral contour enclosed entire elements and followed element borders and the Fintegral contour, a path through the outermost Gaussian sampling points of the integration area border elements [14]. Comparison of J and F with a standard solution [16] for a plane crack in a DENT specimen shows that the stress intensity factor KI can be determined within 0.7% with both routines. This is encouraging since the domain integral method is known in general to yield better accuracy than contour integrals. The eective stress intensity factor KFe for the circular arc crack obtained with the present solution was determined within 4.6% and less (mean 1%) of Cotterell's and Rice's [10] analytical solution. This numerical accuracy is in agreement with results from similar investigations using various methods for plane cracks [14]. Also, the greatest deviations were found for integration paths with some part close to the mesh border, as observed by Bakker [18], who also points out that better accuracy can be obtained with integration paths taken along element boundaries on the ground that a single element is the smallest unit in which equilibrium and other equations are satis®ed. By using entire elements the scaling procedure inherent in the method used [14], is avoided.
438
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
The Cotterell±Rice analytical solution is exact for a linear elastic material and its boundary conditions are ideally simple and applied at in®nity. The F-integral allows treatment of ®nite body, circular arc crack problems with boundary conditions of practical interest, i.e. loading or displacement prescribed on ®nite boundaries. The domain integral method, although strictly applicable to a straight crack only, also oers these features for a linear elastic material, if a vanishingly small integration domain around the crack tip is used. The F-integral on the other hand, allows a remote integration contour enclosing a ®nite part of a curved crack. The F-integral can thus be used to investigate for what portion of a curved crack where the straight crack the domain integral method still yields a reasonable result. Both methods can be used with non-linear materials. The major advantage of the domain integral method is that accurate results are usually obtained even with fairly coarse meshes. This feature was also observed in the present investigation, Table 1. Although de®ned only for a crack that is straight in the direction of crack extension, path independence is fairly well preserved for the innermost contours of the curved crack, inside which the crack is almost straight, Table 2. For the seven innermost and comparable but not identical integration paths, the domain integral J and F dier on the mean 0.4 and 0.7%, respectively, from the Cotterell±Rice solution. For contours enclosing a more curved part of the crack, it was not possible to calculate a J with the present mesh and ®nite element code. It has been reported that calculation of the J-integral in plane problems is less accurate for paths traversing the ®rst rings of element around a crack tip. Such results are usually due to the approximate nature of the numerical solution and are believed to indicate that the mesh should be re®ned in the crack tip region. The values of the F-integral obtained in the present investigation for all integration paths are very close to the analytical solution (24). Thus, the results indicate that the selected mesh and element topology and size were adequate for the present purpose. For a non-linear elastic material (and a deformation theory of plasticity material) loading is most strongly non-proportional close to the crack tip. Path independence is lost for an integration contour/ domain very near the crack tip and requires therefore a remote integration path, that may include such a portion of the crack that the domain integral method no longer applies. A remote integration path is not a problem with the F-integral because it is designed for curved cracks.
6. Conclusion A path independent integral expression for the crack extension force of a curved crack has been derived from the principle of virtual work. The derived integral for the crack extension force can conveniently be implemented into a FEM-post processing routine for numerical calculation. Results from FEM analyses show that the eective stress intensity factor of the circular arc crack can be determined within 4.6% from the exact analytical solution.
References [1] Rice JR. A path-independent integral and approximate analysis of strain concentration by notches and cracks. ASME Journal of Applied mechanics 1968;35:379±86. [2] Parks DM. A stiness derivative ®nite element technique for determination of crack tip stress intensity factors. International Journal of Fracture 1974;10:487±502. [3] Hellen TK. On the method of virtual crack extension. International Journal for Numerical Methods in Engineering 1975;9:187±207.
M. Lorentzon, K. Eriksson / Engineering Fracture Mechanics 66 (2000) 423±439
439
[4] Shih CF, Moran B, Nakamura T. Energy release rate along a three-dimensional crack front in a thermally stressed body. International Journal of Fracture 1986;30:79±102. [5] Moran B, Shih CF. A general treatment of crack tip contour integrals. International Journal of Fracture 1987;35:295±310. [6] Eshelby JD. The continuum theory of lattice defects. In: Seitz F, Turnbull D, editors. Solid State Physics, vol. 3. New York: Academic Press, 1956. [7] Beom HG, Earmme YY, Choi SY. Energy release rates for an arbitrarily curved interface crack. Mechanics of Materials 1994;18:195±204. [8] Choi SY, Beom HG, Earmme YY. A new conservation integral for an arbitrarily kinked interfacial crack. International Journal of Fracture 1996;80:45±57. [9] Eriksson K. The crack extension force of a curved crack derived from the principle of virtual work. International Journal of Fracture, in press. [10] Cotterell B, Rice JR. Slightly curved or kinked cracks. International Journal of Fracture 1980;16:155±69. [11] Flamant. Compt. Rend. 1892;114:1465. [12] Timoshenko SP, Goodier JN. Theory of elasticity. Tokyo: McGraw-Hill Kogakusha, 1970. [13] Kirch G. Verein Deutscher Ingenieure 1898;42. [14] Owen DRJ, Fawkes AJ. Engineering fracture mechanics: numerical methods and applications. Swansea: Pineridge Press Limited, 1983. [15] ABAQUS Users's Manual, Version 5.7. Providence, RI: Hibbit, Karlsson and Sorensen, Inc., 1997. [16] Wu X, Carlsson AJ. Weight functions and stress intensity factor solutions. Oxford: Pergamon Press, 1991. [17] Eriksson K. On the point-wise J-value of axisymmetric plane cracks. International Journal of Fracture 1998;91:L31±L36. [18] Bakker A. On the numerical evaluation of J in three dimensions. In: Sih GC, Sommer E, Dahl W, editors. Proceedings IC Applications of Fracture Mechanics to Materials and Structures. Martinus Nijho Publishers: The Hauge, Boston, Lancaster, 1984.