Influence of the crack propagation rate in the obtaining opening and closing stress intensity factor by finite element method

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 as criterion design. 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. The crack propagation simulation was based on release nodes at the minimum loads to minimize convergence problems. To understand the crack propagation processes under variable amplitude loading, retardation effects are discussed.


INTRODUCTION
he most common technique for predicting the fatigue life of automotive, aircraft, wind turbine and many other 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, 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 T effect generally occurs when the life is dominated by the initiation and growth of small cracks (linear elastic fracture mechanics).Large cyclic strains (elasto-plastic fracture mechanics), which might occur locally at stress raisers due to overload, may predamage 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 of the influence of loading sequences on fatigue crack growth [2].Investigations covering the effects of particular interest, after high overload, in the study of crack growth under variable-amplitude 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 used 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 behavior 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
aris & Erdogan [6] conducted a revision on the crack propagation approach from Head [7] 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, Eq. (2.1), P max min where K max and K min 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, as can be seen in Fig. 1.
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 [9] 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.

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, so, 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 (r) as a correction parameter, which accounts for plasticity effects adjacent to the crack tip, is vital to determine 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 (r).A small plastic zone, (r << a) being a crack length of the structure or specimen, 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, K I (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 should be determined in order to set an approximate limit for both LEFM and EPFM theories [10].Fig. 2 shows schematic plastic zones for plane stress (thin plate) and plane strain (thick plate) conditions Even if in the interior of a plate a condition of plane strain exists, there will always be plane stress at the surface.Stresses perpendicular to the outer surface are nonexistent, and hence σ z = σ 3 = 0 at the surface.If plane strain prevails in the interior of the plate, the stress σ 3 gradually increases from zero (at the surface) to the plane strain value in the interior [11].Consequently, the plastic zone gradually decreases from the plane stress size at the surface, to the plane strain size in the interior of the plate, illustrated schematically in Fig. 2.
Dugdale [12] 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 .In the case of cyclic stress, and as the crack grows, behind the new crack front there is a region with compressive (residual) stresses.The phenomenon of crack closure is caused by these internal (compressive) stresses since they tend to close the crack in the region where a < x < c.In the original Paris crack propagation equation [6] the driving parameters are C, K and m, as shown in Fig. 1.Among other limitations, this equation is valid only in the region (b).So, it does not cover the near threshold region (a) nor the unstable region (c).Some researchers have proposed similar equations that cover one or both extremes of the curve in Fig. 1.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 and Kc, the critical stress intensity factor.Murthy et al. [13] discuss crack growth models for variable amplitude loading and the mechanisms and contribution to overload retardation.There are many authors which have been developing fatigue crack growth models for variable amplitude loading.Tab. 2 presents some authors and the application of their models.

Retardation phenomenon
It has been noted that, under certain conditions, the crack growth presents a slower rate, called retardation, due to several factors.Despite recent large increase in research into the retardation effects in crack propagation there are many aspects of load interaction phenomena that lack adequate explanations.It is presented here some several aspects of the retardation phenomenon by Corbly & Packman [28]. 1. Retardation increases with higher values of peak loading  peak for constant values of lower stress levels [29,30].
2. The number of cycles at the lower stress level required to return to the non-retarded crack growth rate is a function of K peak , K lower , R peak , R lower and number of peak cycles [31].3.If the ratio of the peak stress to lower stress intensity factors is greater than l.5 retardation at the lower stress intensity range is observed.Tests were not continued long enough to see if the crack ever propagated again [31].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 [30,31]. 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 [31].6. Increased percentage delay effects of peak loading given a percent overload are greater at higher baseline stress intensity factors [32].7. Delay is a minimum if compression is applied immediately after tensile overload [33].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 [34].9. Negative peak loads cause up to 50 per cent increase in fatigue crack propagation with R = -1 [33].10.Importance of residual compressive stresses around the tip of crack [35] 11.Low-high sequences cause an initial acceleration of the crack propagation at the higher stress level which rapidly stabilizes [36].

Small Scale Yield Models
While the basic layout of the small scale yield model has been established by Dill & Saff [37], only improvements introduced later by Newman [38] made this approach applicable to general variable amplitude loading.The small scale yield model employs the Dugdale [12] 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. 4. 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 < a afict ), as in the original Dugdale model, and also over some distance in the crack wake (a open  x < a), where the plastic elongations of the strip L(x) exceed the fictitious crack opening displacements V(x).The compressive stress applied in the crack wake to insure L(x)=V(x) are referred to as the contact stresses.Ricardo et al. [39] discuss the importance in the determination of materials properties like crack opening and closing stress intensity factor.The development of crack closure mechanisms, such plasticity, roughness, oxide, corrosion, and fretting product debris, and the use of the effective stress intensity factor range, has provided an engineering tool to predict small and large crack growth rate behavior under service loading conditions.The major links between fatigue and fracture mechanics were done by Elber [21].The crack closure concept put crack propagation theories on a firm foundation and allowed the development of practical life prediction for variable and constant amplitude loading, by such as experienced by modern day commercial aircrafts.Numerical analysis using finite elements has played a major role in the stress analysis crack problems.Swedlow [40] was one of the first to use finite element method to study the elastic-plastic stress field around a crack.
The application of linear elastic fracture mechanics, i.e. the stress intensity factor range, K, to the "small or short" crack growth have been studied for long time to explain the effects of nonlinear crack tip parameters.The key issue for these nonlinear crack tip parameters is crack closure.Analytical models were developed to predict crack growth and crack closure processes like Dugdale [12], or strip yield, using the plasticity induced approach in the models considering normally plane stress or strain effects.Schijve [41], discussing the relation between short and long cracks presented also the significance of crack closure and growth on fatigue cracks under services load histories.The ultimate goal of prediction models is to arrive at quantitative results of fatigue crack growth in terms of millimeters per year or some other service period.Such predictions are required for safety and economy reasons, for example, for aircraft and automotive parts.Sometimes the service load time history (variable amplitude loading) is much similar to constant amplitude loading, including mean load effects.In both cases quantitative knowledge of crack opening stress level S op is essential for crack growth predictions, because Keff is supposed to be the appropriate field parameter for correlating crack growth rates under different cyclic loading conditions.The correlation of crack growth data starts from the similitude approach, based on the Keff, which predicts that same Keff cycles will produce the same crack growth increments.The application of Keff to variable amplitude loading require prediction of the variation of S op , during variable amplitude load history, which for the more advanced prediction models implies a cycle by cycle prediction.The Fig. 5 shows the different K values.The application of Keff is considerably complicated by two problems: (1) small cracks and (2) threshold K values (Kth).Small cracks can be significant because in many cases a relatively large part of the fatigue life is spent in the small crack length regime.The threshold problem is particularly relevant for fatigue under variable amplitude spectrum, if the spectrum includes many "small" cycles, those ones with small stress/load amplitude.
It is important to know whether the small cycles do exceed a threshold K value, and to which extension it will occur.The application of similitude concept in structures can help so much, but the results correlation is not satisfactory and the arguments normally are:  The similarity can be violated because the crack growth mechanism is no longer similar. The crack can be too small for adopting K as a unique field parameter.
 Keff and others conditions being nominally similar, it is possible that other crack tip aspects also affect crack growth, such as crack tip blunting and strain hardening, Schijve 45 .Newman and Armen [42][43][44] and Ohji et al. [45] were the first to conduct the two dimensional analysis of the crack growth process.Their results under plane stress conditions were in quantitative agreement with experimental results by Elber [21], and showed that crack opening stresses were a function of R ratio (Smin/Smax) and the stress level (Smax/0), where  0 is the flow stress i.e: the average between  ys and  u (ultimate stress).
Blom and Holm [46] and Fleck and Newman [47][48] studied crack growth and closure under plane-strain conditions and found that cracks did close but the cracks opening levels were much lower than those under plane stress conditions considering same loading condition.Sehitoglu et al. [49] found later that the residual plastic deformations cause the crack closing.McClung [50][51][52] performed extensive finite element crack closure calculations on small cracks at holes, and various fatigue crack growth models.Solanski et.al [53] found that S max / 0 could correlate the crack opening stresses for different flow stresses ( 0 ).This average value was used as stress level in the plastic zone for the middle crack tension specimen McClung [52] found that K analogy, using K max /K 0 could correlate the crack opening stresses for different crack configurations for small scale yielding conditions where K0=o(a) .(K-analogy assumes that the stress-intensity factor controls the development of closure and crack-opening stresses, and that by matching the K solution among different cracked specimens, an estimate can be made for the crack opening stresses.)

DESCRIPTION OF THE MODEL
compact tension specimen was modeled using a finite element code, MSC/Patran, r1 [54] and ABAQUS Version 68 [55] used as solver.Half of the specimen was modeled and symmetry conditions applied.A plane stress constraint is modeled by the finite element method covering the effects in two dimensional (2D) small scale yielding models of fatigue crack growth under variable spectrum loading, Fig. 7, and the boundary conditions are presented in Fig. 6.The finite element model has triangle and quadrilateral elements with quadratic formulation and spring elements, SPRING1, used to node release in crack surface (this element works only in the y direction).Fatigue Design & Evaluation (FD&E) committee from SAE (Society of Automotive Engineers) has standard fatigue files.The present work used a standard suspension load history.Fig. 7 presents a modified load history, adapted from the FD&E/SAE histogram considering only tractive loads.The maximum load used was scaled to produce a K max  0.6 K IC , using Eq.(4.1), where KIC is the critical stress intensity factor of adopted material in the present study.With the value of Kmax from KIC computed as mentioned above is computed the maximum load using Eq.3.1 to be applied in the specimens as explained in next.

A
In the analysis, fatigue crack growth is simulated by releasing the crack tip node at Pmin, followed by a single loading cycle P min  P max  P min, Fig. 7.The force is divided in nine steps between loads P min -P max and nine steps between the P max -P min , in each cycle.The smallest element size, 0.025 mm, was estimated based on the plastic zone size (rp) ahead of the crack tip and computed by Eq. (3.2).Only the first 20 reverses from load history shown in Fig. 7 were used to identify crack opening/closing and retardation effects.
where: Kmax= maximum stress intensity factor; Pmin= minimum applied load; P max = maximum applied load; B = specimen thickness; a = crack length; W = width of the specimen; a/W = ratio of the crack length to the specimen width; f(a/W) = characteristic function of the specimen geometry .
Antunes & Rodrigues [56] discuss that numerical analysis of plastic induced crack closing (PICC) based on finite element method (FEM) consists of discretising and modeling the cracked body having elastic-plastic behaviour, applying a cyclic load, extending the crack and measuring the crack closure level.The finite element mesh must be highly refined near the crack front, with micron scale, in order to model the forward and reversed crack tip plastic zones.The forward plastic zone is made up of the material near the crack tip undergoing plastic deformation at the maximum load, therefore it is intimately related to Kmax.The reversed plastic zone encompasses the material near the crack tip undergoing compressive yielding at the minimum load and is related to ΔK.Commercial FE software packages offer tools to deal with elasticplastic deformation, crack propagation and contact between crack flanks, and are therefore adequate to model PICC.However, the numerical models have significant simplifications with respect to real fatigue crack propagation, namely: -discrete crack propagations, of the same size as near crack tip elements, which give fatigue crack growth rates significantly higher than real values; -crack propagation is modeled at a constant load when in reality it occurs continuously during the whole load cycle.In numerical simulations, the crack can be incremented at maximum load [57], at minimum load [58,59] or at other positions of the load cycle.Ogura et al. [59] advanced the crack when the crack tip reaction force reached zero during the load cycle.However, none of these approaches truly represents the fatigue process, where, according to slip models of striation formation, crack extension is a progressive process occurring during the entire load cycle.The proposal to increment at minimum load was designed to overcome convergence difficulties caused by propagating the crack at maximum load.This is unrealistic since the crack is not expected to propagate in a compressive stress field.However, several authors [60,61] have already found that the load at which the crack increment occurs does not significantly influence crack closure numerical results.Under constant amplitude loading, crack tip opening load will typically increase monotonically, with increasing crack growth, until a stabilized value is reached.So, it is important to define the minimum crack extension needed to stabilize the opening level.
It is usually sufficient to increase the crack ahead of the monotonic plastic zone resulting from the first load cycle [62,63].The stress level in the crack tip, Fig. 8, must to be positive to characterize the crack opening and negative to characterize the crack closure.Antunes & Rodrigues [56] consider as basic criteria to determine the crack opening or closing: the first contact of the crack flank, which corresponds to the contact of the first node behind the current crack tip.This is the conventional definition proposed by Elber [21] and has been widely used by Jiang et al. [64].
In this work the nodes released in the crack tips were located at the minimum load of a cycle to simulate crack growth and will be considered the first contact of the node behind the crack tip, positive stress (+Syy) to characterize the crack opening and negative stress (-Syy) to characterize the crack closing.

DISCUSSION OF RESULTS
n the present work was identified how difficult is to determine with proper precision the crack opening or closure.It was necessary to use the iterative process in the crack surface step by step during loading and unloading to find the crack opening or closing as shown in details in Ricardo [65].The retard effect is present in some cycles in special cases where there are overloads.In constant amplitude loading, the effective plastic zone increases with the extension of the crack length; the crack propagation rate has no influence in the quality of results, assuming that it is in respect to the Newman [23] recommendation with four elements yielded in the reverse plastic zone.In variable amplitude loading the crack length cannot progress until a new overload occurs or the energy spent during cyclic process creates a new plastic zone and the driving force increases the crack length.The researchers normally work with simple overloads or specific load blocks; this approach can induce some mistakes in terms of results that can be conservative or nonrealistic.Fig. 9 shows the effect of different crack propagation rates in opening stress, op.This graph starts in the second cycle because it was not possible identify the crack opening in all models evaluated when the crack opens, because all stresses in the first cycle were positive.In the beginning there is no representative difference in the four first cycles in all crack propagation models.In the fourth to fifth cycle it is possible identify a difference of crack open stress level from model SAE2 (crack propagation 0.5 mm/cycle) and the others models.The difference of the crack opening stress level from model SAE2 from the others may be related with the overload that the specimen had in the fifth cycle causing the increase of the crack opening stress level to be more representative than in others that suffered the same overload.From the sixth to eight cycles it is possible to identify again little difference in the crack opening stress of the models.The model SAE1 (crack propagation 0.025 mm/cycle) has the lower crack opening stress.In the cycles 8 to 10 there is some difference in the crack opening stress, having as principal cause the different plasticity that the models suffered, due to I different crack propagation rate models.Model SAE2 has the bigger crack opening stress; caused like in the fifth cycle by an overload as in the fifth cycle and again this model had different behavior when compared with others models.The model SAE3 (crack propagation rate 0.75 mm) has no significant difference in the crack opening stress level during all cycles.This could be a good indication that for a first approach in similar conditions the utilization of this crack propagation rate will provide the behavior material faster under similar load history and specimen.Fig. 9 also shows that it is possible to have more different kinds of criteria design.For example for a conservative approach it is possible the utilization of the model SAE1 (crack propagation rate 0.25 mm/cycle).

Surface Crack
Crack Tip y= 103 MPa Fig. 10 presents the results from the crack closing stress against numbers of cycles evaluated for the four different crack propagation models considered.It is possible to observe that in the first four cycles there are no significant difference in the crack closing stress in the models studied.In the other cycles the model SAE1 (crack propagation 0.25 mm/cycle), has no significant difference of crack closing stress during crack propagation.In fact it is the most conservative model from the four evaluated.During the fourth and sixth cycle the models SAE2 (crack propagation model 0.50 mm) and SAE3 (crack propagation model 0.75 mm) have no difference in the crack closing stress.The model SAE4 (crack propagation 1.0 mm/cycle) has representative difference in the crack closing stress when compared with others models in the cycles due to more residual plasticity in the crack tip.The last representative differences between crack closing stress levels in the models happen during propagation in the cycles eight to tenth.An increase of the crack propagation rate will also increase the crack closing stress.Fig. 12 shows that depending on the design criterion it is possible to apply a different crack propagation rate.For example if the criterion is to use a conservative crack closing stress it is recommended utilization of the model SAE1 (crack propagation 0.25 mm).The softest model or that one which allows the bigger crack opening and closing stresses is model SAE4 (crack propagation model 1.0 cycle/mm).
CONCLUSIONS n this work it was possible to identify the crack opening and closure using the finite element method.In the literature there are few works covering crack propagation simulation with random loads like FD&E loads histories from SAE data bank.Normally only a few load blocks are used to reduce the complexity; this should provide conservative answers when used to develop structural components.The use of different crack propagation rate in this work shows that for reproducing the effective plastic zone it is possible to use smaller or larger element sizes compared with the Irwin equation.To improve the correlation between numerical and experimental data it is necessary to increase the crack length to obtain the same qualitative results that is estimated by the Irwin equation.The next step in this work will be to perform some analyses with the same model and load history with different crack propagation rates to identify whether or not the retard effect can be observed.These data will be compared with experimental test and, if necessary, adjustment of the crack propagation model will be done to improve the crack propagation model.
Figs. 11 and 12 present examples of results of post-processing results from model SAE0.50 showing the stress field in the region near where the crack opens and closes.

Table 3 :
4aterial Properties of a Low Alloy Steel.The dimensions of the compact tension specimen were: B= 3.8 mm; W= 50.0 mm; a/W= 0.26.Tab.4shows the estimated and used values of the cyclic plastic zone sizes as well as smaller finite element.Tab. 5 shows the difference crack propagation rates used in the current work.