Microstructure-sensitive critical plastic strain energy density criterion for fatigue life prediction across various loading regimes

In the present work, we postulate that a critical value of the stored plastic strain energy density (SPSED) is associated with fatigue failure in metals and is independent of the applied load. Unlike the classical approach of estimating the (homogenized) SPSED as the cumulative area enclosed within the macroscopic stress–strain hysteresis loops, we use crystal plasticity finite element simulations to compute the (local) SPSED at each material point within polycrystalline aggregates of a nickel-based superalloy. A Bayesian inference method is used to calibrate the critical SPSED, which is subsequently used to predict fatigue lives at nine different strain ranges, including strain ratios of 0.05 and −1, using nine statistically equivalent microstructures. For each strain range, the predicted lives from all simulated microstructures follow a lognormal distribution. Moreover, for a given strain ratio, the predicted scatter is seen to be increasing with decreasing strain amplitude; this is indicative of the scatter observed in the fatigue experiments. Finally, the lognormal mean lives at each strain range are in good agreement with the experimental evidence. Since the critical SPSED captures the experimental data with reasonable accuracy across various loading regimes, it is hypothesized to be a material property and sufficient to predict the fatigue life.

In the present work, we postulate that a critical value of the stored plastic strain energy density (SPSED) is associated with fatigue failure in metals and is independent of the applied load. Unlike the classical approach of estimating the (homogenized) SPSED as the cumulative area enclosed within the macroscopic stress-strain hysteresis loops, we use crystal plasticity finite element simulations to compute the (local) SPSED at each material point within polycrystalline aggregates of a nickel-based superalloy. A Bayesian inference method is used to calibrate the critical SPSED, which is subsequently used to predict fatigue lives at nine different strain ranges, including strain ratios of 0.05 and −1, using nine statistically equivalent microstructures. For each strain range, the predicted lives from all simulated microstructures follow a lognormal distribution. Moreover, for a given strain ratio, the predicted scatter is seen to be increasing with decreasing strain amplitude; this is indicative of the scatter observed in the fatigue experiments. Finally, the lognormal mean lives at each strain range are in good agreement with the experimental evidence. Since the critical SPSED captures the experimental data with reasonable accuracy across 2020 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/ by/4.0/, which permits unrestricted use, provided the original author and source are credited.

Introduction
Engineering components, subjected to cyclic loading, accumulate damage over the period of their operation owing to the irreversibilities associated with microscopic cyclic plasticity. Such degradation is usually referred to as fatigue and is associated with a finite safe-life operation of the component [1,2]. The estimation of the safe life, referred to as fatigue life throughout this article, is one of the critical tasks in the successful design of components. Especially for aircraft structures, where the factor of safety is almost equal to 1, accurate prediction of fatigue life is crucial.
In the past 100 years, scores of researchers offered numerous criteria, including simple phenomenological relationships and elegant physics-based models, to predict fatigue life at a specified loading condition. Basquin [3] proposed a phenomenological model relating applied stress amplitude to the fatigue life via a power-law relationship. Later, Coffin [4] and Manson [5] emphasized the importance of plastic strain, as opposed to stress, to predict life in the low cycle fatigue regime. They related plastic strain amplitude to fatigue life via a similar powerlaw relationship which is widely used in the industry to date. In his classic paper on cyclic plasticity, Morrow [6] commented that at the microscopic level neither the plastic strain (arising from the to-and-fro motion of the dislocations) nor the shear stress (resistance to the motion of the dislocations) was individually responsible for fatigue damage accumulation. He postulated that the plastic work was a composite measure of fatigue damage and thus a unified metric for life prediction. Based on the data from a series of fatigue experiments, Morrow proposed another power-law relationship relating plastic work or energy to fatigue life. Following the footsteps of Morrow, a series of researchers employed a similar energy-based approach to predict fatigue life . In each of these works, plastic work or energy was computed from the area enclosed within macroscopic stress-strain hysteresis loops.
Fatigue failure can be physically considered as a sequence of three successive events; namely, crack initiation, crack growth/propagation and the final fracture. In the second half of the last century, many research efforts were centred on understanding the physics of crack initiation. Based on extensive research of FCC materials, it was discovered that persistent slip bands (PSBs) were primarily responsible for fatigue crack initiation [32]. In 1956, Thompson et al. [33] first described PSBs as specific slip characteristics that would reappear on the surface during cyclic loading even after removing these features by electropolishing. They observed that cracks originated mostly in the interface between the PSB and the matrix. It was also found that the early stage of crack growth was strongly influenced by microstructural attributes (e.g. grain size distribution, texture, etc.) ahead of the crack tip [1,2]. Therefore, microstructural features or characteristics should undoubtedly play a key role in fatigue life prediction. However, in all of the energy-based approaches, as mentioned before, heterogeneities arising in the plastic strain energy at the microstructural length scale were not considered. Instead, a homogenized plastic strain energy density, obtained from the macroscopic stress-strain response, was regarded as a metric for life estimation.
With the advent of high-performance computing, researchers have used crystal plasticity finite element (CPFE) and molecular dynamics (MD) simulations to offer microstructure-sensitive fatigue life prediction frameworks. For example, inspired by the critical plane approach proposed by Fatemi & Socie [34], McDowell and co-workers suggested a fatigue indicative parameter accounting for microstructural heterogeneities in CPFE simulations [35][36][37][38][39][40][41][42]. Sinha & Ghosh [43] proposed a cyclic ratcheting-based model using CPFE simulation results. Further, Ghosh and coworkers used effective traction-based approaches to predict crack initiation during dwell fatigue in titanium alloys [44,45]. On the other hand, using MD simulations, Sangid et al. introduced a PSB-driven failure metric for crack initiation [46] and fatigue life prediction [47]. Building on this approach, Yeratapally et al. [48] used both MD and CPFE simulations to propose another framework combining two different length scales simultaneously (atomic and mesoscopic) and paying particular attention to the role played by twin boundaries in nickel-based superalloys. All of these models capture the physics of crack initiation and satisfactorily predict life. However, these models either involve physical parameters which are very difficult to measure from currently available experiments or include non-physical parameters which are impossible to calibrate using experiments. Thus, calibration of these models poses an unavoidable hindrance towards their applicability in industry.
Almost four decades after Morrow's classic paper on cyclic plasticity [6], researchers have returned to energy-based approaches to address microstructure-sensitive fatigue life performance. For example, building on the work by Korsunsky et al. [49], Wan et al. [50] proposed a Griffith-or Stroh-type critical stored energy density (energy per unit area) as a mesoscale driving force for crack initiation and fatigue life prediction in 2014. Using the same model, Jiang et al. [51] and Chen et al. [52] predicted multiple crack initiations around non-metallic inclusions within nickel-based superalloys, and Wilson et al. [53][54][55] predicted crack growth within polycrystalline microstructures. In a separate work, Cruzado et al. [56] proposed a similar model involving two material parameters. The microstructure-sensitive energy-based approach, introduced by Dunne and co-workers [49][50][51][52][53][54][55], is undoubtedly promising because it (i) is inherently physics based, (ii) involves only one parameter, namely critical stored energy density, (iii) takes care of microstructural heterogeneities, and (iv) offers a path forward for calibration via advanced experimental techniques. However, several aspects associated with such an approach have not been addressed so far in the literature. First, it has not been investigated whether a single critical stored energy density value can reasonably predict fatigue lives at multiple loading regimes, including low and high cycle fatigue and different stress/strain ratios. Second, experimental fatigue life data at a given stress/strain amplitude show a scatter which is typically observed to be (i) lognormally distributed [57] and (ii) increasing with decreasing applied stress/strain amplitude. A microstructure-sensitive framework should be able to capture both trends in the scatter; however, this has not been verified systematically. The primary focus of the current work is to address these two unexplored aspects.
The amount of work done by the externally applied forces during cyclic loading can be equated to the summation of the elastically stored energy, which is recoverable upon unloading, and the internal plastic work, which is non-recoverable upon unloading. The material stores a portion of the internal plastic work by forming dislocation structures and sub-structures, such as PSBs, and the rest is primarily dissipated as heat energy. In the context of fatigue failure, it is the stored portion of the internal plastic work that is of interest. Here, it is postulated that there exists a critical value of the stored plastic strain energy density (SPSED) for a material, which is associated with fatigue failure and is independent of applied loading. In the present work, based on CPFE simulations, the critical plastic strain energy density is calibrated using experimentally available fatigue life data for 718Plus, a nickel-based superalloy. Subsequently, the calibrated energy is used to predict fatigue lives at several loading regimes, including multiple applied strain amplitudes and strain ratios. Moreover, statistically equivalent microstructures (SEMs) are used to generate a scatter in the prediction to test for (i) lognormal distribution at a given applied strain and (ii) increasing spread in the prediction with decreasing applied strain. With this, the paper is structured as follows. The CPFE simulation framework including the materials processing information, the generation of SEMs and the description of the crystal plasticity model is detailed in §2. In §3, the calculation of the plastic strain energy density and the identification of the critical location for failure are discussed. A Bayesian calibration framework is outlined in §4 to estimate the critical plastic strain energy density and quantify the associated uncertainties. In §5, the results for the critical plastic strain energy density and predicted fatigue lives are reported.

Crystal plasticity finite element simulation framework
In the present section, material processing information for 718Plus is summarized. Next, the steps associated with the generation of discretized polycrystalline microstructures are outlined. Subsequently, the kinematic and kinetic equations associated with the crystal plasticity model are presented. Finally, a slip system-based averaging scheme is discussed to regularize the numerical results from the CPFE simulations.

(a) Materials processing
In this study, the material under consideration is 718Plus, a nickel-based superalloy. In particular, this material is processed via additive manufacturing, specifically selective laser melting, using an EOS M280 laser powder-bed machine at Honeywell Aerospace, Phoenix, AZ, USA. The build parameters are optimized based on process modelling, discussed in [58], to result in minimal porosity [59]. The optimized parameter values are 300 W laser power, 0.09 mm hatch spacing, 1.65 m s −1 laser speed and stripe width of 5 mm. The powder used is ATI 718Plus ® [60], which is referred to as 718Plus throughout this paper. Following the build, the material is subjected to a stress relief heat treatment at 1010°C for 1.5 h with a flue gas quench followed by a hot isostatic pressure cycle at 1163°C and 103 MPa for 4 h. Subsequently, the material is solution heat treated at 1010°C for 2 h with a flue gas quench and aged at 788°C for 8 h followed by a ramp down at 39°C h −1 to achieve a temperature of 704°C. At 704°C, the material is held for 8 h then cooled in a furnace. This heat treatment procedure has been shown to provide optimal microstructures and preferred strength characteristics [59]. The resulting microstructure, obtained via electron backscatter diffraction (EBSD), is shown in figure 1a. The bulk material is machined into specimens with a gauge length and diameter of 12.7 mm and 5 mm, respectively. Each specimen is polished to mitigate surface defects due to roughness. A total of 76 samples are tested according to ASTM E606 standards [61]. The fatigue testing is performed on smooth bars at room temperature with a constant amplitude loading frequency of 0.33 Hz. The load values are recorded and, from this information, crack initiation and failure are defined at 10% and 50% load drops, respectively, in the stabilized maximum stress values. Experiments are performed at applied strain ratios of 0.05 and −1 across a range of applied strain amplitudes necessary to achieve a complete strain-life curve, as shown in §5b. From the fractography analyses, the failure locations were typically associated with grain facets; hence, crack initiation was not observed at features such as pores or inclusions. Therefore, microstructural defects, such as pores and inclusions, are not modelled in the CPFE analysis.

(b) Generation of discretized polycrystalline microstructures
The generation of the discretized polycrystalline microstructures is a two-step process. First, a virtual three-dimensional polycrystalline microstructure is created. After that, the microstructure is discretized using a finite element meshing software. The creation of a virtual three-dimensional polycrystalline microstructure is a multi-step iterative process [62]. First, microstructural attributes (grain size distribution, twin area fraction, etc.) are obtained from the EBSD characterization of the material. Second, DREAM.3D software is used to create a synthetic volume with slightly higher grain size distribution [63]. Third, twins are inserted along one of the randomly chosen {111} planes in randomly chosen parent grains within the synthetic volume [64]. The orientation of the inserted twin is assigned to ensure 60°misorientation between the twin and the parent grain. The grain size distribution within the synthetic volume changes after the insertion of the twins. Therefore, the second, third and fourth steps are repeated iteratively until the final grain size distribution and twin area fraction from the virtual three-dimensional polycrystalline microstructure match the same obtained from the EBSD characterization. For the present work, nine virtual three-dimensional polycrystalline microstructures are synthesized (figure 1b). These are cubes of dimension 180 µm and are statistically equivalent to the microstructure features characterized in the EBSD map. Hence, these virtual three-dimensional polycrystalline microstructures are called SEMs representing the strength properties of the material and referred to as SEMs throughout this paper. Microstructural statistics of the SEMs are summarized in table 1.
After performing Laplacian smoothing of the grain boundaries, the surface mesh for each SEM is exported from DREAM.3D. Subsequently, the surface mesh is used to perform volume meshing of the SEMs using Gmsh software [65]. For the present work, computationally efficient linear tetrahedron elements (C3D4) are chosen, and the mean element size is found to be approximately 3.75 µm. The choice of the C3D4 elements often leads to numerical artefacts in the results from an elasto-plastic analysis. A possible mitigation strategy, adopted in the present work, is discussed in §2d. A discretized version of SEM 1 is shown in figure 2a.

(c) Crystal plasticity model and parameter estimation
The total deformation gradient (F) associated with the elasto-plastic deformation of a solid body can be decomposed as [66] F = F e · F p . (2.1) Here, F p and F e are deformation gradients associated with the plastic and elastic parts, respectively. Further, the plastic velocity gradient (L p ) is related to F p as Within a crystal plasticity framework, the plastic velocity gradient (L p ) is written as a sum of the contributions from all active slip systems [67]. Since 718Plus is an FCC material, the sum is performed over all 110 {111} slip systems as follows: Here,γ k , s k and n k are the shearing rate, the unit vector along the slip direction and the unit vector along the slip plane normal associated with the kth slip system. The shearing rate (γ k ) for a slip system is related to the resolved shear stress (τ k ) on the same slip system via a flow rule. In the present work, a Hutchinson-type flow rule [68] is used of the following form: Here,γ 0 is a reference shearing rate; g k and χ k are the reference and back stresses associated with the kth slip system, respectively; and n is the inverse strain rate sensitivity exponent. Both reference and back stresses evolve with plastic deformation. In the present work, the evolutions  are modelled with an Armstrong-Fredrick-type equation [69,70] and are given bẏ Here, H and A are the direct hardening coefficients for the reference and back stresses, respectively; H d and A d are the dynamic recovery coefficients for the reference and back stresses, respectively. Further, q km is the hardening matrix, where diagonal terms represent self-hardening and off-diagonal terms are associated with latent hardening (q km = 1 if k = m, otherwise q km = 1.2 [71]). All the evolution equations, as described above, are discretized and implemented within an ABAQUS user-defined material subroutine (UMAT). For simplicity, the initial values of the reference and back stresses are assumed to be the same for all slip systems. With this, there are eight crystal plasticity parameters (γ 0 , n, g(0), H, H d , χ (0), A, A d ) which are to be calibrated. Starting with the crystal plasticity parameters reported for additively manufactured IN718 in [72], a Matlab-based genetic algorithm framework, as described in [62], is used to identify an optimal set of these parameters such that the homogenized macroscopic response of the model matches the experimentally obtained macroscopic stress-strain response of the material. During simulations, normal displacements on three mutually orthogonal adjacent surfaces of the SEMs are restricted, and normal displacement is specified on another face, as shown in figure 2b, to mimic the uniaxial applied loading condition. The anisotropic elastic constants are adopted from [73]. For the crystal plasticity parameters reported in [72], the macroscopic responses of all SEMs are similar (figure 3a). Therefore, only SEM 1 is used to calibrate the crystal plasticity parameters. The anisotropic elastic constants along with the calibrated crystal plasticity parameters are reported in table 2. The simulated macroscopic stress-strain response with the calibrated crystal plasticity parameters is compared with the experimental stress-strain data for 718Plus in figure 3b. A good agreement between the simulated and experimental stress-strain response is observed in figure 3b. Uncertainty quantification for each parameter of the present crystal plasticity model is discussed in detail by the authors in [62].

(d) Slip system averaging scheme
The mechanical field variables (e.g. stress, plastic strain accumulation) from CPFE simulations often show spurious numerical oscillations within the simulation domain [72,74]. Such behaviour can be related to one or a combination of the following: (i) volumetric locking or shear locking     owing to the choice of the linear tetrahedron (C3D4) elements [74], (ii) lack of mesh refinement near the grain boundaries across which high gradients in field variables are expected, and (iii) the effect of the boundary conditions on the elements near the corresponding surfaces. Under such a situation, one needs to adopt a homogenization scheme that can regularize numerical oscillations, preserving the inherent spatial gradient associated with the field variable. One possibility is a slip system-based averaging [72,[75][76][77][78][79]. To illustrate the averaging scheme, let us consider an integration point as shown in figure 4a. Three directions are identified passing through the integration point; namely, slip direction (s), slip plane normal (n) and a transverse vector (t) on the slip plane such that these directions form a right-hand coordinate system at the integration point. Now, a small rectangular cuboid is considered, with the centroid coinciding with the integration point and the faces parallel to the three directions defined before. The region enclosed by the rectangular cuboid is called the averaging volume. The field variable is averaged over all the elements whose centroids lie within the averaging volume. One can construct 12 such volumes, one for each slip system, around the same integration point. After performing averaging for all 12 volumes, the maximum (by magnitude) of the averaged values is assigned as the slip systemaveraged number to that integration point. If an integration point is near the grain boundary, the averaging volume is not allowed to cross the grain boundary (figure 4b) to preserve the gradient, which is often expected across the grain boundary. The size of the averaging volume is dependent on the size of the finite elements. Presently, the lengths of the averaging volume are defined by 4, 3, 3 elements along the slip direction (s), slip plane normal (n) and transverse direction (t), respectively. Such an averaging scheme is implemented for all relevant field variables under present consideration during the post-processing phase of the CPFE fatigue simulations.

Critical plastic strain energy density criterion for fatigue life estimation
In crystal plasticity analysis, the plastic strain energy density at a material point x is defined as a sum of contributions from all slip systems and its increment over the yth cycle is written as Here, w p y (x) is the increment of the plastic strain energy density at the point x over the yth loading cycle. Now, the total accumulated plastic strain energy density after the yth cycle at the same point, w p y (x), is expressed as In the present work, it is postulated that, if the material fails due to fatigue at x = x * after N f cycles, w p N f (x * ) is a material property and is independent of the applied loading conditions. Further, w p N f (x * ) is referred to as the critical plastic strain energy density and is denoted as

Bayesian calibration of the critical plastic strain energy density
Variability is inherent in the experimental fatigue life data of a metallic material for a given loading condition. Such scatter in fatigue life data primarily comes from the variability in the statistical strength within the material associated with the local microstructure (e.g. grain size distributions, texture, pores, inclusions, precipitates) paired with fatigue crack initiation as a locally driven phenomenon. Therefore, in the end one has multiple W The Bayesian inference method is a powerful tool to calibrate model parameters when variability exists in the input and scatter is inevitable in the prior experimental knowledge of the output [80]. In the present work, the fatigue life prediction model is given by Here, β is a set of parameters which define the loading conditions, such as the applied strain range, strain ratio, temperature; w p N s (β, x * ), w p N s (β, x * ) and N s (β) come from the CPFE simulations; and W p critical is the only model parameter which is to be calibrated. For a given loading situation, i.e. β, variability is expected in w p N s (β, x * ) and w p N s (β, x * ) from different SEMs and scatter is experimentally observed in N expt f (β). Therefore, the Bayesian inference method is an appropriate tool to calibrate W p critical . From traditional fatigue experiments, it is well known that scatter in fatigue life is dependent on the applied strain range showing an increasing trend with decreasing applied strain range [57]. Therefore, it is more appropriate to use a subset of the experimental fatigue life dataset at a higher applied strain range to calibrate W p critical to result in less uncertainty. For the present work, the fatigue experiments are carried out from 0.3% to 1.2% applied strain range. Therefore, only the 1.2% applied strain range dataset is used for calibration, and the rest of the data are kept untouched for subsequent verification.
A broad overview of the Bayesian calibration method is shown in figure 5. Key inputs to the calibration framework are experimental knowledge (i.e. fatigue life data) and prior distribution of the model parameter (W p critical ) based on experience or prior belief. The output from the calibration framework is the posterior distribution, which is the calibrated distribution of the model parameter. Finally, the expected value of the model parameter (W model; e(β) is the error associated with the experimental measurements and uncertainty. The model discrepancy, δ(β), can be assumed to be independently and identically distributed, and hence is considered as normally distributed with zero mean and standard deviation, s 1 (unknown). Similarly, without loss of generality, e(β) can be assumed to be independently and identically distributed, and hence is considered as normally distributed with zero mean and standard deviation, s 2 (unknown). Since the summation of two normal distributions with identical mean is also a normal distribution, the model discrepancy and experimental error terms are combined to form another normal distribution with zero mean and standard deviation, s (unknown). Consequently, in addition to W p critical , the parameter s is also required to be calibrated. Therefore, the new augmented parameter set is given by α = {W p critical , s}. The additional parameter, s, is also known as a hyper-parameter. It is clarified that the fatigue life prediction model, given by equation (4.1), still involves only one parameter. However, because of the choice of the current calibration framework, the model is required to be calibrated by considering equation (4.2), which relates stochastic aspects of the model predictions and experimental observations. Consequently, an additional unknown standard deviation parameter, s, is required to be calibrated. After calibration, s has no role to play in the subsequent fatigue life prediction.
In the present work, the Bayesian calibration framework is composed of three ingredients, namely the Bayesian inference method, Markov chain Monte Carlo (MCMC) sampling and the Metropolis-Hastings algorithm. A brief description of each of these ingredients and their roles in the Bayesian calibration framework are provided below.
Here, π 0 (α) represents the prior distribution of the parameter set α and π (D|α) represents the likelihood of observing experimental data, D, for the given parameter set α; and π (α|D) is the posterior distribution of the parameter set α based on experimental evidence D and likelihood term π (D|α). Since the error terms, as defined in equation (4.2), are assumed to be independently and identically distributed, the likelihood term is expressed as [80] π (D|α) = 1 Here, s is the hyper-parameter as defined before; SSE(α) is the sum of squares of errors for the parameter set α and is expressed as (4.5) In equation (4.5), N of random variables that satisfies the Markov property that the kth term in the sequence depends only on the (k − 1)th term. The kth term is proposed via the Monte Carlo sampling method, which is based on a proposal distribution. The decision whether the new proposed candidate point will be accepted as the kth term in the sequence is taken based on the Metropolis-Hastings algorithm [81,82]. Stepwise computer implementation of the algorithm has been discussed by Yeratapally et al. [83] in relation to the uncertainty quantification of a microstructure-sensitive fatigue life prediction model. However, for the sake of completeness and continuity, important aspects and equations of the Metropolis-Hastings algorithm are outlined below. Let, α k−1 j be the jth parameter (j = 1, 2) in the (k − 1)th term of the Markov chain and let α * j be the proposed jth parameter for the kth term in the sequence. Now, the following number is computed: Here, π (α * j |α k−1 j ) represents the proposal distribution used to generate the proposed candidate point α * j based on the current candidate point α k−1 j . Here, the proposal distribution is assumed to be symmetric, for which Therefore, equation (4.7) can be simplified and rewritten as follows: (4.9) The assumption of a symmetric proposal distribution is a matter of convenience and simplicity. If the posterior distribution is expected to be symmetric, such a choice leads to quick convergence (as discussed later). Otherwise, it may take a longer time to obtain a converged posterior distribution. After obtaining p from equation (4.9), a uniformly distributed random number u is generated within the interval 0 ≤ u ≤ 1. If u ≤ p, the proposed candidate is selected as the kth term in the sequence, i.e. α k j = α * j . Otherwise, the (k − 1)th term is assigned as the kth term, i.e. α k j = α k−1 j . This process is continued until convergence is achieved. To check for convergence, M Markov chains are run in parallel. If B α j and W α j are the variance of the parameter α j between and within chains, respectively, then an estimate of the variance of the parameter α j , V α j , is given by Using V α j and W α j a convergence parameter R α j is defined as With the increase in the number of Monte Carlo sampling, k, in each Markov chain, the parameter R α j approaches unity, resulting in a converged/stationary distribution. Such a converged/stationary distribution is the posterior distribution. The reader is referred to Cross et al. [84] for the expressions of W α j and B α j .    value, as reported in table 3, to predict fatigue lives at each strain range from all SEMs. The probability of failure plots is obtained using the probplot function in Matlab by considering the lognormal distribution as a reference, as shown in figure 8. In figure 8, the vertical axes are depicted in such a way that, if data points are sampled from lognormal distributions, they will be aligned along a straight line. As expected, experimental data points in each plot appear to be aligned along a straight line. Interestingly, all the predicted data points also appear to be aligned along a straight line. Therefore, it can be said that the microstructure-sensitive critical plastic strain energy density criterion consistently predicts lognormally distributed fatigue life across several loading regimes. In each case, predicted lives from all SEMs are fitted to a lognormal distribution. Experimental data, predictions from all SEMs and the lognormal mean of predictions are shown in figure 9. A good agreement between the experimental and predicted lives is seen in figures 8

Discussion
The fatigue life prediction framework, presented in this work, inherently captures the relevant physics or mechanism of the plastic deformation via crystal plasticity constitutive equations. Hence, it can capture the irreversibilities associated with the deformation even if the material is loaded elastically in a macroscopic sense. In other words, the framework can be used for both high cycle fatigue, dominated by microplasticity, and low cycle fatigue, dominated by a large amount of plastic deformation, situations. In this work, dislocation glide is the primary deformation mechanism which has been incorporated in the crystal plasticity framework via equation (2.3).
In other material systems, additional deformation mechanisms, such as deformation twinning, might be dominant as well. In those situations, kinematic equations (equations (2.1)-(2.3)) for the crystal plasticity model should be updated appropriately before using equation (3.1). Further, the present approach can, in principle, be extended to situations where the environment plays an important role (e.g. hydrogen embrittlement, corrosion) in the fatigue performance of the material via incorporating suitable constitutive relations in the crystal plasticity model. The reader is reminded that the CPFE simulation results, i.e. w p 10 (x * ) and w p 10 (x * ), are key inputs to the fatigue life prediction model given by equation (4.1). Therefore, the magnitude of the calibrated W p critical strongly depends on the size of the slip system averaging volume, as discussed in §2c. Since the size of the slip system averaging volume has been consistently kept the same throughout the present work, its effect will not be reflected in the fatigue life predictions. However, the reported value for W p critical should not be considered as a benchmarked number for 718Plus. For more rigorous estimation of the W p critical , the following should be considered. (a) Nonlinear finite elements should be used to avoid volumetric/shear locking and thus eliminating the necessity of the slip system-based averaging. (b) The mesh near the grain boundary should be sufficiently refined to accurately capture the gradients associated with the mechanical field variables in that region.
The fatigue life was experimentally obtained using specimens with a gauge length and diameter of 12.7 mm and 5 mm, respectively, whereas the domain size for the CPFE simulations was restricted to cubes of 180 µm length owing to the higher computational cost. Since simulation domain size was kept identical for calibration of W p critical as well as prediction of fatigue life, its effect is not seen in figures 8 and 9. In the present work, the uncertainty in W p critical as a result of the size of the SEM is implicitly taken care of via the hyper-parameter through equation (4.2). In [48], it is shown that, if the number of grains within the SEM is around or greater than 150, convergence is achieved in the strength properties of the material (in terms of capturing the macroscopic yield, hardening behaviour, Bauschinger effect, etc.). In the present work, the consistent size of the SEMs is necessary to capture the effect of microstructure variability in the fatigue behaviour and thereby the fatigue performance, but each SEM is not designed to predict the statistical minimum fatigue life. The appropriate size of the microstructure instantiation, at the mesoscopic length scale, should be identified more rigorously for fatigue predictions at the length scale of a component, which is a topic for future work.
When the macroscopic stress-strain data are used to calibrate the crystal plasticity parameters, as in the present work, the resulting crystal plasticity parameters may not be unique, which would lead to additional uncertainty in the values of w p 10 (x * ) and w p 10 (x * ). Detailed uncertainty quantification of the crystal plasticity parameters and their influence on the associated fatigue metrics (e.g. the SPSED) have been performed by the authors, as reported in [62]. The results indicate that the uncertainty in the crystal plasticity parameters may propagate to approximately 25% variability in the values of w p 10 (x * ) and w p 10 (x * ) [62]. In the present context, it would imply that the calibrated distribution of W p critical might result in a slightly higher standard deviation, which would not affect the mean of W p critical and hence the predicted fatigue life.
As discussed in §1, the material stores a portion of its internal plastic work by forming dislocation structures and sub-structures, and the rest is primarily dissipated as heat energy. The stored portion of the internal plastic work is of interest in the context of fatigue failure. It is often argued that 5% of the total plastic work is stored in the material and the remaining 95% is lost as heat energy [50]. However, the number 95% is not universally accepted (see, for example, [85,86]); in some cases, it is found to be increasing with increasing plastic strain [85]. Hence, there exists a critical need to carry out rigorous studies to characterize the percentage of plastic work lost as heat energy and, in turn, the percentage of plastic work stored within the material as a function of applied strain and temperature. In the absence of a lack of agreement between the reported values of the partitioning coefficient in the literature, researchers often neglect the partitioning effect in their numerical calculations (see, for example, ). With this, while computing W p critical , we have not considered the partitioning of the plastic work. The entire plastic work has been assumed to be stored by the material. In the present context, such an assumption only affects the calibrated value of the W p critical , not the predicted fatigue lives. In other words, a more realistic calibration of W p critical would require incorporation of a scalar constant in equations (3.1)-(4.1) to incorporate the partitioning aspect, which would change only the calibrated numerical value of W p critical but not the subsequent fatigue life prediction.
One can observe from figure 9 that the predicted life at lower strain amplitude, especially at R ε = 0.05, spans from approximately 10 5 cycles to approximately 10 6 cycles. Therefore, if experimental data at these strain amplitudes are used to calibrate W p critical , the resulting posterior distributions will involve a higher level of uncertainty, characterized by a higher standard deviation in W p critical and higher mean and standard deviation for the hyper-parameter s. In other words, just like the model discrepancy δ(β) and experimental uncertainty e(β) terms in equation  (e) ( f ) (g) (h) (i) Figure 10. The location of maximum w p 10 within SEM 1 at (a) ε = 1.2%, R ε = −1, (b) ε = 0.95%, R ε = −1, (c) ε = 0.75%, R ε = −1, (d) ε = 0.675%, R ε = −1, (e) ε = 0.625%, R ε = −1, (f ) ε = 0.5%, R ε = −1, (g) ε = 0.5%, R ε = 0.05, (h) ε = 0.43%, R ε = 0.05, and (i) ε = 0.36%, R ε = 0.05. The yellow arrows (arrows in the top row of images) indicate the critical location near or at the free surface and the white arrows (arrows in the middle and bottom rows of images) indicate hotspots at the twin boundaries. Each colour map is a cross-sectional view of the SEM 1 on the x-z plane, normal to the loading direction y. (Online version in colour.) (4.2), the hyper-parameter s is also dependent on the loading conditions defined by β. Since the life predictions from different SEMs involve least scatter at the 1.2% applied strain range, the experimental data corresponding to the same loading regime are used in the present work to calibrate W p critical . If the entire experimental and simulation datasets for R ε = −1 are used for the calibration, 2 the resulting mean and standard deviation of W p critical are found to be approximately 4% lower and approximately 4 times higher, respectively, whereas the order of magnitude of the hyper-parameter is found to be approximately 100 times higher than the values reported in table 3. These observations further support the foregoing discussion.
It has already been mentioned in §5b that an increasing spread in the fatigue life predictions, from an identical set of SEMs, is observed with a decreasing applied strain. Such progressive uncertainty in predictions, from the same set of SEMs, with decreasing applied strain is certainly a manifestation of the change in x * , i.e. failure location, with applied loading situation. To confirm, x * is tracked within SEM 1 with the change in applied loading. A shift of the failure location from the free surface to an interior twin boundary is observed with the decrease in applied strain (figure 10). At ε = 1.2%, 0.95% and 0.675%, the failure location is observed near or at the free surface ( figure 10a-c). Subsequently, at lower ε, the failure location is identified at interior twin boundaries (figure 10d-i). The ability of the microstructure-sensitive plastic strain energy density criterion to capture (i) the shift in failure location, often seen in real experiments, and (ii) real fatigue life data with reasonable accuracy across various loading regimes certainly demands attention of the industry as well as academia to reduce the extensive time and cost associated with experimental fatigue life characterization.
It is well known that properties associated with the mechanical behaviour of materials are inherently microstructure dependent. In the present work, we have hypothesized the critical plastic strain energy density as a material property and have characterized its statistical distribution. While predicting fatigue lives, we have used the mean value, which is applicable for the current microstructure. By altering the statistical descriptors of the microstructure, a different mean value of the critical plastic strain energy density is expected. However, such dissimilarity has to be verified in future work. Further, for a given statistical description of the microstructure, the applicability of a single critical plastic strain energy density is also required to be tested and verified in (i) the presence of defects, such as pores and inclusions, and (ii) other loading conditions, including the role of environment, various operating temperatures and multi-axial applied loading.

Conclusion
A microstructure-sensitive critical plastic strain energy density is postulated to be the driving mechanism of fatigue crack initiation and is applicable to predict failure across several loading regimes. The model is inherently physics based, and requires only one parameter, namely the critical plastic strain energy density. This critical plastic strain energy density is calibrated using experimental fatigue life data under fully reversed type loading at 1.2% applied strain range via the Bayesian inference method. Subsequently, the calibrated critical energy value is used to predict fatigue lives at eight strain ranges including strain ratios −1 and 0.05. The key findings are summarized below.
-Predicted lives from nine statistically equivalent microstructures for each applied strain range are found to follow a lognormal distribution which is consistent with observations in traditional fatigue experiments. -A good agreement is observed between the experimental fatigue life and the lognormal mean of the predictions at eight strain ranges including strain ratios −1 and 0.05. -For a given strain ratio, increasing spread in the predictions from the identical set of statistically equivalent microstructures is observed with decreasing applied strain amplitude.