- Email: [email protected]
A new finite element fatigue modeling approach for life scatter in tensile steel specimens
International Journal of Fatigue 32 (2010) 685–697
Contents lists available at ScienceDirect
International Journal of Fatigue journal homepage: www.elsevier.com/locate/ijfatigue
A new finite element fatigue modeling approach for life scatter in tensile steel specimens Anurag Warhadpande, Behrooz Jalalahmadi, Trevor Slack, Farshid Sadeghi * School of Mechanical Engineering, Rm #303, Purdue University, West Lafayette, IN 47907, United States
a r t i c l e
i n f o
Article history: Received 15 April 2009 Received in revised form 7 September 2009 Accepted 7 October 2009 Available online 7 November 2009 Keywords: Damage mechanics Fatigue Voronoi finite element fatigue modeling
a b s t r a c t This paper presents the results of a finite element model developed to investigate the effects of material microstructure topology on fatigue damage evolution in dog-bone tensile specimens. The Voronoi finite element model (VFEM) developed was used to study fatigue damage evolution in AISI 4142 steel under various loading conditions. An environmental scanning electron microscope (ESEM) coupled with an Instron fatigue testing machine was used to determine the strain and fatigue life of the dog-bone 4142 steel specimens at different mean stress levels. The damage variable was obtained using the method of variation of elasticity modulus. The experimental damage versus cycle curve was then used to determine the material properties (i.e. resistance stress rR and exponent m) needed for the VFEM fatigue damage model. The VFEM model was used to corroborate the analytical and experimental fatigue results. A comparison of the results obtained from the VFEM model and ESEM indicate that they are in good agreement. Ó 2009 Elsevier Ltd. All rights reserved.
1. Introduction Fatigue is primarily a process of damage accumulation which is manifested by progressive internal deterioration of material properties due to nucleation and growth of micro-cracks, debonding, voids, etc. The mechanism of fatigue failure in general consists of three stages: (i) crack initiation, (ii) crack propagation and (iii) final catastrophic failure. The percentage of the total life spent in each stage is related to material properties and loading conditions. Designing against fatigue failures is a time consuming and involves extensive experimentation to account for the statistical nature of the fatigue process. For high cycle fatigue much of the variability in fatigue life data can be accounted for in the crack initiation stage [1]. Crack initiation is a highly localized event which is strongly affected by the random distribution of material properties and defects that exist at the microstructural level [2,3]. Furthermore, engineering materials are polycrystalline, consisting of individual grains with random geometrical properties separated by grain boundaries that can act as physical discontinuities. Modeling approaches which include this topological and material property randomness can be used to explain fatigue scatter and to reduce the number of experiments required to quantify it. Damage mechanics approach [4] for modeling material deterioration and fatigue is increasing in popularity. Damage mechanics
* Corresponding author. Tel.: +1 765 494 5719; fax: +1 765 494 0539. E-mail addresses: [email protected] (A. Warhadpande), bjalalah@ purdue.edu (B. Jalalahmadi), [email protected] (T. Slack), [email protected] (F. Sadeghi). 0142-1123/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijfatigue.2009.10.003
introduces a damage variable, either scalar or tensorial, into the material constitutive model to account for the degradation in material properties that occurs as damage grows. The evolution of the damage variable is governed either by a phenomenological relationship or by a relationship developed based on micromechanical principals. Fatigue models based on damage mechanics concepts have been implemented in both continuum and discrete frameworks. Rinaldi et al. [5] proposed a lattice model to simulate polycrystalline microstructures and studied the resulting scatter in fatigue response at the macro level. In their model the fatigue behavior of each discrete spring was described using the Basquin law and damage accumulated according to a Palmgren–Miner law. By introducing a normal distribution of static strengths for the springs fatigue scatter was found to follow a lognormal distribution. Bolotin et al. [2] developed a model for the nucleation and early growth of fatigue cracks using the finite element method. In their model each grain was approximated by a square finite element and was randomly assigned a resistance stress from a Weibull distribution. Here resistance stress refers to a material parameter used in the fatigue damage evolution law. Their results exhibited many features such as a non-monotonic change in growth rate and scatter in the crack dimensions and growth rates that have been experimentally observed in the early propagation of fatigue cracks. Recently, Raje et al. [6,7] proposed a model for the fatigue of bearing contacts using the discrete element method. In their model the material was represented as an assemblage of discrete, rigid Voronoi polygons connected by normal and tangential springs. Damage accumulated according to a continuum damage model and was used to degrade the tangential spring stiffnesses.
686
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697
Nomenclature ~ A A b D Dcr e E E e Ec e Et F h L1 L2 m N R
effective cross-sectional area supporting the load overall cross-sectional area modification co-efficient damage variable critical value of damage variable damaged modulus of elasticity modulus of elasticity damaged elastic moduli in compression damaged elastic moduli in tension force acting on area A crack closure factor distance between the spots before application of load distance between the spots after application of load power exponent number of cycles minimum cyclic stress divided by maximum cyclic stress
By including randomness in the material topology and material properties they obtained Weibull slopes of the simulated life data and final spall crack patterns that compared well with experimental results. In this investigation a fatigue life model is developed using a Voronoi finite element method (VFEM) with a damage coupled constitutive relationship. The randomness of the material microstructure is simulated using a Voronoi tessellation process. The Voronoi cells are discretized using constant strain triangular finite elements and damage accumulates in each triangle according to a continuum damage model. Uniaxial in situ fatigue tests were conducted using an environmental scanning electron microscope (ESEM) equipped with an Instron fatigue testing machine to measure both damage as a function of loading cycle and life as a function of applied stress. The damage as a function of loading cycle data is used to identify the empirical constants used in the damage evolution law by the method of variation of elasticity modulus. The VFEM model was used to predict the stress versus life results which was then compared against experimental results. Initial imperfections in the form of material flaws and voids were then introduced into the model to investigate their effects on fatigue life scatter. Experimental and numerical VFEM results demonstrate good agreement over a wide range of loading conditions.
DD DN Dr
e m~ m
r~ r rmax rmean rR rR0 ru
increment in damage increment in number of cycles stress range axial strain damaged Poisson’s ratio Poisson’s ratio effective stress stress acting on area A peak value of stress mean value of stress resistance stress resistance stress corresponding to zero mean stress level ultimate strength
investigation. Care is taken that during the tessellation process; the element sizes are relatively uniform so that there is no significant variation in element size across the domain. This is achieved by setting an upper and lower bounds on the distance between the nucleation points. 2.2. Finite element modeling Previous investigations have in general used the classical finite element methods (for example using quadrilateral elements for the domain), which employs the homogenization techniques to simulate the material as an equivalent continuum model on the macroscale. However, using standard finite element software that contains a limited variety of elements having predefined shapes is not suitable because they cannot be used to properly simulate randomly shaped topology of material microstructure. Ghosh and coworkers [10–12] introduced the Voronoi Cell Finite Element Method (VCFEM) for deformation and stress analysis of domains with arbitrary microstructural distributions. Other investigators
2. Voronoi finite element method based fatigue modeling In the following section, the approach used to develop the Voronoi finite element model to investigate fatigue of dog-bone shaped steel specimens is presented. We will briefly describe the Voronoi tessellation process, the new damage mechanics based finite element model developed and proceed to describe the test rig used in this investigation. Then the experimental and analytical results of this study are presented with some concluding remarks. 2.1. Voronoi tessellation In the present study, each simulation domain is generated using the Voronoi tessellation method. A detailed description of the Voronoi tessellation process can be found in [8,9]. The set of points used to generate the Voronoi diagram are called the nucleation points. These nucleation points are distributed randomly before the start of the Voronoi tessellation process resulting in distinct simulation domains. The density of the nucleation points is decided in accordance with the grain size of the material under
Fig. 1. (a) Discretizing a domain with Voronoi elements. (b) Dividing the Voronoi cell into triangular elements. (c) Linear triangular finite element.
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697
[13–15] have also applied VCFEM to study various problems. Voronoi polygons generated for this investigation have random size and shape and are used to simulate material topological disorder and its corresponding effects on stress distribution and fatigue life. In this investigation, a new Voronoi finite element model [16] was developed to investigate fatigue in steel tensile specimens. For this purpose, the domain is divided into Voronoi cells using the Voronoi tessellation and then the geometrical center of each Voronoi element is determined. Because Voronoi elements are convex polygons, it is assured that the center of each element is located within the element. This property is used to divide the Voronoi element into finer triangular elements by connecting the center of the cell to each node. Fig. 1 illustrates the process of discretizing a domain into the Voronoi cells and then dividing the cell into finer triangular elements. A two-step assembly process is used to obtain the global stiffness matrix (K). In the first step, assuming linear shape functions for each triangle, its stiffness matrix is obtained and then the stiffness matrix for each Voronoi element was set up using the stiffness matrices for each individual triangle. In the second step, the global K matrix is assembled using the stiffness matrices of the Voronoi elements. 2.3. Damage variable Fatigue damage leads to progressive degradation of the material due to formation and growth of micro-cracks and voids. Using damage mechanics, these microscopic failure mechanisms can be treated in an empirical fashion by using the internal damage variable D. The damage variable D in general is a tensor, however, under isotropic conditions the damage variable D reduces to a scalar. It is also assumed that the Poisson’s ratio remains unaffected by damage [4]. The damage variable D accounts for the surface density of the micro-cracks and cavities in any plane of a representative volume element of the material. Kachanov [17] defined D as,
e AA D¼ A
ð1Þ
e is the effective area which carwhere A is the apparent area and A ries the load. Robotnov [18] introduced the concept of an effective stress to describe the effect of damage on strain behavior. The effec~ is defined as, tive stress r
F
r~ ¼ e ¼ A
F A r ¼ e A A 1D
ð2Þ
Applying the Hooke’s law Eq. (2) can be written in terms of the effective stress as,
e¼
r~ E
¼
r Eð1 DÞ
ð3Þ
687
where E is the modulus of elasticity and e is the elastic strain of the undamaged material. Defining e E as the elastic modulus of the damaged material,
e E ¼ Eð1 DÞ
ð4Þ
The reduction in elastic modulus with the increase in fatigue cycles is shown in Fig. 2. The damage variable is thus given by,
D¼1
e E E
ð5Þ
Eq. (5) can be used to indirectly measure the damage variable D in a process known as the method of variation of elasticity modulus. 2.4. Fatigue damage modeling Damage was incorporated into the VFEM model using the damage coupled constitutive relationship, Eq. (4), for each triangular element. In this analysis the evolution of the damage variable is described by the following model proposed by Xiao et al. [19] for high cycle fatigue
dD ¼ dN
Dr rR ð1 DÞ
m ð6Þ
Here, Dr is the stress range, N is the cycle number, and rR and m are material parameters. The parameter rR is sometimes referred to as the resistance stress [20] since it controls the ability of a material to resist damage accumulation. As the resistance stress increases the rate of damage accumulation decreases and consequently the fatigue life is longer. Xiao et al. [19] indicate that the resistance stress (rR) is a function of the mean stress and the relationship between rR and rmean is given by,
rR ¼ rR0 1 b
rmean ru
ð7Þ
where b is the modification co-efficient, rR0 is the resistance stress under completely reversed loading conditions and ru is the ultimate strength of the material. In order to capture the effects of micro-crack closure in compression, an additional term h is introduced in the constitutive model. The closure parameter h = 1 in tension and 0 6 h < 1 in compression. A value of h = 0 corresponds to complete micro-crack closure. In the present analysis, a value of h = 0.2 as proposed by Lemaitre [4] is used. The damaged coupled constitutive model becomes:
e E c ¼ Eð1 hDÞ e E t ¼ Eð1 DÞ
ð8Þ
m~ ¼ m
Fig. 2. Stiffness degradation with increase in fatigue cycles.
where E is the initial elastic modulus under pristine conditions and e E t are the damaged elastic moduli in compression and tenE c and e sion respectively. Fatigue modeling involves the simultaneous solution of Eqs. (6)–(8) to account for the coupling between damage and the material constitutive equations for each triangular element. Thus, the stress–strain relationships need to be solved and the triangular element elastic moduli need to be updated after every stress cycle to incorporate the current state of damage in the triangular element. However, it is impractical to perform the large number of increments needed for the millions of stress cycles arising in high cycle fatigue. Therefore, a method that involves the use of the ‘jump-incycles’ procedure [4] was used for achieving computational efficiency.
688
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697
The method assumes a piecewise periodic constant amplitude loading, i.e., the stress histories for the triangular elements are assumed to remain unaltered over a finite number of cycles DN constituting a block i. The increment in damage DD is assumed to be constant over a block of cycles. Thus, the damage evolution is assumed to be piecewise linear with respect to the number of cycles. It is to be noted that the shape of the damage evolution curve is not predetermined but it is an outcome of the numerical simulations due to the stress–damage coupling. The value of DD is a tradeoff between accuracy and computation time. For higher resolution, DD needs to be chosen as small as possible, however it directly affects the computation time. Stress calculations are performed once for every block of cycles and not for each individual cycle N. This saves a considerable amount of computation time and effort since the major part of the computational expense is spent in solving for the stress fields in the domain. The number of cycles in a block is variable and it is obtained using the following procedure: (i) Initial damage states for each triangular element are assigned. For a pristine material domain containing no initial flaws, the initial damage in each triangular element is assumed to be 0. Thus,
D0j ¼ 0;
j ¼ 1 . . . ntriangular elements
ð9Þ
The number of cycles elapsed is initialized to 0. Thus,
N¼0
ð10Þ i min Þj
(ii) Stress histories ðrmax r for each triangular element are computed for the current block i. (iii) The damage evolution rate in each triangular element is evaluated knowing the stress history for the present block and the current state of damage in the triangular element using,
dD dN
"
i ¼ j
ðDrÞij
#m
rR ð1 Dij Þ
ð11Þ
(iv) The triangular element with the maximum damage evolution rate is chosen as the critical triangular element for the current block of cycles. Thus,
dD dN
i crit
dD i ¼ Max dN j
ð12Þ
(v) The number of cycles in the current block of cycles is computed as,
DD DNi ¼ i dD
ð13Þ
dN crit
(vi) The number of cycles elapsed is updated to,
N ¼ N þ DN i
ð14Þ
(vii) The increment in damage for each triangular element during the current block of cycles is then given by,
DDij ¼
i dD DN i dN j
ð15Þ
(vii) The damage states for each triangular element at the start of the next block of cycles are updated to,
Djiþ1 ¼ Dij þ DDij
ð16Þ
(ix) The elastic moduli are modified at the start of the next block of cycles according to,
i ðEc Þij ¼ ðEÞ 1 hDj
ðEt Þij ¼ ðEÞ 1 Dij
ð17Þ
Steps (i) through (ix) are successively repeated for each block of cycles. A micro-crack is assumed to be initiated at a triangular element when the accumulated damage in the triangular element equals Dcr i.e. the critical value of the damage variable D. 3. Experimental approach to determine the material dependent constants The damage evolution law introduces three empirical constants
rR0 , m and b that have to be experimentally evaluated. The following sections describe the test apparatus, fatigue test specimen geometry and the experimental scheme developed to estimate the material constants. 3.1. Experimental setup The fatigue experiments were conducted using an ESEM equipped with an Instron load frame [21]. The Instron load frame has a servo-hydraulic actuator and controller which can be used to induce the load cycles in the test specimen. Fig. 3 illustrates the ESEM and Instron load frame with the test specimen installed between the grippers. One end of the specimen is firmly held in the fixed jaws while the other end is actuated using the servo-hydraulic system. Fig. 4 depicts the load frame installed inside the ESEM which allows for in situ monitoring of the test specimen during fatigue testing. The tests were conducted with water vapor as the gaseous medium inside the ESEM vacuum chamber at a partial pressure of 2–4 Torr and at ambient room temperature. 3.2. Specimen geometry The specimens used in this investigation were made of AISI 4142 steel. Fig. 5 depicts the two different test specimen geometries evaluated for fatigue testing. One of the test specimens had a uniform gauge cross section, while the other had a non-uniform gauge cross-sectional area with the minimum cross section occurring at the mid-span. When the test specimen with uniform gauge cross-sectional area (Fig. 5) was subjected to cyclic loading, failure occurred at various locations with the majority occurring near the shoulder of the fillet radius. This was attributed to manufacturing errors associated with the fillet radius; and the effects of clamping leading to fracture as shown in Fig. 5. It is to be noted that the ESEM and load frame are designed such that the electron detector can focus only in the central region of specimen gauge section. Hence the specimen geometry was modified to one with a nonuniform cross section leading to a region of maximum stress concentration occurring at the central section of the specimen. In order to obtain a continuously decreasing cross-sectional area an elliptical profile was machined on the specimen with the major and minor axis of ellipse as indicated in Fig. 6. With this specimen geometry the electron detector could be focused on the region of maximum stress concentration where the failure occurred. During the course of the fatigue testing, a crack typically appeared at the central region of the specimen where the cross-sectional area is minimum. The crack initiated at one of the edges and propagated along the thickness and width of the specimen as shown in Fig. 7. Once the crack extended through the thickness of the specimen, instantaneous fracture occurred.
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697
689
Fig. 3. ESEM and Instron fatigue test machine.
Fig. 4. Instron load frame installed inside the ESEM vacuum chamber.
3.3. Strain measurement In this section, a brief description of the strain measurement technique used in this study is presented. The digital images obtained from ESEM during in situ fatigue testing were used to determine the strain data [22–24] from which the damage variable is measured. The method adopted here utilizes features present on the surface of the test specimen (such as microstructure and material inhomogeneity) which provide useful contrast. This approach in strain measurement requires monitoring of the distance between two dark spots on the grayscale image. The change in length divided by the original length between the two spots provides the strain under the applied load. While performing the in situ fatigue testing, the electron detector is focused along the specimen thickness near the minimum cross-sectional area of the specimen gauge section. The region is scanned and two spots along the loading direction are selected.
The spots are chosen such that they span a large-enough number of pixels so that their shapes and positions are well defined and the area of interest is fully covered [25]. Fig. 8 depicts ESEM micrographs of two spots monitored during the loading and unloading of a specimen. While generating these images, the contrast and brightness were kept constant. Digital images were produced at a resolution of 2048 2048 pixels. From the ESEM images, grayscale intensity plots were generated along the two selected spots as shown in Fig. 8. The two deep valleys in the figures correspond to the two dark spots on the test specimen. Both valleys are associated with particular pixel numbers. The difference in pixels between the valleys is converted to the distance between the dark spots using a scaling factor. The scaling factor is a function of image magnification and provides the relationship between pixel distances on the digital image and the actual distance on the specimen. It is obtained from the scale on the digital image. Once the distances are obtained the strain can be calculated using,
690
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697
Fig. 5. Specimen geometry.
3.4. Identification of material dependent constants The empirical constants rR0 , m and b introduced by the damage evolution law are evaluated as follows. Assuming that the stress range is constant in the damage evolution equation, Eq. (6), is integrated to obtain the following relationship between the damage variable D and fatigue cycle N,
( D¼
Fig. 6. Drawing of the dog-bone shaped specimen (all dimensions are in mm).
e¼
L2 L1 L1
ð18Þ
where L1 and L2 are the distances between the dark spots before and after the application of load, respectively.
1 ) m mþ1 Dr 1 1 ðm þ 1Þ N
rR
ð19Þ
For a given fatigue test, the damage variable D was measured as a function of fatigue cycle N, generating an experimental D–N curve, using the following procedure. During a fatigue test, after DN cycles, the fatigue test was stopped. The specimen was loaded to a certain known value of stress and the corresponding strain induced in the specimen was measure using the method described in the previous section (Eq. (18)). The modulus of elasticity of the specimen at that instant of time was determined by dividing the applied value of stress by the measured value of strain. Since the test specimen has undergone fatigue cycling, the material softens
Fig. 7. Typical crack propagation pattern.
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697
691
Fig. 8. ESEM micrographs and grayscale intensity plots used to obtain strain data before and after load application.
and the modulus measured is the damaged modulus of elasticity e E (refer to Eq. (5)). Thus, using Eq. (5) the damage variable D was determined. In Eq. (5) E is the undamaged modulus of elasticity of the specimen measured at the start of each fatigue test. The fatigue test was continued for another DN cycles and the above mentioned procedure was repeated until the specimen fractured. Measuring D as a function of N during the course of fatigue test resulted in the D–N curve. Fig. 9 depicts the degradation in modulus of elasticity E and corresponding increase in damage variable D for a sample test (rmax = 300 MPa, rmean = 0 MPa, R = 1). Eq. (19) was used to curve fit the D–N data and obtain the empirical constants rR and m. Least square regression was used to fit the equation and the R-square value for the fit was 96%. It is to be noted that the value of the damage variable D is not equal to 1 when the final rupture occurs but was found to be D = 0.11. This is the critical value of the damage variable (Dcr) at which instantaneous rupture of the specimen occurs during uniaxial fatigue loading.
D–N curves were obtained as illustrated in Fig. 10. Table 1 contains the parameters rR and m obtained from these results. Similar tests were conducted at rmean = 50 MPa and for different values of rmax. The results are shown in Fig. 11 and the corresponding material parameters are listed in Table 2. With reference to Eq. (7) it can be seen that the tests conducted at zero mean stress yield directly the material parameter rR0 . The value used in the simulations was taken as the average of the four measured values. The power law exponent in Eq. (6) is not a function of the mean stress so the value used in the simulations was taken as the average of all the values measured under each mean stress condition. In order to determine the constant b in Eq. (7), tests were conducted at different mean stress levels. The value of resistance stress (rR) was calculated for each mean stress level (rmean). The results are illustrated in Fig. 12. A linear fit between rR and rmean results in b = 0.41. While
Fig. 9. Decrease in modulus of elasticity and increase in damage variable.
Fig. 10. D–N curves for rmean = 0 MPa.
Fatigue experiments were carried out for different values of
rmax while keeping rmean = 0 MPa (R = 1) and the corresponding
692
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697
Table 1 Material dependent parameters at rmean = 0 MPa.
rmean (MPa)
rmax (MPa)
R
rR (MPa)
m
Dcr
0 0 0 0
300 250 200 150
1 1 1 1
6535 6425 6503 6483
4.048 3.894 3.697 3.812
0.1039 0.1163 0.0913 0.1219
Fig. 13. AISI 4142 steel microstructure observed under optical microscope.
Fig. 11. D–N curves for rmean = 50 MPa.
Table 2 Material dependent parameters at rmean = 50 MPa.
rmean (MPa)
rmax (MPa)
R
rR (MPa)
m
Dcr
50 50 50 50
300 250 200 150
0.67 0.6 0.5 0.33
6490 6521 6475 6506
3.989 4.114 4.412 3.793
0.1340 0.1166 0.1294 0.1419
applying a linear fit, the ultimate stress (ru) for AISI 4142 steel was assumed to be 1250 MPa. 4. Application of VFEM model to tensile fatigue specimens The VFEM fatigue model described in the previous section was used to investigate the fatigue behavior of the tensile dog-bone shaped specimens. In order to identify the average grain size of steel, the optical micrograph of the material microstructure of AISI
Fig. 12. Relationship between rmean and rR.
4142 steel was used, as shown in Fig. 13. Note that the discrete material grains have been marked in the image for illustration purposes. Microstructural examination revealed grain sizes ranging from 30 lm to 50 lm which is consistent with the data reported by Sarıog˘lu [26]. Hence, while modeling the dog-bone shaped specimen in the VFEM based fatigue model, the mean Voronoi cell size was kept at 40 lm. Considering the specimen dimensions and the grain size, modeling of the entire specimen with this grain size would require a large number of elements which is computationally very expensive. Therefore, a representative volume element (RVE) technique was adopted to model only the critical region of the specimen near the mid-span, as shown in Fig. 14. Using the VFEM fatigue model, the RVE was discretized into Voronoi polygons which then were further divided into triangular elements. The finite element developed for this investigation is a 2dimensional plane stress model. The boundary conditions imposed on the RVE were such that the x displacements on one end are fixed while the x displacements on the other end were specified such that the desired stress range is obtained in the RVE. Note, that the stress range is specified for the average axial stress acting on the specimen mid-span. As mentioned earlier, the triangular elements in the RVE undergo stiffness degradation due to fatigue loading (Fig. 2). This stiffness degradation is governed by the damage evolution law which was experimentally determined and incorporated into the fatigue model. It has been recognized that in the case of high cycle fatigue and uniaxial loading much of the life is spent in the crack initiation stage [27–30]. This is due to the uniform stress field under the tensile loading, hence once an element fails (crack initiates), the neighboring elements fail shortly after (the crack grows fast). The same observation was made while performing the in situ testing (i.e., once a surface crack appeared), the specimen failed almost instantaneously. Consequently, in the current investigation initiation of the first crack in the model was chosen as the life criteria of our numerical analysis. Note that in the fatigue model, a crack is formed when the damage in a triangular element reaches the critical value of damage variable Dcr. The value of Dcr for AISI 4142 steel was measured to be 0.12 with a small degree of scatter as indicated in Tables 1 and 2. Also, a damage increment of DD = 0.01 was chosen for the numerical simulations. Fig. 15 shows the size, shape, location and orientation of the first triangle that failed during the simulation for a sample domain. The triangle failing first is marked in red color. The location of the first failed triangle was found to be consistent with experimental observation.
693
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697 Table 3 Simulation parameters. Elastic modulus of material, E Poisson’s ratio, t Ultimate strength Mean element size Damage increment, DD Resistance stress, rR0 Power exponent, m Modification co-efficient, b
210 GPa 0.3 1250 MPa 40 lm 0.01 6486.5 MPa 3.97 0.41
Fig. 14. RVE selected for analysis.
First, the effect of the topology of material microstructure on fatigue life was studied. For this purpose, forty domains with the different topographical microstructural distributions were generated. The material properties for all 40 domains were assumed to be homogenous and isotropic so that the only difference between the specimens was the difference in topology. These domains were subjected to alternating axial tensile loading (Fig. 14) with rmax = 300 MPa and rmean = 0 MPa (R = 1). Table 3 contains the simulation parameters used in this investigation. The scatter of the fatigue lives obtained are shown on a Weibull probability graph in Fig. 16. This indicates that the results follows a twoparameter Weibull distribution with a slope of 14.36 and Weibull strength (i.e. number of cycles for which the probability of failure is 63.2%) of 29,794 cycles. The high value of Weibull slope indicates that there is little scatter in fatigue life for axial loaded members and homogenous material. This can be attributed to the fact that in this part of investigation the material was assumed to be pristine with uniform material properties and no voids or flaws (micro-cracks) in the material (domain). Second, the effect of applied load on fatigue life was studied and compared with the experimental results. In this part of investigation, ten of the previously generated domains were subjected
Fig. 16. Weibull life plot for 40 different domains with pristine conditions.
Fig. 17. Comparison of experimental and numerical results.
Fig. 15. First triangle failing during simulation for a sample domain.
694
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697
to a load of rmean = 0 and rmax ranging from 350 MPa to 100 MPa. Fatigue experiments were also carried out at the same loading condition. Three specimens were tested at each of these loading conditions. Fig. 17 provides a comparison of the fatigue lives obtained from the experimental and VFEM model. Fig. 17 indicates
that the experimental and analytical VFEM are in good agreement. The degree of scatter seen in experimental data is not observed in numerical results because while performing the simulations the material was assumed to be pristine without voids or flaws.
Fig. 18. (a) Close-up view of the initial flaw. (b) Tensile stress distribution in the RVE with one initial flaw for stress amplitude of 300 MPa. (c) Zoomed view of stress distribution around the initial flaw.
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697
5. Effects of internal voids and flaws on fatigue life Voids and flaws are inherently present in material microstructure. It has been shown experimentally that these discontinuities lead to local areas of stress concentrations initiating early crack formation [28]. Internal material voids and flaws were modeled using the VFEM fatigue model in order to study their effects on fatigue life. For this purpose, a random initial flaw and an internal void were introduced in the 40 domains previously used and then the domains were subjected to rmean = 0 MPa and rmax = 300 MPa (R = 1) until a new crack was initiated. 5.1. Initial flaws An initial flaw (micro-crack) was randomly placed in the domain before loading was initiated. In order to create this initial flaw (micro-crack), the nodes of the two neighboring elements were uncoupled by defining two different nodes at the same location. Considering the average size of element sides in the RVE (40 lm), the average length for initial flaws is 40 lm. However, the exact size, location, and orientation of these initial flaws are random which cause more randomness in specimens creating more scatter in the obtained fatigue lives. Fig. 18a depicts an initial flaw in one of the specimens. The tensile stress (rx) distribution influenced by this flaw is shown in Fig. 18b and c. As expected,
695
the initial flaw acts as the internal source for stress concentration. The crack initiation lives obtained for 40 specimens having one initial flaw are illustrated in Fig. 19 along with the results obtained for the pristine specimens. The results indicate that the Weibull slope reduces from 14.36 for the pristine specimens to 1.54 and the Weibull strength decreases from 29,794 cycles to 13,140 cycles. As expected, the introduction of initial flaws decreases the fatigue lives and increases their scatter. This can be explained by considering that when the initial flaws are applied to the domain, a new randomness is introduced in addition to the randomness caused by the microstructure topology. This consequently reduces the fatigue lives and increases their scatter compared with the fatigue lives of the pristine domains. In the previous paragraph, all of the simulations were carried out at a single stress level in order to highlight the stochastic behavior of materials under fatigue loading in the presence of initial flaws. Now, we consider the effect of flaws on the stress – life ( S–N) results. The S–N data for the domains having one initial flaw was generated using the same 10 domains used for generating S–N curve of the pristine domains (Fig. 17). Fig. 20 provides a comparison between the stress – life results obtained for the 10 domains with and without one initial flaw. As seen, presence of the initial flaws in the domains causes the reduction in fatigue lives and increase in their scatter. 5.2. Internal void
Fig. 19. Weibull life plot for crack initiation in 40 different domains with initial material imperfections.
The effects of initial internal void on fatigue life are discussed in this section. A void was incorporated into the RVE by removing one randomly selected Voronoi polygon. In order to compare the results obtained for the pristine specimens and the ones having an internal void, the same 40 domains previously generated were used again. One internal void was distributed randomly inside the domains at the beginning of simulations to investigate the effects of initial internal voids on fatigue life. Fig. 21a depicts an initial internal void introduced in a domain. The tensile stress (rx) distribution influenced by this internal void is shown in Fig. 21b and c. As expected, the void acts as the internal source for stress concentration. The initiation lives obtained for the domains with one internal void are also shown in Fig. 19. The Weibull slope obtained for the domains having one internal void is 1.88 and the Weibull strength is 5419 cycles. As expected, introduction of an internal initial void causes more reduction in the fatigue lives than the introduction of initial flaws does. This can be explained by considering that when a void is introduced in the domain, the effective area resisting the load reduces. This increases the stress and consequently reduces the life. Table 4 contains comparison of Weibull slopes and strengths obtained for all three cases considered. In order to investigate the effects of internal voids on S–N results, the same 10 domains (previously used for generating S–N results for pristine condition) were used and one internal void was incorporated in them. Fig. 22 compares the S–N data obtained for domains with and without an internal void. From the results it is clear that internal voids reduces the fatigue lives and increases their scatter. 6. Summary and conclusions
Fig. 20. S–N curve for crack initiation in a tensile specimen with one initial flaw.
In this study, fatigue damage evolution in AISI 4142 steel under uniaxial variable loading was studied using a Voronoi finite element method (VFEM) based fatigue model that incorporates an empirical damage evolution law. In order to determine the empirical constants needed for the damage evolution law, in situ fatigue tests were conducted using an Instron fatigue test machine inside an environmental scanning electron microscope (ESEM). Strain
696
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697
Fig. 21. (a) Close-up view of the internal void. (b) Tensile stress distribution in the RVE with internal void for stress amplitude of 300 MPa. (c) Zoomed view of stress distribution around the internal void.
data obtained from ESEM images were used to determine the variation of the elastic modulus during fatigue loading. This variation was used to obtain the damage variable. Fatigue tests were conducted at various mean stress levels to acquire the D–N data which later was used to determine the empirical constants involved in the fatigue model. The empirical constants were incorporated into the VFEM fatigue model and simulations were performed on a representative volume element. Forty different domains with different microstructural distributions were generated for numerical analysis and initiation of first crack was selected as the life criteria. The number of cycles to failure predicted by the model at various loading conditions was corroborated with the experimental results.
There was a good agreement between the experimental and numerical findings. The model was also used to study the effects of initial flaws and internal voids on the fatigue life to demonstrate
Table 4 Weibull probability statistics. Material condition
Weibull slope, e
Weibull strength (cycles)
Pristine One initial flaw Initial internal void
14.36 1.54 1.88
29,794 13,140 5419
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697
Fig. 22. S–N curve for crack initiation in a tensile specimen with initial internal void.
the ability of model to capture material imperfections. The presence of one internal void was found to be more detrimental than one initial flaw whereas more scatter in fatigue life data was observed in the presence of an initial flaw. Overall, the required data for the VFEM fatigue model is the material dependent parameters which are used in the damage evolution law. Once these parameters are obtained from a simple tension test, the VFEM fatigue model can then be used to study the behavior of different components subject to diverse types of loading. The model can also be used to investigate material flaws and voids and obtain a fatigue life distribution for the component. Acknowledgement The authors would like to express their deepest appreciations to the sponsors of the Mechanical Engineering Tribology Laboratory (METL), Purdue University, West Lafayette, IN for their continued interest and support of METL. References [1] Bolotin VV, Belousov IL. Early fatigue crack growth as the damage accumulation process. Probab Eng Mech 2001;60:279–87. [2] Bolotin VV, Babkin AA, Belousov IL. Probabilistic Model of Early Fatigue Crack Growth 1998;13:227–32. [3] Suresh S. Fatigue of materials. Cambridge (UK): Cambridge University Press; 1998.
697
[4] Lemaitre J. A course on damage mechanics. Berlin, Germany: Springer-Verlag; 1992. [5] Rinaldi A, Peralta P, Krajcinovic D, Lai YC. Prediction of scatter in fatigue properties using discrete damage mechanics. Int J Fatigue 2006;28:1069–80. [6] Raje N, Sadeghi F, Rateick RG. A statistical damage mechanics model for subsurface initiated spalling in rolling contacts. J Trib 2008;130:042201. [7] Raje N, Slack T, Sadeghi F. A discrete damage mechanics, model for high cycle fatigue in polycrystalline materials subject to rolling contact. Int J Fatigue 2009;31:346–60. [8] Okabe A, Boots B. Spatial tessellations: concepts and applications of Voronoi diagrams. New York: John Wiley & Sons; 1992. [9] Moller J. Lectures notes on random Voronoi tessellations. Berlin: Springer; 1994. [10] Ghosh S, Mallett RL. Voronoi cell finite elements. Comp Struct 1994;50(1994):33–46. [11] Ghosh S, Moorthy S. Elastic–plastic analysis of arbitrary heterogeneous materials with the Voronoi finite element method. Comp Methods Appl Mech Eng 1995;121:373–409. [12] Moorthy S, Ghosh S. A model for analysis of arbitrary composite, porous microstructures with Voronoi cell finite elements. Int J Numer Methods Eng 1996;39:2363–98. [13] Grujicic M, Zhang Y. Determination of effective elastic properties of functionally graded materials using Voronoi cell finite element method. Mater Sci Eng A – Struct Mater Proper Microstruct Proc 1998;251:64–76. [14] Guo R, Shi HJ, Yao ZH. Modeling of interfacial debonding crack in particle reinforced composites using Voronoi cell finite element method. Comp Mech 2003;32:52–9. [15] Vena P, Gastaldi D. A Voronoi cell finite element model for the indentation of graded ceramic composites. Composites B – Eng 2005;36:115–26. [16] Jalalahmadi B, Sadeghi F. A Voronoi finite element study of fatigue life scatter in rolling contacts. ASME J Trib 2009;131(2). [17] Kachanov LM. Time of the rupture process under creep conditions. Izv Akad Nauk SSR, Otd Tekh Nauk 1958;8:26–31. [18] Robotnov YN. Creep problems in structural mechanics. North-Holland; 1969. [19] Xiao YC, Li S, Gao Z. A continuum damage mechanics model for high cycle fatigue. Int J Fatigue 1998;20:503–8. [20] Bolotin VV. Mechanics of fatigue. Boca Raton, Florida: CRC Press; 1999. [21] Biallas G, Maier HJ. In-situ fatigue in an environmental scanning microscopepotential and current limitation. Int J Fatigue 2007;29:1413–25. [22] Wissuchek DJ, Mackin TJ, DeGraef M, Lucas GE, Evans AG. A simple method for measuring surface strains around cracks. Exp Mech 1996;36:173–9. [23] Hung PC, Voloshin AS. In-plane strain measurement by digital image correlation. J Braz Soc Mech Sci Eng 2003;25:215–21. [24] James MR, Morris WL, Cox BN. A high accuracy automated strain-field mapper. Exp Mech 1990;30:60–7. [25] Biery N, De Graef M, Pollock TM. A method for measuring microstructuralscale strains using a scanning electron microscope: application to c-titanium aluminides. Metall Mater Trans 2003;34 A:2301–13. [26] Sarıog˘lu F. The effect of tempering on susceptibility to stress corrosion cracking of AISI 4140 steel in 33% sodium hydroxide at 80 °C. Mater Sci Eng A 2001;315:98–102. [27] Andersson J. The influence of grain size variation on metal fatigue. Int J Fatigue 2005;27(8):847–52. [28] Grandt Jr AF. Fundamentals of structural integrity. NJ (USA): John Wiley & Sons, Inc.; 2004. [29] Miller KJ, Zachariah KP. Cumulative damage laws for fatigue crack initiation and stage I propagation. J Strain Anal 1977;12:262–70. [30] Lin TH, Chen QY. Micromechanics of high-cycle initiation of fatigue cracking intermetallic compounds. Acta Metall Mater 1992;40:3255–63.
Contents lists available at ScienceDirect
International Journal of Fatigue journal homepage: www.elsevier.com/locate/ijfatigue
A new finite element fatigue modeling approach for life scatter in tensile steel specimens Anurag Warhadpande, Behrooz Jalalahmadi, Trevor Slack, Farshid Sadeghi * School of Mechanical Engineering, Rm #303, Purdue University, West Lafayette, IN 47907, United States
a r t i c l e
i n f o
Article history: Received 15 April 2009 Received in revised form 7 September 2009 Accepted 7 October 2009 Available online 7 November 2009 Keywords: Damage mechanics Fatigue Voronoi finite element fatigue modeling
a b s t r a c t This paper presents the results of a finite element model developed to investigate the effects of material microstructure topology on fatigue damage evolution in dog-bone tensile specimens. The Voronoi finite element model (VFEM) developed was used to study fatigue damage evolution in AISI 4142 steel under various loading conditions. An environmental scanning electron microscope (ESEM) coupled with an Instron fatigue testing machine was used to determine the strain and fatigue life of the dog-bone 4142 steel specimens at different mean stress levels. The damage variable was obtained using the method of variation of elasticity modulus. The experimental damage versus cycle curve was then used to determine the material properties (i.e. resistance stress rR and exponent m) needed for the VFEM fatigue damage model. The VFEM model was used to corroborate the analytical and experimental fatigue results. A comparison of the results obtained from the VFEM model and ESEM indicate that they are in good agreement. Ó 2009 Elsevier Ltd. All rights reserved.
1. Introduction Fatigue is primarily a process of damage accumulation which is manifested by progressive internal deterioration of material properties due to nucleation and growth of micro-cracks, debonding, voids, etc. The mechanism of fatigue failure in general consists of three stages: (i) crack initiation, (ii) crack propagation and (iii) final catastrophic failure. The percentage of the total life spent in each stage is related to material properties and loading conditions. Designing against fatigue failures is a time consuming and involves extensive experimentation to account for the statistical nature of the fatigue process. For high cycle fatigue much of the variability in fatigue life data can be accounted for in the crack initiation stage [1]. Crack initiation is a highly localized event which is strongly affected by the random distribution of material properties and defects that exist at the microstructural level [2,3]. Furthermore, engineering materials are polycrystalline, consisting of individual grains with random geometrical properties separated by grain boundaries that can act as physical discontinuities. Modeling approaches which include this topological and material property randomness can be used to explain fatigue scatter and to reduce the number of experiments required to quantify it. Damage mechanics approach [4] for modeling material deterioration and fatigue is increasing in popularity. Damage mechanics
* Corresponding author. Tel.: +1 765 494 5719; fax: +1 765 494 0539. E-mail addresses: [email protected] (A. Warhadpande), bjalalah@ purdue.edu (B. Jalalahmadi), [email protected] (T. Slack), [email protected] (F. Sadeghi). 0142-1123/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijfatigue.2009.10.003
introduces a damage variable, either scalar or tensorial, into the material constitutive model to account for the degradation in material properties that occurs as damage grows. The evolution of the damage variable is governed either by a phenomenological relationship or by a relationship developed based on micromechanical principals. Fatigue models based on damage mechanics concepts have been implemented in both continuum and discrete frameworks. Rinaldi et al. [5] proposed a lattice model to simulate polycrystalline microstructures and studied the resulting scatter in fatigue response at the macro level. In their model the fatigue behavior of each discrete spring was described using the Basquin law and damage accumulated according to a Palmgren–Miner law. By introducing a normal distribution of static strengths for the springs fatigue scatter was found to follow a lognormal distribution. Bolotin et al. [2] developed a model for the nucleation and early growth of fatigue cracks using the finite element method. In their model each grain was approximated by a square finite element and was randomly assigned a resistance stress from a Weibull distribution. Here resistance stress refers to a material parameter used in the fatigue damage evolution law. Their results exhibited many features such as a non-monotonic change in growth rate and scatter in the crack dimensions and growth rates that have been experimentally observed in the early propagation of fatigue cracks. Recently, Raje et al. [6,7] proposed a model for the fatigue of bearing contacts using the discrete element method. In their model the material was represented as an assemblage of discrete, rigid Voronoi polygons connected by normal and tangential springs. Damage accumulated according to a continuum damage model and was used to degrade the tangential spring stiffnesses.
686
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697
Nomenclature ~ A A b D Dcr e E E e Ec e Et F h L1 L2 m N R
effective cross-sectional area supporting the load overall cross-sectional area modification co-efficient damage variable critical value of damage variable damaged modulus of elasticity modulus of elasticity damaged elastic moduli in compression damaged elastic moduli in tension force acting on area A crack closure factor distance between the spots before application of load distance between the spots after application of load power exponent number of cycles minimum cyclic stress divided by maximum cyclic stress
By including randomness in the material topology and material properties they obtained Weibull slopes of the simulated life data and final spall crack patterns that compared well with experimental results. In this investigation a fatigue life model is developed using a Voronoi finite element method (VFEM) with a damage coupled constitutive relationship. The randomness of the material microstructure is simulated using a Voronoi tessellation process. The Voronoi cells are discretized using constant strain triangular finite elements and damage accumulates in each triangle according to a continuum damage model. Uniaxial in situ fatigue tests were conducted using an environmental scanning electron microscope (ESEM) equipped with an Instron fatigue testing machine to measure both damage as a function of loading cycle and life as a function of applied stress. The damage as a function of loading cycle data is used to identify the empirical constants used in the damage evolution law by the method of variation of elasticity modulus. The VFEM model was used to predict the stress versus life results which was then compared against experimental results. Initial imperfections in the form of material flaws and voids were then introduced into the model to investigate their effects on fatigue life scatter. Experimental and numerical VFEM results demonstrate good agreement over a wide range of loading conditions.
DD DN Dr
e m~ m
r~ r rmax rmean rR rR0 ru
increment in damage increment in number of cycles stress range axial strain damaged Poisson’s ratio Poisson’s ratio effective stress stress acting on area A peak value of stress mean value of stress resistance stress resistance stress corresponding to zero mean stress level ultimate strength
investigation. Care is taken that during the tessellation process; the element sizes are relatively uniform so that there is no significant variation in element size across the domain. This is achieved by setting an upper and lower bounds on the distance between the nucleation points. 2.2. Finite element modeling Previous investigations have in general used the classical finite element methods (for example using quadrilateral elements for the domain), which employs the homogenization techniques to simulate the material as an equivalent continuum model on the macroscale. However, using standard finite element software that contains a limited variety of elements having predefined shapes is not suitable because they cannot be used to properly simulate randomly shaped topology of material microstructure. Ghosh and coworkers [10–12] introduced the Voronoi Cell Finite Element Method (VCFEM) for deformation and stress analysis of domains with arbitrary microstructural distributions. Other investigators
2. Voronoi finite element method based fatigue modeling In the following section, the approach used to develop the Voronoi finite element model to investigate fatigue of dog-bone shaped steel specimens is presented. We will briefly describe the Voronoi tessellation process, the new damage mechanics based finite element model developed and proceed to describe the test rig used in this investigation. Then the experimental and analytical results of this study are presented with some concluding remarks. 2.1. Voronoi tessellation In the present study, each simulation domain is generated using the Voronoi tessellation method. A detailed description of the Voronoi tessellation process can be found in [8,9]. The set of points used to generate the Voronoi diagram are called the nucleation points. These nucleation points are distributed randomly before the start of the Voronoi tessellation process resulting in distinct simulation domains. The density of the nucleation points is decided in accordance with the grain size of the material under
Fig. 1. (a) Discretizing a domain with Voronoi elements. (b) Dividing the Voronoi cell into triangular elements. (c) Linear triangular finite element.
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697
[13–15] have also applied VCFEM to study various problems. Voronoi polygons generated for this investigation have random size and shape and are used to simulate material topological disorder and its corresponding effects on stress distribution and fatigue life. In this investigation, a new Voronoi finite element model [16] was developed to investigate fatigue in steel tensile specimens. For this purpose, the domain is divided into Voronoi cells using the Voronoi tessellation and then the geometrical center of each Voronoi element is determined. Because Voronoi elements are convex polygons, it is assured that the center of each element is located within the element. This property is used to divide the Voronoi element into finer triangular elements by connecting the center of the cell to each node. Fig. 1 illustrates the process of discretizing a domain into the Voronoi cells and then dividing the cell into finer triangular elements. A two-step assembly process is used to obtain the global stiffness matrix (K). In the first step, assuming linear shape functions for each triangle, its stiffness matrix is obtained and then the stiffness matrix for each Voronoi element was set up using the stiffness matrices for each individual triangle. In the second step, the global K matrix is assembled using the stiffness matrices of the Voronoi elements. 2.3. Damage variable Fatigue damage leads to progressive degradation of the material due to formation and growth of micro-cracks and voids. Using damage mechanics, these microscopic failure mechanisms can be treated in an empirical fashion by using the internal damage variable D. The damage variable D in general is a tensor, however, under isotropic conditions the damage variable D reduces to a scalar. It is also assumed that the Poisson’s ratio remains unaffected by damage [4]. The damage variable D accounts for the surface density of the micro-cracks and cavities in any plane of a representative volume element of the material. Kachanov [17] defined D as,
e AA D¼ A
ð1Þ
e is the effective area which carwhere A is the apparent area and A ries the load. Robotnov [18] introduced the concept of an effective stress to describe the effect of damage on strain behavior. The effec~ is defined as, tive stress r
F
r~ ¼ e ¼ A
F A r ¼ e A A 1D
ð2Þ
Applying the Hooke’s law Eq. (2) can be written in terms of the effective stress as,
e¼
r~ E
¼
r Eð1 DÞ
ð3Þ
687
where E is the modulus of elasticity and e is the elastic strain of the undamaged material. Defining e E as the elastic modulus of the damaged material,
e E ¼ Eð1 DÞ
ð4Þ
The reduction in elastic modulus with the increase in fatigue cycles is shown in Fig. 2. The damage variable is thus given by,
D¼1
e E E
ð5Þ
Eq. (5) can be used to indirectly measure the damage variable D in a process known as the method of variation of elasticity modulus. 2.4. Fatigue damage modeling Damage was incorporated into the VFEM model using the damage coupled constitutive relationship, Eq. (4), for each triangular element. In this analysis the evolution of the damage variable is described by the following model proposed by Xiao et al. [19] for high cycle fatigue
dD ¼ dN
Dr rR ð1 DÞ
m ð6Þ
Here, Dr is the stress range, N is the cycle number, and rR and m are material parameters. The parameter rR is sometimes referred to as the resistance stress [20] since it controls the ability of a material to resist damage accumulation. As the resistance stress increases the rate of damage accumulation decreases and consequently the fatigue life is longer. Xiao et al. [19] indicate that the resistance stress (rR) is a function of the mean stress and the relationship between rR and rmean is given by,
rR ¼ rR0 1 b
rmean ru
ð7Þ
where b is the modification co-efficient, rR0 is the resistance stress under completely reversed loading conditions and ru is the ultimate strength of the material. In order to capture the effects of micro-crack closure in compression, an additional term h is introduced in the constitutive model. The closure parameter h = 1 in tension and 0 6 h < 1 in compression. A value of h = 0 corresponds to complete micro-crack closure. In the present analysis, a value of h = 0.2 as proposed by Lemaitre [4] is used. The damaged coupled constitutive model becomes:
e E c ¼ Eð1 hDÞ e E t ¼ Eð1 DÞ
ð8Þ
m~ ¼ m
Fig. 2. Stiffness degradation with increase in fatigue cycles.
where E is the initial elastic modulus under pristine conditions and e E t are the damaged elastic moduli in compression and tenE c and e sion respectively. Fatigue modeling involves the simultaneous solution of Eqs. (6)–(8) to account for the coupling between damage and the material constitutive equations for each triangular element. Thus, the stress–strain relationships need to be solved and the triangular element elastic moduli need to be updated after every stress cycle to incorporate the current state of damage in the triangular element. However, it is impractical to perform the large number of increments needed for the millions of stress cycles arising in high cycle fatigue. Therefore, a method that involves the use of the ‘jump-incycles’ procedure [4] was used for achieving computational efficiency.
688
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697
The method assumes a piecewise periodic constant amplitude loading, i.e., the stress histories for the triangular elements are assumed to remain unaltered over a finite number of cycles DN constituting a block i. The increment in damage DD is assumed to be constant over a block of cycles. Thus, the damage evolution is assumed to be piecewise linear with respect to the number of cycles. It is to be noted that the shape of the damage evolution curve is not predetermined but it is an outcome of the numerical simulations due to the stress–damage coupling. The value of DD is a tradeoff between accuracy and computation time. For higher resolution, DD needs to be chosen as small as possible, however it directly affects the computation time. Stress calculations are performed once for every block of cycles and not for each individual cycle N. This saves a considerable amount of computation time and effort since the major part of the computational expense is spent in solving for the stress fields in the domain. The number of cycles in a block is variable and it is obtained using the following procedure: (i) Initial damage states for each triangular element are assigned. For a pristine material domain containing no initial flaws, the initial damage in each triangular element is assumed to be 0. Thus,
D0j ¼ 0;
j ¼ 1 . . . ntriangular elements
ð9Þ
The number of cycles elapsed is initialized to 0. Thus,
N¼0
ð10Þ i min Þj
(ii) Stress histories ðrmax r for each triangular element are computed for the current block i. (iii) The damage evolution rate in each triangular element is evaluated knowing the stress history for the present block and the current state of damage in the triangular element using,
dD dN
"
i ¼ j
ðDrÞij
#m
rR ð1 Dij Þ
ð11Þ
(iv) The triangular element with the maximum damage evolution rate is chosen as the critical triangular element for the current block of cycles. Thus,
dD dN
i crit
dD i ¼ Max dN j
ð12Þ
(v) The number of cycles in the current block of cycles is computed as,
DD DNi ¼ i dD
ð13Þ
dN crit
(vi) The number of cycles elapsed is updated to,
N ¼ N þ DN i
ð14Þ
(vii) The increment in damage for each triangular element during the current block of cycles is then given by,
DDij ¼
i dD DN i dN j
ð15Þ
(vii) The damage states for each triangular element at the start of the next block of cycles are updated to,
Djiþ1 ¼ Dij þ DDij
ð16Þ
(ix) The elastic moduli are modified at the start of the next block of cycles according to,
i ðEc Þij ¼ ðEÞ 1 hDj
ðEt Þij ¼ ðEÞ 1 Dij
ð17Þ
Steps (i) through (ix) are successively repeated for each block of cycles. A micro-crack is assumed to be initiated at a triangular element when the accumulated damage in the triangular element equals Dcr i.e. the critical value of the damage variable D. 3. Experimental approach to determine the material dependent constants The damage evolution law introduces three empirical constants
rR0 , m and b that have to be experimentally evaluated. The following sections describe the test apparatus, fatigue test specimen geometry and the experimental scheme developed to estimate the material constants. 3.1. Experimental setup The fatigue experiments were conducted using an ESEM equipped with an Instron load frame [21]. The Instron load frame has a servo-hydraulic actuator and controller which can be used to induce the load cycles in the test specimen. Fig. 3 illustrates the ESEM and Instron load frame with the test specimen installed between the grippers. One end of the specimen is firmly held in the fixed jaws while the other end is actuated using the servo-hydraulic system. Fig. 4 depicts the load frame installed inside the ESEM which allows for in situ monitoring of the test specimen during fatigue testing. The tests were conducted with water vapor as the gaseous medium inside the ESEM vacuum chamber at a partial pressure of 2–4 Torr and at ambient room temperature. 3.2. Specimen geometry The specimens used in this investigation were made of AISI 4142 steel. Fig. 5 depicts the two different test specimen geometries evaluated for fatigue testing. One of the test specimens had a uniform gauge cross section, while the other had a non-uniform gauge cross-sectional area with the minimum cross section occurring at the mid-span. When the test specimen with uniform gauge cross-sectional area (Fig. 5) was subjected to cyclic loading, failure occurred at various locations with the majority occurring near the shoulder of the fillet radius. This was attributed to manufacturing errors associated with the fillet radius; and the effects of clamping leading to fracture as shown in Fig. 5. It is to be noted that the ESEM and load frame are designed such that the electron detector can focus only in the central region of specimen gauge section. Hence the specimen geometry was modified to one with a nonuniform cross section leading to a region of maximum stress concentration occurring at the central section of the specimen. In order to obtain a continuously decreasing cross-sectional area an elliptical profile was machined on the specimen with the major and minor axis of ellipse as indicated in Fig. 6. With this specimen geometry the electron detector could be focused on the region of maximum stress concentration where the failure occurred. During the course of the fatigue testing, a crack typically appeared at the central region of the specimen where the cross-sectional area is minimum. The crack initiated at one of the edges and propagated along the thickness and width of the specimen as shown in Fig. 7. Once the crack extended through the thickness of the specimen, instantaneous fracture occurred.
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697
689
Fig. 3. ESEM and Instron fatigue test machine.
Fig. 4. Instron load frame installed inside the ESEM vacuum chamber.
3.3. Strain measurement In this section, a brief description of the strain measurement technique used in this study is presented. The digital images obtained from ESEM during in situ fatigue testing were used to determine the strain data [22–24] from which the damage variable is measured. The method adopted here utilizes features present on the surface of the test specimen (such as microstructure and material inhomogeneity) which provide useful contrast. This approach in strain measurement requires monitoring of the distance between two dark spots on the grayscale image. The change in length divided by the original length between the two spots provides the strain under the applied load. While performing the in situ fatigue testing, the electron detector is focused along the specimen thickness near the minimum cross-sectional area of the specimen gauge section. The region is scanned and two spots along the loading direction are selected.
The spots are chosen such that they span a large-enough number of pixels so that their shapes and positions are well defined and the area of interest is fully covered [25]. Fig. 8 depicts ESEM micrographs of two spots monitored during the loading and unloading of a specimen. While generating these images, the contrast and brightness were kept constant. Digital images were produced at a resolution of 2048 2048 pixels. From the ESEM images, grayscale intensity plots were generated along the two selected spots as shown in Fig. 8. The two deep valleys in the figures correspond to the two dark spots on the test specimen. Both valleys are associated with particular pixel numbers. The difference in pixels between the valleys is converted to the distance between the dark spots using a scaling factor. The scaling factor is a function of image magnification and provides the relationship between pixel distances on the digital image and the actual distance on the specimen. It is obtained from the scale on the digital image. Once the distances are obtained the strain can be calculated using,
690
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697
Fig. 5. Specimen geometry.
3.4. Identification of material dependent constants The empirical constants rR0 , m and b introduced by the damage evolution law are evaluated as follows. Assuming that the stress range is constant in the damage evolution equation, Eq. (6), is integrated to obtain the following relationship between the damage variable D and fatigue cycle N,
( D¼
Fig. 6. Drawing of the dog-bone shaped specimen (all dimensions are in mm).
e¼
L2 L1 L1
ð18Þ
where L1 and L2 are the distances between the dark spots before and after the application of load, respectively.
1 ) m mþ1 Dr 1 1 ðm þ 1Þ N
rR
ð19Þ
For a given fatigue test, the damage variable D was measured as a function of fatigue cycle N, generating an experimental D–N curve, using the following procedure. During a fatigue test, after DN cycles, the fatigue test was stopped. The specimen was loaded to a certain known value of stress and the corresponding strain induced in the specimen was measure using the method described in the previous section (Eq. (18)). The modulus of elasticity of the specimen at that instant of time was determined by dividing the applied value of stress by the measured value of strain. Since the test specimen has undergone fatigue cycling, the material softens
Fig. 7. Typical crack propagation pattern.
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697
691
Fig. 8. ESEM micrographs and grayscale intensity plots used to obtain strain data before and after load application.
and the modulus measured is the damaged modulus of elasticity e E (refer to Eq. (5)). Thus, using Eq. (5) the damage variable D was determined. In Eq. (5) E is the undamaged modulus of elasticity of the specimen measured at the start of each fatigue test. The fatigue test was continued for another DN cycles and the above mentioned procedure was repeated until the specimen fractured. Measuring D as a function of N during the course of fatigue test resulted in the D–N curve. Fig. 9 depicts the degradation in modulus of elasticity E and corresponding increase in damage variable D for a sample test (rmax = 300 MPa, rmean = 0 MPa, R = 1). Eq. (19) was used to curve fit the D–N data and obtain the empirical constants rR and m. Least square regression was used to fit the equation and the R-square value for the fit was 96%. It is to be noted that the value of the damage variable D is not equal to 1 when the final rupture occurs but was found to be D = 0.11. This is the critical value of the damage variable (Dcr) at which instantaneous rupture of the specimen occurs during uniaxial fatigue loading.
D–N curves were obtained as illustrated in Fig. 10. Table 1 contains the parameters rR and m obtained from these results. Similar tests were conducted at rmean = 50 MPa and for different values of rmax. The results are shown in Fig. 11 and the corresponding material parameters are listed in Table 2. With reference to Eq. (7) it can be seen that the tests conducted at zero mean stress yield directly the material parameter rR0 . The value used in the simulations was taken as the average of the four measured values. The power law exponent in Eq. (6) is not a function of the mean stress so the value used in the simulations was taken as the average of all the values measured under each mean stress condition. In order to determine the constant b in Eq. (7), tests were conducted at different mean stress levels. The value of resistance stress (rR) was calculated for each mean stress level (rmean). The results are illustrated in Fig. 12. A linear fit between rR and rmean results in b = 0.41. While
Fig. 9. Decrease in modulus of elasticity and increase in damage variable.
Fig. 10. D–N curves for rmean = 0 MPa.
Fatigue experiments were carried out for different values of
rmax while keeping rmean = 0 MPa (R = 1) and the corresponding
692
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697
Table 1 Material dependent parameters at rmean = 0 MPa.
rmean (MPa)
rmax (MPa)
R
rR (MPa)
m
Dcr
0 0 0 0
300 250 200 150
1 1 1 1
6535 6425 6503 6483
4.048 3.894 3.697 3.812
0.1039 0.1163 0.0913 0.1219
Fig. 13. AISI 4142 steel microstructure observed under optical microscope.
Fig. 11. D–N curves for rmean = 50 MPa.
Table 2 Material dependent parameters at rmean = 50 MPa.
rmean (MPa)
rmax (MPa)
R
rR (MPa)
m
Dcr
50 50 50 50
300 250 200 150
0.67 0.6 0.5 0.33
6490 6521 6475 6506
3.989 4.114 4.412 3.793
0.1340 0.1166 0.1294 0.1419
applying a linear fit, the ultimate stress (ru) for AISI 4142 steel was assumed to be 1250 MPa. 4. Application of VFEM model to tensile fatigue specimens The VFEM fatigue model described in the previous section was used to investigate the fatigue behavior of the tensile dog-bone shaped specimens. In order to identify the average grain size of steel, the optical micrograph of the material microstructure of AISI
Fig. 12. Relationship between rmean and rR.
4142 steel was used, as shown in Fig. 13. Note that the discrete material grains have been marked in the image for illustration purposes. Microstructural examination revealed grain sizes ranging from 30 lm to 50 lm which is consistent with the data reported by Sarıog˘lu [26]. Hence, while modeling the dog-bone shaped specimen in the VFEM based fatigue model, the mean Voronoi cell size was kept at 40 lm. Considering the specimen dimensions and the grain size, modeling of the entire specimen with this grain size would require a large number of elements which is computationally very expensive. Therefore, a representative volume element (RVE) technique was adopted to model only the critical region of the specimen near the mid-span, as shown in Fig. 14. Using the VFEM fatigue model, the RVE was discretized into Voronoi polygons which then were further divided into triangular elements. The finite element developed for this investigation is a 2dimensional plane stress model. The boundary conditions imposed on the RVE were such that the x displacements on one end are fixed while the x displacements on the other end were specified such that the desired stress range is obtained in the RVE. Note, that the stress range is specified for the average axial stress acting on the specimen mid-span. As mentioned earlier, the triangular elements in the RVE undergo stiffness degradation due to fatigue loading (Fig. 2). This stiffness degradation is governed by the damage evolution law which was experimentally determined and incorporated into the fatigue model. It has been recognized that in the case of high cycle fatigue and uniaxial loading much of the life is spent in the crack initiation stage [27–30]. This is due to the uniform stress field under the tensile loading, hence once an element fails (crack initiates), the neighboring elements fail shortly after (the crack grows fast). The same observation was made while performing the in situ testing (i.e., once a surface crack appeared), the specimen failed almost instantaneously. Consequently, in the current investigation initiation of the first crack in the model was chosen as the life criteria of our numerical analysis. Note that in the fatigue model, a crack is formed when the damage in a triangular element reaches the critical value of damage variable Dcr. The value of Dcr for AISI 4142 steel was measured to be 0.12 with a small degree of scatter as indicated in Tables 1 and 2. Also, a damage increment of DD = 0.01 was chosen for the numerical simulations. Fig. 15 shows the size, shape, location and orientation of the first triangle that failed during the simulation for a sample domain. The triangle failing first is marked in red color. The location of the first failed triangle was found to be consistent with experimental observation.
693
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697 Table 3 Simulation parameters. Elastic modulus of material, E Poisson’s ratio, t Ultimate strength Mean element size Damage increment, DD Resistance stress, rR0 Power exponent, m Modification co-efficient, b
210 GPa 0.3 1250 MPa 40 lm 0.01 6486.5 MPa 3.97 0.41
Fig. 14. RVE selected for analysis.
First, the effect of the topology of material microstructure on fatigue life was studied. For this purpose, forty domains with the different topographical microstructural distributions were generated. The material properties for all 40 domains were assumed to be homogenous and isotropic so that the only difference between the specimens was the difference in topology. These domains were subjected to alternating axial tensile loading (Fig. 14) with rmax = 300 MPa and rmean = 0 MPa (R = 1). Table 3 contains the simulation parameters used in this investigation. The scatter of the fatigue lives obtained are shown on a Weibull probability graph in Fig. 16. This indicates that the results follows a twoparameter Weibull distribution with a slope of 14.36 and Weibull strength (i.e. number of cycles for which the probability of failure is 63.2%) of 29,794 cycles. The high value of Weibull slope indicates that there is little scatter in fatigue life for axial loaded members and homogenous material. This can be attributed to the fact that in this part of investigation the material was assumed to be pristine with uniform material properties and no voids or flaws (micro-cracks) in the material (domain). Second, the effect of applied load on fatigue life was studied and compared with the experimental results. In this part of investigation, ten of the previously generated domains were subjected
Fig. 16. Weibull life plot for 40 different domains with pristine conditions.
Fig. 17. Comparison of experimental and numerical results.
Fig. 15. First triangle failing during simulation for a sample domain.
694
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697
to a load of rmean = 0 and rmax ranging from 350 MPa to 100 MPa. Fatigue experiments were also carried out at the same loading condition. Three specimens were tested at each of these loading conditions. Fig. 17 provides a comparison of the fatigue lives obtained from the experimental and VFEM model. Fig. 17 indicates
that the experimental and analytical VFEM are in good agreement. The degree of scatter seen in experimental data is not observed in numerical results because while performing the simulations the material was assumed to be pristine without voids or flaws.
Fig. 18. (a) Close-up view of the initial flaw. (b) Tensile stress distribution in the RVE with one initial flaw for stress amplitude of 300 MPa. (c) Zoomed view of stress distribution around the initial flaw.
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697
5. Effects of internal voids and flaws on fatigue life Voids and flaws are inherently present in material microstructure. It has been shown experimentally that these discontinuities lead to local areas of stress concentrations initiating early crack formation [28]. Internal material voids and flaws were modeled using the VFEM fatigue model in order to study their effects on fatigue life. For this purpose, a random initial flaw and an internal void were introduced in the 40 domains previously used and then the domains were subjected to rmean = 0 MPa and rmax = 300 MPa (R = 1) until a new crack was initiated. 5.1. Initial flaws An initial flaw (micro-crack) was randomly placed in the domain before loading was initiated. In order to create this initial flaw (micro-crack), the nodes of the two neighboring elements were uncoupled by defining two different nodes at the same location. Considering the average size of element sides in the RVE (40 lm), the average length for initial flaws is 40 lm. However, the exact size, location, and orientation of these initial flaws are random which cause more randomness in specimens creating more scatter in the obtained fatigue lives. Fig. 18a depicts an initial flaw in one of the specimens. The tensile stress (rx) distribution influenced by this flaw is shown in Fig. 18b and c. As expected,
695
the initial flaw acts as the internal source for stress concentration. The crack initiation lives obtained for 40 specimens having one initial flaw are illustrated in Fig. 19 along with the results obtained for the pristine specimens. The results indicate that the Weibull slope reduces from 14.36 for the pristine specimens to 1.54 and the Weibull strength decreases from 29,794 cycles to 13,140 cycles. As expected, the introduction of initial flaws decreases the fatigue lives and increases their scatter. This can be explained by considering that when the initial flaws are applied to the domain, a new randomness is introduced in addition to the randomness caused by the microstructure topology. This consequently reduces the fatigue lives and increases their scatter compared with the fatigue lives of the pristine domains. In the previous paragraph, all of the simulations were carried out at a single stress level in order to highlight the stochastic behavior of materials under fatigue loading in the presence of initial flaws. Now, we consider the effect of flaws on the stress – life ( S–N) results. The S–N data for the domains having one initial flaw was generated using the same 10 domains used for generating S–N curve of the pristine domains (Fig. 17). Fig. 20 provides a comparison between the stress – life results obtained for the 10 domains with and without one initial flaw. As seen, presence of the initial flaws in the domains causes the reduction in fatigue lives and increase in their scatter. 5.2. Internal void
Fig. 19. Weibull life plot for crack initiation in 40 different domains with initial material imperfections.
The effects of initial internal void on fatigue life are discussed in this section. A void was incorporated into the RVE by removing one randomly selected Voronoi polygon. In order to compare the results obtained for the pristine specimens and the ones having an internal void, the same 40 domains previously generated were used again. One internal void was distributed randomly inside the domains at the beginning of simulations to investigate the effects of initial internal voids on fatigue life. Fig. 21a depicts an initial internal void introduced in a domain. The tensile stress (rx) distribution influenced by this internal void is shown in Fig. 21b and c. As expected, the void acts as the internal source for stress concentration. The initiation lives obtained for the domains with one internal void are also shown in Fig. 19. The Weibull slope obtained for the domains having one internal void is 1.88 and the Weibull strength is 5419 cycles. As expected, introduction of an internal initial void causes more reduction in the fatigue lives than the introduction of initial flaws does. This can be explained by considering that when a void is introduced in the domain, the effective area resisting the load reduces. This increases the stress and consequently reduces the life. Table 4 contains comparison of Weibull slopes and strengths obtained for all three cases considered. In order to investigate the effects of internal voids on S–N results, the same 10 domains (previously used for generating S–N results for pristine condition) were used and one internal void was incorporated in them. Fig. 22 compares the S–N data obtained for domains with and without an internal void. From the results it is clear that internal voids reduces the fatigue lives and increases their scatter. 6. Summary and conclusions
Fig. 20. S–N curve for crack initiation in a tensile specimen with one initial flaw.
In this study, fatigue damage evolution in AISI 4142 steel under uniaxial variable loading was studied using a Voronoi finite element method (VFEM) based fatigue model that incorporates an empirical damage evolution law. In order to determine the empirical constants needed for the damage evolution law, in situ fatigue tests were conducted using an Instron fatigue test machine inside an environmental scanning electron microscope (ESEM). Strain
696
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697
Fig. 21. (a) Close-up view of the internal void. (b) Tensile stress distribution in the RVE with internal void for stress amplitude of 300 MPa. (c) Zoomed view of stress distribution around the internal void.
data obtained from ESEM images were used to determine the variation of the elastic modulus during fatigue loading. This variation was used to obtain the damage variable. Fatigue tests were conducted at various mean stress levels to acquire the D–N data which later was used to determine the empirical constants involved in the fatigue model. The empirical constants were incorporated into the VFEM fatigue model and simulations were performed on a representative volume element. Forty different domains with different microstructural distributions were generated for numerical analysis and initiation of first crack was selected as the life criteria. The number of cycles to failure predicted by the model at various loading conditions was corroborated with the experimental results.
There was a good agreement between the experimental and numerical findings. The model was also used to study the effects of initial flaws and internal voids on the fatigue life to demonstrate
Table 4 Weibull probability statistics. Material condition
Weibull slope, e
Weibull strength (cycles)
Pristine One initial flaw Initial internal void
14.36 1.54 1.88
29,794 13,140 5419
A. Warhadpande et al. / International Journal of Fatigue 32 (2010) 685–697
Fig. 22. S–N curve for crack initiation in a tensile specimen with initial internal void.
the ability of model to capture material imperfections. The presence of one internal void was found to be more detrimental than one initial flaw whereas more scatter in fatigue life data was observed in the presence of an initial flaw. Overall, the required data for the VFEM fatigue model is the material dependent parameters which are used in the damage evolution law. Once these parameters are obtained from a simple tension test, the VFEM fatigue model can then be used to study the behavior of different components subject to diverse types of loading. The model can also be used to investigate material flaws and voids and obtain a fatigue life distribution for the component. Acknowledgement The authors would like to express their deepest appreciations to the sponsors of the Mechanical Engineering Tribology Laboratory (METL), Purdue University, West Lafayette, IN for their continued interest and support of METL. References [1] Bolotin VV, Belousov IL. Early fatigue crack growth as the damage accumulation process. Probab Eng Mech 2001;60:279–87. [2] Bolotin VV, Babkin AA, Belousov IL. Probabilistic Model of Early Fatigue Crack Growth 1998;13:227–32. [3] Suresh S. Fatigue of materials. Cambridge (UK): Cambridge University Press; 1998.
697
[4] Lemaitre J. A course on damage mechanics. Berlin, Germany: Springer-Verlag; 1992. [5] Rinaldi A, Peralta P, Krajcinovic D, Lai YC. Prediction of scatter in fatigue properties using discrete damage mechanics. Int J Fatigue 2006;28:1069–80. [6] Raje N, Sadeghi F, Rateick RG. A statistical damage mechanics model for subsurface initiated spalling in rolling contacts. J Trib 2008;130:042201. [7] Raje N, Slack T, Sadeghi F. A discrete damage mechanics, model for high cycle fatigue in polycrystalline materials subject to rolling contact. Int J Fatigue 2009;31:346–60. [8] Okabe A, Boots B. Spatial tessellations: concepts and applications of Voronoi diagrams. New York: John Wiley & Sons; 1992. [9] Moller J. Lectures notes on random Voronoi tessellations. Berlin: Springer; 1994. [10] Ghosh S, Mallett RL. Voronoi cell finite elements. Comp Struct 1994;50(1994):33–46. [11] Ghosh S, Moorthy S. Elastic–plastic analysis of arbitrary heterogeneous materials with the Voronoi finite element method. Comp Methods Appl Mech Eng 1995;121:373–409. [12] Moorthy S, Ghosh S. A model for analysis of arbitrary composite, porous microstructures with Voronoi cell finite elements. Int J Numer Methods Eng 1996;39:2363–98. [13] Grujicic M, Zhang Y. Determination of effective elastic properties of functionally graded materials using Voronoi cell finite element method. Mater Sci Eng A – Struct Mater Proper Microstruct Proc 1998;251:64–76. [14] Guo R, Shi HJ, Yao ZH. Modeling of interfacial debonding crack in particle reinforced composites using Voronoi cell finite element method. Comp Mech 2003;32:52–9. [15] Vena P, Gastaldi D. A Voronoi cell finite element model for the indentation of graded ceramic composites. Composites B – Eng 2005;36:115–26. [16] Jalalahmadi B, Sadeghi F. A Voronoi finite element study of fatigue life scatter in rolling contacts. ASME J Trib 2009;131(2). [17] Kachanov LM. Time of the rupture process under creep conditions. Izv Akad Nauk SSR, Otd Tekh Nauk 1958;8:26–31. [18] Robotnov YN. Creep problems in structural mechanics. North-Holland; 1969. [19] Xiao YC, Li S, Gao Z. A continuum damage mechanics model for high cycle fatigue. Int J Fatigue 1998;20:503–8. [20] Bolotin VV. Mechanics of fatigue. Boca Raton, Florida: CRC Press; 1999. [21] Biallas G, Maier HJ. In-situ fatigue in an environmental scanning microscopepotential and current limitation. Int J Fatigue 2007;29:1413–25. [22] Wissuchek DJ, Mackin TJ, DeGraef M, Lucas GE, Evans AG. A simple method for measuring surface strains around cracks. Exp Mech 1996;36:173–9. [23] Hung PC, Voloshin AS. In-plane strain measurement by digital image correlation. J Braz Soc Mech Sci Eng 2003;25:215–21. [24] James MR, Morris WL, Cox BN. A high accuracy automated strain-field mapper. Exp Mech 1990;30:60–7. [25] Biery N, De Graef M, Pollock TM. A method for measuring microstructuralscale strains using a scanning electron microscope: application to c-titanium aluminides. Metall Mater Trans 2003;34 A:2301–13. [26] Sarıog˘lu F. The effect of tempering on susceptibility to stress corrosion cracking of AISI 4140 steel in 33% sodium hydroxide at 80 °C. Mater Sci Eng A 2001;315:98–102. [27] Andersson J. The influence of grain size variation on metal fatigue. Int J Fatigue 2005;27(8):847–52. [28] Grandt Jr AF. Fundamentals of structural integrity. NJ (USA): John Wiley & Sons, Inc.; 2004. [29] Miller KJ, Zachariah KP. Cumulative damage laws for fatigue crack initiation and stage I propagation. J Strain Anal 1977;12:262–70. [30] Lin TH, Chen QY. Micromechanics of high-cycle initiation of fatigue cracking intermetallic compounds. Acta Metall Mater 1992;40:3255–63.