Numerical Simulation and Experimental Studies on Elastic-Plastic Fatigue Crack Growth

Jie Wang,Wei Jiang, and Qi Wang

Abstract: A elastic-plastic fatigue crack growth (FCG) finite element model was developed for predicting crack growth rate under cyclic load.The propagation criterion for this model was established based on plastically dissipated energy.The crack growth simulation under cyclic computation was implemented through the ABAQUS scripting interface.The predictions of this model are in good agreement with the results of crack propagation experiment of compact tension specimen made of 304 stainless steel.Based on the proposed model,the single peak overload retardation effect of elastic-plastic fatigue crack was analyzed.The results shows that the single peak overload will reduce the accumulation rate of plastic energy dissipation of elements at crack tip plastic zone,so that crack growth will be arrested.The crack growth rate will not recover until the crack tip exceed the affected region.Meanwhile,the crack growth rate is mainly determined by the amplitude rather than the mean load under the condition of small scale yielding.The proposed model would be helpful for predicting the growth rate of mode I elastic-plastic fatigue crack.

Keywords:Plastic dissipation energy,elastic-plastic fatigue crack,finite element model,crack growth rate.

1 Introduction

Fatigue failure of components usually involves the initiation and propagation of cracks.Accurate prediction of crack growth rate has a great significance for evaluating the service life of equipment and components.Currently,there are many approaches which base on corresponding failure criteria has been proposed to predict fatigue crack growth rate,including damage mechanics,stress intensity factors and energy criteria [Meneghetti and Ricotta (2016);Noroozi,Glinka and Lambert (2007);Kujawski (2001);Paris (1963);Huffman (2016)].The formulas based on stress intensity factors have been widely used for calculating fatigue crack growth rate,which are mainly derived based on the theory of linear elastic fracture mechanics or empirical formula obtained by fitting experimental results.The application scope of these formulas is limited.Especially for the metal materials with high toughness,the calculation accuracy is relatively low [Huseyin,Gall and Garcia (1996)].However,with the continuous deepening and development of relevant researches,plastic deformation near the crack tip usually occurs during crack propagation.The existence of plastic zone will introduce residual compressive stress near the tip which can rationally explain the single peak overload retardation effect of fatigue cracks.Many scholars have carried out a lot of experimental and theoretical studies on the plastic zone near the crack tip.Gao et al.[Gao,Wang,Kang et al.(2010)] presented an analytical solution to calculate the crack tip plastic zone in plane stress states based on stress field equations and Hill yield criterion.Besel et al.[Besel and Breitbarth (2016)] studied the types of plastic zones at crack tip of a commercially available aluminium alloy under cyclic loading and found that there are five kinds of plastic zones (PZ) including primary plastic zone (PPZ),persistent part of primary plastic zone (PPPZ),cyclic plastic zone (CPZ),backward cyclic plastic zone (BCPZ) and forward cyclic plastic zone (FCPZ).Rice [Rice (1967)] pointed out that the plastic deformation of metallic materials is mainly due to dislocation movement caused by fatigue.The accumulation of plastic strain can directly lead to the increase of plastically dissipated energy.Microplastic strain hysteresis energy is considered to be an index for fatigue damage.The fracture failure would occur when the plastic dissipation energy of the material reaches a critical value [Feltenr and Morrow (1961)].Therefore,the propagation of cracks can be regarded as the failure of materials near the crack tip with the accumulation of plastic dissipation energy under cyclic loading.The application of critical plastic dissipation energy criterion for fatigue crack extension was first suggested by Rice [Rice (1967)],and large number of numerical and experimental studies were followed.Klingbeil [Klingbeil (2003)] proposed a new theory for predicting fatigue crack growth based on the total plastic energy dissipation per circle ahead of the crack.Dang [Dang (1993)] corroborated that dissipated energy criteria were versatile both at the microscopic and macroscopic analysis of fatigue through a macro-micro approach in high-cycle multiaxial fatigue.The experimental carried out by Ranganathan et al.[Ranganathan,Chalon and Meo (2008)] also show that the plastic dissipation energy can be used to determine crack propagation rates under both constant amplitude and variable amplitude cyclic loading.Pandey [Pandey (2003)] developed a FCG model under constant amplitude loading by considering energy balance during growth of the crack.Zuo et al.[Zuo,Kermanidis and Pantelakis (2002)] proposed a strain energy density crack growth model which can predict the lifetime of fatigue crack growth.

Numerical methods have been increasingly used to solve complex engineering problems since the development of computers and numerical method,especially,since the significant growth of computing power in the last decades [Ding,Karlsson and Santare (2017);Breitbarth and Besel (2017);Camas,Lopez,Gonzalez et al.(2017);Paul (2016);Paul and Tarafder (2013)].The plastic zone at crack tip is relatively small.The numerical model with refining grid can be used to simulate the change rule of crack tip plastic zone under periodic load,and the plastic dissipation energy of the elements can also be determined easily.Nittur et al.[Nittur,Karlsson and Carlsson (2013)] investigated the crack growth rate during cyclic loading via numerical simulation,and the crack extension is governed by a propagation criterion that relates the increment in plastically dissipated energy ahead of the crack tip to a critical value.Ding et al.[Ding,Karlsson and Santare (2017)] simulated the paris-regime fatigue crack growth in polymers via a numerical procedure.The propagation is calculated based on the assumption that crack extension is controlled by the plastically dissipated energy in the plastic zone around crack tip.Compared with the approach based on stress,strain or displacement,the prediction method of crack growth rate based on plastic dissipation energy is attractive for the following two reasons:(1) even with classical sharp crack modeling,the plastic strain energy is nonsingular,which makes for straightforward interpretation of numerical results [Klingbeil (2003)];(2) predicting the crack growth rate based on the total dissipated energy is in good agreement with broad trends in crack growth data for a variety of ductile metals,such as several aluminum,low carbon steel,titanium and nickel based alloys used in aerospace applications.

In this study,a two-dimensional finite element model of precracked compact tension (CT) specimen was established to simulate crack growth under cyclic computation.Meanwhile,the critical plastic dissipation energy has been developed as the failure criterion of elements to study the crack growth rate under cyclic load,including constant amplitude loads,single peak overload,and random loads.The two-dimensional plane strain finite element model is suitable for simulating stress state of thick plate,which can also save computing cost.It was found that this model can predict the crack growth rate accurately by comparing the calculation results and experimental results.Based on the finite element model proposed in this paper,the distribution of plastic dissipation energy in the vicinity of crack tip under cyclic loading and crack propagation rate under different loads were studied.

2 Modeling framework

2.1 Material and determination of critical plastic dissipative energy

The plastic dissipation energy is caused by plastic deformation and the accumulation law of which is determined by the constitutive relation of materials under cyclic loading.The 304 stainless steel (304SS) has good mechanical properties and toughness,which has been widely used in many fields of manufacturing industry.This paper will take 304SS as the research object and the mechanical properties of this material are shown in Tab.1 [Syed,Jiang,Wang et al.(2015)].The cyclic stress-strain curve of 304SS is shown in Fig.1.

Table 1:Mechanical properties of 304SS

According to Ramberg-Osgood [Nies?ony,Dsoki and Kaufmann (2008)] cyclic stress strain curve,the correlation between stress and strain under uniaxial load can be described with the following expression.

wherek' is the cyclic strength coefficient;n' is the cyclic strain hardening exponent;εandσare the strain and stress;Eis the modulus of elasticity.

The correlation between strain amplitude and cycle numbers can be expressed by the Manson-Coffin-Basquin strain-life fatigue curve [Ricotta (2015)],as shown below.

where Δε2 is the strain amplitudes;2Nfis the number of reversals to strain failure;is the fatigue strength coefficient,calculated as the value of the stress amplitude for 2Nf= 1;bis the fatigue strength exponent;ε'fis the fatigue ductility coefficient,defined as the value of the plastic strain component for 2Nf= 1;cis the fatigue ductility exponent andE' is the material dynamic elastic modulus,which is usually considered to be equal to elastic modulusE.

The Ramberg-Osgood stress-strain equation of 304SS can be obtained by least squares fitting method,the calculation process is shown as follows.

The Eq.(1) can be rearranged in the form ofy Ax B= + is shown as follows [Liu,Zhou,Xia et al.(2016)]:

The fatigue life of the material under a strain amplitudeaεcan be obtained by conducting fatigue test,the total strain amplitudeaεis consists of elastic strain amplitudeeεand plastic strain amplitudepε.The fatigue life under corresponding elastic strain amplitude and plastic strain amplitude can be obtained according to the strain-life data of stainless steel 304 [Wang,Kan and Zhang (2003)] and Eq.(2).The undetermined coefficient of Eq.(6) and Eq.(7) isandc,which can be got by linear regression.Then,the Manson-Coffin-Basquin strain-life fatigue equation can be determined ultimately.The calculation process can be completed through the Curve Fitting tool in the commercial software Matlab.

Figure 2:Fitting curve of strain-life

Pandey [Pandey (2003)] presents the concept that energycWabsorbed till fracture is the area below the cyclic stress strain curve.The formula for calculating critical plastic energy is shown as follows.

wherecWis critical plastic dissipation energy;is fatigue strength coefficient;is fatigue ductility coefficient;k' is cyclic strength coefficient;is cyclic strain hardening exponent.

2.2 Calculation of crack propagation rate

According to the Neuber micro-support concept,the engineering metallic materials can be modeled as a medium made of element blocks of dimensionδ[Mikheevskiy,Bogdanov and Glinka (2015)].The propagation of cracks can be considered as the fracture failure of elements at the front of crack tip.As shown in Fig.3a,there are eight elements arranged in the front of the crack tip.As they are symmetrical about ligament,the plastic dissipative energy of corresponding elements are equal because of their same plastic deformation,such as element 1 and element 5.For mode I crack,it can be assumed that the crack propagates along the ligament.The plastic deformation of the elements in plastic zone would be produced under cyclic loads.Therefore,the plastic dissipation energy in the elements will continue to accumulate with the increase of load cycles.The numerical results show that the plastic dissipation energy of element 1 and element 5 accumulate rapidly and much higher than that of other elements.The plastic dissipation energy of these two elements will reach the critical value first.When the plastic dissipation energy of element 1 and element 5 reaches the critical value,crack propagation will be accomplished by releasing the constraints on the current crack front nodes at the end of the load cycle,i.e.at minimum load,as shown in Fig.3b.The increment of crack propagation is increased with the binding of nodes,then,the current crack growth rate can be expressed as Eq.(10).For standard CT specimen,the 1/2 model can be adapted for numerical simulation,and symmetry boundary is the ligament of crack.

Figure 3:Schematic of crack propagation

where dadNrepresent fatigue crack propagation rate,Ncis the cycles required for plastic dissipation energy of a element reach to the critical value.

3 Parameterized finite element modeling and calculation process description

In this paper,the standard CT specimen has been taken as the research object,and dimensions are shown in Fig.4.The specimen with a 2 mm precrack is made by wire cutting machine and the total lengthaof the initial crack is 18 mm.The numerical simulation is implemented in the commercially available finite element simulation package ABAQUS and the element type is CPE4.The linear kinematic hardening model was adopted to consider the Basinger effect under cyclic loading.The vicinity of crack tip is refined by structural mesh with the size ofδ.A reference point is set in the center of pin hole,and the kinematic type coupling constraint is established between the reference point and the inner wall of pin hole.The load in the direction of the X-axis is applied on the reference point.

Figure 4:Dimensions and finite element model of CT specimen

In order to achieve computational scheme for cyclic crack propagation,the parameterized modeling was complete through ABAQUS script file.The script file of ABAQUS is written in an object-oriented programming language based on Python.The contents of script file include creation of model,setting the parameters of material and element type,specifications of loads and boundary conditions,finite element mesh generation and so on.Then the Job file was automatically submitted to the ABAQUS to complete the calculation.Once the computation was completed,the plastic dissipative energy of specific element was extracted from the database file.Then compare the values with the critical dissipation energy to judge whether the crack can growth,otherwise the number of loads would be further increased.If the extension criterion was satisfied,the crack would propagate a length of dawhich is equal toδ.When the crack extends to a specified lengthca,the program would automatically exit.The data transfer function of ABAQUS can achieve the goals of setting previous database as the predefined field of current model,which can effectively avoid repeated computation and save computing time.The specific computing process can refer to the flow chart of Fig.5.

Figure 5:Flow chart for prediction of crack growth life based on ABAQUS

4 Crack propagation experiment of 304CT specimen

In order to verify the accuracy of numerical model established in this paper,the crack propagation experiment of CT specimen made of 304 stainless steel was carried out by MTS 810 material test system.The CT specimens with 2 mm prefabricated cracks were prepared in accordance with ASTM E647-12 standards and the dimensions of which is shown in Fig.4.As the precrack was created by a wire electric discharge machining system,the irregular crack front may affect the generation of initial cracks.Before the experiment was started,the fatigue crack with a length about 0.2 mm was prefabricated with low amplitude load.Then the distance between the newly generated crack tip and load action line (The connection between two centers of pin holes) is taken as the initial length of the crack.The sinusoidal constant amplitude periodic load was applied in the fatigue test,the frequency is 20 Hz,the load ratio is 0.1,and the maximum load is 8 kN.Three samples were used in this experiment and the basic parameters are shown in Tab.2.

Table 2:Basic parameters of specimen

Figure 6:Experimental equipment and specimen.a.Material test system:1.MTS 810 material test system;2.CT specimen;3.measuring microscope.b.The photograph of crack tip prefabricated by wire cutting machine

In general,the Paris formula is wildly used to evaluate the crack growth rate under the condition of small scale yielding [Mikheevskiy,Bogdanov and Glinka (2015);Toshihisa (2005)].Technical manuals [Zhao (2015)] proposed a formula based on Paris' law for calculating the crack growth rate of 304SS,in which the material constants arec=7.45 ×10-14andm= 3.1,as shown below.

The range of stress intensity factor in above formular can refer to the empirical formula ofKΔ presented in ASTM E647-11 standards,as shown below.

in whichα=a W,ΔPis load range,Bis the thickness of CT specimen.

The experimental results are compared with the results of numerical simulation and theoretical formula,as shown in Fig.7.The results show that the simulation results are in good agreement with the results of experiment and theoretical formula.However,the crack growth rate evaluated by theoretical formula gradually exceeds the value with the increase of crack length.The reason maybe that the plastic region in the vicinity of crack tip is increased with the increase of crack length and the Paris formula is established on the basis of linear elastic fracture mechanics theory.Therefore,the accuracy of the theoretical formula will be reduced with the increase of the plastic zone at the crack tip.This result also shows the superiority of the finite element model in describing the propagation of elastic-plastic fatigue cracks.The curve of simulation results also show that crack growth rate of CT specimen under constant amplitude cyclic loads increased with the extension of crack.The main reason contributing to this phenomenon is the increase of both stress intensity factor and stress intensity factor amplitude with the increase of crack length.The load cycles required for the same length of crack increment is reduced.It is proved that this numerical model is applicable for predicting fatigue crack growth rate.

Table 3:Experimental results of crack propagation

Figure 7:Experimental results and numerical simulation results of crack propagation

5 Study on the law of crack propagation

5.1 The retardation effect of single overloading

Numerous studies show that the single peak overload will restrain the crack growth rate.In this section,the retardation effect will be explained from the change of plastic dissipation energy of elements ahead of crack tip after the single peak overloading.The load ratio of constant amplitude is 0.1,and the maximum load is 4 kN.As shown in Fig.8,the plastic zone near the crack tip usually presents butterfly shape (Only 1/2 model is shown).The elements located on the surface of crack ligament are numbered from 1 to 6,which is located in the plastic zone.The simulation results in Fig.8 (detailed data are listed in Tab.4) show that the plastic dissipation energy of elements mentioned above increased with the increase of load cycles.However,the plastic dissipation energy of element 1 adjacent to the crack tip increases fastest and its value is much higher than that of other elements.The simulation results also show that the area of plastic zone will increase with the increase of load cycles.There is no detail for such phenomenon,because the size of plastic zone is growing slowly.The plastic dissipation energy of these elements in newly generated plastic zone is relatively low,which can be ignored.The growth rate of plastic dissipation energy in element 1 determines the crack growth rate under current computing step.Of course,the accumulation of plastic dissipative energy in other elements also affects subsequent crack propagation.However,the basic values and amplitude are relatively low when compared with critical plastic dissipation energy,which has a little effect on the crack growth rate.The overloading is set at 16th load step,which is 1.25 times of the maximum load of constant amplitude cyclic loads.The load spectrum is shown in Fig.9.The results show that the plastic dissipation energy of both element 1 and element 2 suddenly increases after overloading.If the overload is too large and make the plastic dissipation energy of element 1 exceeds the critical value,and the crack growth rate will increase slightly.But due to the effects of overloading,the average accumulative rate of strain dissipative energy was reduced to some degree.This trend shown in Fig.10 is not evident because the total number of cycles is too small.Actual computation shows that the number of cycles required for the plastic dissipation of element to reach the critical value is about 7000 under current load condition.Therefore,the number of cycles required for crack to extend an incremental length will increase.It also means that the growth rate of crack is reduced.The plastic dissipation energy of element 2 increased a lot after single overloading,but the magnitude is relative low.The accumulative rate of the plastic dissipative energy is almost fell to 0.

Figure 8:The plastic dissipative energy of elements varies with the number of cycles

Figure 9:Load spectrum

Figure 10:The effect of single peak overload on the plastic dissipative energy of elements.The left ordinate corresponds to element 1,and the right ordinate corresponds to element 2

Figure 11:Effect of single peak overload on crack propagation

Figure 12:Effect of single peak overload on crack growth rate

When the single peak overloading occurs,the length of crack are 20 mm.Fig.11 shows that the number of cycles required to extend the same length of crack are increased due to the retardation effect of single peak overload.The result in Fig.12 also shows that the crack growth rate is reduced after single overloading.With the crack continues to extend,the rate gradually returns to normal level.In order to find out the fundamental cause of this phenomenon,the stress of four nodes were extracted.The stress of four nodes in the direction of X-axis has been extracted for research (shown in Fig.14a and Fig.14b).The node 1 is directly located at the tip of the crack.The node 2 is located about 0.2 mm ahead of the crack tip.Node 3 is located at the edge of plastic zone,which is about 0.4 mm ahead of the crack tip.Node 4 is located well outside the plastic zone with a distance of 0.7 mm to the crack tip.The location of these four points is shown in Fig.13.First,we can find that the plastic zone increases significantly after overloading.The single peak overloading has different effects on the stress of the four nodes.

For node 1,during the first load step the local stress reaches about 600 MPa.During the first unloading phase the local stress values change to -390 MPa.During the next fifteen load cycles the local stress under maximum and minimum load change to 580 MPa and -480 MPa respectively.However,the maximum stress and minimum stress reaches to 660 MPa and -520 MPa after overloading.Consequently,the effective local stress range under cyclic loading changes to about 1100 MPa which is essentially equal to the value without overloading.But the mean stress is relatively low.

For node 2,during the first load step the local stress reaches about 480 MPa.During the first unloading phase the local stress values change to 40 MPa.During the next fifteen load cycles this local value under maximum and minimum loading has a little change and reaches to 470 MPa and 50 MPa,respectively.However,the overloading has a great influence on the local stress,and the maximum stress and minimum stress reaches to 300 MPa and -80 MPa after overloading.

Different from the other three nodes,the overloading has a little effect on the local stress of node 4.The maximum stress and minimum stress of node 4 nearly keep a constant of 340 MPa and 80 MPa,respectively.Due to single overloading,the maximum stress change to 440 MPa.However,the value returns to normal level during next load cycle.The stress range is basically maintained at 520 MPa.

Comparing with node 1 and node 2,the average stress of which is decreased by overloading.The average stress of node 4 has a small increase after overloading.The maximum stress increased from 240 MPa to 340 MPa after overloading and then back to 270 MPa during next following load cycle.

As we can see from the variation of local stress of four key nodes mentioned above that large compressive residual stress in plastic zone ahead of crack tip has been induced by the overloading.Meanwhile,a relative small residual tensile stress also has been induced in some part of area ahead of plastic zone.As the compressive residual stress would decrease plastic deformation of the elements,thus the crack growth rate will be reduced caused by the decrease of accumulation rate of plastic dissipation energy.The domains influenced by single peak overloading are limited.When the increment of crack exceeds the affected region,the crack growth rate will return to normal level.

Figure 13:Plastic zone and location of four key nodes

Figure 14:The cyclic stress ranges of nodes;a.node 1 and node 2;b.node 3 and node 4

5.2 Effect of Load amplitude and mean load on crack propagation

In order to study the effect of load amplitude and mean load on the crack growth rate,the crack growth rate under different load ratios were simulated based on above numerical model under the condition of the mean load and load amplitude are constant,respectively.When the effect of load amplitude on crack growth rate was studied,the mean load is 4.4 kN,the range of load ratio is 0.05 to 0.4,the increment of crack propagation is 5.8 mm,and the results are shown in the Fig.15.When studying the effect of mean load on crack growth rate,the amplitude of the load is 7.2 kN,the range of load ratio is 0.05 to 0.4,the increment of crack propagation is 5.8 mm,and the results are shown in the Fig.16.

As can be seen from the results shown in Fig.15,the crack propagation rate seems to increase exponentially with the load amplitude.Meanwhile,the results of Fig.16 show that the mean load has little effect on the crack growth rate.Therefore,crack growth rate are mostly affected by load amplitude.That is why the crack growth rate based on linear elastic fracture mechanics is expressed by stress intensity factor amplitude [Paris (1963);Noroozi,Glinka and Lambert (2005)].With the discovery of crack closure effect,more and more crack growth rate expression based on effective stress intensity factors are proposed [Antunes,Correia and Ramalho (2015);Jiang,Feng and Ding (2005);Fukumura,Suzuki and Hamada (2015)].According to above numerical model,it also can be speculated that the plastic deformation in the vicinity of crack tip would not be induced by a certain low external load.The crack will stop propagating,as the plastic dissipation energy cannot be accumulated.The specific effect of crack closure needs to be further explored.

Figure 15:The effect of load amplitude on crack growth rate

Figure 16:The effect of mean load on crack growth rate

6 Conclusion

In this work,a numerical model for predicting crack growth rate under cyclic loads has been established based on plastic dissipative energy.The critical plastic dissipation energy is taken as the propagation criterion of elements,and the crack growth simulation under cyclic computation is implemented through the scripting interface of commercial finite element software ABAQUS.The critical plastic dissipation energy is mainly determined by the constitutive relation of metallic materials under uniaxial cyclic loading.The crack propagation experiment of 304SS CT specimen under cyclic load has been carried out through material testing system.The predictions of this model are in good agreement with the results of experiment.The retardation effect of single overloading and the influence of mean load and load amplitude on crack growth rate are studied based on the numerical model established in this paper.The analysis results show that crack propagation rate is mainly affected by the load amplitude,while the effect of mean load on crack propagation rate is relatively small.The residual compressive stress induced by single peak overloading will reduce the accumulation rate of plastic dissipation,thus restraining the crack growth to a certain extent.This model will help to further study the propagation of elastic-plastic fatigue crack.

Acknowledgment:Project No.51575076 supported by the National Natural Science Foundation of China.