Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 45, 2, pp. 225-237, Warsaw 2007 INFLUENCE OF THE CRACK TIP MODEL ON RESULTS OF THE FINITE ELEMENTS METHOD Marcin Graba Jarosław Gałkiewicz Faculty of Mechatronics and Machine Design, Kielce University of Technology e-mail: mgraba@eden.tu.kielce.pl; jgalka@eden.tu.kielce.pl The paper contains recommendations of finite element models of the cracktipneighborhoodtoobtain results independentof thefinite element mesh. The recommendations are valid for elastic-plastic problems and finite strains. As an example, analysis of single edge notched specimens under bending is presented. Key words: fracture mechanics, finite deformation, stress distribution, J-integral, CTOD, finite elements method 1. Introduction The Finite ElementsMethod (FEM) requires propermodeling of a given pro- blem. The modeling concerns specimen geometry, description of a material, selection of the proper Finite Elements (FE) mesh. The Hutchinson, Rice and Rosengren (HRR) solution (Hutchinson, 1968; Rice and Rosengren, 1968) of stress, strain and displacement fields near the crack tip for non-linear (elastic-plastic) materials concerns small strains and consists of one singular term only. The amplitude of the singular stress field is the J-integral. The J-integral is path independent when the strain energy is a unique function of strains. However, the stress distribution for a plane strain model is in most cases different from theHRR solution (see Fig.1). This difference between theFEM stress distribution and the HRR solution was called by O’Dowd and Shih the ”Q-parameter” (O’Dowd and Shih, 1991, 1992). In fact, the Qσ0 term (where σ0 is yield stress), when added to the HRR singular term, replaces all neglected terms in the asymptotic expansions of the stress field in front of the crack. 226 M. Graba, J. Gałkiewicz Fig. 1. (a) The stress distribution near the crack tip. Curves obtained using FEM for small and finite strain and HRR formula; (b) the stress-strain curve of the material used in FEM analysis In a real structural element with a crack, stresses near the crack tip are finite. The stress infinity is a result of the assumption that the crack tip is perfectly sharp and it remains sharp during loading of the crack tip. When the assumption of small strains is relaxed, the crack tip blunts and the stres- ses in front of the crack become finite. The opening stress reaches maximum at the distance equal to r = 0.5J/σ0 to r = 2J/σ0, and its value depends on material properties, specimen geometry and external loading. This fe- ature was first noticed by Rice and Johnson (1970), McMeeking and Parks (1979). FEanalysis of the stress and strainfield in front of the crackwhen thefinite strains are used is not a trivial problem. The level of the stress maximumand its localization depends on FE mesh details when it is not properly selected. For the small strain option, such a problem is not observed (ODowd and Shih, 1992; Al-Ani and Hancock, 1991; O’Dowd, 1995; O’Dowd et al., 1995). 2. Numerical model Numerical results are presented for single edge notched bend speci- mens (SEN(B)). Dimensions of the specimens satisfy requirements of the ASTM E 1820-05 standard. Computations were performed for a plane strain using the finite strain option. The relative crack lengthwas a/W =0.5, where a is the crack length, and the width of specimens W was equal to 40mm. Influence of the crack tip model... 227 The computationswere performedusingADINASYSTEM8.3.Due to the symmetry, only a half of the specimen was modeled. The finite element mesh was filled with 9-node plane strain elements. The size of finite elements in the radial direction was decreasing towards the crack tip, while in the angular direction the size of each elementwas kept constant. It varied from ∆θ = π/13 to ∆θ = π/23 for various cases tested. The crack tip region was modeled using 6, 16, 35 and 75 semicircles. The crack tip wasmodeled as an arc whose radius varied from 0.0001mm to 0.01mm, (10−6(W − a) to 10−4(W − a), respectively). The Poisson ratio which was ν = 0.3, the yield stress σ0 = 315MPa, Young’s modulus E = 206000MPa and the work hardening exponent n = 5 definedmechanical properties of thematerial.The true stress-strain curveused in FEM analysis is presented in Fig.1b. In the model, the stress-strain curve was approximated by relation ε ε0 =      σ σ0 for σ ¬ σ0 α ( σ σ0 )n for σ > σ0 (2.1) where ε0 denotes the reference strain (ε0 = σ0/E) and α is a hardening constant, which was assumed 1 in this case. 3. Results of finite elements analysis The numerical analysis of the stress distribution in the domain next to the crack tip revealed that results dependon the details of FEmodeling, when the finite strain option is adopted. It is shown in Fig.2 andFig.3. Onemay notice that when the number of the FEs between the crack tip and the maximum location of theopening stress is not large enough for the results to converge toa single curve. It is recommendedtouseat least 20FEsbetween thecrack tipand the stress maximum location. Thus, the FE size in the radial direction should be smaller than 0.1δT , where δT is the Crack Tip Opening Displacement (CTOD). In contrast to the FE size in the radial direction, the size in the angular direction does not effect the stress distribution in front of the crack at θ =0◦. Several cases for the FE modeled as semicircular rings were tested (from 13 to 23 segments). The results are shown in Fig.4. In the literature, the crack tip is modeled in two different ways shown in Fig.5c and Fig.5d. 228 M. Graba, J. Gałkiewicz Fig. 2. Influence of the number of FE between the crack tip and the maximum location of the opening stress σ22 along the θ =0 ◦ direction on the stress distribution: (a) r ¬ 6J/σ0; (b) r ¬ 2J/σ0 Brocks et al. (2003) andBrocks andScheider (2003) suggests that the crack tip shouldbemodeled in theway shown inFig.5d.Computations confirmthat when this model is used, the radius of the crack tip can be smaller than for the model shown in Fig.5c to conduct the FE analysis for larger external loads. O’Dowd and Shih (1991, 1992), O’Dowd (1995), O’Dowd et al. (1995) suggests that the crack tip radius should be smaller than the half of the crack tip opening displacement δTc for the critical moment rw = 1 2 δTc = 1 2 dn Jc σ0 (3.1) where Jc is the critical value of the J-integral (for the material used in the FEM analysis: Jc = 40kN/m) and dn is a parameter introduced by Shih, which connects the J-integral, yield stress and crack tip opening displacement (in our case dn =0.297). The crack tip radius calculated from equation (3.1) is equal to rw =1.89 ·10 −5mand is greater than the crack tip radii tested in the FEM analysis which seemed to be advisable. However O’Dowd was not interested in the opening stress distribution at the distance r < J/σ0 which is of vital interest in our analysis. He used the small strain option during computation. The FEM analysis shows that when the crack tip radius ρ decreases, the level of the maximum of the opening stress in front of the crack increases Influence of the crack tip model... 229 Fig. 3. Influence of the number of FE between the crack tip and the maximum location of the opening stress σ22 along the θ =0 ◦ direction: (a) the level of the maximum opening stress as a function of the external loading (JIC =40000N/m for this material); (b) normalized location of the opening stress maximum as a function of the external load; (c) location of the opening stress maximum as a function of the external loading (P0 is the limit load) and appears closer to the crack tip. However, for a sufficiently small value of ρ, both the opening stress maximum and its location in front of the crack become independent of the crack tip radius. For increasing external load, the saturation of the ξ0 = ξ0(J) and ψ0 = ψ0(J) (where ξ0 = (σ22)max/σ0 and ψ0 =(r22)maxσ0/J) curves was observed (see Fig.6). The curves level off, and the level is independentof the crack tip radius.Convergence of theFEMresults is observed,when ρ is about 2.5·10−6m (δTc/15 for the criticalmomentwhen the J-integral value is equal to Jc =40kN/m, respectively). 230 M. Graba, J. Gałkiewicz Fig. 4. Influence of the FE size in the tangential direction on FE results: (a) stress distributions in front of the crack tip for different levels of external loads (results for the θ =0◦ direction); (b) normalized distance of the maximum opening stress to the crack tip; (c) distance of the maximum opening stress to the crack tip Thus the crack tip model shown in Fig.5d is recommended. 4. Estimation of the J-integral In this section, the influence of size of contour of integration on the J-integral is investigated. In ADINA SYSTEM8.3, the J-integral may be calculated by making use of two methods. The first method is based on the J-integral definition J = ∫ C ( w dx2− t ∂u ∂x1 ) ds (4.1) Influence of the crack tip model... 231 Fig. 5. (a) SEN(B) specimen; (b) FEMmesh for SEN(B) specimen; (c), (d) two alternative crack tip models 232 M. Graba, J. Gałkiewicz Fig. 6. Effect of size of the crack tip radius on: (a) level of the maximum opening stress; (b) normalized location of the maximum opening stress where w is the strain energy density, t is the stress vector acting on the contour C drawn around the crack tip, u denotes the displacement vector and ds is an infinitesimal segment of the contour C. The secondmethod, called the ”virtual shift method”, uses the concept of virtual crack growth to compute the virtual energy change. The J-integral is path independent for small strain formulation only. Ho- wever, onemayalso calculate the J-integral using the large strain formulation, but the contour of integration should be sufficiently distant from the crack tip and not too close to the specimen borders. Influence of the crack tip model... 233 InTable 1, the contours of integration used in this investigation are shown. Table 1.Definition of J-integral contours Number of Distance form contour crack tip [mm] 1 0.02 2 0.15 3 0.36 4 0.67 5 1.08 6 1.59 7 2.25 8 4.64 9 9.39 10 18.0 11 18.0 Fig. 7. Effect of size of the integral contour and the mesh size on the numerically found J-integral Figure 7 presents J-integral values obtained for different contours and different meshes. The smallest FE mesh (the first mesh in Table 2) used in the analysis was extremely dense. Subsequent meshes were gradually sparse. Every mesh but the last (the fifth mesh in Table 2) were filled with 9-nodes elements. In the fifth mesh, 4-nodes finite elements were used. The calculations confirm that the J-integral is almost independent of the finite elements used in FEM analysis. The most important is the way the contour of integration is drawn. It should not lie too close to the crack tip nor 234 M. Graba, J. Gałkiewicz to the edge of the specimen.Thebest recommendation is to use a fewdifferent integral contours in FEM analysis and compare the results. Table 2.Definition of FE meshes used for calculation of the J-integral Number of elements between crack tip and Nodes per element Mesh maximun opening stress for critical moment Jc =40kN/m mesh 1 72 9 mesh 2 35 9 mesh 3 16 9 mesh 4 6 9 mesh 5 36 4 5. Determination of the crack tip opening displacement The Crack Tip Opening Displacement (CTOD) is determined using the me- thod proposed by Shih (1981) (see Fig.8). Themost important problem here is the choice of the crack tipmodel. In our calculations, we used threemodels of the crack tip. Fig. 8. Shih’s method proposed for determination of the crack opening displacement In the most common model called hereafter ”sharp”, the crack tip radius is equal to zero. In the secondmodel, the crack tip is modeled as a half of the circle (see Fig.5c) and in the last one, the crack tip is in a form of the quarter of the circle (see Fig.5d). In the third model, the initial value of CTOD is equal to the radius of the crack tip circle. The results obtained show (see Fig.9) that the best results, which are close to the theoretical plane strain value of CTOD (Shih, 1981), are obtained for the sharp crack tip model, and the worst results – with a half of the circle. Unfortunately, two of the threemodels are unable to determine theCTOD for small external loads, which inFig.9 are represented by the J-integral level. To Influence of the crack tip model... 235 determine the CTOD for the complete spectrum of external loads, we suggest to use the thirdmodel (see Fig.5d). Fig. 9. Changes of the CTOD for growing external load 6. Conclusions In the paper, the influence of finite elementmesh size and the crack tipmodel on FEM results was discussed. In particular, the stress distribution in front of the crack for the finite strain option was analyzed. It is important that the size of elements in front of the blunted crack should be less than 0.1 ·CTOD for the critical moment. The crack tip radius should be of the order of 10−6-10−7m.The size of the crack tip radius is not important if the stresses are measured at the distance r > J/σ0 andwhen it increases to the value three times greater than the initial one. We propose to use the crack tip model shown in Fig.5d when the large strain formulation is used. The calculations show that the J-integral can be determined using quite large finite elements if ”the virtual shift method” is used. The contour of integration, which is used to determine the J-integral, must not lie too close to the crack tip and to the edge of the specimen. It shouldnot cross the border between the plastic and elastic zone if it is possible. To calculate the crack tip openingdisplacement,we suggest using the crack tip model shown in Fig.5d. This model allows one to determine the CTOD value for the whole spectrum of external loads. 236 M. Graba, J. Gałkiewicz Acknowledgements We are pleased to acknowledge helpful interactions and discussions with prof. A. Neimitz fromKielce University of Technology. The work presented in this paper was conducted with support of Polish State Committee for Scientific Research, grant No. 5T07C00425. References 1. Al-Ani A.M., Hancock J.W., 1991, J-dominance of short cracks in tension and bending, J. Mech. Phys. Solids, 39, 1, 22-43 2. Brocks W., Cornec A., Scheider I., 2003, Computational aspects of nonlinear fracture mechanics,Bruchmechanik, GKSS-Forschungszentrum, Ge- esthacht, Germany, Elsevier, 127-209 3. Brocks W., Scheider I., 2003, Reliable J-values. Numerical aspects of the path-dependence of the J-integral in incremental plasticity, Bruchmechanik, GKSS-Forschungszentrum, Geesthacht, Germany, Elsevier, 127-209 4. Graba M., Gałkiewicz J., 2005, Wpływ modelu wierzchołka pęknięcia na wyniki uzyskane metodą elementów skończonych, Materiały Konferencyjne X Konferencji Mechaniki Pękania, 333-342,Opole,Wisła 5. Hutchinson J.W., 1968, Singular behaviour at the end of a tensile crack in a hardening material, Journal of the Mechanics and Physics of Solids, 16, 13-31 6. McMeeking R.M., ParksD.M., 1979,On criteria for J-dominance of crack tip fields in large scale yielding,ASTMSTP 668, American Society for Testing and Materials, Philadelphia, 175-194 7. O’Dowd N.P., 1995, Applications of two parameter approaches in elastic- plastic fracturemechanics,Engineering Fracture Mechanics, 52, 3, 445-465 8. O’Dowd N.P., Shih C.F., 1991, Family of crack-tip fields characterized by triaxiality parameter. I – Structure of fields, Journal of Mechanics and Physics of Solids, 39, 8, 989-1015 9. O’Dowd N.P., Shih C.F., 1992, Family of crack-tip fields characterized by triaxiality parameter. II – Fracture applications, Journal of Mechanics and Physics of Solids, 40, 5, 939-963 10. O’Dowd N., Shih C.F., Dodds R.H. Jr., 1995, The role of geometry and crack growth on constraint and implications of ductile/brittle fracture,ASTM STP 1244, American Society for Testing and Materials, Philadelphia 11. Rice J.R., Johnson M.A., 1970, The role of large crack tip geometry chan- ges in plane strain fracture, In: Inelastic Behaviour of Solids, M.F. Kanninen, McGraw-Hill, 641-672 Influence of the crack tip model... 237 12. Rice J.R., RosengrenG.F., 1968,Plane strain deformation near a crack tip in a power-law hardening material, Journal of the Mechanics and Physics of Solids, 16, 1-12 13. ShihC.F., 1981,Relationbetween theJ-integral and the crackopening displa- cement for stationary and extending cracks, Journal of Mechanics and Physics of Solids, 29, 305-329 Wpływ modelu wierzchołka pęknięcia na wyniki uzyskane metodą elementów skończonych Streszczenie WpracyprzeprowadzonoanalizęwpływumodeluMESnawartośćnaprężeńprzed frontem pęknięcia w materiałach sprężysto-plastycznych, wyznaczaną w sposób nu- meryczny wartość całki J oraz rozwarcie wierzchołka pęknięcia (RWP). Obliczenia prowadzono dla płaskiego stanu odkształcenia przy założeniu dużych odkształceń. Manuscript received June 14, 2006; accepted for print February 21, 2007