- Email: [email protected]
On the calculation of energy release rate for curved cracks of rubbery materials
FINITE ELEMENTS IN ANALYSIS A N D DESIGN
ELSEVIER
Finite Elements in Analysis and Design 17 (1994) 323-338
On the calculation of energy release rate for curved cracks of rubbery materials J.H. C h a n g Department of Civil Engineering, National Central University, Chungli, Taiwan, ROC Received January 1994; revised May 1994
Abstract
The energy release rate G for 2-D rubbery material problems with curved cracks is calculated with finite element solutions in this paper. Two approaches, a generalized domain integral method and the virtual crack extension method, are investigated. The generalized domain integral method is demonstrated to be "patch independent" and therefore a complicated finite element mesh adjacent to the crack tip area is not required.
I. Introduction
The overall strength of engineering structures relies on integrity of their components, which is essentially limited by the growth of cracks. For those materials with mechanical behavior remaining (nonlinearly) elastic at large strains, we often call them rubbery. Rubbery materials find applications in many structure components. Thus the prediction of crack initiation and growth of rubbery materials is of practical importance. Typically, two criteria for crack growth are used in fracture mechanics. One is critical stress criterion, the other is critical energy criterion. For problems involving nonlinear responses, the stress state around the crack tip area becomes rather complicated such that the use of stress criterion appears to be impractical. Therefore, the Energy Release Rate (ERR), usually denoted G, is often employed as a critical energy criterion for predicting the initiation and growth of cracks in nonlinear problems. The parameter, ERR, is defined by the concept of Griffith energy criterion for fracture states. For a mechanical system with a crack of length l, crack growth occurs only if the energy delivered by the system reaches the value required to form an additional crack of length dl, i.e. when the corresponding energy is up to G dl. The critical value of G is dependent upon a number of factors e.g. loading state, material property, and especially, crack geometry. Investigations and modifications on computation methods of ERR calculation for problems with cracks of general curved shape is hence necessary. 0168-874X/94/$07.00 © 1994 Elsevier Science B.V. All rights reserved SSDI 0 1 6 8 - 8 7 4 X ( 9 4 ) 0 0 0 5 0 - P
J.H. Chang/Finite Elements in Ana(vsis and Design 17 (1994) 323 338
324
Analytically, calculation of G for 2-D problems can be performed by a line integral, denoted J, which was first introduced by Eshelby [-1] for the problem of a cracked elastic body. By definition, J is an integral evaluated along a counterclockwise contour enclosing and shrinking onto the tip of a crack. For elastic media with traction-free crack surfaces, J is path independent, i.e. the value of integration remains unchanged along any shape of contour encircling the tip as long as the crack is straight inside the contour. This result is valid for both problems with small and large deformations. In the context of large deformation, path independence is valid when the crack remains straight inside the contour in its reference configuration. From the experimental results (e.g. [2] ), rubbery materials at large deformation state show substantially reversible elasticity except in the immediate neighborhood of the crack tip. The validity of the J-integral for these materials is thus presumed. However, for the case of curved cracks, the property of path independence no longer holds. Modifications on the formulation of J-integral need to be applied. In finite element calculations, two other approaches the virtual crack extension method (VCEM) and the domain integral method (DIM) - have been proposed as alternatives to direct application of J-integral. The first technique, VCEM, was introduced for problems of linear elasticity [3, 4] and deformation theory of plasticity [5], and later used for rubbery material problems [6 8]. The other method, DIM, performs a domain integral to determine ERR for deformation theory of plastic [9], thermoelastic [10], and rubbery material problems [11]. It should be noted that both methods were developed and considered for problems containing straight crack advancing along its original direction. Under this situation, they are equivalent being alternative interpretations of J-integral, as described in [-10]. One common advantage of both methods is patch independence, a counterpart of path independence of J. This characteristic eliminates the need for a complicated finite element modeling around the crack tip. The object of this paper is to develop a generalized formulation, based on the domain integral method, along with finite element solutions, for calculation of energy release rate of rubbery materials with curved cracks. The applicability of VCEM on curved cracks will be investigated, both theoretically and numerically. 2-D problems, including plane stress, and plane strain cases, are considered.
2. Curved cracks
Consider an elastic body (linear or nonlinear) in a 2-D field, containing a curved crack in its reference configuration (Fig. 1). We introduce a local coordinate system originating at the crack tip O, with the crack lying tangentially along the negative xl-direction. The body is subjected to a system of loads (not shown in the figure) and the crack will spread by an infinitesimal amount dl only if the corresponding value of G dl exceeds the critical surface free energy. The energy release rate can be evaluated by J-integral defined as G = J = lira
F+O
Wnl
--
aijn
j
ds,
(1)
where W is the strain energy density of the material, 6ij and u~ are the Cartesian components of the stress tensor and displacement vector, nj are the Cartesian components of outward unit vector
J.H. Chang~FiniteElementsin Analysisand Design17 (1994) 323 338
325
n
×2 :ii
x I
Fig. 1. A curved crack in its reference configuration with a local coordinate system.
Fig. 2. Path-independence of J for an elastic body with a straight crack.
normal to F, nl is the xa-component of nj, and s is the arc length along the contour. By definition, F is a counterclockwise contour encircling and shrinking onto the tip of the crack (this limiting case is not shown in Fig. 1). For the problems of large deformation, all the state variables are interpreted with respect to this specific reference configuration. When the crack surfaces are traction free and the material is h o m o g e n e o u s in the xa-direction, it can be proved (e.g. [6]) that the line integral can be executed along any other counterclockwise contour, Fi, which encloses a straight portion of the crack (Fig. 2), with the result remaining unchanged. This is termed "path independence". Note that the complicated stress behavior around the crack tip can be ignored in evaluating J-integral due to path independence. Modifications of J should apply in the case of curved cracks, when the straight portion can hardly be defined. If we choose an arbitrary contour Fi, as shown in Fig. 3 and from Eq. (1), J then becomes c3u~
J = fr, [ Wna - a~jnj(-~--x~x~)ldS + fc, +c2Wna ds.
(2)
In Eq. (2), two additional line integrals along Ca and C2, are included. A schematic diagram of the geometry is shown in Fig. 3. Although the path Fi can still be chosen arbitrarily, the extra curves C~ and C:, should both be terminated at the crack tip O. With this portion of line integrals, the singular behavior in the region near the crack tip is thus involved in the calculation.
3. Virtual crack extension method
V C E M was developed to calculate the energy difference between two configurations of finite element mesh with slightly altered crack length in their straight portions, as shown in Fig. 4.
J.H. Chang~Finite Elements in Analysis and Design 17 (1994) 323 338
326
//t--..\
!
Fig. 3. An arbitrary contour F~, along with two additional curves C1 and C2 along the crack, is used for the calculation.
(a)
(b)
Fi
I
---~
Fig. 4. (a) Finite element mesh before the crack advances. A band of elements are delimited for VCEM calculation. (b) Finite element mesh with a band of deformed elements.
Consider the mesh in Fig. 4(a) under a system of mechanical loads, with its discretized potential energy denoted by rc. The energy release rate is the decrease of x due to a unit advance of crack, i.e. the derivative of rc with respect to the crack length l, d~
~c~u
dl -
Ou ~I
+
c~x~x Ox
8l'
(3)
J.H. Chang~Finite Elements in Analysis and Design 17 (1994) 323-338
327
where the position and displacement vectors at the nodal points, x and u, are functions of I. Note that x is defined in the undeformed reference configuration. Under the state of equilibrium -
-
c~u
= 0,
(4)
Eq. (3) reduces to 0
-
drc dl
•rc t3x #x 81"
(5)
To compute the derivative of potential energy with respect to l, an imaginary extension of crack length is accommodated by rigidly translating all nodes on and within a contour Fi around crack tip O by an amount of Al along direction of the crack, with all other nodes remaining in their original positions (Fig. 4(b)). Thus the elements outside Fo retain both their initial positions and shapes, while those inside/'i remain unchanged in shape with a translation of Al. Due to this perturbation of mesh, Eq. (5) can then be evaluated. As verified in I-3] that the contour £o, which delimits the translating patch, can be chosen arbitrarily except for the requirements to be inside the body, enclose the crack tip and the straight portion of the crack. We then identify patch independence in this method as a counterpart of path independence in J-integral. This implies that results can be achieved by adopting a contour remote from the crack tip which will not depend directly on the finite element solutions in the immediate vicinity of the singularity. Checking feasibility of VCEM for curved cracks is necessary. For clarity, here we consider a special case in which a kinked crack is involved. For the configurations of finite element mesh in Figs. 5(a) and (b), all the nodes on and within Fi subject to a rigid translation by an amount Al along
(a)
(b)
ro
Fig. 5. (a) Original configuration of the finite element mesh for the virtual crack extension method. (b) New configuration after the curved crack extends along the xl direction.
J.H. Chanq/Finite Elements in Analysis and Design 17 (1994) 323 338
328
the x~ direction. Note that the altered mesh in Fig. 5(b) does not represent the proper geometry of the problem after crack extension along the x~-direction. It appears that further investigations on proper methods of ERR calculations for curved cracks is required.
4. Generalized domain integral method The domain integral m e t h o d was introduced as an alternative to calculate energy release rate for problems of straight cracks. Basically, an area integral is performed for 2-D problems, while a volume integral is calculated for 3-D problems. Again, we consider an elastic body in a twodimensional field containing a straight crack, with homogeneous material in the xl-direction, and no crack surface tractions (Fig. 2). The area integral, which evaluates the energy released from the mechanical system as the straight crack advances a unit length along the x 1-direction, is derived from the following path independent J-integral
( 0u; '~]ds.
(6)
The contour F; is usually chosen to be far from the crack tip so that a more accurate result can be ensured. In this integral method, we select an arbitrary exterior contour Fo and introduce a smooth weight function, q, which is unity on F i and vanishes on F o. By applying divergence theorem, Eq. (6) is then transformed into the following area integral: J =
-
W~)lj-
Gij~:XI~--
~ x i ] (](A '
(7)
where da is the infinitesimal integration area, 6aj is Kronecker delta with its first c o m p o n e n t as 1. As shown in Fig. 6, the integral is carried out in the shaded region A. This region A is often chosen to coincide with the original finite element mesh so that no extra attention is needed for the calculation. In general, the area integral appears to be superior to line integral in finite element applications due to its convenience in programming and better accuracy. As in VCEM, the area A can be chosen arbitrarily except for the requirement to be inside the body, which is also termed patch independence. It should be noted that D I M is equivalent to V C E M under a specific choice of integration domain and weighting function. To this end, the two contours Y~ and Fo are defined similar to those in the approximate finite element domain shown in Fig. 4(a), and a function for q is chosen to be unity on Fi, zero on Fo, and continuous between F~ and Fo. Although q is only C o continuous over the discretized domain, it is C ~' continuous inside each element and regarded as sufficiently smooth, With this special choice of q, only the elements between Fi and Fo make contribution to the domain integral in Eq. (7). It can be shown (e.g. [6]) that the domain integral in Eq. (71 is equivalent to the right-hand side of Eq. (5), which is the released potential energy. For the problems of curved cracks, a generalized formulation of D I M should be derived. In this case, we consider the shaded region A around the crack tip in Fig. 7 with Fi replaced by F, which shrinks onto the crack tip (this is not shown in the plot). The region A is bounded by a closed curve
J.H. Chang~Finite Elements in Analysis and Design 17 (1994) 323-338
329
t0
Fig. 6. The domain integral is performed in the region A enclosed by the contours Fi and Fo.
Fig. 7. Generalized DIM with a curved crack is performed on A and along C1 and C2.
C, where C = C 1 - F + C2 + Fo and m is the outward unit normal vector on C. Again, by introducing a weighting function q to Eq. (1), we have
J ~--
I/~?ui'~lqmjds fC [~V(~lj-- ffij~-~x1)
- f, +c2WmlqdS,
(8'
where q is a sufficiently smooth function in A, taken as unity on F (i.e. on O) and vanishing on Fo. Applying divergence theorem to the line integral on C in Eq. (8) and considering the state of equilibrium result in
J = -
~g(~lj-
i - x j xjjua-
Wmlqds.
(9)
I +C2
As observed, the additional term on the right-hand side of Eq. (9) - the line integral on C1 and C2, both terminated at the crack tip 0 - accounts for the energy release due to effect of the curved crack. With this extra term, the adjacent region surrounding the crack tip, where singular behavior is involved, is included in the formulation. However, the singular region can be avoided with appropriate choices of weighting function q, as described in the next section. In the finite element formulation, the area A is approximated by a finite element mesh Ah. Over each element, the displacement vectors u, and the weighting function q are interpolated with the local shape functions and their nodal point values. The energy release rate is then evaluated by Eq. (9) in the postprocessor of the finite element code once the nodal displacements under deformation are solved.
330
J.H. Chang/Finite Elements in Analysis and Design 17 (1994) 323 338
5. Numerical examples A n u m b e r of numerical investigations are carried out to test our formulations in this section. Three sets of problems are presented in the following subsections. The first considers the case of symmetric cracks emanating from a circular hole, and the next is a plane strain specimen with a central symmetric 'kinked' crack, followed by a problem of the same specimen with a curved crack. The calculation is generated by T E X P A C - N L finite element program. Results from both V C E M and the generalized D I M will be presented. Calculations under a number of patches, each covering different areas of the same finite element mesh, are performed in order to investigate the patch independence property. Note that the validity of V C E M calculations for problems with straight crack under both small and large deformation has been verified in the author's previous paper [6]. The nonlinear rubbery material behavior can be characterized in terms of the isotropic strain energy density function W as (10)
I/V= W ( I 1 , 1 2 , 1 3 ) ,
where I1, I2, 13 are the three principal invariants of the C a u c h y - G r e e n strain tensor V2. When the material is assumed to be incompressible, we then require 13 = 1,
(11)
i.e. square of the ratio of current to initial volume of the material remains unchanged under the deformation. In this case, the n u m b e r of independent variables reduces from three to two, i.e. 11 and 12. A particular isotropic incompressible model is used to simulate the behavior of a specific rubbery material here. This model is calibrated from a set of experimental data fit within a certain range of deformation [12]. The corresponding strain energy density is expressed as a function of the first two principal invariants, as suggested in [13] W = 390(I1 - 3) + 105(I2 - 3)
(unit: kPa).
{12)
In the following studies, quadratic finite elements are used to interpolate the displacement vector. No particular singular elements are used t h r o u g h o u t the calculation. 5.1. Crack emanating from a circular hole Problem 1 This test case considers a plane strain problem under uniform tensile loading a with a central circular hole, from which two straight cracks emanate symmetrically (Fig. 8). For isotropic linear elastic material, an approximate solution for this problem is expressed in [13] as
1 1 - v2
- 2
E
2 le. Sec
/'~leff"~
T ),
t131
where o- is the applied boundary tensile stress,/eff is the effective crack length and the secant is an approximate term to deal with the effect of finite width of the specimen. The value of/eff is given empirically in [14] as a function of the ratio of physical crack length I and diameter of the hole d.
J.H. Chang~Finite Elements in Analysis and Design 17 (1994) 323 338
331
A I I I I A
I 71.1 cm
d = 10.2 cm
i::~'+: d l
l = 1.61 cm
I
I
101.6 cm
I
I
Y
Y
V
Fig. 8. A plain strain specimen with symmetric cracks emanating from a circular hole.
Fig. 9. FE mesh of right lower quarter of the specimen in Fig. 8.
The factor ½ in Eq. (13) accounts for the definition of ERR as the surface free energy "per unit area", which is one-half in its value of the usual definition "unit length of extension in 2-D field". Because of symmetry, only the right lower quarter is considered in the finite element calculation, as shown in Fig. 9. Three patches are used to evaluate ERR. The exterior contours Fo'S of the patches are depicted in Fig. 10, while the associated interior contours Fi's are chosen one layer of elements inward (not shown here). The weighting function q of the generalized D I M in Eq. (9) is taken as zero on the outer contour of each patch, and unity on and within an interior contour inside the patch. With the choice of q, only the elements between the outer and interior contours, along with those containing the curved portion of crack surface, make contribution to the calculation. The singular area around the crack tip is thus not involved in this case. For this problem, it is interesting to note that singularity also happens around the edge of the hole. However, by choosing patches which are "large" enough, we observe that the area around the edge can be avoided in the calculation. We then have the following results given in Table 1. Although evaluation by Eq. (13) is only empirically approximated, it does appear to be very close to the results from our calculation. Also, both VCEM and D I M appear to be patch independent and yield almost the same values of ERR. This is due to very little contribution from the line integral on the right-hand side of Eq. (9) along portions of the traction-free circular hole. Note that the calculation is analytically patch independent. Nevertheless, small differences are observed from the above numerical results. This is because the equilibrium equations and the traction boundary conditions are satisfied only weakly in the finite element model.
332
J.H. Chang~Finite Elements in Analysis and Design 17 (1994) 323 338
(1)
(2)
/
Fig. 10. Three patches are used for calculation of ERR for both DIM and V C E M .
Table l Eq. (13) 1.07
VCEM
DIM Note: G (10 -2 Pam),/eft
=
Patch 1
Patch 2
Patch 3
1.16 1.16
1.14 1.14
1.17 1.16
3.5 cm, a = 1.38 M P a , E = 207 GPa,
v = 0.25.
5.2. Kinked crack
Here, we consider a plane strain specimen of width w in Fig. 11, with a symmetric central kinked crack. The problem is chosen to verify our DIM formulation in that a patch enclosing a straight portion along the crack surface can be appropriately defined. The corresponding ERR evaluated by VCEM with this particular patch is then used as a reference for comparison of subsequent calculations by the generalized DIM with different patches.
J,H. Chang~Finite Elements in Analysis and Design 17 (1994) 323-338
333
i
A I
J J
I
I I
ll.6cm
~ - -~
J,0.6 cm I
Fig. 11. A plain strain specimen with kinked crack.
Fig. 12. FE mesh of right half of the specimen in Fig. 11.
Table 2
VCEM DIM
Patch 1
Patch 2
Patch 3
Patch 4
99.9 99.9
99.4 99.4
96.4 99.2
97.1 99.1
Note: G (10 -2 Pam), a = 1.38 MPa, E = 207 GPa, v = 0.25.
Problem 2 In this test case, the specimen is considered of isotropic linear elastic material, under infinitesimal deformation. Only half of the piece is modeled with finite elements due to its symmetry, as shown in Fig. 12. The specimen is subjected to tensile loading o- uniformly along its top and bottom boundaries. Four patches (Fig. 13) are used to calculate the energy release rate. Each of the first two patches encloses a straight portion of the crack surface, with which validity of the calculation by VCEM is thus ensured. The numerical results are listed in Table 2. As expected, the results from VCEM and D I M match pretty well in the first two patches. The additional line integral on the right-hand side of Eq. (9) accounts for the discrepancy of the values between the two approaches in Patches 3 and 4. While D I M is proved to be patch independent, deviation of the results by VCEM does not appear to be substantial due to very small contribution from the line integral, in which the value of strain energy along the traction-free crack surface is rather small. Again, singularity at both the crack tip and the kinks is not explicitly involved in our calculation.
334
J.H. Chang~Finite Elements in Analysis and Design 17 (1994) 323 338
-q I
I
/
-
-
i (2)
(I)
/---
(3)
(4)
Fig. 13. Four patches are used for calculation of ERR for both DIM and VCEM.
A study on local mesh refinement around the crack tip area has been performed. The results show that our calculation is very insensitive to the crack tip models due to the choice of weighting function q. Problem 3 The same specimen as in the previous example, yet of rubbery material modeled in terms of the strain energy density function Eq.(12), is considered. Instead of the tensile loading, a uniform tensile stretch AL is applied along the top and bottom boundaries. In the calculation, non-conforming elements with reduced integration are used for the incompressible constraint in order to maintain
J.H. Chang / Finite Elements in Analysis and Design 17 (1994) 323-338
335
I 1
I
ll.6cm
I
~ - --'~
' V 10.6 cm I
~1__ _ 30.5cm
-- - - )
I
Fig. 14. The deformed mesh of the specimen in Fig. 1 l. The original configuration is depicted in dashed lines.
Fig. 15. A plain strain specimen with curved crack.
numerical stability. The deformed mesh, along with the original undeformed configuration, is shown in Fig. 14. Still, the four patches in Fig. 13 are used to evaluate the energy release rate, and the results are shown in Table 3. We have similar conclusions from the results. Patch independence holds for DIM well, while rather small deviation is observed for VCEM.
5.3. Curved crack Problem 4 In this case, we have a plane strain rubbery piece of width w (Fig. 15), with a symmetric central curved crack, subjected to uniform unilateral boundary stretch AL. A corresponding finite element mesh in Fig. 16 is defined to model half of the specimen due to its symmetry. The deformed mesh is shown in Fig. 17. As observed in this case, delimiting a patch containing a "straight portion" of the crack is possible only when a very complex finite element model is used around the crack tip area. It is therefore essential to investigate the results from the two approaches with solutions from the current finite element model. The same weighting function q is used here, with which those elements containing the crack tip are not involved in the calculations. Similarly, with the four patches in Fig. 18, we have the results, given in Table 4. As expected, DIM appears to be patch independent, and so does VCEM with only very small deviation. As a matter of fact, the results are very close to those in the previous problem, in which only the crack configuration is different.
336
J.H. Chang/Finite Elements in Analysis and Design 17 (1994) 323-338
I,, 111111
I I
//
IIII Illl
rlL
IIII Fig. 16. FE mesh of right half of the specimen in Fig. 15.
Fig. 17. The deformed FE mesh of the specimen in Fig. 15. The original configuration is depicted in dashed lines.
Table 3
VCEM DIM
Patch 1
Patch 2
Patch 3
Patch 4
43.5 43.5
43.8 43.8
42.7 43.3
42.7 43.0
Note: G (kPam), AL = 10.16 cm.
Table 4
VCEM DIM
Patch 1
Patch 2
Patch 3
Patch 4
44.0 43.8
44.0 43.8
43.3 43.5
43.5 43.8
Note: G (kPam), AL = 10.16 cm.
6. Conclusion T w o c o m m o n l y used a p p r o a c h e s for c a l c u l a t i o n of e n e r g y release rate for straight cracks, V C E M a n d D I M , are investigated for c o n s i d e r a t i o n of c u r v e d cracks u n d e r a specific a p p l i c a t i o n of 2-D r u b b e r y material problems. F o r straight cracks, these t w o m e t h o d s are equivalent with
J.H. Chang~Finite Elements in Analysis and Design 17 (1994) 323-338
337
f
J
(2)
(1)
f
(3)
(4)
Fig. 18. Four patches are used for calculation of ERR for both DIM and VCEM.
special choices of patches and weighting functions. For the case of curved cracks, VCEM does not appear to be feasible in its formulation with perturbation of two configurations, while an additional line integral along the crack surfaces should be included in derivation of the generalized DIM. Although VCEM gives good results in the presented examples due to small contribution from the line integral, the general DIM is still recommended, especially for those problems with more complex crack configurations and with crack surface tractions. It is demonstrated that our formulation, with appropriate choices of weighting function, is patch independent. This characteristic allows us to choose a patch that is remote from the crack tip and eliminate the need for complicated finite element modeling around it.
338
J.H. Chan#/Finite Elements in Analysis and Design 17 (1994) 323 338
Acknowledgments This work has been supported by National Science Council Grant No. NSC 82-0410-E-008-141 to National Central University.
References [1] J.D. Eshelby, "Calculation of energy release rate", in: Sih, Van Elst, Broek (ed.), Posp. Frac. Mech. Noordhoff, pp. 69 84, 1974. [2] R.S. Rivlin and A.G. Thomas, "'Rupture of rubber, I. Characteristic energy for tearing", J. Poly. Sci. X, pp. 291 318, 1953. [3] D.M. Parks, "A stiffness-derivative finite element technique for determination of elastic crack tip stress intensity factors", Int. J. Fraet. 10, pp. 487 502, 1974. [4] T.K. Helen, "On the method of virtual crack extensions", Int. J. Numer. Methods En{t. 9, pp. 187 207, 1975. [5] D.M. Parks, "The virtual crack extension method for nonlinear material behavior", Comput. Methods Appl. Mech. Eng. 12, pp. 353 364, 1977. [6] J.H. Chang and E.B. Becker, "'Finite element calculations of energy release rate of 2-D rubbery material problems with nonconservative crack surface tractions", Int. J. Numer. Methods Enq. 33, pp. 907 927, 1992. [7] K.M. Liechti, E.B. Becker, C. Lin and T.H. Miller, "A fracture analysis of cathodic delamination in rubber to metal bonds", Int. J. Fract. 39, pp. 217, 1989. [8] R.M.V. Pidaparti, T.Y. Yang and W. Soedel, "A plane stress finite element method for the prediction of rubber fracture", Int. J. Fract. 39, pp. 255 268, 1989. [9] F.Z. Li, C.F. Shih and A. Needleman, "A comparison of methods for calculating energy release rate", En~t. Fract. Mech. 21, pp. 405 421, 1985. [10] C.F. Shih, B., Moran and T. Nakamura, "'Energy release rate along a three-dimensional crack front in a thermally stressed body", Int. J. Fract. 30, pp. 79 102, 1986. [1 l] J.H. Chang, "Calculation of energy release rate for rubbery material problems with FE solutions", Chinese J. Mech., 1994, to appear. [12] E. Jankovich, G. Evrard, F. Leblanc and J.P. Nottin, Mixed finite element for the calculation of rubber articles on the basis of nonlinear elasticity laws, Final Study Report of a Research Project Financed by the DGRST, Tire Research Office and Laboratory, Colombes and Bezons, Information Research and Processing, Scientific Group, Paris, 1978. [13] R.W. Ogden, Non-linear Elastic" Dc#~)rmations, Wiley, New York, 1984. [14] D. Broek, Elementary En~ineeriny Fracture Mechanics, Martinus Nijhoff, 4th edn., 1986.
ELSEVIER
Finite Elements in Analysis and Design 17 (1994) 323-338
On the calculation of energy release rate for curved cracks of rubbery materials J.H. C h a n g Department of Civil Engineering, National Central University, Chungli, Taiwan, ROC Received January 1994; revised May 1994
Abstract
The energy release rate G for 2-D rubbery material problems with curved cracks is calculated with finite element solutions in this paper. Two approaches, a generalized domain integral method and the virtual crack extension method, are investigated. The generalized domain integral method is demonstrated to be "patch independent" and therefore a complicated finite element mesh adjacent to the crack tip area is not required.
I. Introduction
The overall strength of engineering structures relies on integrity of their components, which is essentially limited by the growth of cracks. For those materials with mechanical behavior remaining (nonlinearly) elastic at large strains, we often call them rubbery. Rubbery materials find applications in many structure components. Thus the prediction of crack initiation and growth of rubbery materials is of practical importance. Typically, two criteria for crack growth are used in fracture mechanics. One is critical stress criterion, the other is critical energy criterion. For problems involving nonlinear responses, the stress state around the crack tip area becomes rather complicated such that the use of stress criterion appears to be impractical. Therefore, the Energy Release Rate (ERR), usually denoted G, is often employed as a critical energy criterion for predicting the initiation and growth of cracks in nonlinear problems. The parameter, ERR, is defined by the concept of Griffith energy criterion for fracture states. For a mechanical system with a crack of length l, crack growth occurs only if the energy delivered by the system reaches the value required to form an additional crack of length dl, i.e. when the corresponding energy is up to G dl. The critical value of G is dependent upon a number of factors e.g. loading state, material property, and especially, crack geometry. Investigations and modifications on computation methods of ERR calculation for problems with cracks of general curved shape is hence necessary. 0168-874X/94/$07.00 © 1994 Elsevier Science B.V. All rights reserved SSDI 0 1 6 8 - 8 7 4 X ( 9 4 ) 0 0 0 5 0 - P
J.H. Chang/Finite Elements in Ana(vsis and Design 17 (1994) 323 338
324
Analytically, calculation of G for 2-D problems can be performed by a line integral, denoted J, which was first introduced by Eshelby [-1] for the problem of a cracked elastic body. By definition, J is an integral evaluated along a counterclockwise contour enclosing and shrinking onto the tip of a crack. For elastic media with traction-free crack surfaces, J is path independent, i.e. the value of integration remains unchanged along any shape of contour encircling the tip as long as the crack is straight inside the contour. This result is valid for both problems with small and large deformations. In the context of large deformation, path independence is valid when the crack remains straight inside the contour in its reference configuration. From the experimental results (e.g. [2] ), rubbery materials at large deformation state show substantially reversible elasticity except in the immediate neighborhood of the crack tip. The validity of the J-integral for these materials is thus presumed. However, for the case of curved cracks, the property of path independence no longer holds. Modifications on the formulation of J-integral need to be applied. In finite element calculations, two other approaches the virtual crack extension method (VCEM) and the domain integral method (DIM) - have been proposed as alternatives to direct application of J-integral. The first technique, VCEM, was introduced for problems of linear elasticity [3, 4] and deformation theory of plasticity [5], and later used for rubbery material problems [6 8]. The other method, DIM, performs a domain integral to determine ERR for deformation theory of plastic [9], thermoelastic [10], and rubbery material problems [11]. It should be noted that both methods were developed and considered for problems containing straight crack advancing along its original direction. Under this situation, they are equivalent being alternative interpretations of J-integral, as described in [-10]. One common advantage of both methods is patch independence, a counterpart of path independence of J. This characteristic eliminates the need for a complicated finite element modeling around the crack tip. The object of this paper is to develop a generalized formulation, based on the domain integral method, along with finite element solutions, for calculation of energy release rate of rubbery materials with curved cracks. The applicability of VCEM on curved cracks will be investigated, both theoretically and numerically. 2-D problems, including plane stress, and plane strain cases, are considered.
2. Curved cracks
Consider an elastic body (linear or nonlinear) in a 2-D field, containing a curved crack in its reference configuration (Fig. 1). We introduce a local coordinate system originating at the crack tip O, with the crack lying tangentially along the negative xl-direction. The body is subjected to a system of loads (not shown in the figure) and the crack will spread by an infinitesimal amount dl only if the corresponding value of G dl exceeds the critical surface free energy. The energy release rate can be evaluated by J-integral defined as G = J = lira
F+O
Wnl
--
aijn
j
ds,
(1)
where W is the strain energy density of the material, 6ij and u~ are the Cartesian components of the stress tensor and displacement vector, nj are the Cartesian components of outward unit vector
J.H. Chang~FiniteElementsin Analysisand Design17 (1994) 323 338
325
n
×2 :ii
x I
Fig. 1. A curved crack in its reference configuration with a local coordinate system.
Fig. 2. Path-independence of J for an elastic body with a straight crack.
normal to F, nl is the xa-component of nj, and s is the arc length along the contour. By definition, F is a counterclockwise contour encircling and shrinking onto the tip of the crack (this limiting case is not shown in Fig. 1). For the problems of large deformation, all the state variables are interpreted with respect to this specific reference configuration. When the crack surfaces are traction free and the material is h o m o g e n e o u s in the xa-direction, it can be proved (e.g. [6]) that the line integral can be executed along any other counterclockwise contour, Fi, which encloses a straight portion of the crack (Fig. 2), with the result remaining unchanged. This is termed "path independence". Note that the complicated stress behavior around the crack tip can be ignored in evaluating J-integral due to path independence. Modifications of J should apply in the case of curved cracks, when the straight portion can hardly be defined. If we choose an arbitrary contour Fi, as shown in Fig. 3 and from Eq. (1), J then becomes c3u~
J = fr, [ Wna - a~jnj(-~--x~x~)ldS + fc, +c2Wna ds.
(2)
In Eq. (2), two additional line integrals along Ca and C2, are included. A schematic diagram of the geometry is shown in Fig. 3. Although the path Fi can still be chosen arbitrarily, the extra curves C~ and C:, should both be terminated at the crack tip O. With this portion of line integrals, the singular behavior in the region near the crack tip is thus involved in the calculation.
3. Virtual crack extension method
V C E M was developed to calculate the energy difference between two configurations of finite element mesh with slightly altered crack length in their straight portions, as shown in Fig. 4.
J.H. Chang~Finite Elements in Analysis and Design 17 (1994) 323 338
326
//t--..\
!
Fig. 3. An arbitrary contour F~, along with two additional curves C1 and C2 along the crack, is used for the calculation.
(a)
(b)
Fi
I
---~
Fig. 4. (a) Finite element mesh before the crack advances. A band of elements are delimited for VCEM calculation. (b) Finite element mesh with a band of deformed elements.
Consider the mesh in Fig. 4(a) under a system of mechanical loads, with its discretized potential energy denoted by rc. The energy release rate is the decrease of x due to a unit advance of crack, i.e. the derivative of rc with respect to the crack length l, d~
~c~u
dl -
Ou ~I
+
c~x~x Ox
8l'
(3)
J.H. Chang~Finite Elements in Analysis and Design 17 (1994) 323-338
327
where the position and displacement vectors at the nodal points, x and u, are functions of I. Note that x is defined in the undeformed reference configuration. Under the state of equilibrium -
-
c~u
= 0,
(4)
Eq. (3) reduces to 0
-
drc dl
•rc t3x #x 81"
(5)
To compute the derivative of potential energy with respect to l, an imaginary extension of crack length is accommodated by rigidly translating all nodes on and within a contour Fi around crack tip O by an amount of Al along direction of the crack, with all other nodes remaining in their original positions (Fig. 4(b)). Thus the elements outside Fo retain both their initial positions and shapes, while those inside/'i remain unchanged in shape with a translation of Al. Due to this perturbation of mesh, Eq. (5) can then be evaluated. As verified in I-3] that the contour £o, which delimits the translating patch, can be chosen arbitrarily except for the requirements to be inside the body, enclose the crack tip and the straight portion of the crack. We then identify patch independence in this method as a counterpart of path independence in J-integral. This implies that results can be achieved by adopting a contour remote from the crack tip which will not depend directly on the finite element solutions in the immediate vicinity of the singularity. Checking feasibility of VCEM for curved cracks is necessary. For clarity, here we consider a special case in which a kinked crack is involved. For the configurations of finite element mesh in Figs. 5(a) and (b), all the nodes on and within Fi subject to a rigid translation by an amount Al along
(a)
(b)
ro
Fig. 5. (a) Original configuration of the finite element mesh for the virtual crack extension method. (b) New configuration after the curved crack extends along the xl direction.
J.H. Chanq/Finite Elements in Analysis and Design 17 (1994) 323 338
328
the x~ direction. Note that the altered mesh in Fig. 5(b) does not represent the proper geometry of the problem after crack extension along the x~-direction. It appears that further investigations on proper methods of ERR calculations for curved cracks is required.
4. Generalized domain integral method The domain integral m e t h o d was introduced as an alternative to calculate energy release rate for problems of straight cracks. Basically, an area integral is performed for 2-D problems, while a volume integral is calculated for 3-D problems. Again, we consider an elastic body in a twodimensional field containing a straight crack, with homogeneous material in the xl-direction, and no crack surface tractions (Fig. 2). The area integral, which evaluates the energy released from the mechanical system as the straight crack advances a unit length along the x 1-direction, is derived from the following path independent J-integral
( 0u; '~]ds.
(6)
The contour F; is usually chosen to be far from the crack tip so that a more accurate result can be ensured. In this integral method, we select an arbitrary exterior contour Fo and introduce a smooth weight function, q, which is unity on F i and vanishes on F o. By applying divergence theorem, Eq. (6) is then transformed into the following area integral: J =
-
W~)lj-
Gij~:XI~--
~ x i ] (](A '
(7)
where da is the infinitesimal integration area, 6aj is Kronecker delta with its first c o m p o n e n t as 1. As shown in Fig. 6, the integral is carried out in the shaded region A. This region A is often chosen to coincide with the original finite element mesh so that no extra attention is needed for the calculation. In general, the area integral appears to be superior to line integral in finite element applications due to its convenience in programming and better accuracy. As in VCEM, the area A can be chosen arbitrarily except for the requirement to be inside the body, which is also termed patch independence. It should be noted that D I M is equivalent to V C E M under a specific choice of integration domain and weighting function. To this end, the two contours Y~ and Fo are defined similar to those in the approximate finite element domain shown in Fig. 4(a), and a function for q is chosen to be unity on Fi, zero on Fo, and continuous between F~ and Fo. Although q is only C o continuous over the discretized domain, it is C ~' continuous inside each element and regarded as sufficiently smooth, With this special choice of q, only the elements between Fi and Fo make contribution to the domain integral in Eq. (7). It can be shown (e.g. [6]) that the domain integral in Eq. (71 is equivalent to the right-hand side of Eq. (5), which is the released potential energy. For the problems of curved cracks, a generalized formulation of D I M should be derived. In this case, we consider the shaded region A around the crack tip in Fig. 7 with Fi replaced by F, which shrinks onto the crack tip (this is not shown in the plot). The region A is bounded by a closed curve
J.H. Chang~Finite Elements in Analysis and Design 17 (1994) 323-338
329
t0
Fig. 6. The domain integral is performed in the region A enclosed by the contours Fi and Fo.
Fig. 7. Generalized DIM with a curved crack is performed on A and along C1 and C2.
C, where C = C 1 - F + C2 + Fo and m is the outward unit normal vector on C. Again, by introducing a weighting function q to Eq. (1), we have
J ~--
I/~?ui'~lqmjds fC [~V(~lj-- ffij~-~x1)
- f, +c2WmlqdS,
(8'
where q is a sufficiently smooth function in A, taken as unity on F (i.e. on O) and vanishing on Fo. Applying divergence theorem to the line integral on C in Eq. (8) and considering the state of equilibrium result in
J = -
~g(~lj-
i - x j xjjua-
Wmlqds.
(9)
I +C2
As observed, the additional term on the right-hand side of Eq. (9) - the line integral on C1 and C2, both terminated at the crack tip 0 - accounts for the energy release due to effect of the curved crack. With this extra term, the adjacent region surrounding the crack tip, where singular behavior is involved, is included in the formulation. However, the singular region can be avoided with appropriate choices of weighting function q, as described in the next section. In the finite element formulation, the area A is approximated by a finite element mesh Ah. Over each element, the displacement vectors u, and the weighting function q are interpolated with the local shape functions and their nodal point values. The energy release rate is then evaluated by Eq. (9) in the postprocessor of the finite element code once the nodal displacements under deformation are solved.
330
J.H. Chang/Finite Elements in Analysis and Design 17 (1994) 323 338
5. Numerical examples A n u m b e r of numerical investigations are carried out to test our formulations in this section. Three sets of problems are presented in the following subsections. The first considers the case of symmetric cracks emanating from a circular hole, and the next is a plane strain specimen with a central symmetric 'kinked' crack, followed by a problem of the same specimen with a curved crack. The calculation is generated by T E X P A C - N L finite element program. Results from both V C E M and the generalized D I M will be presented. Calculations under a number of patches, each covering different areas of the same finite element mesh, are performed in order to investigate the patch independence property. Note that the validity of V C E M calculations for problems with straight crack under both small and large deformation has been verified in the author's previous paper [6]. The nonlinear rubbery material behavior can be characterized in terms of the isotropic strain energy density function W as (10)
I/V= W ( I 1 , 1 2 , 1 3 ) ,
where I1, I2, 13 are the three principal invariants of the C a u c h y - G r e e n strain tensor V2. When the material is assumed to be incompressible, we then require 13 = 1,
(11)
i.e. square of the ratio of current to initial volume of the material remains unchanged under the deformation. In this case, the n u m b e r of independent variables reduces from three to two, i.e. 11 and 12. A particular isotropic incompressible model is used to simulate the behavior of a specific rubbery material here. This model is calibrated from a set of experimental data fit within a certain range of deformation [12]. The corresponding strain energy density is expressed as a function of the first two principal invariants, as suggested in [13] W = 390(I1 - 3) + 105(I2 - 3)
(unit: kPa).
{12)
In the following studies, quadratic finite elements are used to interpolate the displacement vector. No particular singular elements are used t h r o u g h o u t the calculation. 5.1. Crack emanating from a circular hole Problem 1 This test case considers a plane strain problem under uniform tensile loading a with a central circular hole, from which two straight cracks emanate symmetrically (Fig. 8). For isotropic linear elastic material, an approximate solution for this problem is expressed in [13] as
1 1 - v2
- 2
E
2 le. Sec
/'~leff"~
T ),
t131
where o- is the applied boundary tensile stress,/eff is the effective crack length and the secant is an approximate term to deal with the effect of finite width of the specimen. The value of/eff is given empirically in [14] as a function of the ratio of physical crack length I and diameter of the hole d.
J.H. Chang~Finite Elements in Analysis and Design 17 (1994) 323 338
331
A I I I I A
I 71.1 cm
d = 10.2 cm
i::~'+: d l
l = 1.61 cm
I
I
101.6 cm
I
I
Y
Y
V
Fig. 8. A plain strain specimen with symmetric cracks emanating from a circular hole.
Fig. 9. FE mesh of right lower quarter of the specimen in Fig. 8.
The factor ½ in Eq. (13) accounts for the definition of ERR as the surface free energy "per unit area", which is one-half in its value of the usual definition "unit length of extension in 2-D field". Because of symmetry, only the right lower quarter is considered in the finite element calculation, as shown in Fig. 9. Three patches are used to evaluate ERR. The exterior contours Fo'S of the patches are depicted in Fig. 10, while the associated interior contours Fi's are chosen one layer of elements inward (not shown here). The weighting function q of the generalized D I M in Eq. (9) is taken as zero on the outer contour of each patch, and unity on and within an interior contour inside the patch. With the choice of q, only the elements between the outer and interior contours, along with those containing the curved portion of crack surface, make contribution to the calculation. The singular area around the crack tip is thus not involved in this case. For this problem, it is interesting to note that singularity also happens around the edge of the hole. However, by choosing patches which are "large" enough, we observe that the area around the edge can be avoided in the calculation. We then have the following results given in Table 1. Although evaluation by Eq. (13) is only empirically approximated, it does appear to be very close to the results from our calculation. Also, both VCEM and D I M appear to be patch independent and yield almost the same values of ERR. This is due to very little contribution from the line integral on the right-hand side of Eq. (9) along portions of the traction-free circular hole. Note that the calculation is analytically patch independent. Nevertheless, small differences are observed from the above numerical results. This is because the equilibrium equations and the traction boundary conditions are satisfied only weakly in the finite element model.
332
J.H. Chang~Finite Elements in Analysis and Design 17 (1994) 323 338
(1)
(2)
/
Fig. 10. Three patches are used for calculation of ERR for both DIM and V C E M .
Table l Eq. (13) 1.07
VCEM
DIM Note: G (10 -2 Pam),/eft
=
Patch 1
Patch 2
Patch 3
1.16 1.16
1.14 1.14
1.17 1.16
3.5 cm, a = 1.38 M P a , E = 207 GPa,
v = 0.25.
5.2. Kinked crack
Here, we consider a plane strain specimen of width w in Fig. 11, with a symmetric central kinked crack. The problem is chosen to verify our DIM formulation in that a patch enclosing a straight portion along the crack surface can be appropriately defined. The corresponding ERR evaluated by VCEM with this particular patch is then used as a reference for comparison of subsequent calculations by the generalized DIM with different patches.
J,H. Chang~Finite Elements in Analysis and Design 17 (1994) 323-338
333
i
A I
J J
I
I I
ll.6cm
~ - -~
J,0.6 cm I
Fig. 11. A plain strain specimen with kinked crack.
Fig. 12. FE mesh of right half of the specimen in Fig. 11.
Table 2
VCEM DIM
Patch 1
Patch 2
Patch 3
Patch 4
99.9 99.9
99.4 99.4
96.4 99.2
97.1 99.1
Note: G (10 -2 Pam), a = 1.38 MPa, E = 207 GPa, v = 0.25.
Problem 2 In this test case, the specimen is considered of isotropic linear elastic material, under infinitesimal deformation. Only half of the piece is modeled with finite elements due to its symmetry, as shown in Fig. 12. The specimen is subjected to tensile loading o- uniformly along its top and bottom boundaries. Four patches (Fig. 13) are used to calculate the energy release rate. Each of the first two patches encloses a straight portion of the crack surface, with which validity of the calculation by VCEM is thus ensured. The numerical results are listed in Table 2. As expected, the results from VCEM and D I M match pretty well in the first two patches. The additional line integral on the right-hand side of Eq. (9) accounts for the discrepancy of the values between the two approaches in Patches 3 and 4. While D I M is proved to be patch independent, deviation of the results by VCEM does not appear to be substantial due to very small contribution from the line integral, in which the value of strain energy along the traction-free crack surface is rather small. Again, singularity at both the crack tip and the kinks is not explicitly involved in our calculation.
334
J.H. Chang~Finite Elements in Analysis and Design 17 (1994) 323 338
-q I
I
/
-
-
i (2)
(I)
/---
(3)
(4)
Fig. 13. Four patches are used for calculation of ERR for both DIM and VCEM.
A study on local mesh refinement around the crack tip area has been performed. The results show that our calculation is very insensitive to the crack tip models due to the choice of weighting function q. Problem 3 The same specimen as in the previous example, yet of rubbery material modeled in terms of the strain energy density function Eq.(12), is considered. Instead of the tensile loading, a uniform tensile stretch AL is applied along the top and bottom boundaries. In the calculation, non-conforming elements with reduced integration are used for the incompressible constraint in order to maintain
J.H. Chang / Finite Elements in Analysis and Design 17 (1994) 323-338
335
I 1
I
ll.6cm
I
~ - --'~
' V 10.6 cm I
~1__ _ 30.5cm
-- - - )
I
Fig. 14. The deformed mesh of the specimen in Fig. 1 l. The original configuration is depicted in dashed lines.
Fig. 15. A plain strain specimen with curved crack.
numerical stability. The deformed mesh, along with the original undeformed configuration, is shown in Fig. 14. Still, the four patches in Fig. 13 are used to evaluate the energy release rate, and the results are shown in Table 3. We have similar conclusions from the results. Patch independence holds for DIM well, while rather small deviation is observed for VCEM.
5.3. Curved crack Problem 4 In this case, we have a plane strain rubbery piece of width w (Fig. 15), with a symmetric central curved crack, subjected to uniform unilateral boundary stretch AL. A corresponding finite element mesh in Fig. 16 is defined to model half of the specimen due to its symmetry. The deformed mesh is shown in Fig. 17. As observed in this case, delimiting a patch containing a "straight portion" of the crack is possible only when a very complex finite element model is used around the crack tip area. It is therefore essential to investigate the results from the two approaches with solutions from the current finite element model. The same weighting function q is used here, with which those elements containing the crack tip are not involved in the calculations. Similarly, with the four patches in Fig. 18, we have the results, given in Table 4. As expected, DIM appears to be patch independent, and so does VCEM with only very small deviation. As a matter of fact, the results are very close to those in the previous problem, in which only the crack configuration is different.
336
J.H. Chang/Finite Elements in Analysis and Design 17 (1994) 323-338
I,, 111111
I I
//
IIII Illl
rlL
IIII Fig. 16. FE mesh of right half of the specimen in Fig. 15.
Fig. 17. The deformed FE mesh of the specimen in Fig. 15. The original configuration is depicted in dashed lines.
Table 3
VCEM DIM
Patch 1
Patch 2
Patch 3
Patch 4
43.5 43.5
43.8 43.8
42.7 43.3
42.7 43.0
Note: G (kPam), AL = 10.16 cm.
Table 4
VCEM DIM
Patch 1
Patch 2
Patch 3
Patch 4
44.0 43.8
44.0 43.8
43.3 43.5
43.5 43.8
Note: G (kPam), AL = 10.16 cm.
6. Conclusion T w o c o m m o n l y used a p p r o a c h e s for c a l c u l a t i o n of e n e r g y release rate for straight cracks, V C E M a n d D I M , are investigated for c o n s i d e r a t i o n of c u r v e d cracks u n d e r a specific a p p l i c a t i o n of 2-D r u b b e r y material problems. F o r straight cracks, these t w o m e t h o d s are equivalent with
J.H. Chang~Finite Elements in Analysis and Design 17 (1994) 323-338
337
f
J
(2)
(1)
f
(3)
(4)
Fig. 18. Four patches are used for calculation of ERR for both DIM and VCEM.
special choices of patches and weighting functions. For the case of curved cracks, VCEM does not appear to be feasible in its formulation with perturbation of two configurations, while an additional line integral along the crack surfaces should be included in derivation of the generalized DIM. Although VCEM gives good results in the presented examples due to small contribution from the line integral, the general DIM is still recommended, especially for those problems with more complex crack configurations and with crack surface tractions. It is demonstrated that our formulation, with appropriate choices of weighting function, is patch independent. This characteristic allows us to choose a patch that is remote from the crack tip and eliminate the need for complicated finite element modeling around it.
338
J.H. Chan#/Finite Elements in Analysis and Design 17 (1994) 323 338
Acknowledgments This work has been supported by National Science Council Grant No. NSC 82-0410-E-008-141 to National Central University.
References [1] J.D. Eshelby, "Calculation of energy release rate", in: Sih, Van Elst, Broek (ed.), Posp. Frac. Mech. Noordhoff, pp. 69 84, 1974. [2] R.S. Rivlin and A.G. Thomas, "'Rupture of rubber, I. Characteristic energy for tearing", J. Poly. Sci. X, pp. 291 318, 1953. [3] D.M. Parks, "A stiffness-derivative finite element technique for determination of elastic crack tip stress intensity factors", Int. J. Fraet. 10, pp. 487 502, 1974. [4] T.K. Helen, "On the method of virtual crack extensions", Int. J. Numer. Methods En{t. 9, pp. 187 207, 1975. [5] D.M. Parks, "The virtual crack extension method for nonlinear material behavior", Comput. Methods Appl. Mech. Eng. 12, pp. 353 364, 1977. [6] J.H. Chang and E.B. Becker, "'Finite element calculations of energy release rate of 2-D rubbery material problems with nonconservative crack surface tractions", Int. J. Numer. Methods Enq. 33, pp. 907 927, 1992. [7] K.M. Liechti, E.B. Becker, C. Lin and T.H. Miller, "A fracture analysis of cathodic delamination in rubber to metal bonds", Int. J. Fract. 39, pp. 217, 1989. [8] R.M.V. Pidaparti, T.Y. Yang and W. Soedel, "A plane stress finite element method for the prediction of rubber fracture", Int. J. Fract. 39, pp. 255 268, 1989. [9] F.Z. Li, C.F. Shih and A. Needleman, "A comparison of methods for calculating energy release rate", En~t. Fract. Mech. 21, pp. 405 421, 1985. [10] C.F. Shih, B., Moran and T. Nakamura, "'Energy release rate along a three-dimensional crack front in a thermally stressed body", Int. J. Fract. 30, pp. 79 102, 1986. [1 l] J.H. Chang, "Calculation of energy release rate for rubbery material problems with FE solutions", Chinese J. Mech., 1994, to appear. [12] E. Jankovich, G. Evrard, F. Leblanc and J.P. Nottin, Mixed finite element for the calculation of rubber articles on the basis of nonlinear elasticity laws, Final Study Report of a Research Project Financed by the DGRST, Tire Research Office and Laboratory, Colombes and Bezons, Information Research and Processing, Scientific Group, Paris, 1978. [13] R.W. Ogden, Non-linear Elastic" Dc#~)rmations, Wiley, New York, 1984. [14] D. Broek, Elementary En~ineeriny Fracture Mechanics, Martinus Nijhoff, 4th edn., 1986.