Engineering Fracture Mechanics Vol. 51, No. 6, pp. 7 1S-720, 1997 Q 1997 Elsevier Science Ltd. All rights reserved
Pergamon
PII: soo13-7944(97)ooo6o-x
Printed in Great Britain 0013-7944/97 $17.00 + 0.00
TECHNICAL NOTE FIRST AND SECOND ORDER STRESS EXTRAPOLATION FOR MODE I AND MODE II EDGE CRACKS BYEONG S. KIM and ALAN W. EBERHARDT Department of Materials and Mechanical Engineering, University of Alabama at Birmingham, Birmingham, Alabama 35294, U.S.A. Abstract-The use of linear and second order stress extrapolation to obtain Kt and &t in two-dimensional finite element models of a thick plate containing an edge crack was examined. Three loading cases were studied, including classical Mode I and Mode II problems and a problem of tribological contact. Linear extrapolation was observed to yield accurate predictions of 4 in cases of dominant Mode I loading. In Mode II situations, notably where the crack faces experienced compressive normal stresses, second order extrapolation was observed to improve estimates of 4,. 0 1997 Elsevier Science Ltd Keywords-finite
elements, stress extrapolation, Mode II, tribology.
1. INTRODUCTION MANY DIFFERENTmethods may be used for calculating stress intensity factors (SIFs) at crack tips using finite elements, including extrapolation, nodal force, J-integral, Griffith energy and stiffness derivative techniques (see ref. [l] for a complete review). Each of these methods has its benefits and shortcomings, and to date, no consensus exists on the “best” method. Although it is considered less accurate than other methods, extrapolation is still commonly used (e.g. ref. [2]), recommended for its simplicity and ability to yield “reasonable results”[3]. In extrapolation techniques, finite elements are used to estimate stresses or displacements at nodes in the vicinity of the crack tip. The “apparent SIFs” are calculated and plotted as a function of node location. The crack tip SIFs are traditionally determined using linear extrapolation, either by collocation (fit through straight line sections of the curves)[4], or by linear regression including nonlinear regions of the curves [3,5]. Studies have reported accurate KI predictions for problems of symmetric, Mode I loading, however, the method has been less successful in predicting KIr in Mode II or mixed mode loading[3,5]. Surface cracking and subsequent pitting is one source of tribological failure in gears, cams and other contacting systems. The mechanisms by which such surface cracks propagate are not completely understood, as contact pressures near a surface crack cause complex states of tensile, compressive and shear stresses at the crack tip [6]. In a recent investigation using stress extrapolation to study surface cracks in tribological coatings the authors observed that, in cases where Mode II dominates, the apparent K,t curves may be strongly nonlinear [7]. In such cases, traditional linear regression yielded erroneous and ambiguous results, depending on the number of nodes used in the regression. The purpose of this technical note is to present findings which demonstrate that second order stress extrapolation may improve 4, predictions.
2. METHODS The present methodology includes the finite element model construction, calculation of apparent stress intensities at nodes ahead of the crack tip from nodal stresses, linear and second order extrapolation of stress intensities and comparison with closed-form analytical results. Three models were constructed using COSMOS/M finite element software (Structural Research and Analysis Corporation, Los Angeles, CA) as a plate in plane strain containing a vertical edge crack on the upper surface and subject to tensile, shear or contact loading [Fig. l(a)-(c)]. The majority of elements used were two-dimensional QM6 elements, with triangular elements used to refine the mesh near the crack, where convergence was confirmed. The QM6 element is a nonconforming element which introduces nodeless variables which generate additional degrees of freedom, providing improved accuracy over traditional four-node quadrilateral elements [8]. Material properties of steel (E = 200 GPa and v = 0.3) were prescribed to all models. In Model 1, a uniform tensile (Mode I) pressure, u = 20 psi (138 kPa), was applied to the right model edge in the xdirection, while elements along the left edge were fully constrained. The uniform mesh, containing 4860 elements, is shown in Fig. 2(a). The stress intensity factor for this Mode I loading is given by the well-known formula[9]: Kl = o(m~)“~F,(a/b),
(1)
where a is the crack length, b is the plate width and F&/b) is obtained as FI(n/b) = [(2b/rru)tan(nu/2b)]‘/2{0.752 + 2.02(u/b) + 0.37[1 - sin(no/2b)‘])/[cos(rru/2@].
(2)
This model was constructed with a = 16 in (40.6 cm) and b = 40 in (102 cm), therefore, u/b = 0.4 and F,(u/b) was obtained as 2.108. Model 2, that of an edge cracked structure loaded in shear, is dominated by Mode II. The mesh was similar to Model 1 with reduced crack length. In this case, a line load, Q = 100 lb/in (175 N/cm), was applied as direc715
716
Technical Note
(4
b
Q
Q
b) L a v b
2a
z
Fig. 1. Schematic diagrams of geometries and loadings for (a) Model 1: a plate with a vertical edge crack loaded in tension; (b) Model 2: a plate with a vertical edge crack loaded in shear; and (c) elastic half-space with a vertical surface crack loaded by rolling and sliding Hertz contact. Note that in (a) and (b) the crack tip openings are exaggerated for illustration.
tionally opposite forces at the nodes on opposite faces at the crack opening. In this case[9] ~11 =
Q[2/(na)"*l41(alb).
Ftr(a/b) = [1.3 - 0.65(a/b) + 0.37(a/b)* + 0.28(a/b)3]/(1 - a/b)“*.
(3)
(4)
In Model 2, (I = 2 in (5.08 cm) and b = 10 in (25.4 cm), therefore, a/b = 0.2 and F&/b) was obtained as 1.325. A third model was developed of normal and sliding Hertzian pressure on a uniform elastic media containing a vertical surface crack. This problem was solved analytically by Keer ef ul.[6], who presented stress intensity values graphically for a range of loads, crack lengths and frictional loads. In Model 3, the coefficient of surface friction was f = 0.25. The ratio of crack length to contact radius, c/a = 0.5. Load locations, g/a = 1 and g/a = 5 were examined, where g is the distance from crack to load center, as shown in Fig. l(c). At g/a = 1, the crack experiences high shear stresses and compressive Mode 1 stresses, while at g/a = 5, the load is primarily tensile, due to the surface friction and the direction of sliding. Figure 2(b) illustrates the mesh used in Model 3.
(b)
Fig. 2. (a) Finite element mesh of Model 1, showing full constraint along the left edge and uniform xdirected stress applied at right. The crack is the bold line in the center of the mesh. (b) Finite element mesh of Model 3, showing the use of triangular elements to refine the mesh in the vicinity of the crack tip. The model is fully constrained along the bottom, and constrained horizontally at the left and right edges. Normal and horizontal surface tractions are illustrated and the crack is the bold line in the center of the mesh.
718
Technical Note
664x”2
77.38 + 0.96786
RA2 = 0.998
270
260 y = 287.02 -
5.0560X
R”2
2
3
= 0.970
250 0
1
4
5
6
Distancefrom crack tip [in] Fig. 3. First and second order extrapolations
of Ki for Model 1 using seven nodes ahead of the crack tip.
In Model 1, the crack faces were free of constraint. In Models 2 and 3, frictionless contact between elements on the opposing crack faces was obtained using node-to-line gap elements[8]. Crack mouth openings between 0.01-0.1 pm produced consistent and accurate results. For mode I crack opening, the well-known equation for the SIF is K] =
s/%&r),
(5)
where ON(r), are nodal stresses directly ahead of, and normal to, the crack tip. Apparent K,* values were calculated at each node approaching the crack tip and plotted versus node location. The crack tip K, was extrapolated at the intersec-
120 110
)’= 105.78 - 54.471x + 10.515x”2 R*2 = 1.000
100
60 50
0.0
0.5
1.0
1.5
Distance from the crack tip [in] Fig. 4. First and second order extrapolations of 41 for Model 2 using seven nodes ahead of the crack tip.
Technical Note
719
280
200 206.59x
120
3 + 679.09x - 792.37xA2
R”2 = 0.862
R”2 = 0.997 I
80 0.0
0.1
0.2
0.3
0.4
0.5
Distance from the crack tip [mm] Fig. 5. First and second order extrapolations of KI for Model 3 using seven nodes ahead of the crack tip.
tion of the vertical (r = 0) axis. K II* curves were calculated similarly for shear stresses, crs, where KII = &us(r).
(6)
Linear and second order extrapolations were performed for each of the models using the method of least squares in Cricket Graph (Cricket Software, Malvern, PA). In each case, K,* and K I,* values from the two nodes nearest the crack tip were not included in the regression, as recommended in ref.[4].
= - 526.96 - 3587.4x + 2831.9xA2
$ B
RA2 = 0.995
-1000 -1200
3 -1400
-1600 b
-1800 0.0
0.1
0.2
0.3
0.4
0.5
0.6
Distance from crack tip [mm] Fig. 6. First and second order extrapolations of 41 for Model 3 using eight nodes ahead of the crack tip.
720
Technical Note
3. RESULTS AND DISCUSSION For Model 1, eqs (1) and (2) yield the analytical solution, Ki = 299 psi& (329 kPa&i). Figure 3 illustrates the seven-ooint linear extrapolation of the FE stresses. for which Kr = 287 nsi& (316 kPa./m). This value deviates from the analytical solution by 3.3%. The second order extrapolation yields 7.4% error, where ki = 277 psi& (305 kPa,/Z). The linear extrapolation was improved as more nodes were included and the fit approached that of the coilocation technique[l]. The second order extrapolation did not improve with the inclusion of more nodes. For Model 2 the analytical solution is K 11 = 106 psi& (117 kPa,/iii). Linear extrapolation of nodal stresses ahead of the crack tip yielded Kn = 97 psi& (107 kPa,/iii=8% error), while the second order fit yielded Ku = 106 psi& (117 kPafi; Fig. 4). Inclusion of more data points in the extrapolation did not improve the linear estimate. In this case, therefore, where Mode II dominates, the second order extrapolation provided more accurate and consistent predictions of Ku and the linear extrapolations underestimate the analytical solution. The theoretical solution for KI in Model 3, as visually extracted from graphs in ref. [6], for g/a = 5, lies between 168 MPa,/?ii (153 psi&) and 174 MPaJiir (158 psi&). The results of linear extrapolation (Fig. 5), using seven nodes from the present FE models, yielded the value of K, = 173 MPa,/iYi (158 psi& . This falls within the range extracted from ref. [6], while the second order prediction of Ki = 114 MPaJiii (103 psi > in) is well below the cited range. The value of Ku, as visually extracted from the graphs in ref. [6], for g/a = 1, lies between -515 MPaJiii (469 psi&) and -540 MPa,/iii (492 psi&). The results of the linear and second order extrapolations, using eight nodes from the present FE models, are shown in Fig. 6. The figure illustrates that the value of Ku = - 712 MPafi (-647 psi&) predicted by the linear extra olation falls well outside the cited range, while the second order prediction of Ku = - 527 MPafi (480 PSIY- m) is consistent with the results from ref. [6]. It was further observed that, in this case, at least seven nodes were necessary to achieve accurate regression, and that the corresponding correlation coefficient for the second order extrapolation approached r * = 1 as more nodes were included.
4. CONCLUSIONS The use of linear and second order stress extrapolation of K, and Kn for a few select edge crack problems was examined. For the limited cases studied, linear extrapolation best predicted K, for cases dominated by Mode I loading. In cases where Mode II dominated, second order extrapolation was observed to yield better estimates of Ku. It should be noted that singular elements were not used in any of the present models. Their use may affect the present trends and conclusions. Acknowledgements-The authors gratefully acknowledge support of the UAB Graduate Faculty Research Grants Programs.
Research Assistantship
and
REFERENCES 1. Anderson, T. L., Fracture Mechanics: Fundamental and Applications, 2nd edn. CRC Press, Boca Raton, 1995. 2. Menandro, F. C. N., Moyer, E. T. and Liebowitz, H., A near optimal crack tip mesh. Engineering Fracture Mechanics, 1995, !!0(5/6), 703-71
I.
3. Banks-Sills. L. and Sherman, D., Comparison of methods for calculating stress intensity factors with quarter point elements. International Journal of Fracture, 1986, 32, 127-140. 4. Chan, S. K., Tuba, I. S. and Wilson, W. K., On the finite element method in linear fracture mechanics. Engineering Fracture Mechanics, 1970, 2, I-17.
5. Fawkes, A. J., Owen, D. R. J. and Luxmoore, A. R., An assessment of crack tip singularity models for use with isoparametric elements. Engineering Fracture Mechanics, 1979, 13, 143-159. 6. Keer, L. M., Bryant, M. D. and Haritos, G. K., Subsurface and surface cracking due to Hertzian contact. Journal of Lubrication Te&nology, 1982, 104, 347-351. 7. Eberhardt. A. W. and Kim. B. S.. Crack face friction effects on Mode II stress intensities for a surface cracked coating in two-dimensional rolling contact, to appear in Trib. Trans., 1997. 8. COSMOS/M Finite Element Analysis System, Version 1.75, User Guide, Vol. I, 1995. 9. Tada, H., The Stress Analysis of Cracks Handbook. Del Research Corporation, 1973. (Received 15 October 1996, infinalform
30 April 1997, accepted 4 May 1997)