Crack Propagation 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. 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.


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 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), 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, 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.

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 K I (plastic) > K I (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 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 zdirection 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.
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 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 P s and P ys due to yielding, as indicated in Fig. 2, are equal; that is AP s =AP ys when the plastic zone size is r << a Mathematically, these loads are the equilibrium forces per unit length defined by [16].  where B=thickness λ = 1 for plane stress λ = 3 for plane Irwin's yielding factor for plane strain [7] For equilibrium conditions, the force balance Inserting eq. (2.12.) into (12.3) gives 2r 1 = r 1 +r 2 which implies that r 1 =r 2 and from Fig. 2, r 1 = r 1 +r 2. Hence, a e = a+r is the virtual crack length proposed by Irwin [6]. Obviously, eq. (2.14) provides the effective stress intensity factor 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.
Furthermore, the plastic zone size for plane conditions can easily be determined by combining eqs. (2.4) and (2.12). Thus, 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].

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 σ = -K I, where K σ is the applied stress intensity factor and K I is due to yielding ahead of the crack tip [6]. Hence, the stress intensity factors due to wedge internal forces are defined by According to the principle of superposition, the total stress intensity factor is 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: Substituting eq. (2.28) into (2.14) gives the corrected stress intensity factor due plasticity at the crack tip and crack geometry 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]. 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].
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].

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 <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. The fatigue crack growth is simulated using the strip material as shown schematically in Fig. 5. Ricardo et al. [43] 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 Christensen [44] and Elber [45]. 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 [35] 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 [17], or strip yield, using the plasticity induced approach in the models considering normally plane stress or strain effects. Schijve [36], 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 Sop 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 K eff 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. 6 shows the different K values. The application of K eff is considerably complicated by two problems: (1) small cracks and (2) threshold K values (K th ).
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. 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 are 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 [36]. Newman and Armen [46][47][48] and Ohji et al. [49] 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 [45], and showed that crack opening stresses were a function of R ratio (S min /S max ) and the stress level (S max / 0 ), where  0 is theflow stress i.e: the average between σ ys and σ u . Blom and Holm [50] and Fleck and Newman [51][52] 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. [53] found later that the residual plastic deformations that cause closure came from the crack. McClung [54][55][56] performed extensive finite element crack closure calculations on small cracks at holes, and various fatigue crack growth models. Newman [57] 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 [56] found that K analogy, using Kmax/K0 could correlate the crack opening stresses for different crack configurations for small scale yielding conditions where K 0 = 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.)

Crack opening stresses
The applied stress level at which the crack surfaces are fully open (no surfaces contact), denoted as S 0 , was calculated from the contact stresses at Smin. To have no surface contact, the stress-intensity factor due to applied stress increment (S0 -S min ) is set equal to the stress intensity factor due to the contact stresses. Solving for S 0 gives : and c 0 is the current crack length minus c*, where c* is the crack growth increment over which S 0 is held constant. The analytical closure model provides no information about the amount of crack growth per cycle. Crack growth is simulated by extending the crack an incremental value and the moment of the maximum applied stress. The amount of crack extension (c*), was arbitrarily defined as c*,=0.05max, where max is the plastic zone size caused by the maximum applied stress occurring during the c* growth increment. The increment in width of element n, and its significance is discussed in the next section. If j = 0 for j = 11 to n -1 at the minimum applied stress, then the crack is already open, and S 0 cannot be determined from Eq. (14). The stress  j at the crack tip changes from compression to tension when the applied stress level reaches S0, [56]. The crack growth equation proposed by Elber [45] states that the crack growth rate is a power function of the effective stress intensity factor range only. Later, Hardraht et al. [58] showed that the power law was inadequate at high growth rates approaching to fracture. The results presented herein show that it is also inadequate at low growth rates approaching threshold. To account for these effects, the power law was modified to: where: The constants C1 to C5 are determined by experimental test under constant amplitude loading. The factor F is the boundary correction factor on stress intensity. The analytical closure model provides no information about the amount of crack growth per cycle. Extending the crack an incremental value at the maximum applied stress simulates crack growth.

Modified Dugdale model (crack closing)
Various modified Dugdale models were proposed [59][60]. After Elber [45] defined crack closure, the research community developed analytical or numerical models to simulate fatigue crack growth and closure. These models were designed to calculate the growth and closure behavior instead of assuming such behavior as in the empirical models. Seeger [61] and Newman [48]

CRACK PROPAGATION BY FINITE ELEMENT METHOD
lber´s [45] experiments of crack closure with constant amplitude loading proposed the following equation for fatigue crack propagation rates: where C and n are constants of the material and Keff is the effective stress intensity factor range that can be calculated by: where: c -half length of the crack, F -boundary correction factor S eff -effective stress range Fig. 7 shows the center crack panel that was used to evaluate the crack propagation. Fig. 8 shows the panel idealized to finite element method. The panel material was assumed to be elastic perfect plastic with a tensile (and compressive) yield stress,  0 , of 350.0 MPa and modulus of elasticity of 70000 MPa. These properties are for an aluminum alloy. The released nodes range from node A to node F. The accuracy of the calculated crack opening stresses would be affected by the mesh size chosen to model the crack tip region. A finer element mesh size would give more accurate results. Newman [29] evaluated three types of mesh, as shows in Tab. 3. E Figure 7: Center Crack Panel, Newman [29]. Figure 8:     Fig. 9 displays how is obtained the stabilization of opening stress during crack propagation. The best agreement with experimental results was mesh II, although mesh III provides also good results. Nowadays, with the computer facilities, the time to evaluate this mesh is reduced, and the size of the element and the increment are mostly used today to evaluate crack closure analysis or propagation. Matos & Nowell [66] present a literature review of the phenomenon of plasticityinduced fatigue crack closure under plane strain conditions and mention that there are controversial topics concerning the mechanics of crack propagation. In general there is no consensus in the scientific community. Fleck [67] used finite elements to simulate plasticity induced crack closure under plane strain conditions and predicted that the nature of the closure process changes from continuous to discontinuous after a sufficient increment of crack growth. He suggested that closure involves only a few elements relatively distant from the current crack tip and the closure levels decay steadily as the crack grows beyond its initial length. In the limit, the closure would not occur at all. In Singh et al. [67] the authors provide a review of some crack propagation issues. The paper covers the transients and single overload effects as well as the plasticity induced crack closure. In this topic Singh et al [68] presented a discussion regarding how the researchers normally work in crack propagation simulation considering overload-induced crack closure. Lei [69] determine the crack closure by finite element method in a compact specimen. In the work Lei [69] use ABAQUS [70] to perform the crack propagation simulation using the crack face method was good agreement with experimental data.

DESCRIPTION OF MODEL
compact tension specimen was modeled using a finite element code, MSC/Patran, r1 [71] and ABAQUS Version 68 [70] 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.10, and the boundary conditions are presented in Fig. 11. The finite element models 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. 10 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 Kmax 0.6 K IC , using eq. (4.1), where K IC is the critical stress intensity factor of adopted material in the present study. With the value of K max from K IC computed as mentioned above is computed the maximum load using Eq. 3.1 to be applied in the specimens as explained in next.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. 10.
The force is divided in nine steps between loads Pmin-Pmax and nine steps between the Pmax-Pmin,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. (4.2). Only the first 20 reverses from load history shown in Fig. 7 were used to identify crack opening/closing and retardation effects.
where: K max = 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 [72] 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 [73], at minimum load [74,75] or at other positions of the load cycle. Ogura et al. [49] 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 [76,77] 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 [78,79]. The stress level in the crack tip, Fig.12, must to be positive to characterize the crack opening and negative to characterize the crack closure. Antunes & Rodrigues [72] 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 [45] and has been widely used by Jiang et al. [73]. 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 [80]. 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 [47] recommendation with four elements yielded in the reverse plastic zone. In variable amplitude loading the crack length can not 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. 13 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 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. 14 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). Fig. 15 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), I

Surface Crack
Crack Tip  y= -85 MPa 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. 16 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.