Microsoft Word - numero_43_art_4 L.C.H. Ricardo, Frattura ed Integrità Strutturale, 43 (2018) 57-78; DOI: 10.3221/IGF-ESIS.43.04 57 Crack Propagation by Finite Element Method Luiz Carlos H. Ricardo Materials Technology Department, IPEN, University of São Paulo, Brazil, Instituto de Pesquisas Energéticas e Nucleares Av. Lineu Prestes 2242 - Cidade Universitária - São Paulo - SP BRASIL- CEP: 05508-000. lricardo@ipen.br,https://orcid.org/0000-0002-1712-1437 ABSTRACT. Crack propagation simulation began with the development of the finite element method; the analyses were conducted to obtain a basic understanding of the crack growth. Today structural and materials engineers develop structures and materials properties using this technique. The aim of this paper is to verify the effect of different crack propagation rates in determination of crack opening and closing stress of an ASTM specimen under a standard suspension spectrum loading from FD&E SAE Keyhole Specimen Test Load Histories by finite element analysis. To understand the crack propagation processes under variable amplitude loading, retardation effects are observed. KEYWORDS. Fatigue; Crack propagation simulation; Finite element method; Retardation. Citation: L.C.H. Ricardo, Crack Propagation by Finite Element Method, Frattura ed Integrità Strutturale, 43 (2018) 57-78. Received: 02.10.2017 Accepted: 01.11.2017 Published: 01.01.2018 Copyright: © 2018 This is an open access article under the terms of the CC-BY 4.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. INTRODUCTION he most common technique for predicting the fatigue life of automotive, aircraft and wind turbine structures is Miner’s rule [1]. Despite the known deviations, inaccuracies and proven conservatism of Miner’s cumulative damage law, it is even nowadays being used in the design of many advanced structures. Fracture mechanics techniques for fatigue life predictions remain as a back up in design procedures. The most important and difficult problem in using fracture mechanics concepts in design seems to be the use of crack growth data to predict fatigue life. The experimentally obtained data is used to derive a relationship between stress intensity range (K) and crack growth per cycle (da/dN). In cases of fatigue loaded parts containing a flaw under constant stress amplitude fatigue, the crack growth can be calculated by simple integration of the relation between da/dN and K. However, for complex spectrum loadings, simple addition of the crack growth occurring in each portion of the loading sequence produces results that, very often, are more erroneous than the results obtained using Miner’s rule with an S-N curve. Retardation tends to cause conservative results using Miner’s rule when the fatigue life is dominated by the crack growth. However, the opposite effect generally occurs when the life is dominated by the initiation and growth of small cracks. In these cases, large cyclic strains, which might occur locally at stress raisers due to overload, may pre-damage the material and lower its resistance to fatigue. The experimentally derived crack growth equations are independent of the loading sequence and depend only on the stress intensity range and the number of cycles for that portion of the loading sequence. The central problem in the successful utilization of fracture mechanic techniques applied to the fatigue spectrum is to obtain a clear understanding T L.C.H. Ricardo, Frattura ed Integrità Strutturale, 43 (2018) 57-78; DOI: 10.3221/IGF-ESIS.43.04 58 of the influence of loading sequences on fatigue crack growth [2]. Investigations covering the effects of particular interest, after high overload, loading in the growth rate region, called crack growth retardation, seem to have little interest nowadays. Stouffer & Williams [3] and other researchers show a number of attempts to model this phenomenon through manipulation of the constants and stress intensity factors in the Paris-Erdogan equation however little appears to have been done in the effort to develop a completely rational analysis of the problem. Probably, the only one reason that the existing models of retarded crack growth are not satisfactory is that these models are deterministic whereas the fatigue crack growth phenomenon shows strong random features. In addition, most of the reported theoretical descriptions of the retardation are based on data fitting techniques, which tend to hide the behaviour of the phenomenon. If the retarding effect of a peak overload on the crack growth is neglected, the prediction of the material lifetime is usually very conservative [4]. Accurate predictions of the fatigue life will hardly become possible before the physics of the peak overload mechanisms is better clarified. According to the existing findings, the retardation is a physically very complicated phenomenon which is affected by a wide range of variables associated with loading, metallurgical properties, environment, etc., and it is difficult to separate the contribution of each of these variables [5]. CRACK PROPAGATION CONCEPTS rwin [6,7] defines in his work a release energy rate G, which is a measure of the available energy, dП-potential of energy and A-crack area, to provoke crack propagation as shown in Eq. (2.1). The term rate as employed is not related to a derivate in relation to the time but is referred to a change in the potential energy rate in the crack area. Later, this quantity has been called K, and is used to characterize the stress state ("stress intensity") near a crack tip caused by a remote load or residual stress in isotropic and elastic bodies. The stress field in the crack tip is given by Eq. (2.2),   d G dA (2.1)         1/2 1/2 2 3(2 ) ( ) ( ) ( ) ......ij ij ij ijK r f A g A h r (2.2) where K is the stress intensity factor; r and  are the distance from the crack tip and the angle between the crack tip and the plane of the crack, respectively; Ai is a constant of the material; fij (), gij ()and hij() are functions of . After years, the stress-intensity factors for a large number of crack configurations have been generated; and these have been collated into several handbooks (see, for example, Refs [8,9]). The use of K is meaningful only when small-scale yielding conditions exist. Plasticity and nonlinear effects will be covered in the next section.Because fatigue-crack initiation is, in general, a surface phenomenon, the stress-intensity factors for a surface- or corner-crack in a plate or at a hole, such as those developed by Raju and Newman [10,11], are solutions that are needed to analyze small-crack growth. Some of these solutions are used later to predict fatigue-crack growth and fatigue lives for notched specimens made of a variety of materials [12]. Paris & Erdogan [13] conducted a revision on the crack propagation approach from Head [14] and others and discussed the similarity of these theories and the differences of results between them, isolated and in group tests. Paris suggested that, for a cyclical load variation, the stress field in the crack tip for a cycle can be characterized by a variation of the stress intensity factor,   max minK K K (2.3) where Kmax and Kmin are the maximum and the minimum stress intensity factors, respectively. In the crack propagation curve, the linear part represents the Paris - Erdogan law, when plotting the values of K vs da/dN in logarithmic scale. Fatigue crack initiation and growth under cyclic loading conditions is controlled by the plastic zones that result from the applied stresses and exist in the vicinity (ahead) of a propagating crack and in its wake or flanks of the adjoining surfaces. For example, the fatigue characteristics of a cracked specimen or component under a single overload or variable amplitude loading situations are significantly influenced by these plastic zones. In modelling the fatigue crack growth rate this is accounted by the incorporation of accumulative damage cycle after cycle and should include plasticity effects. During the crack propagation the plastic zone should grown and the plastic wake will have compressive plastic zones that can help to keep the crack close. Hairman & Provan [15] discuss the problems pertaining to fatigue loading of engineering structures under single overload and variable amplitude loading involving the estimation of plasticity affected zones ahead of the crack tip. I L.C.H. Ricardo, Frattura ed Integrità Strutturale, 43 (2018) 57-78; DOI: 10.3221/IGF-ESIS.43.04 59 Crack tip plasticity Most solid materials develop plastic strains when the yield strength is exceeded in the region near a crack tip. Thus, the amount of plastic deformation is restricted by the surrounding material, which remains elastic during loading. Theoretically, linear elastic stress analysis of sharp cracks predicts infinite stresses at the crack tip. In fact, inelastic deformation, such as plasticity in metals and crazing in polymers, leads to relaxation of crack tip stresses caused by the yielding phenomenon at the crack tip. As a result, a plastic zone is formed containing microstructural defects such dislocations and voids. Consequently, the local stresses are limited to the yield strength of the material. This implies that the elastic stress analysis becomes increasingly inaccurate as the inelastic region at the crack tip becomes sufficiently large and linear elastic fracture mechanics (LEFM) is no longer useful for predicting the field equations. The size of the plastic zone can be estimated when moderate crack tip yielding occurs. Thus, the introduction of the plastic zone size as a correction parameter that accounts for plasticity effects adjacent to the crack tip is vital in determining the effective stress intensity factor (Keff) or a corrected stress intensity factor. The plastic zone is also determined for plane conditions; that is, plane strain for maximum constraint on relatively thick components and plane stress for variable constraint due to thickness effects of thin solid bodies. Moreover, the plastic zone develops in most common in materials subjected to an increase in the tensile stress that causes local yielding at the crack tip. Most engineering metallic materials are subjected to an irreversible plastic deformation. If plastic deformation occurs, then the elastic stresses are limited by yielding since stress singularity cannot occur, but stress relaxation takes place within the plastic zone. This plastic deformation occurs in a small region and it is called the crack-tip plastic zone. A small plastic zone, (r << a) is referred to as small-scale yielding. On the other hand, a large-scale yielding corresponds to a large plastic zone, which occurs in ductile materials in which r >> a. This suggests that the stress intensity factors within and outside the boundary of the plastic zone are different in magnitude so that KI (plastic) > KI (elastic). In fact, KI (plastic) must be defined in terms of plastic stresses and displacements in order to characterize crack growth, and subsequently ductile fracture. As a consequence of plastic deformation ahead of the crack tip, the linear elastic fracture mechanics (LEFM) theory is limited to r << a; otherwise, elastic-plastic fracture mechanics (EPFM) theory controls the fracture process due to a large plastic zone size (r ≥ a). This argument implies that r may be determined in order to set an approximate limit for both LEFM and EPFM theories. Fig. 1.b shows schematic plastic zones for plane stress (thin plate) and plane strain (thick plate) conditions [16]. Plane strain: 1. Large thickness B, and εz ~ 0 on in an internal region and. σz = υ(σ x + σ y). This means that the material is constrained in the z-direction due to a sufficiently large thickness and the absence of strain in this axis. In fact, the stress in the z- direction develops due to the Poisson’s effect as explicitly included in the equation that defines σz. 2. Yielding is suppressed due to the kinematics constrain from the surrounding elastic material. 3. Plastic deformation is associated with the hinge mechanism (internal necking) Fig. 1.a) 4. The plastic zone size is small in the midsection of the plate (Fig. 1.a).This condition implies that the plastic zone must be smaller than the crack length Plane stress: 1. The thickness B is small, σz = 0 and εz ≠ 0 on the surface (external region) and through the whole thickness. This means that the stresses normal to the free surface are absent and therefore, σz = 0 through the thickness. Consequently, a biaxial state of stress results. 2. If σ y ≥ σ x>0 (Tresca Criterion), then yielding occurs by a cumulative slip mechanism (Fig. 1.b). 3. The height of the yielded zone is limited due to the slip mechanism. 4. The total motion has a necking effect in front of the crack as it opens.           2 1 2 I p ys K r Plane Stress B ≤ 2.5( KIC/σ ys)2 (2.4) Irwin [6] has shown that the effect on the plastic zone is to artificially extend the crack by a distance r1 (Fig. 2) known as Irwin’s plastic zone correction. The elastic stress distribution shown in Fig. 2 indicates that as σy → ∞. Actually, σy is limited to σys as shown by the elastic-plastic stress distribution. This means that σy → ∞ occurs mathematically, not physically. In order to account for the changes due to the artificial crack extension or virtual crack length and to visualize the plastic zone as r → 0 a cylinder, the crack length a can be replaced by ae in eqs. (2.4 and 2.5). Moreover, the virtual L.C.H. Ricardo, Frattura ed Integrità Strutturale, 43 (2018) 57-78; DOI: 10.3221/IGF-ESIS.43.04 60 crack length defined by ae is referred to as the effective crack length in the literature. The conditions of equilibrium for an immobile crack tip include internal and external forces per unit length [15,16]. In such a case, the areas related to the shedding loads Ps and Pys due to yielding, as indicated in Fig. 2, are equal; that is APs=APys when the plastic zone size is r << a Mathematically, these loads are the equilibrium forces per unit length defined by [16]. Figure 1: Yielding Mechanism of a Plate [6, 7]. Figure 2: Crack Tip Plastic Zone [16].     1 0 r s ysP B dx (2.5) L.C.H. Ricardo, Frattura ed Integrità Strutturale, 43 (2018) 57-78; DOI: 10.3221/IGF-ESIS.43.04 61   2 0 r ys ysP B dx (2.6) where B=thickness λ = 1 for plane stress λ = 3 for plane Irwin’s yielding factor for plane strain [7] For equilibrium conditions, the force balance    0s ysP P leads to the determination of the of the plastic zone size Hence,        1 2 0 0 0 r r ys ysdx dx (2.7) Considering:           1 0 lim 2 2yy I yy r K r f r for     , 0yy yy r (2.8) Inserting (2.2) into (2.8) and integrating yields            1 2 0 0 0 2 r r I ys ys K dx dx r (2.9)      1 1 2 2 0 2 I ys r K r r r (2.10)     1 1 22 0y ysr r r (2.11) The elastic stress can be defined by  y ys (2.12) Inserting eq. (2.12.) into (12.3) gives 2r1= r1+r2 which implies that r1=r2 and from Fig. 2, r1= r1+r2. Hence, ae= a+r is the virtual crack length proposed by Irwin [6]. Obviously, eq. (2.14) provides the effective stress intensity factor  IK a (2.13)       I eK a r a (2.14) The plastic zone size can be calculated by eqs. 2.4 and 2.5. This KI equation is the corrected stress intensity factor due to finite specimen size and plasticity. Now, inserting eqs. (2.4) into (2.14) yields.           2 2 ys a r (2.15) L.C.H. Ricardo, Frattura ed Integrità Strutturale, 43 (2018) 57-78; DOI: 10.3221/IGF-ESIS.43.04 62 where σ = Applied stress (MPa) σys= Yield strength (MPa) a = Crack length (m)                1 1 2I ys K a (2.16) Furthermore, the plastic zone size for plane conditions can easily be determined by combining eqs. (2.4) and (2.12). Thus,                     2 2 1 2 2 I ys ys K a r (2.17) In plane strain condition, yielding is suppressed by the triaxial state of stress and the plastic zone size is smaller than that for plane stress as predicted by the α parameter in eq. (2.17). The same reasoning can be used for mode III. Thus, the plastic zone becomes [12].           2 1 2 III ys K r (2.18) Dugdale’s approximation Dugdale [17] proposed a strip yield model for the plastic zone under plane stress conditions. Consider Fig. 3 which shows the plastic zones in the form of narrow strips extending a distance r each, and carrying the yield stress σys The phenomenon of crack closure is caused by internal stresses since they tend to close the crack in the region where a < x < c. Furthermore, assume that stress singularities disappear when the following equality is true Kσ = - KI, where Kσ is the applied stress intensity factor and KI is due to yielding ahead of the crack tip [6]. Hence, the stress intensity factors due to wedge internal forces are defined by      a r A a P a x K dx a xa (2.19)      a r B a P a x K dx a xa (2.20) According to the principle of superposition, the total stress intensity factor is KI = KA+ KB so that             a r I a P a x a x K dx a x a xa (2.21)   2 cosI a x K P ar a (2.22) The plastic zone correction can be accomplished by replacing the crack length a for the virtual crack length ( a+r ), and P for σys Thus, the stress intensity factor are: L.C.H. Ricardo, Frattura ed Integrità Strutturale, 43 (2018) 57-78; DOI: 10.3221/IGF-ESIS.43.04 63 Figure 3: Dugdale Plastic Zone Strip Model [17]          2 cos ys I x K a r ar a r (2.23)     K a r (2.24) But, Kσ = KI and the simplified equation takes the form     cos 2 ys x ar a r (2.25) Let    2 ys y so that   cos x y a r (2.26)  (sec 1)r a y (2.27) Expanding the trigonometric function eq. (2.27) yields and neglecting the high order terms from (2.27) becomes:            2 2 2 2 2 ys ay a r (2.28) Substituting eq. (2.28) into (2.14) gives the corrected stress intensity factor due plasticity at the crack tip and crack geometry L.C.H. Ricardo, Frattura ed Integrità Strutturale, 43 (2018) 57-78; DOI: 10.3221/IGF-ESIS.43.04 64                2 1 2I ys K a (2.29) Expression (2.29) is similar to Irwin’s expression, eq. (2.16). In addition, if r << a, plasticity corrections are not necessary. Fig. 4 compares the normalized stress intensity factors as per Irwin’s and Dugdale’s approximations. The curves significantly differ as σ/σys → 1; however, similarities occur at r < σ/σys ≤ 0.2. This strongly suggests that both Irwin’s and Dugdale’s approximation methods should be used very carefully because their differences in normalized stress intensity factor . Figure 4: Normalized Stress Intensity factor as function of stress ratio [16].          min max ( ) 1 m c C Kda dN K K K K    max ( )m c C Kda dN K K  1 max( ) ( )m mda C K K dN Table 1: Empirical Crack Growth Equations for Constant Amplitude Loading [18]. The models of Irwin [6,7] and Dugdale [17] give an idea of the size of the plastic zone but not of its shape. The size, in general, is estimated as a circle of certain diameter (ryor rp) obtained on the basis of reasoning given in the above models for crack-tip-plasticity. In these models the effect of the shape of the plasticity affected zones is not taken into account. In the original Paris crack propagation equation [18] the driving parameters are C, K and m. In Tab. 1 it is possible to see some other crack propagation equations for constant amplitude loading, which are modifications of the Paris equation, relating the mentioned parameters. Murthy et al. [19] discuss crack growth models for variable amplitude loading and the mechanisms and contribution to overload retardation. Tab. 2 presents some authors and the application of their models. Retardation phenomenon Corbly & Packman [34] present some aspects of the retardation phenomenon some of which are presented below. 1. Retardation increases with higher values of peak loading peak for constant values of lower stress levels [35,36]. 2. The number of cycles at the lower stress level required to return to the non-retarded crack growth rate is a function of Kpeak, Klower, R peak, Rlower and number of peak cycles [37]. 3. If the ratio of the peak stress to lower stress intensity factors is greater than l.5 complete retardation at the lower stress intensity range is observed. Tests were not continued long enough to see if the crack ever propagated again [37]. 4. With a constant ratio of peak to lower stress intensity the number of cycles to return to non-retarded growth rates increases with increasing peak stress intensity [36,37]. L.C.H. Ricardo, Frattura ed Integrità Strutturale, 43 (2018) 57-78; DOI: 10.3221/IGF-ESIS.43.04 65 5. Given a ratio of peak stress to lower stress, the number of cycles required to return to non-retarded growth rates decreases with increased time at zero load before cycling at the lower level [37]. 6. Increased percentage delay effects of peak loading given a percent overload are greater at higher baseline stress intensity factors [38]. 7. Delay is a minimum if compression is applied immediately after tensile overload [39]. 8. Negative peak loads cause no substantial influence of crack growth rates at lower stress levels if the values of R > 0 for the lower stress [40]. 9. Negative peak loads cause up to 50 per cent increase in fatigue crack propagation with R = - 1 [39]. 10. Importance of residual compressive stresses around the tip of crack [41] 11. Low-high sequences cause an initial acceleration of the crack propagation at the higher stress level which rapidly stabilizes [42]. Yield Zone Concept Crack Closure Concept Wheeler [20] Elber [27] Willenborg, Engle, Wood [21] Bell and Creager (Generalized Closure) [28] Porter [22] Newman (Finite Element Method) [29] Gray (Generalized Wheeler) [23] Dill and Staff (Contact Stress ) [30] Gallagher and Hughes [24] Kanninen, Fedderson, Atkinson [31] Johnson [25] Budiansky and Hutchinson [32] Chang et al. [26] de Koning [33] Table 2: Fatigue Crack Growth Models [19]. Small Scale Yield Models While the basic layout of the small scale yield model has been established by Newman [29] and this approach was applicable to general variable amplitude loading. The small scale yield model employs the Dugdale [17] theory of crack tip plasticity modified to leave a wedge of plastically stretched material on the fatigue crack surfaces. The fatigue crack growth is simulated by severing the strip material over a distance corresponding to the fatigue crack growth increment as shown Fig. 5. In order to satisfy the compatibility between the elastic plate and the plastically deformed strip material, a traction must be applied on the fictitious crack surfaces in the plastic zone (a  x > /ColorImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000ColorACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000ColorImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasGrayImages false /CropGrayImages true /GrayImageMinResolution 300 /GrayImageMinResolutionPolicy /OK /DownsampleGrayImages true /GrayImageDownsampleType /Bicubic /GrayImageResolution 300 /GrayImageDepth -1 /GrayImageMinDownsampleDepth 2 /GrayImageDownsampleThreshold 1.50000 /EncodeGrayImages true /GrayImageFilter /DCTEncode /AutoFilterGrayImages true /GrayImageAutoFilterStrategy /JPEG /GrayACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /GrayImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000GrayACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000GrayImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasMonoImages false /CropMonoImages true /MonoImageMinResolution 1200 /MonoImageMinResolutionPolicy /OK /DownsampleMonoImages true /MonoImageDownsampleType /Bicubic /MonoImageResolution 1200 /MonoImageDepth -1 /MonoImageDownsampleThreshold 1.50000 /EncodeMonoImages true /MonoImageFilter /CCITTFaxEncode /MonoImageDict << /K -1 >> /AllowPSXObjects false /CheckCompliance [ /None ] /PDFX1aCheck false /PDFX3Check false /PDFXCompliantPDFOnly false /PDFXNoTrimBoxError true /PDFXTrimBoxToMediaBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXSetBleedBoxToMediaBox true /PDFXBleedBoxToTrimBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXOutputIntentProfile () /PDFXOutputConditionIdentifier () /PDFXOutputCondition () /PDFXRegistryName () /PDFXTrapped /False /CreateJDFFile false /Description << /ARA /BGR /CHS /CHT /CZE /DAN /DEU /ESP /ETI /FRA /GRE /HEB /HRV (Za stvaranje Adobe PDF dokumenata najpogodnijih za visokokvalitetni ispis prije tiskanja koristite ove postavke. Stvoreni PDF dokumenti mogu se otvoriti Acrobat i Adobe Reader 5.0 i kasnijim verzijama.) /HUN /ITA /JPN /KOR /LTH /LVI /NLD (Gebruik deze instellingen om Adobe PDF-documenten te maken die zijn geoptimaliseerd voor prepress-afdrukken van hoge kwaliteit. De gemaakte PDF-documenten kunnen worden geopend met Acrobat en Adobe Reader 5.0 en hoger.) /NOR /POL /PTB /RUM /RUS /SKY /SLV /SUO /SVE /TUR /UKR /ENU (Use these settings to create Adobe PDF documents best suited for high-quality prepress printing. Created PDF documents can be opened with Acrobat and Adobe Reader 5.0 and later.) >> /Namespace [ (Adobe) (Common) (1.0) ] /OtherNamespaces [ << /AsReaderSpreads false /CropImagesToFrames true /ErrorControl /WarnAndContinue /FlattenerIgnoreSpreadOverrides false /IncludeGuidesGrids false /IncludeNonPrinting false /IncludeSlug false /Namespace [ (Adobe) (InDesign) (4.0) ] /OmitPlacedBitmaps false /OmitPlacedEPS false /OmitPlacedPDF false /SimulateOverprint /Legacy >> << /AddBleedMarks false /AddColorBars false /AddCropMarks false /AddPageInfo false /AddRegMarks false /ConvertColors /ConvertToCMYK /DestinationProfileName () /DestinationProfileSelector /DocumentCMYK /Downsample16BitImages true /FlattenerPreset << /PresetSelector /MediumResolution >> /FormElements false /GenerateStructure false /IncludeBookmarks false /IncludeHyperlinks false /IncludeInteractive false /IncludeLayers false /IncludeProfiles false /MultimediaHandling /UseObjectSettings /Namespace [ (Adobe) (CreativeSuite) (2.0) ] /PDFXOutputIntentProfileSelector /DocumentCMYK /PreserveEditing true /UntaggedCMYKHandling /LeaveUntagged /UntaggedRGBHandling /UseDocumentProfile /UseDocumentBleed false >> ] >> setdistillerparams << /HWResolution [2400 2400] /PageSize [612.000 792.000] >> setpagedevice