A Probabilistic Fatigue Framework to Enable Location-Specific Lifing for Critical Thermo-mechanical Engineering Applications

The present paper describes a probabilistic framework to predict the fatigue life and failure mode under various thermo-mechanical loading conditions. Specifically, inclusion- and matrix-driven competing failure modes are examined within nickel-based superalloys. The critical accumulated plastic strain energy density (APSED) is employed as a unified metric to predict fatigue crack initiation in metals, which is favorable due to the usage of a single unknown parameter and its capability to predict failure across loading conditions and failure modes. In this research, we characterize the temperature-dependent variation of the critical APSED using a Bayesian inference framework and predict the competing failure modes in a coarse grain variant of RR1000 with varying strain range and temperature. The critical APSED appears to decrease along a vertically reflected sigmoidal curve with increasing temperature. Further, (a) the prediction of a failure mode, (b) failure mode associated with the minimum life, and (c) the change in the location associated with the matrix-driven failure mode with increasing temperature and decreasing strain range are consistent with the experimentally observed trends in RR1000, as well as other Nickel-based superalloys, documented in the literature. Finally, for each simulated loading condition, the uncertainty in the fatigue life is quantified as a prediction interval computed based on a 95%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$95\%$$\end{document} confidence level of the critical APSED and the computed APSED from simulations. The overall framework provides a promising step towards microstructural-based fatigue life determination of components and enables a location-specific lifing approach.


Background
The fatigue performance of metallic materials is closely tied to their underlying microstructure, which in turn is dependent on the material processing route. Nickel (Ni)based superalloys, being an essential material for the aerospace propulsion industry, have been the focus of extensive research and development for the past several decades to seek new processing routes to achieve tailored microstructure and hence improved mechanical performance [1][2][3]. As an example, one may consider turbine blades. In the 1940s, wrought materials were used for such applications, whereas later, cast materials were preferred due to their better creep performance. The turbine blade materials went through an evolution from the equiaxed to columnar (directionally solidified) to single-crystal microstructure, enabling modern turbine blades to withstand almost 60% higher temperature [1]. In turbine disks, a goal is to tailor the microstructure to achieve targeted location-specific properties and fatigue performance. Powder metallurgy (PM) processing routes, along with dual microstructure heat treatment (DMHT) strategy, are used to produce such a tailored microstructure in turbine disks [4,5]. DMHT processing can result in a fine grain microstructure near the disk's bore region, whereas a coarse grain microstructure is formed near the rim region. Typically, a fine grain microstructure leads to higher strength and fatigue life, whereas a coarse grain microstructure offers better creep resistance. From the DMHT turbine disk example, it is evident that a shift in the lifing methodology from classical approaches of lifing the entire component, as a monolithic structure, to new strategies of lifing the 1 3 component at specific locations, based on the microstructure, are necessary, in order to achieve targeted location-specific performance [6]. With this in mind, the present work's objective is to predict fatigue life and associated mechanisms in a PM produced Ni-base superalloy, RR1000, at various mechanical and thermal loading conditions using a microstructure-sensitive energy-based approach. We start with a brief review of the existing literature on the microstructuresensitive fatigue life prediction and discuss the advantages of the energy-based approaches.
The energy-based approaches consider the fact that the motion of dislocation, i.e., plastic strain, and resistance to the motion of dislocation, i.e., shear stress, are not individually responsible for fatigue damage. It is the combination of these two local quantities, i.e., plastic work or energy, which is the driver for fatigue crack nucleation [35]. Although a significant portion of the plastic work is lost, primarily, as heat energy, the remaining portion is stored by dislocation structures and substructures. This stored energy is associated with the fatigue crack nucleation in metals. Dunne and colleagues [22][23][24][25][26][27][28][29][30] postulate and rationalize that the magnitude of the stored energy, as well as the length-scale over which it is stored, are collectively responsible for fatigue crack initiation. Accordingly, they have used a Griffith-or Stroh-type critical stored energy density (energy per unit area) to predict fatigue life under uniaxial [22] and multiaxial [23] loading conditions, to predict crack initiation around non-metallic inclusions [24][25][26], to predict crack growth within polycrystalline aggregates [27][28][29], and to assess fatigue performance under non-proportional loading situations [30]. Cruzado et al. [31,32], Bandyopadhyay et al. [33], and Prithivirajan and Sangid [34] have used a volumetric critical accumulated plastic strain energy density (APSED), as opposed to energy per unit area used by Dunne et al. [22][23][24][25][26][27][28][29][30], to predict cycles to crack initiation. The researchers, in Refs. [31][32][33][34], have taken into account the length-scale effect implicitly via nonlocal averaging of the APSED values. In addition to the critical APSED, Cruzado et al. [31,32] used an additional parameter to predict fatigue life (e.g., an exponent to account for nonlinear behavior). Bandyopadhyay et al. [33] have shown that a single critical APSED value is sufficient to predict the median fatigue life and associated scatter successfully at multiple applied strain ranges and strain ratios under uniaxial loading conditions. Additionally, Prithivirajan and Sangid [34] have demonstrated that the critical APSED is well-suited to predict matrix-and pore-driven competing failure modes in additively manufactured IN718. From Refs. [23,33], it is at least indicative, if not apparent, that the critical energy value, associated with crack initiation, for a given material is constant at a given temperature. Bandyopadhyay et al. [33] used a Bayesian framework that considers the uncertainty in the simulation and error associated with experimental fatigue life measurement to identify the critical APSED value. Thus, the energy-based fatigue life prediction approaches in Refs. [21][22][23][24][25][26][27][28][29][30][31][32][33][34], along with the calibration strategy proposed by the present authors [33], offer a path forward for validation, verification, and uncertainty quantification, which are crucial for the qualification of a new/improved material and in general, building trust in the microstructure-sensitive fatigue life prediction.
In this paper, we extend the work in [33]. First, using crystal plasticity finite element (CPFE) simulations, experimental fatigue life, and the Bayesian framework, we calibrate the critical value of APSED for RR1000 as a function of temperature. We also quantify the uncertainty in the critical APSED value originating due to the imprecise calibration of the crystal plasticity (CP) parameters and variability associated with the experimental measurement of fatigue life. Second, we use the calibrated APSED values to predict fatigue life and associated failure modes, namely inclusionand matrix-driven failures, in RR1000 as a function of strain range and temperature. Third, based on the calibrated uncertainty in the critical APSED values, we quantify uncertainty in the predicted fatigue life. Prediction of the competing failure modes, namely inclusion-driven failure (IDF) and matrix-driven failure (MDF), as a function of strain range and temperature, requires careful consideration of microstructure and defects structure that influence the emergence of these failure modes. Therefore, we briefly discuss these parameters and their implications on the competing failure modes.

Parameters Influencing Competing Failure Modes in Ni-Base Superalloys
Many researchers have experimentally investigated crack initiation from inclusions and competing failure modes in Ni-base superalloys. For example, Hyzak and Bernstein [36], and Huron and Roth [37] reported a higher percentage of IDF at a higher temperature in AF-115 and AF2-1DA, and René 88 DT, respectively. However, counter-evidence was also reported by Bhowal and Wusatowska-Sarnek [38] in their notched low-cycle fatigue experiments with IN718. Additionally, Hyzak and Bernstein [36] observed a change in the location of crack initiation from the surface to the interior of the specimens with decreasing strain range. Caton et al. [39], Qiu and Wu [40], and Texier et al. [41,42] reported that surface and sub-surface inclusions were mostly associated with IDF in René 88 DT, RR1000, and IN718DA, respectively. Texier et al. [42] also noted that all surface inclusions in IN718DA were not necessarily associated with crack initiation. Caton et al. [39] observed that inclusions did not result in significant debit in life for stress amplitude greater than 1000 MPa, whereas the debit was significantly higher at lower stress amplitudes leading to a bifurcation in the stress-life plot. However, Shi et al. [43] showed counter-evidence in their research with FGH96. In Refs. [40,41], one may notice a single isolated particle, i.e., Type 1, or an agglomerate of small particles, i.e., Type 2, can be associated with IDF. For example, Qiu and Wu [40] observed IDF from Type 1 alumina inclusions and Type 2 hafnia inclusions in RR1000, whereas Texier et al. [41] reported IDF due to Type 1 and 2 titanium nitride inclusions and Type 2 niobium carbide inclusions in IN718DA. Alexandre et al. [44] showed that particle-induced failure is more prevalent below a critical grain size of IN718. Based on an extensive experimental investigation with René 88 DT, Stinville et al. [45,46] reported that surface crack initiation due to stress-assisted oxidation was dominant in the low-cycle fatigue regime. In the high-cycle fatigue regime, internal crack initiation, due to inclusions and crystallographic facets, was observed, leading to almost 51% MDF and 49% IDF. They noticed failure due to crystallographic facets in the very high-cycle fatigue regime was more dominant, resulting in 95% MDF and 5% IDF. Zhang et al. [47][48][49], Naragani et al. [50], Jiang et al. [24,51], and Chen et al. [25] carried out coupled experimental and numerical modeling work to investigate the micromechanical basis of crack initiation in RR1000 [47][48][49][50] and FGH96 [24,25,51] due to the presence of inclusions. From Refs. [24,25,[47][48][49][50][51], we observe that significant strain-gradient, higher stored energy density, debonding of the inclusion-matrix interface, and residual stress around the inclusion are critical factors that influence slip-mediated IDF in these materials.
From the literature review, as presented above, it is evident that a multitude of parameters influences the emergence of IDF and MDF in Ni-base superalloys. We can summarize these parameters as (a) the type of the inclusion (Type 1 or 2), (b) the fracture strength of the inclusion particles, (c) the inclusion-matrix interfacial strengths, (d) the presence of a partially debonded region in the inclusionmatrix interface due to the thermo-mechanical processing of the material, (e) microstructural variability around the inclusions in terms of grain size and orientation distributions, (f) dissimilar coefficient of thermal expansion (CTE) of the inclusion and the matrix material, (g) distance of the inclusion from the free surface, (h) residual stress in the material due to the thermo-mechanical processing, and (i) the operating loading conditions including stress/ strain range, stress/strain ratio, and temperature. The type of the inclusion, i.e., Type 1 or 2, determines the effective load it carries, thereby its structural stiffness. The fracture strength of the inclusion particles and the inclusion-matrix interfacial strength determine if an IDF would originate due to particle cracking or interfacial decohesion, respectively. The remaining parameters, including inclusion type, determine if crack initiation would occur near the inclusion or at a location sufficiently far from the inclusion due to slip accumulation in the matrix material. Accordingly, one would observe either an IDF or MDF, respectively.
Recently, Bandyopadhyay and Sangid [52] and Bergsmo and Dunne [26] have investigated the effects of the abovelisted parameters on the emergence of these two failure modes. Specifically, Bergsmo and Dunne [26] reported that the oxide particles' interfacial strength in RR1000 was 2050 MPa, and the fracture strength of the same particles was 2300 MPa. They showed that the particle-cracking and interfacial decohesion would be more prevalent above 1230 MPa applied peak stress, whereas slip-induced IDFs or MDFs would be more dominant at lower stress levels. In IN718DA, Texier et al. [53] also reported particle cracking occurring at similar macroscopic stress levels (somewhere between 1147 MPa and 1292 MPa) under monotonic loading. For slip-induced competing failure modes, the parametric study by Bandyopadhyay and Sangid [52] shows that none of the parameters, as mentioned above, are individually responsible for the emergence of the either of these failure modes if the inclusion size is on the order of the grain size; it is the complex interaction between these parameters, that leads to either of these failure modes. The results of [52] suggest an inherent stochastic aspect in the emergence of these failure modes, which is supported by experimental evidence in Refs. [24,25,[36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][51]. Realistic prediction of these two competing failure modes requires careful modeling of these parameters' effect in the simulation framework, as discussed in section "Polycrystalline Microstructures with Inclusions and Debonded Region".
So far, we have motivated the need for location-specific lifing (section "Background"), described the advantages of the microstructure-sensitive energy-based approaches for location-specific lifing (section "Microstructure-Sensitive Fatigue Life Prediction"), and reviewed parameters influencing the competing failure modes in Ni-base superalloys (section "Parameters Influencing Competing Failure Modesin Ni-Base Superalloys"). Now, we summarize the location-specific lifing framework for the present work and layout the paper's structure.

Summary of the Location-Specific Lifing Framework
For illustration, consider the schematic of a turbine disk in Fig. 1. We choose a location A in this component and demonstrate how the probability of failure plots and strainlife plots can be obtained for location A. The steps are as follows.  Fig. 1, we demonstrate how to obtain the probability of failure plots using steps (i)-(v). Strain-life plots can also be generated using the same steps. For example, simulations carried out at a constant strain ratio, temperature, and varying strain ranges lead to a strain-life plot. Hence, simulations can be carried out at multiple strain ratios and temperatures to generate strain-life plots for a spectrum of loading conditions.
In section "Crystal Plasticity Finite Element Simulation Framework", steps (i)-(iii) are described in detail. The fatigue life model, i.e., steps (iv) and (v), is described in section "Fatigue Life Prediction Model". In section "Bayesian Calibration of the Critical Plastic Strain Energy Density", the Bayesian calibration framework of the fatigue life model parameter, namely critical APSED, is outlined. The results are reported and discussed in section "Results". Finally, conclusions are drawn in section "Conclusions".

Material Characterization
A PM produced Ni-base superalloy, RR1000, is the focus of this work. Specifically, we consider the coarse grain variant of the material, which is supersolvus heat-treated and, hence, does not contain the primary ′ phase. The processing of this material is described in Ref. [54]. A representative electron backscatter diffraction (EBSD) characterization of the material is shown in Fig. 2a. From the EBSD characterization, the mean grain size is 25.24 μm, and a random texture is observed in the material (see Fig. 2b).
Micro-computed tomography (µ-CT) technique is used to characterize the inclusions in this material. Specifically, inclusions in five RR1000 samples, obtained from work in Ref. [55], are characterized using µ-CT. During µ-CT, a 433 μm × 300 μm × 1250 μm rectangular cuboid volume of the material is probed with the resolution, i.e., pixel size, of 1.5 μm . Afterward, Avizo software is used to process the raw data to obtain equivalent diameter and sphericity of the inclusions (see Fig. 2c, d). The sphericity [56] is computed as where A o and V o are the surface area and volume of the inclusion, respectively, such that the sphericity represents a measure of the deviation in shape from a sphere (i.e., for a sphere, sphericity = 1). The material characterization results, as shown in Fig. 2, are used to obtain 3D synthetic polycrystalline microstructures of the material, as described below.

Polycrystalline Microstructures with Inclusions and Debonded Region
A flowchart, outlining the procedure to obtain a 3D synthetic polycrystalline microstructure with an embedded inclusion, local mesh refinement, and a debonded region, is shown in Fig. 3. The process can be broadly divided into two steps: (a) creation of a 3D synthetic polycrystalline microstructure without an inclusion using statistical descriptors (grain size,

Fig. 1
Schematic of location-specific lifing of an aerospace component shape, and orientation distribution, along with the twin area fraction) of the microstructure from EBSD characterization, and (b) insertion of inclusion within these 3D synthetic polycrystalline aggregates using µ-CT characterization data of the inclusion population. Each of these steps is discussed below.
Equiaxed grains are observed in RR1000. For such a case, a framework to obtain 3D synthetic polycrystalline microstructures, starting with the 2D EBSD characterization of the material, has been previously described by the present authors [20,34,52,[57][58][59]. These 3D synthetic polycrystalline microstructures are statistically equivalent to the characterized microstructural features and are sufficiently large to reproduce the microstructural strength characteristics of the material. Hence, these 3D synthetic polycrystalline microstructures are called statistically equivalent microstructures and referred to as SEMs in the remainder of this paper. Previously, Bandyopadhyay and Sangid [52] described the procedure to geometrically model an inclusion within an SEM as a solid ellipsoid. Specifically, an SEM is first discretized using finite elements (the mesh details are described in section "Discretization of the Statistically Equivalent Microstructures"). Subsequently, a set of elements is defined as an inclusion.
In order to predict fatigue life and associated failure mode, we need to create SEMs with realistic attributes, including (a) the number and size of the SEMs, (b) the mesh size within the discretized SEMs, and the degree of mesh refinement around the inclusion, (c) thermo-mechanical properties of the inclusion, (d) size and morphology of the inclusions, and (e) the position of the inclusions within the SEMs. Each of these items is described below.

Number and Size of Statistically Equivalent Microstructures
Ideally, a large number of SEMs is desirable to predict a stationary distribution of fatigue life. From a component design perspective, engineers are often interested in the B0.1 life or fatigue life associated with a 0.001 probability of failure [60]. The number of SEMs, necessary to obtain a stationary B0.1 life, is still an outstanding question and is not the focus of the present work. Here, we extend the framework this material, c equivalent diameter, and d sphericity obtained from the µ-CT characterization of five RR1000 specimens in Ref. [55] described by Bandyopadhyay et al. [33] to predict the competing failure modes in RR1000 over a range of operating temperatures and applied strain range. With this, a total of 20 SEMs are created, which is consistent with previous microstructure-sensitive fatigue life predictions using SEMs in Refs. [20,31,33,34]. In previous experimental research with more than 500 specimens of René 88 DT, Stinville et al. [45,46] observed crack initiation from alumina and hafnia inclusions in nearly 49% of the specimens. In the present work, 10 SEMs (out of a total of 20 SEMs) are created with inclusions, with each SEM containing one inclusion. With this, a maximum of 50% IDF can be predicted from this study.
CPFE simulations, carried out at the specimen lengthscale, would be ideal for comparative assessment between simulated and experimental fatigue life. However, such a task is impractical because of computationally expensive CPFE simulations. Recently, Prithivirajan and Sangid [34] conducted a sensitivity analysis for the SEM size compared to the full gauge volume of the microtensile fatigue specimens. Their results suggest an SEM containing ∼ 200 grains is reasonable for fatigue life prediction. Given the grain size distribution of RR1000, cubes with dimensions of 160 μm contain ∼ 200 grains and are therefore used as SEMs in this work.

Discretization of the Statistically Equivalent Microstructures
In microstructure-sensitive fatigue life predictions, it is essential to preserve the morphology of defects, such as the inclusions, during meshing to ensure the gradients in the micromechanical field variables can be captured faithfully. Tetrahedral elements in 3D are appropriate for capturing microstructural morphologies with fewer nodes than the hexahedron elements. Hence, in this work, linear tetrahedron or C3D4 elements are chosen to discretize the SEMs.
Due to numerical issues, many researchers prefer nonlinear elements over C3D4 elements. However, such a choice is infeasible in the present work involving 240 CPFE simulations (see section "Results") with each of the SEMs having ∼ 200 grains. Such a large number of 3D simulations with nonlinear elements is currently impractical on supercomputer clusters in a reasonable time. Moreover, many finite element solvers cannot scale down simulation time appropriately (due to non-optimal parallelization schemes) with the increase in the number of computing nodes. Hence, the choice of the C3D4 element is a trade-off between the accurate geometric representation of the polycrystalline aggregates, computational expenses, and numerical accuracy. There exist multiple ways to tackle the numerical problems that might arise due to the choice of the C3D4 elements. In this work, we adopt a regularization scheme in the post-processing phase of the analysis described in section "Nonlocal Averaging Scheme". Based on the mesh sensitivity study performed by Prithivirajan and Sangid [58], the element size is chosen as ∼ 3 μm . Further, to faithfully capture gradients in the mechanical field variables across the inclusion-matrix interfaces, the mesh is refined around the inclusion. Specifically, a spherical region with its center coinciding with the centroid of the inclusion and diameter approximately twice the equivalent diameter of the inclusion is identified and discretized using ∼ 1.2 μm C3D4 elements.

Thermo-Mechanical Properties of Inclusions
In this work, the inclusions are modeled as homogeneous isotropic linear elastic materials. Such an idealization is consistent with the fact that the non-metallic inclusions are brittle and often polycrystalline. Moreover, these inclusions are often pulverized into granules [61]. Hence, as explained below, we adopt a structural stiffness description of the homogenous inclusion. With this, the thermo-mechanical properties that are of interest are structural stiffness and CTE.
Depending on the chemistry (e.g., oxides, carbides) and type (Type 1 or 2), structural stiffness of the inclusions is expected to span a wide range resulting in lower or higher stiffness compared to the matrix. Here, the matrix's stiffness corresponds to Young's modulus 1 of RR1000; however, the inclusion's structural stiffness is different from its Young's modulus. For example, a fully dense alumina particle has a well-defined Young's modulus value. However, depending on the type of the alumina inclusion, i.e., Type 1 or 2, the structural stiffness, associated with the load carried by the inclusion, would be lower or higher than its Young's modulus. Prior work by the authors [52] showed that for cases when the inclusion's structural stiffness is higher than the Young's modulus of the matrix, IDFs are more likely, irrespective of the applied strain range. However, if the inclusion's structural stiffness is lower than Young's modulus of the matrix, a higher probability of IDFs is observed at a lower applied strain range compared to the higher applied strain ranges. The fatigue life dataset of RR1000 in this research shows IDF evidence for applied strain range ≤ 0.4% . Hence, in this study, inclusions are predominantly modeled as structurally soft materials. However, to avoid any bias towards the stiffness value choice, the inclusion stiffness values are sampled from the shaded region of the probability distribution shown in Fig. 4a. It is to be noted that the probability distribution in Fig. 4a is hypothetical and does not represent measured structural stiffness values of inclusions in RR1000. However, the lower and upper bounds of the shaded region are motivated by the fact that the Young's modulus of porous alumina can lie somewhere between 150 and 400 GPa , depending on the pore volume fraction [62].
Due to the varying chemical composition and crystallographic structure, one may expect a wide range of CTE values for the inclusions as well, resulting in a relatively lower or higher CTE than the matrix material. In previous research with RR1000 [40] and similar PM Ni-base superalloys [45,46], IDFs have been reported to be due to alumina and hafnia inclusions. Additionally, from chemical analysis, it has been observed that the hafnia inclusions are surrounded by alumina particles [47]. Although the CTE of bulk alumina and hafnia is lower than RR1000, depending on the monoclinic or tetragonal structure, the anisotropic CTE of hafnia [63] could be higher than the CTE of RR1000 [64]. With this, the inclusion CTEs in this work are primarily considered lower than the matrix material. However, to consider possible variabilities and avoid any bias, the CTE of the inclusions in this research is sampled from the shaded region of the probability distribution shown in Fig. 4b. The probability distribution in Fig. 4b is hypothetical and does not represent the measured CTE of inclusions in RR1000. However, in the shaded region in Fig. 4b, the lower bound is based on the CTE of isotropic polycrystalline alumina and hafnia, whereas the upper bound is based on the anisotropic CTE of hafnia reported in Ref. [63].

Inclusion Size and Morphology
Based on prior work by the authors [52], the damage associated with the presence of an inclusion during cyclic loading is not directly correlated with the size and curvature of the inclusion, if the inclusion size is of the order of the mean grain size. As seen in section "Material Characterization", in the present work, the inclusion equivalent diameter is of the order of the mean grain size. To avoid any bias, the inclusion equivalent diameter and sphericity are sampled from uniform distributions in Fig. 4c,   The Young's modulus and CTE values of RR1000 are indicated for reference and obtained from Refs. [64,65] at room temperature, respectively the interval [0.5, 0.95] . Such a choice, especially on the lower bound of the inclusion equivalent diameter, is due to the mesh size, which is ∼ 1.2 μm surrounding the inclusion.
The presence of a debonded region may not necessarily lead to an IDF [52]. However, in previous research with RR1000, IDFs have been reported to be due to the cluster of alumina and hafnia particles, i.e., Type 2 inclusions [40]. In the present experimental dataset, all IDFs originate from Type 2 inclusions, i.e., structurally soft inclusions, as well. Therefore, in this study, we introduce a partially debonded region in the inclusion-matrix interface to implicitly capture some of the effects of a Type 2 inclusion in the material, and the number of SEMs containing a partially debonded region is randomly sampled. The inclusion-matrix interface in all SEMs is assumed to be perfectly bonded, and a partially debonded region is realized by identifying and deleting a set of elements based on the study in Ref. [52]. While introducing a partially debonded region, a void of ∼ 5 μm thickness is created between the inclusion and matrix (see Fig. 3), which is sufficiently large to ensure no contact between inclusion and matrix surfaces for up to 3% applied strain. Since, in the present research, the simulated applied strain range is ≤ 0.75% , the necessity to add friction or contact surfaces between the inclusion and matrix is avoided. We note that although all inclusions in this work are geometrically modeled as ellipsoids, i.e., Type 1 inclusions, a structurally soft inclusion with a partially debonded region would implicitly capture some of the effects of a Type 2 inclusion in the material.

Inclusion Position Within the Polycrystalline Aggregate
The emergence of IDF or MDF is a consequence of complex interaction between the inclusion attributes and the microstructure surrounding the inclusion [52]. In other words, a specific location of the inclusion within the microstructure does not solely lead to a preferred failure mode. Although all surface-connected or subsurface inclusions do not lead to fatigue failure, in the present experimental dataset for RR1000, most of the IDFs originate from the surface and subsurface inclusions. Hence, in this study, the inclusions are placed at randomly chosen points on the grain boundaries, and it is ensured that inclusions in (a) three SEMs are surface-connected, (b) five SEMs are at the subsurface, and (c) one SEM is at the interior of the microstructure (see Fig. 4e). As the name suggests, the surface-connected inclusions touch one of the two free surfaces (see Fig. 5 and accompanying discussions in section "Boundary Conditions for Crystal Plasticity Finite Element Simulations)" of the SEMs. The effect of the inclusion influences the mechanical response of the surrounding anisotropic matrix material up to a spatial distance of its equivalent diameter [52] (which is indicated for each SEM based on the sampling described in section "Inclusion Size and Morphology"). Hence, the subsurface inclusions are placed within a distance of the equivalent diameter of the inclusion, d eq , from the free surface, whereas the interior inclusions are placed at a distance greater than d eq from the free surface.

Crystal Plasticity Constitutive Model
As mentioned before, in this work, the inclusions are modeled as homogeneous isotropic linear elastic materials, and crystal plasticity constitutive equations are adopted for the matrix material. Such an idealization is reasonable because the non-metallic inclusions are brittle, whereas plasticity is the consequence of dislocation motion in the matrix. At each integration point within a grain, the matrix material is homogenized to capture both the phase and ′ precipitates (secondary and tertiary). A rate-dependent crystal plasticity framework, incorporating the temperature effect, is adopted to model such a homogenized response at the integration points. With this, we decompose the total deformation gradient F at a material point into a plastic, thermal, and elastic parts as follows: Here, F e is associated with elastic stretch and rotation of the crystal, F is associated with the thermal expansion, and F p captures the plastic deformation. Moreover, F does not evolve with deformation, and is written as where I is the identity tensor, is the CTE, and ΔT is the temperature change. The plastic velocity gradient, L p , is written as [66] where, ̇j , ŝ j , and n j are the shear strain rate, a unit vector along slip direction, and a unit vector along the slip plane normal associated with the jth slip system, respectively; N s is the number of slip systems, which corresponds to twelve ⟨110⟩{111} octahedral slip systems in RR1000. The shear strain rate ( ̇j ) and the resolved shear stress ( j ) on the jth slip system is related via a power-law flow rule [67] as follows: where ̇0 is the reference shearing rate, g j is the reference shear stress or slip system resistance, j is the back stress, and m is the strain-rate sensitivity exponent. In this work, g j is described via a Taylor-type hardening law where g 0 is the initial slip resistance, j is the total dislocation density on the jth slip system, b is the Burger's vector, μ is the shear modulus, and h n is a scaling parameter. We adopt a Kocks-Mecking's approach [68][69][70] to describe the evolution of the total dislocation density j as Here, k 1 and k 2 (̇, T) are material constants and correspond to the dislocation storage and annihilation rates, respectively; T is temperature, and ̇ is the applied strain rate. Further, the constants k 1 and k 2 (̇, T) can be explicitly related to temperature and strain-rate [71,72] as follows: Here, Γ act is the activation energy, D is a scaling constant, k B is Boltzmann constant, and ̇0 is the reference strain rate. Equation (8) is primarily applicable for dislocation glide mediated plasticity if 10 −5 s −1 ≤̇≤ 10 3 s −1 . Finally, the evolution of back stress is defined via an Armstrong-Fredrick type [73,74]

description as
Since the present CPFE framework is temperaturedependent, we need to define the temperature-dependent variations for the elastic constants ( C 11 , C 12 , and C 44 ) to accurately compute stress at a material point as a function of temperature. Since our working range of temperature is 20−750 • C , the elastic constants can be assumed to be linearly decreasing with temperature [75]. Therefore, we define: where p 1 , p 2 , and p 3 are positive and non-dimensional constants, and ascertained by regressions through experimentally determined Young's modulus values as a function temperature; T 0 is a reference temperature and set to 20 • C . The elastic constants at the reference temperature, i.e., C 11 T 0 , C 12 T 0 , and C 44 T 0 , are adopted from Ref. [48]. The remaining CP parameters are directly used from the previous work in Ref. [52].

Boundary Conditions for Crystal Plasticity Finite Element Simulations
As indicated in Fig. 5, during CPFE simulations, the normal displacement is set to zero on three mutually orthogonal and adjacent surfaces of the SEMs, and a nonzero normal displacement is specified on one of the remaining three faces, which is orthogonal to the y-axis. Such a choice of the boundary conditions ensures only the yy component of the Cauchy stress tensor to be dominant when is homogenized over the entire simulation domain. Therefore, nominally, the stress state within the simulation domain is uniaxial, which would allow us to compare the simulated fatigue life data to experimental fatigue life obtained from uniaxial tests (see section "Results"). Additionally, two surfaces of the SEMs are unconstrained via the choice of

Residual Stress
The effect of residual stress on the fatigue performance of metallic materials is well-recognized [76]. Since nonmetallic inclusions and matrix material in Ni-base superalloys have dissimilar CTEs, residual stress is induced in the material during the thermo-mechanical processing. Previously, Naragani et al. [50] have quantified the residual stress distribution in the vicinity of a large-seeded inclusion within RR1000 using high energy X-ray diffraction microscopy (HEDM) experiments. A temperature of 760 • C (which is a common aging temperature [54]) was selected for the stress-free condition or ΔT = −740 • C in Eq. (3), which produces a type 2 residual stress distribution similar to the HEDM measurements in Ref. [50]. Therefore, in this work, we carry out an elastic cooling simulation considering ΔT = −740 • C in Eq. (3) to initialize the residual stresses in the SEMs containing inclusion. Subsequently, fatigue load is applied to the SEMs utilizing the full crystal plasticity simulations.

Nonlocal Averaging Scheme
The micromechanical field outputs from the CPFE simulations often show spurious numerical oscillations [58,77]. Such a trend can be related to one or a combination of the following: (a) volumetric locking or shear locking due to the choice of the linear tetrahedron (C3D4) elements [77], (b) lack of mesh refinement near the grain boundaries across that may result in high gradients in field variables, (c) the effect of the boundary conditions on the elements near the corresponding surfaces, and (d) poor quality of C3D4 elements, especially near the grain boundaries [34]. To further explore the last point, i.e., the quality of the C3D4 elements, we consider two metrics, namely and , where and In Eq. (13), CR and IR represent the radii of the circumscribed and inscribed spheres to a tetrahedron element. In Eq. (14), S rms represents the root-mean-square value of the element's edge lengths, and V is the element's volume. A . For a randomly chosen SEM, cumulative probability plots for and are shown in Fig. 6a, b, respectively. Most of the elements ( ∼ 99% ) are good and fall within the shaded region in Fig. 6a, b. However, ~ 1% of the elements are poor, which might lead to spurious results from the CPFE simulations. In the literature, multiple strategies exist to address this situation. For example, Cheng et al.
[77] explored three methods that are implemented during the finite element simulation to mitigate this issue, namely node-based uniform strain method, B-bar-based element method, and F-bar-based element methods. The use of an appropriate regularization scheme during the post-processing phase presents another opportunity to address implications due to poor elements.
Here, we choose the second route, which also involves a significantly lower overall computation cost than the previously mentioned methods. In recent work, Prithivirajan and Sangid have made a comparative assessment of three different averaging schemes, namely (a) nonlocal averaging, (b) band-averaging, and (c) grain-averaging [34]. Based on their analysis, it is apparent that nonlocal averaging is appropriate to predict the location of fatigue failure, fatigue life, and associated failure mode. Hence, we adopt the nonlocal averaging scheme in this work. Such an averaging scheme has previously been rationalized and used by other researchers (see, for example, [33,52,58,[79][80][81][82]). To illustrate the nonlocal averaging scheme, we consider an integration point as shown in Fig. 6c and identify three mutually orthogonal unit vectors passing through the integration point, namely the slip direction ( ŝ j ), slip plane normal ( n j ), and a transverse vector ( t j ) on the slip plane corresponding to the jth slip system, such that these directions form a right-hand coordinate system at the integration point.
Now, a small rectangular cuboid of volume V j is considered having its center coinciding with the integration point and faces parallel to the three orthogonal directions defined in Fig. 6c. The micromechanical field variable is averaged over all the elements whose centroids lie within the averaging volume. One can construct twelve such volumes around the same integration point corresponding to each of the twelve octahedral slip systems, i.e., for j = 1 to 12 . After performing averaging for all twelve such volumes, the maximum (by magnitude) of the averaged values is assigned as the slip system averaged number to that integration point.
To illustrate the nonlocal averaging scheme mathematically, let ϕ CPFE (x) and ϕ NLA (x) be the solution of a micromechanical field variable ϕ(x) from the CPFE simulation and a nonlocally averaged value of the CPFE solution, respectively. Hence, ϕ NLA (x) is related to ϕ CPFE (x) as For integration points near the grain boundaries, the averaging volume is not allowed to cross the grain boundary (see Fig. 6d) 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. Based on a sensitivity analysis, the lengths of the averaging volume are defined by the size of 5, 3, 3 elements along the slip direction ( ŝ j ), slip plane normal ( n j ), and transverse direction ( t j ), respectively.

Fatigue Life Prediction Model
From the CPFE analysis, the increment of plastic strain energy density at a material point x over the ith cycle, Δw p i (x;B) , can be written as where B is a set of parameters that define loading condition, such as strain range ( Δ ), strain ratio ( R ), and temperature ( T ). Accordingly, the APSED after the ith cycle at the same point, w p i (x;B) , is expressed as From classical fatigue experiments, the hysteresis behavior saturates upon cyclic loading when stable dislocation structures and substructures are formed within the material [76]. If cycle N s (B) is associated with the onset of saturation at the loading condition B , subsequent accumulation of plastic strain energy density can be considered to be steady, at the rate of Bandyopadhyay et al. [33] showed that if the material fails due to fatigue at x = x * , w p N f (x * ;B) is a material property and independent of Δ and R at room temperature. They referred to w p N f (x * ;B) as the critical APSED of the material and denoted as W p critical . In this work, we hypothesize that the critical APSED, just like other material properties, such as the stiffness and strength properties of a material, is a function of temperature, i.e., W p critical = W p critical (T) . One objective of the present work is to ascertain the temperature-dependent variation of W p critical (T) . Subsequently, W p critical (T) can be used to predict fatigue life using a slightly rearranged version of Eq. (18) as

Bayesian Calibration of the Critical Plastic Strain Energy Density
If  [33] constructed a Bayesian inference-based framework that takes into account the variabilities associated with the simulated quantities and experimental fatigue life and produces a distribution of W p critical (T) at a given temperature, T . Subsequently, the mean or expected value of the distribution can be treated as a representative value of the quantity and used in Eq. (19) to predict fatigue life at temperature T . In Ref. [33], the experimental fatigue life data was used for a given Δ and R to calibrate the distribution of W p critical (T) at room temperature and subsequently, fatigue life was predicted at eight different loading conditions, including multiple Δ and R . In this work, we use all the fatigue life data that is available at a given temperature T to obtain the distribution of (18) W p critical (T) and repeat the process at multiple temperatures to characterize the temperature-dependent behavior of W p critical (T). A schematic of the Bayesian calibration approach for W p critical (T) is shown in Fig. 7. Starting from the uniform prior distribution, the posterior distribution for W p critical (T) at room temperature results in a normal distribution [33]. Therefore, in this work, we start with normal prior distributions of W p critical (T) at a given temperature, T , and, subsequently, update the mean and the standard deviation of the distribution. With available w p N s (x * ;B) and Δw

Comparison of the Experimental and Predicted Fatigue Life
The experimental ( N , are assumed to be IID random variables, the likelihood term in Eq. (22) can be expressed as [83]  Typically, it is challenging to evaluate the integral in the denominator of Eq. (22) analytically or using numerical quadrature. The MCMC sampling method is preferred in such situations to obtain (A|D(B)) . A Markov chain is a sequence of random variables satisfying the property that the hth term in the sequence depends only on the (h − 1)th term. Based on a proposal distribution, the hth term is identified via Monte Carlo sampling and the decision whether the new proposed candidate point will be accepted as the hth term in the sequence is taken based on the MH algorithm [84,85]. Now, the MH algorithm is invoked to decide if A * r can be accepted as the hth candidate in the sequence. Particularly, a uniformly distributed random number is gener- This process is continued until convergence is achieved. M Markov chains are run in parallel to check convergence. Let, B A r and W A r be the variance of the parameter A r between and within chains, respectively. An estimate of the variance of the parameter A r , V A r , can be computed as Now, the convergence parameter R A r is defined as With the increase in Monte Carlo sampling, h , in each Markov chain, the parameter R A r approaches unity resulting in a converged/stationary distribution, also known as the posterior distribution. The readers may refer to Cross et al. [86] for the expressions of W A r and B A r .

Results
To

Critical Accumulated Plastic Strain Energy Density as a Function of Temperature
As representative examples, the evolution of Δw p k (x * ;B) with the number of cycles, k , at the critical location ( x * ) within one of the 20 SEMs, is shown in Fig. 8. 3 It is evident from Fig. 8 that the saturation of Δw p k (x * ;B) is inherently dependent on Δ and T , i.e., B . For the current set of 20 SEMs, each containing ∼ 2.5 M elements, it is impractical to carry

3
out CPFE simulations for hundreds of cycles. Under such situations, researchers either continue simulations up to a point, usually three to 20 cycles, where the saturation of the fatigue indicative metric is acceptable [9,11,13,14,22,23,26,[31][32][33][34] or adopt accelerated simulation procedures based on the method of multiple time-scales [87][88][89][90] or perform reduced-order CPFE simulations [91]. The first approach has been widely used in the literature and is considered in this work. In such an approach, we define the degree of saturation at the kth cycle, ξ k (x * ;B) , as The degree of saturation is considered to be acceptable if ξ k (x * ;B) ≥ 99.9% . With this, for Δ = 0.6% and 0.75% , the simulations are continued up to 10 cycles; for Δ = 0.4% , simulations are extended up to 20 cycles. Consequently, N s = 10 for Δ = 0.6% and 0.75% ; N s = 20 for Δ = 0.4%.
For the Bayesian calibration of W Fig. 10 The probability of failure plots considering log-normal distribution as reference. In each subplot, fatigue life is normalized with respect to the median predicted life. Each row corresponds to a constant temperature with T = 20 • C for the first row (from top), T = 300 • C for the second row, T = 600 • C for the third row, and T = 750 • C for the fourth row. Each column corresponds to a constant applied strain range with Δ = 0.4% for the first column (from left), Δ = 0.6% for the second column, and Δ = 0.75% for the third column. MDF stands for matrix-driven failure, and IDF stands for inclusion-driven failure. 'free surface' represents failure originating from the free surface of the microstructure, and 'interior' represents cases where failure originates within the bulk of the microstructure ◂

Prediction of Fatigue Life and Associated Failure Modes
The simulated and nonlocally averaged w  Fig. 9, to predict fatigue lives and the associated failure modes. Specifically, for SEMs with inclusions, if x * is located within the distance, d eq , from the centroid of the inclusion, where d eq is the equivalent diameter of the inclusion, such a situation would indicate an IDF; otherwise, the simulation predicts an MDF. The probability of failure plots are obtained using the probplot function in MATLAB, considering a log-normal distribution as a reference and shown in Fig. 10. The normalized fatigue life is plotted in Fig. 10, where the normalized factor is the ratio of the actual life and median life, whereas a life of 10 on the scale represents a factor of 10 multiplied to the median life.
As previously shown [57], this framework, using the W p critical (T) , is capable of predicting log-normally distributed fatigue life at room temperature, and the predicted scatter increases with decreasing Δ . Both characteristics of the predicted life distributions are emergent behavior from the model and indicate the physical nature of fatigue behavior. In this work, a similar observation is also made in Fig. 10. For example, (a) the predicted data points in each subplot in Fig. 10 appear to be aligned along a straight line, confirming that the predicted fatigue life follows a log-normal distribution and (b) for each temperature, scatter in the predicted fatigue life increases with decreasing Δ .
In section "Critical Accumulated Plastic Strain Energy Density as a Function of Temperature", we quantified the variability in W p critical (T) , i.e., (B) , using Eqs. (32)- (35). We also mentioned that the variability in the simulated w p N s (x * ;B) and Δw Although inclusions are present in 10 SEMs, only one IDF is predicted at Δ = 0.4% , consistently, at all temperatures (see Fig. 10). Further, the fatigue life associated with the IDF does not correspond to the minimum life predicted from the simulations. In the present experimental dataset, evidence of IDF is only observed for Δ ≤ 0.4% , and the life associated with IDF does not correspond to the minimum life. Although Δ < 0.4% regime is not explored in this study, however, based on the results of the prior parametric study and accompanied discussion therein [52], a higher probability of IDF is expected in the Δ < 0.4% regime.
As mentioned in section "Background", a turbine disk produced with a DMHT shows fine grains near the bore region and coarse grain microstructure near the rim region. In experimental research with a turbine disk made of LSHR, Gabb et al. [93] showed that the fine-grained microstructure leads to higher fatigue life than the coarse grain microstructure. In their research, Type 2 IDF is observed in the fine-grained region, whereas crystallographic surface and subsurface facets are responsible for failure in the coarse grain region. The crystallographic facet resulted in the minimum fatigue life observed in their samples. The experimentally observed minimum life for the coarse grain variant of RR1000 is associated with MDF as well. From Fig. 10, one can verify that the predicted minimum lives in all simulated loading conditions correspond to an MDF. Further, the MDF corresponding to the minimum life originated from a microstructure without inclusion, thereby confirming that the minimum life is not necessarily always influenced by the presence of inclusions.
Gabb et al. [93] reported that the location of the failure, associated with MDF, shifted more towards the interior of the samples with increasing temperature as opposed to originating from the free or near free-surface. Further, failure originating at the free surface was associated with lower life than failure arising from the interior. Additionally, Hyzak and Bernstein [36] observed a change in the crack initiation location from the surface to the interior of the AF-115 and AF2-1DA specimens with decreasing strain range. The present fatigue life prediction framework captures these trends in temperature-and strain range-dependent behavior of failure location. For example, it can be seen in Fig. 10 that (a) the minimum life is always associated with the MDF originating from the free surface, (b) failure arising at the interior of the microstructures are, in general, associated with relatively higher fatigue life, (c) failure originates at the interior of the SEMs more frequently with decreasing strain range. Further, from Fig. 11, it is observed that the percentage of free surface failures decreases with increasing temperatures and thereby indicating the location of the failure, associated with MDF mode, shifting more towards the interior of the SEMs. From Fig. 12a, we can verify that such a trend is in good qualitative agreement with the experimental observation in RR1000.
From the simulation results in Fig. 11, the percentage of IDF remains unchanged with changing temperature. Such an observation is perhaps because the same set of SEMs is simulated across all temperatures. Overall, the modeling framework predicts 2% IDFs 5 (Fig. 12b). These failures originate in an SEM having a structurally soft, partially debonded, and subsurface inclusion. On the other hand, the present Fig. 12 a Comparison of the percentage of free surface failure from simulation and experiments. Percentage of inclusionand matrix-driven failures from b simulation and c experiment. MDF stands for matrix-driven failure, and IDF stands for inclusion-driven failure experimental dataset show 6% IDFs 6 (Fig. 12c), where all such failures originate from subsurface Type 2 inclusions. Although all inclusions in this work are geometrically modeled as ellipsoids, i.e., Type 1, a structurally soft inclusion with a partially debonded region would implicitly capture some of the effects of a Type 2 inclusion in the material. Thus, the predictions regarding the inclusion type associated with IDF are in good agreement with the experimental observations. Further, a recognized underestimation of the percentage of IDFs from the present modeling framework is attributed to the presence of a structurally soft inclusion with a partially debonded region, e.g., the worst-case scenario, in only one SEM. From the results, as presented above, it is clear that the microstructure-sensitive energybased approach [33] is capable of capturing the median and scatter trend in fatigue life data, as well as the multiple failure modes and realistic shift in the location of failure with varying applied strain range and temperature.

Conclusions
In this work, a microstructure-sensitive energy-based approach is used to predict the fatigue life, associated failure mode, and uncertainty in the fatigue life prediction as a function of applied strain range and temperature. The fatigue life prediction model used in this work involves only one material constant: the critical accumulated plastic strain energy density (APSED). Using a Bayesian inference method, the critical APSED is calibrated for a powder metallurgy produced Ni-base superalloy, namely RR1000, as a function of temperature. Subsequently, the calibrated critical APSED is used to predict inclusion-and matrix-driven competing failure modes in RR1000 as a function of applied strain range ( Δ = 0.4% , 0.6% , and 0.75% ) and temperature ( T = 20 • C , 300 • C , 600 • C , and 750 • C ) for R = 0 loading regime. Specifically, for fatigue life prediction, 20 statistically equivalent microstructures, 10 of these containing inclusions and the remaining 10 without inclusions, are used to conduct simulations using crystal plasticity finite element (CPFE). The key findings are summarized below.
• The microstructure-sensitive critical APSED decreases along a vertically reflected sigmoidal curve with increasing temperature, which is consistent with the fatigue strength of the material decreasing at the elevated temperature.
• The predicted fatigue life from the statistically equivalent microstructures follows a log-normal distribution at each simulated loading condition, which is consistent with the traditional fatigue experimental results. Moreover, the predicted scatter in fatigue life also increases with decreasing strain range for a given temperature. • For each simulated loading condition, the uncertainty in fatigue life is quantified as a prediction interval computed based on the 95% confidence level of the critical APSED and the simulated APSED from CPFE. Experimental data for each loading condition is seen to fall within such a prediction interval. • From CPFE simulations, inclusion-driven failure is observed only at Δ = 0.4% across all temperatures. Experimentally, inclusion-driven failures are only observed for Δ ≤ 0.4% . Although Δ < 0.4% regime is not explored explicitly in this study, it has been previously shown in a parametric study [52] that there exists a higher probability of inclusion-driven failure in the Δ < 0.4% regime. • From the simulations, the inclusion-driven failure mode is not (necessarily) associated with the minimum life, which is consistent with experimental observation. Albeit, this observation may change with microstructure inputs that include larger inclusions (far greater than the average grain size of the material). Experimentally, the minimum life corresponds to matrix-driven failure. In simulations, minimum life also corresponds to matrixdriven failure and originates from a microstructure that does not contain an inclusion. • In the literature, it is seen that the location of the failure, associated with matrix-driven failure, shifts towards the interior of the specimen from the free surface with increasing temperature, and the failures originating from the interior are typically associated with higher life. In simulations, the percentage of failure originating from the free surface is seen to be decreasing with increasing temperature as well. Further, life, associated with failure from the interior, is predicted to be consistently higher than the free surface failure life.
The above findings extend the framework proposed by Bandyopadhyay et al. [33] across various temperature regimes and competing failure modes while offering uncertainty quantification and propagation. Thus, this work provides a valuable step towards the adoption of microstructural-based modeling assessments for fatigue lifing.
We also thank Dr. Prithivirajan Veerappan from Purdue University for sharing the scripts to quantify the quality of C3D4 elements.

Compliance with Ethical Standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.