Introduction

Ni-based superalloys, such as René 88DT, are widely used in high-temperature applications especially in the aerospace propulsion industry, on account of their excellent thermomechanical and chemical properties. Despite significant advantages, the fatigue life of these materials is compromised due to the presence of heterogeneities in the microstructure. While very few metallurgical defects like inclusions or voids are found in their polycrystalline microstructures, typically manufactured by powder metallurgy techniques, heterogeneities manifest as spatial variations in the crystallographic and grain morphology distributions, as well as intragranular compositions1. Additionally, annealing twins are critical heterogeneities that form in the polycrystalline microstructures during processing. These microstructural variations are major contributing factors to their fatigue failure under operating conditions2.

A number of studies, e.g. in3,4,5,6,7 have observed nucleating cracks near coherent Σ3 twin boundaries under fatigue loading conditions. For the polycrystalline René 88DT, fatigue cracks are often found to nucleate at these twin boundaries, especially where the adjacent grains also have a high Schmid factor1,8. Crack nucleation is preceded by the formation of slip bands that are parallel to the coherent twin boundaries, yet slightly offset from the boundary itself9,10. In addition to the high Schmid factor, high elastic modulus mismatch across the boundary and long twin boundaries have also been associated with the probable locations of crack nucleation in11,12,13. A hypothesis behind such observations of crack nucleation is that these geometric and crystallographic metrics can be correlated with high values of the resolved shear stress on shared slip systems between twin and parent grains, parallel to the twin boundary14,15,16,17,18,19. The studies conjecture that rather than high Schmid factors, the high resolved shear stresses promote dislocation glide along the shared plane leading to nucleation. Out-of-plane forest dislocations do not disrupt slip on the shared plane, and hence dislocation glide is unimpeded along the boundary. With the progression of cyclic loading, highly localized slip parallel to the twin boundary trigger the formation of persistent slip bands (PSB), which precede crack nucleation and propagation20. Experiments on Ni-based superalloys have shown that PSBs can form at strains below 0.35% within a single cycle of deformation9.

These observations necessitate the development of robust mechanism-driven predictive models that can relate evolving state variables in the underlying polycrystalline microstructure to the nucleation of cracks. The last two decades have witnessed significant advances in the understanding and prediction of the connection between microstructure and local mechanical fields, along with critical failure indicators under fatigue loading. This has been enabled by a combination of advanced experimental characterization techniques like electron backscatter diffraction (EBSD), high-energy diffraction microscopy (HEDM), and X-ray tomography/diffraction techniques, and computational models like adaptive finite element methods (FEM) and fast Fourier transform (FFT) methods integrated with crystal plasticity constitutive models. In the computational modeling approaches, virtual 2D and 3D polycrystalline microstructures are generated from materials characterization data, and subsequently image-based micromechanical simulations are conducted to interrogate the locations of crack evolution21,22,23. Comprehension of material states near the critical cracked regions has guided researchers in constructing metrics of fatigue and failure indicators. Both physics-based and data-driven approaches have been pursued to create the appropriate fatigue metrics, such as dislocation density and plastic slip-based criteria combined with crack opening stresses for susceptible grain configurations24,25, energetics and stability of PSBs in26, accumulated plastic slip and elastic anisotropy in14, and cyclic stored plastic energy density21,23,27,28. Important fatigue indicator parameters have thereby been derived through the use of correlation metrics, variance-based sensitivity, or physics-derived arguments in combination with machine learning and Bayesian techniques21,22,29,30,31.

Despite significant advances, shortcomings, e.g. in experimental throughput rates and computational efficacy, have in some cases mitigated the effectiveness of these models. Limitations in the collection of experimental data sets linking microstructural configurations to crack nucleation sites have often led to the development of models from very limited microstructural data sets, typically encompassing a single grain configuration. Alternatively, the microstructural geometry has been abandoned entirely in favor of higher-scale statistical analysis. The lack of realistic experimental data leads to inherent uncertainty in the model output. The resulting fatigue indicators are not likely to retain fidelity with the introduction of new data. Furthermore, probabilistic methods, required to handle multiple experiments in a high-dimensional input data space, suffer from various issues resulting from the problem dimensionality. From a modeling point of view, the computational methods must ensure accurate determination of state variables like stresses, plastic strains, dislocation densities, etc. at potential crack nucleation sites for effective development of mechanics-based models of crack nucleation. Boundary conditions play a very important role in the prediction of localized variables in the image-based micromechanical models32,33. Conventionally used periodicity boundary conditions can lead to gross inaccuracies when modeling fatigue behavior.

This paper introduces a data-driven Bayesian model framework for predicting fatigue crack nucleation in polycrystalline microstructures of the Ni-based superalloy René 88DT, overcoming some of the limitations mentioned above. The overall approach leverages large data sets obtained from crystal plasticity FEM-based micromechanical simulations with statistics and machine learning algorithms. Within this framework, the probability of observing a crack nucleation event is derived based on computationally generated micromechanical data on state variables at a potential nucleation point. High fidelity, image-based crystal plasticity finite element (CPFE) models developed in34,35,36,37 are adopted for simulating evolving state variables in microstructural representative volume elements (RVE) of polycrystalline Ni-based superalloys. The deficiencies with improper boundary conditions on the RVE is alleviated in this work through a concurrent multiscale model, where the microstructural RVE is embedded in a homogenized exterior material domain. The exterior domain is modeled using an upscaled, anisotropic continuum plasticity constitutive model that is obtained from homogenization of crystal plasticity RVE simulations38,39,40,41,42,43,44. This model is discussed in the Supplementary Methods Section.

The paper introduces a robust Bayesian model framework for extracting fatigue nucleation criterion using a combination of simulation and experiments. The primary contributions of this framework are delineated below.

  • The probabilistic crack nucleation model accounts for multiple experimental microstructures that are simulated using crystal plasticity FEM with one-to-one virtual representations of EBSD scans. A major factor for model calibration and validation is the generation of statistically significant results over a set of multiple experimental microstructures, covering a range of microstructural variabilities. By incorporating multiple one-to-one experimental data sets directly in the probabilistic modeling framework, the approach mitigates uncertainties arising from a single microstructure dependency and insufficient validation data. An important additional aspect is the accounting for material states of both cracked and uncracked regions in the development, through efficient algorithms for handling large data sets.

  • The Bayesian crack nucleation model is formulated using a customized Nataf transformation to represent multivariate joint posteriors. The resulting probabilistic framework yields several advantages over other comparable methods, viz. (i) semi-analytical expressions for high-dimensional joint probability distributions, (ii) easy sampling of arbitrary points in the data-space with no reliance on specialized sampling techniques, (iii) simple calculation of moments and measures on the joint probability distribution, and (iv) straightforward regularization of discrete data that reduces a generalized model error. These advantages enable efficient use of the final posterior joint distribution in predicting crack nucleation.

  • The automated simulation framework creates pipelines, which take a set of EBSD scans delineating grain configurations and crack nucleation sites as input, and returns a reduced order state-variable-based Bayesian model for probabilistic nucleation prediction as shown in Fig. 1. The framework integrates multiple computational procedures into a continuous workflow.

  • The model introduces an algorithm for selecting important state variables and isolating critical metrics that indicate locations of crack nucleation. This algorithm assesses successive multi-dimensional joint distributions and accounts for highly correlated sets of variables to be incrementally optimal in its predictions.

  • The automated procedure with minimal intervention uses experimental data to inform the probabilistic model of the most informative parameters. These parameters yield mechanistic insights into the drivers of crack nucleation for this material. The resulting probabilistic model demonstrates robust performance across variations in the microstructure and local grain configuration.

  • A scalar crack nucleation indicator is derived from the simulation platform.

The process of developing a probabilistic crack nucleation modeling framework integrates experimental testing and characterization with a pipeline of simulation tools for mechanical modeling and data-driven, supervised learning techniques. The sequence of four major modules in this pipeline is schematically depicted in Fig. 1.

Fig. 1: Four modules comprising the pipeline for developing the data-driven crack nucleation model.
figure 1

(I) experimental testing and characterization data, (II) development of the concurrent multiscale model with embedded microstructural RVE, (III) creation of a database of microstructural state variables, and (IV) development of the Bayesian crack nucleation model.

Experimental testing and characterization

In this module, specimens of the Ni-based superalloy René 88DT are subjected to interrupted fatigue loading at various stages of deformation, to track the precursors of crack nucleation. Scanning electron microscopy (SEM) maps are acquired and utilized for identifying microstructural locations of crack nucleation. Simultaneously, EBSD maps are acquired to determine the crystallographic orientations and grain morphologies at locations of interest in the microstructure.

Concurrent multiscale model development with embedded microstructural RVE

This module builds a two-scale concurrent model, consisting of an embedded polycrystalline microstructural RVE in a homogenized outer domain. Important microstructural deformation variables are extracted from simulations of the concurrent model of the experimental specimen. The two-scale concurrent representation is constructed to avoid direct application of boundary conditions on the RVE, thereby averting boundary condition-related inaccuracies in micromechanical fields. The embedded microstructural region is generated by extruding 2D EBSD scans in the out-of-plane direction. It is modeled using a crystal plasticity constitutive model for individual grains in the polycrystalline microstructure, while the outer domain is modeled by a homogenized, rate-dependent elasto-plastic constitutive model that is calibrated to the average mechanical response of the embedded RVE.

Creating a database of microstructural state variables

The concurrent multiscale model is simulated with fatigue loading conditions applied on the external boundary, from which the evolving state variables in the embedded microstructural RVE region are extracted. The state variables are then processed and archived in a crystal plasticity simulation database for subsequent stages of the data-driven modeling.

Development of a Bayesian crack nucleation model

The crystal plasticity simulation database provides the input data for the modeling framework to determine the probability of crack nucleation corresponding to a mechanical state in the material volume. The Bayesian model selects the most informative state variables for predicting crack nucleation by navigating the covariance structure of candidate predictor variables. The results are used to create a simple model for identifying potential locations of crack nucleation in the microstructure.

Each module is standardized and automated so that the model building process can be rapidly and systematically repeated for all experimental microstructures. Furthermore, this approach allows for the extension to other critical events as well as the inclusion of new data sets.

Results

Material description

The material investigated in this study is a polycrystalline Ni-based superalloy René 88DT. Generally produced by powder metallurgy techniques45,46, René 88DT has a nominal chemical composition (weight percent) of 13% Co, 16% Cr, 4% Mo, 4% W, 3.7% Ti, 0.7% Nb 0.03% C, 0.015% B, and the balance in Ni47. The intragranular microstructure of René 88DT consists of an FCC-structured Ni matrix (γ phase) and an embedded population of L12-structured Ni3Al precipitates (γ' phase) with a bimodal size distribution of large (~100 nm) secondary and small (<10 nm) tertiary \(\gamma ^{\prime}\) precipitates. At the polycrystalline scale, the material has an average grain size of about 26 μm, very weak crystallographic texture, and a significant population of Σ3 twin boundaries. Details on the René 88DT microstructure are given in8,9.

Experiments on fatigue tests and data acquisition

An experimental data set consisting of 12 specimens with independent grain structure configurations is utilized for developing the Bayesian model in the present study. Room temperature, low-cycle fatigue tests are conducted to evaluate the fatigue life of the Ni-based superalloy René 88DT, for which the mechanical testing procedures are detailed in9. The tests are performed in air on an electromechanical machine in a uniaxial, symmetric push-pull mode. Fully reversed loading tests with an R-ratio equal to −1 at a frequency of 1 Hz, are carried out in a stress-controlled mode at a maximum stress of 758 MPa. Cylindrical specimens with a gauge diameter of 5 mm and a gauge length of 16 mm are used in these tests. The strains in the specimen are measured by a mechanical extensometer that is positioned on the specimen gauge. The average fatigue lifetime of the specimens is observed to be 65340 cycles with a scatter (standard deviation) of 5530 cycles. A maximum and minimum total strain of 0.35% ± 0.04% and −0.35% ± 0.03% are measured for the first cycle. The material exhibits cyclic evolution of the strain, with the maximum and minimum strain stabilizing at 0.40% ± 0.05% and −0.30% ± 0.04% respectively, after half of the lifetime i.e. around 30–40,000 cycles. The strain range remains constant at 0.70% ± 0.03% during the entire fatigue test. Interrupted tests are performed to enable specimen surface examination using SEM and EBSD measurements at about 45% of the average fatigue lifetime (~30,000 cycles). The SEM scans identify regions exhibiting evidence of micro-crack formation, while the EBSD measurements interrogate the configuration of the surface grains surrounding the micro-cracks. The interrupted tests capture two types of cracks, viz. those confined to one grain48 (type I) and those extended between one and three-grain diameters, crossing one or two grain boundaries on the surface (type II). This approach facilitates the identification of the low-cycle crack nucleation sites that are characteristic for each microstructure. A few typical microstructures, showing nucleated cracks under the fatigue loading conditions, are shown in Fig. 2.

Fig. 2: Two sample experimental data sets used in the development of the Bayesian crack nucleation model.
figure 2

The SEM images display fatigue crack nucleation near and parallel to the twin boundaries. The associated EBSD maps present the grain structure around the crack nucleation sites.

Experimentally generated polycrystalline microstructures

Subsequent to EBSD acquisition, an automated data processing pipeline is applied to each of the 12 investigated grain structures in the microstructure generator software DREAM.3D49 to remove any measurement defects. Low image quality voxels are identified and automatically filled. Grains are determined by grouping voxels within a 5 misorientation angle threshold. Grains containing an area <20 μm2 are merged with their surrounding grains using a nearest-neighbor identification method. These processing stages result in a collection of 2D segmented images, free of artifacts that can cause difficulties with generating robust finite element meshes.

Key experimental observations on fatigue crack nucleation

The primary sources of heterogeneity in René 88DT are derived from variations in microstructural features since minimal metallurgical defects like inclusions and voids are generated with powder metallurgy processing. The most prominent microstructural features are annealing twins, whose boundaries are uniformly dispersed throughout the grain structure during heat treatment. These boundaries have unique crystallographic and morphological properties and thereby affect dislocation motion during material deformation.

Prior research, e.g. in1,9,50,51, suggests a number of experimental correlations between the material grain structure and physical location of crack nucleation sites. A few key observations are delineated below:

  1. 1.

    Fatigue crack nucleation and the highest degree of strain localization are systematically reported to occur near and parallel to Σ3 twin boundaries in the low, high and very high cycle regimes9;

  2. 2.

    Crack nucleation always occurs between pairs of grains, which are relatively oriented in a parallel slip configuration;

  3. 3.

    Additional geometric, crystallographic, and mechanical factors, such as high Schmid factor, large difference in elastic modulus between adjacent grains, and long straight twin boundary length, have been shown to be associated with crack nucleation sites;

  4. 4.

    Slip bands at crack nucleation sites can form as early as the first cycle, suggesting state variables at early stages of deformation are predictive of the crack nucleation locations;

  5. 5.

    For both high and low-cycle fatigue regimes, the previous studies have reported that all crack nucleation sites are associated with twin boundaries that display parallel slip configurations15.

The above observations are taken into account in this study for predicting the microstructure and grain configurations that generate crack nucleation. Next the pipeline requires the development of a concurrent multiscale model for fatigue analysis, followed by the creation of a microstructural state variable database related to fatigue nucleation. These are respectively discussed in Sections I and II of the Methods Section.

A Bayesian framework for predicting crack nucleation

A stochastic machine learning (ML) method is developed for assessing the probability of crack nucleation in the René 88DT superalloy microstructure. The Bayesian framework, outlined in the Methods Section III, is chosen for generating the probability of observing a crack nucleation event, as a function of the micromechanical state in the microstructure. This is a supervised learning classification approach that is effective with a limited set of high-dimensional data. It has model parameters that are physically interpretable functions of the micromechanical state variables and offers a level of uncertainty in its predictions. The Bayesian modeling framework naturally yields a probabilistic outcome as opposed to a binary (0 or 1) response, which is advantageous over other machine learning classification techniques. The finite element simulation data, generated by the concurrent model in the Methods Section I, is aggregated over all microstructures and used as training data for this stochastic model. The results of these simulations are collected into a segmented database containing over 20 potential crack nucleation predictors in Methods Section II.

A Bayesian framework for predicting crack nucleation

Constructing a Bayesian model using all of these parameters would result in over-fitting, providing little useful information. To remedy this, a model selection technique is developed to systematically reduce the number of state variables to only those that most influence the posterior distribution p(Cx) as a combined set. The objective is to select a reduced order model with only the most significant state variables from the possible k = 20 state variables in the data set X listed in Table 2 as well as their regularized difference terms.

A recursive build-up approach is executed to add state variables one by one to the model and assess the incremental value of each new addition. In order to compare across multiple models, an objective function is designed to quantify the overall effectiveness of each posterior distribution. Let the objective function ρC be defined as:

$${\rho }_{C}=\exp \left(\mathop{\sum }\limits_{i}^{{N}_{C}}{w}_{C}^{(i)}\log (p(C\,| \,{{{{\bf{x}}}}}^{(i)}))\right)$$
(1)

where NC is the total number of cracked material volumes across all microstructures and \({w}_{C}^{(i)}\) are the normalized volumes for each, i.e. \(\mathop{\sum }\nolimits_{i}^{{N}_{C}}{w}_{C}^{(i)}=1\). This objective function ρC represents the geometric mean of the probability that the model correctly labels each known cracked volume of the observed data set. In practice, \(\log {\rho }_{C}\) is used for numerical computations.

The reduced order posterior model p(Cx*) is determined by choosing a subset of state variables X*X that incrementally maximizes ρC for a fixed number of state variables k* ≤ k. This optimization procedure is described in Algorithm 1 in the Supplementary Methods Section. The fundamental idea is to train the posterior PDF with each of the state variables individually, and assess the effectiveness of each model by comparing their respective ρC. After determining the most informative state variable with the maximal objective value, the model is trained again with each of the remaining state variables individually, as well as the previously chosen most informative state variable, in order to select the second most informative state variable. This incremental process is repeated until a priori chosen k* ≤ k state variables are established. The incremental selection process algorithm is given in the Supplementary Methods Section under Algorithm 1.

The incremental selection process is executed with the present data set X comprised of k = 20 state variables, using k* = 8. In each iteration, the method optimally adds the most informative state variable. With each additional state variable included in the probabilistic model, the overall performance characterized by ρC increases. Figure 3a plots the incremental effectiveness of the posterior with additional state variables. The slope of the plot gradually saturates as less new and independent crack-delineating information is provided to the joint posterior PDF from the additional simulation output data.

Fig. 3: Most predictive state variables for crack nucleation locations.
figure 3

a Increasing geometric mean probability of prediction for known cracked material volumes with incremental addition of state variables into the posterior PDF; Cracked and uncracked marginal likelihood distributions for b the von Mises stress σVM, c the maximum plastic slip-rate \(\dot{\gamma }\), and d the plastic defect energy ψp.

The above procedure yields a systematic ranking of state variables from the data set in terms of effectiveness, into a set X*. It provides a clear path to reduce the degrees of freedom in the final model. The top three state variables determined by this process are:

  1. 1.

    The von Mises stress σVM.

  2. 2.

    The maximum plastic slip-rate \(\dot{\gamma }\).

  3. 3.

    The plastic defect energy ψp.

These three quantities contain the most combined information that is capable of distinguishing between cracked and uncracked material volumes within the regions undergoing parallel slip. The marginal distributions for these state variables are shown in Fig. 3a–c, illustrating the overall trends for the crack data separation. Generally, low σVM, high \(\dot{\gamma }\), and high ψp all correlate with an increased likelihood of cracked nucleation at the material volumes in a parallel slip configuration.

A reduced-order posterior is generated from these three state variables by defining a set \({{{{\bf{x}}}}}_{{{\mbox{r}}}}=[{\sigma }_{{{\mbox{VM}}}},\dot{\gamma },{\psi }_{{{\mbox{p}}}}]\). The three-dimensional likelihood PDFs p(xrC, ξ) and \(p(\,{{{{\bf{x}}}}}_{{{\mbox{r}}}}\,| \,\bar{C},\xi )\) are constructed, and subsequently the resultant posterior PDF p(Cxr) is calculated. Two-dimensional projections of the likelihoods are visualized in Fig. 4a, which shows the tendency of the cracked data to separate according to the values of σVM and \(\dot{\gamma }\). Incorporating the priors yields the final projected posterior as shown in Fig. 4b, where the highest probability of crack nucleation commences for low σVM and high \(\dot{\gamma }\) values.

Fig. 4: Illustration of crack nucleation model performance.
figure 4

a 2D projection of the iso-probablistic contours of the final uncracked and cracked likelihood PDFs, and b 2D projection of the final posterior joint PDF indicating high crack probability at low σVM and high \(\dot{\gamma }\) values. c Log-odds of the geometric mean of the probability of correct crack nucleation prediction for increasingly model complexity. The included models are \({{{{\mathcal{M}}}}}_{0}=p(C)\), \({{{{\mathcal{M}}}}}_{\xi }=p(C\,| \,\xi )\), \({{{{\mathcal{M}}}}}_{{{\mbox{r}}}}=p(C\,| \,\xi ,{{{{\bf{x}}}}}_{{{\mbox{r}}}})\), and \({{{{\mathcal{M}}}}}_{* }=p(C\,| \,\xi ,{{{{\bf{x}}}}}_{* })\).

The posterior PDF, when compared with other existing crack nucleation models of varying complexity, demonstrates excellent improvement. The original prior \({{{{\mathcal{M}}}}}_{0}=p(C)\), the prior after conditioning on parallel slip configuration regions \({{{{\mathcal{M}}}}}_{\xi }=p(C\,| \,\xi )\), and additional posteriors including an increasing number of state variables, are ranked according to their predictive capacity on the training data set. In Fig. 4c, the log-odds of each PDF’s geometric mean probability ρC (in bits) is used as the measure of comparison and is calculated as \({{{\mbox{log}}}}_{2}\left(\frac{{\rho }_{C}}{1-{\rho }_{C}}\right)\). The log-odds gives a natural interpretation of the strength of evidence gained by adding model complexity to a probabilistic model52. A unit increase in log-odds corresponds to an additional bit of evidence. This is roughly interpreted as observing 2:1 evidence in support of the model’s hypothesis, i.e. that it correctly predicts crack nucleation sites. In this context, p(Cξ) adds about 6.5 bits of evidence over the naive p(C), by requiring that crack nucleation occurs at material volumes in a parallel slip configuration. Furthermore, \({{{{\mathcal{M}}}}}_{{{\mbox{r}}}}=p(C\,| \,\xi ,{{{{\bf{x}}}}}_{{{\mbox{r}}}})\) adds about 1.5 bits of evidence on top of p(Cξ), by incorporating the top three state variables in the posterior PDF. The model therefore greatly improves the probability of locating a rare crack nucleation event, by about 8 bits relative to the original prior from the data.

From a mechanistic point of view, an interpretation is offered for the model’s emphasis on σVM, \(\dot{\gamma }\), and ψp in terms of the experimental evidence supporting long, thin cracks in PSB regions at parallel slip configurations. The von Mises stress, which is a norm of the deviatoric stress, typically represents a driving force for the volume-preserving shear deformation such as plasticity and dislocation motion. However, since the material volume is already in a parallel slip configuration, a high directionally independent measure of shear stress, such as σVM, indicates that dislocation slip is likely to occur on multiple slip systems, rather than focusing on only the parallel slip plane. This importance of directionality is further evidenced by a large \(\dot{\gamma }\) measure, also acting as a crack-delineating state variable. The combination of the low σVM and high \(\dot{\gamma }\) is an environment that promotes dislocation motion along the primary slip plane, which is the parallel slip plane in a parallel slip configuration. It represses slip along with other secondary slip systems that interferes with the main deformation mode. The plastic defect energy as the third major indicative variable of crack nucleation is linked to the build-up of dislocation substructures that are typical of persistent slip band (PSB) regions. Regions of large ψp store more external work into complex and organized arrays of tangled dislocations, rather than dissipating it entirely as plastic work.

After constructing the full posterior distribution p(Cxr), probability fields are overlaid on the microstructures to determine the most likely crack nucleation sites. A sample polycrystalline microstructure is shown in Fig. 5 with the corresponding orientation map from the EBSD scan. The probability contours by the Bayesian model concur with the experimental crack nucleation sites, as probable locations for observing a crack.

Fig. 5: Comparison of the crack nucleation model predictions with the observed crack locations of an experimental microstructure.
figure 5

The probability of crack nucleation predicted from the posterior distribution in an embedded microstructure of the concurrent multiscale model (red corresponding to the most probable), compared with the corresponding EBSD map and experimental nucleation site.

Cross-validation

The predictive capacity of the probabilistic model is assessed in this section. Delineation of the training and test data allows for a comparison of the performance of the methodology in terms of both training error and generalization error. Model over-fitting is a concern when the only efficacy measure of a model is on the data it is trained on. Therefore, a “leave-two-out” cross-validation approach is used to ensure the crack nucleation model generalizes onto unseen data. The leave-two-out method is a standard, exhaustive cross-validation approach, in which the model is trained on every available set of data except for two microstructures. The model is then tested/validated on the two microstructures that are left out from the training set. This process is repeated, leaving two different data sets out each time, and retraining the model on the remaining ten microstructures. This approach is one of the best for small data sets and for training a computationally expensive model, as it increases the amount of training data available for each validation evaluation.

For the case of M = 12 microstructures, there are M(M − 1)/2 = 66 ways to leave two unique data sets out and train on the remaining ten. Therefore, the model is trained and tested 66 times for all combinations. The geometric mean of correct crack nucleation prediction ρC is used as a measure of the effectiveness of the model in both training and validation. An additional related measure is defined as:

$${\rho }_{\bar{C}}=\exp \left(\mathop{\sum }\limits_{i}^{{N}_{\bar{C}}}{w}_{\bar{C}}^{(i)}\log \left(1-p(C\,| \,{{{{\bf{x}}}}}^{(i)})\right)\right)$$
(2)

where \({N}_{\bar{C}}\) is the total number of uncracked material volumes across all microstructures of interest, and \({w}_{\bar{C}}^{(i)}\) are the normalized volumes for each. This measure is the geometric mean of the probability that the model correctly labels every known uncracked material volume of the observed data set.

The arithmetic mean and geometric mean of ρC and \({\rho }_{\bar{C}}\) over all 66 cases are given in Table 1 for both the training and testing data sets, by employing the leave-two-out method. The resulting evaluation measures indicate that the probabilistic model yields similar abilities for both the training and testing sets. Only the geometric mean for the validation cracked cases is slightly lower than that of the test case. This deviation is due to the geometric mean’s sensitivity to variance and bias towards low values, which are higher in the testing examples. Overall, these results demonstrate that the model generalizes well, and is predictive on new data.

Table 1 Mean probabilities of nucleation prediction for all “leave-two-out” cases.

The predictive behavior of the model is further emphasized by plotting the posterior distribution over a grain structure for one of the 66 leave-two-out cases. Fig. 1 in the Supplementary Methods Section shows the probability field and the corresponding EBSD map over one of the two cases that are not included in the training. The probability field illustrates likely locations for crack nucleation, and in particular, highlights the locations near the experimental crack nucleation site. These regions are comprised of material volumes in a parallel slip configuration with high propensity for primary slip and have already undergone substantial defect build-up.

Discussion

The probabilistic crack nucleation model consisting of a validated three-dimensional posterior PDF can be used to create a simple, heuristic scalar indicator model of crack nucleation probability. A one-dimensional functional form is constructed to become a separator of the cracked and uncracked regions. The dimensionless crack nucleation indicator is defined at material volumes in a parallel slip configuration (Ξ = ξ) as:

$${{{\Psi }}}_{{{\mbox{c}}}}={\kappa }_{1}\left(\frac{\log \left(\frac{\dot{\gamma }}{\dot{{\gamma }_{0}}}\right)}{{\sigma }_{\,{{\mbox{VM}}}}^{2}}\right)+{\kappa }_{2}{\psi }_{{{\mbox{p}}}}$$
(3)

where \(\dot{\gamma }\) is the slip-rate, \(\dot{{\gamma }_{0}}=5\times 1{0}^{7}\) is the reference slip-rate, ψp is the plastic defect energy, and κ1, κ2 are parameters to be determined. The crack nucleation indicator incorporates the basic mechanistic ideas discussed in the previous section. The first term of Ψc corresponds to a state of material which heavily favors parallel slip along a boundary rather than homogeneous plastic deformation across multiple slip systems. The numerator \(\log \left(\frac{\dot{\gamma }}{\dot{{\gamma }_{0}}}\right)\) is approximately proportional to the effective RSS driving slip through thermal activation, and the denominator \({\sigma }_{\,{{\mbox{VM}}}\,}^{2}\) is proportional to the elastic distortion energy. The second term of Ψ, often used in the Helmholtz energy representation of phase-field crack models53, indicates that the material has a tendency to build up dislocation substructures to promote the onset of persistent slip band formation. The proposed crack nucleation indicator encompasses and combines features of fatigue indicator parameters proposed in the literature, such as the emphasis on the plane of maximum plastic slip in24 and the significance of the stored plastic defect energy in23.

The two parameters κ1 and κ2 of Ψc are linear coefficients with respect to the inputs of the function. Correspondingly they can be determined using linear discriminant analysis (LDA). LDA is a classification method that finds the linear combination of inputs that best separates two distinct distributions, to find a reduced order model. The application of LDA to the crack nucleation data yields explicit values for the linear coefficients of Ψc, i.e. κ1 = − 0.28  10−18 Pa−2 and \({\kappa }_{2}=0.26\cdot 1{0}^{-6}{\left(\frac{{{\mbox{J}}}}{{{{\mbox{m}}}}^{3}}\right)}^{-1}\). This process results in a single scalar value for crack nucleation prediction and is used as a sole indicator in the probabilistic model p(Cx), where x = [Ψc].

Visual analysis of the likelihood distributions of Ψc in Fig. 6a establishes it as a useful separator of the cracked and uncracked experimental data. Furthermore, the model in Eq. 3 is assessed similar to the previous posteriors by calculating ρC. Applying the cracked data to the model p(C Ψc) yields ρC = 11.8%, which is superior to building a model with any of the other previously evaluated state variables individually. Consequently, the crack nucleation indicator, which embodies the mechanistic motivation of parallel slip within a persistent slip band, is a consistent measure for identifying cracked regions of the microstructures.

Fig. 6: Performance and robustness of the scalar crack nucleation indicator.
figure 6

a The likelihood distributions of the scalar crack nucleation indicator for cracked and uncracked states, b aggregated posterior distribution of the 66 models, conditioned on only the scalar crack nucleation indicator, for different combinations of ten microstructure in the training data sets. The average of all models is depicted by the dark solid line and the standard deviation bounds are depicted by the dotted lines.

The contour plots of the crack nucleation indicator field are shown for four polycrystalline microstructures in Fig. 2 of the Supplementary Methods Section. This enables a visual comparison to assess its performance with respect to predicting crack nucleation sites. Qualitatively, large values of Ψc highlight areas associated with the experimentally identified nucleation sites for each case.

A final study is performed to quantify the sensitivity of the scalar-valued probabilistic model development to the inclusion of new data. The 66 cross-validation cases are revisited for this analysis. For each case, the posterior model is constructed by training with data from ten microstructures, using only the data from the crack nucleation indicator Ψc. The posterior p(CΨc) is generated and plotted as a function of Ψc for each case. In Fig. 6b, the 66 models are aggregated into an average model shown with the solid dark line, and standard deviation bounds depicted with dotted lines. This result demonstrates that most models lie within about 5–10% of the average model, and therefore are robust and relatively insensitive to the choice of training data set.

In summary, a robust Bayesian inference-based probabilistic predictive model for crack nucleation in the Ni-based superalloy René 88DT under fatigue loading is developed in this paper. The process fuses several statistical and data-driven methods operating on micromechanical and experimental data with important mechanistic observations. The model systematically reduces a large set of experimentally acquired microstructural data and corresponding crystal plasticity-based micromechanical simulation results to a 3D posterior probability density function (PDF) for crack nucleation. Three hierarchically important state variables in the microstructural material volume, relevant to the crack nucleation process, are identified by the framework. They are: (i) the von Mises stress, (ii) the maximum plastic slip rate, and (ii) the plastic defect energy. This data-driven selection process is performed in an automated fashion without the need for external intervention. It follows a statistically unbiased approach that makes an informed decision while accounting for various experimental microstructures and mechanical simulations.

Given a set of EBSD scans along with the associated pre-identified crack nucleation sites as input, the procedure ranging from the processing of experimental EBSD and SEM data, multiscale simulations for creating a reference database, to the development of the probabilistic crack nucleation model is implemented through an automated computational pipeline. It comprises four sequential modules given below.

  1. 1.

    Collection of post-fatigue loading local microstructural data, and subsequent cleaning and processing of EBSD data;

  2. 2.

    Development of a concurrent multiscale model with the image-based 3D microstructure embedded in a homogenized material modeled by an upscaled, anisotropic continuum plasticity constitutive model;

  3. 3.

    Multiscale simulations with different experimental microstructures to create a database of micromechanical state variables relevant to the crack nucleation process;

  4. 4.

    Formulation of a Bayesian crack nucleation model that identifies crack nucleation sites from the database of local state variable fields.

The probabilistic models for indicating crack nucleation are based on the core concept of parallel slip at susceptible coherent twin boundaries. The material near these boundaries possesses particular crystallographic configurations that promote an environment for recurring dislocation slip parallel to the boundary. State variables, such as the von Mises stress, the maximum plastic slip rate, and the plastic defect energy in this region, have a strong effect on this mechanism. Slip is driven on the primary plane and builds up dislocation substructures like PSBs that further promote future slip along the same path. Simultaneously, little driving force exists for forest dislocations to impede obstacles on the primary plane. The combination of these events creates a probable setting for the incidence of crack nucleation.

The multi-dimensional Bayesian inference model is constructed to capture this sequence of mechanistic events. It is validated for various microstructures and probabilistic measures, using the “leave-two-out” cross-validation process. Finally, a simple crack nucleation indicator is formulated to integrate the parallel slip concepts into a single scalar model. This composite indicator outperforms other individual state variables in differentiating between cracked and uncracked material volumes. In conclusion, this paper provides a robust framework for employing sophisticated data analysis techniques in conjunction with multiscale simulations to investigate the mechanisms underlying complex physical phenomena, and produces interpretable results. The Bayesian crack nucleation model should be valid for different microstructural compositions and loading conditions, as long as the main mechanisms of deformation such as plastic slip, dislocation glide, etc., remain the same.

The overall modeling framework developed in this paper is not limited to only superalloys considered in this study. Most aspects of the framework are transferable to other polycrystalline alloys for which, state variable evolution corresponding to parallel slip and PSB formation govern their fatigue behavior. A majority of steps in the automated framework, including the embedding methodology, finite element meshing, state variable calculation, probabilistic model procedure, and the reduced model selection technique, are material agnostic within the realm of metallic materials. Transferring the microstructure to deformation behavior and fatigue platform to other metallic materials will require the corresponding experimental data, viz. EBSD scans and crack nucleation sites, a calibrated crystal plasticity constitutive model, and the associated upscaled anisotropic continuum plasticity model. Additionally, the probabilistic classification framework has many other potential applications in data-rich and knowledge-poor circumstances.

For component-scale analysis, it is necessary to introduce a detailed parametrically upscaled constitutive model (PUCM) and associated parametrically upscaled crack nucleation model (PUCNM) developed in39,40,41,42,43. The thermodynamically consistent PUCMs, incorporate a parametric representation of lower-scale microstructural descriptors in higher-scale constitutive coefficients. The PUCM coefficients are expressed as functions of representative aggregated microstructural parameters (RAMPs), representing lower-scale descriptors of microstructural morphology and crystallography. While the upscaled anisotropic continuum plasticity model developed for the concurrent model is a step in this direction, it does not have the parametric forms in terms of the RAMPs. The PUCM, in conjunction with the PUCNM, readily performs component-scale analysis with direct connection to the microstructure. The models developed in this study can be used for SERVE analysis to develop the PUCM and PUCNM for the superalloy considered for this study.

Methods

Concurrent multiscale model for fatigue analysis

The concurrent multiscale model in Fig. 7 consists of a polycrystalline microstructural RVE that is embedded in an exterior computational domain. The RVE and the exterior domain are modeled by a crystal plasticity constitutive model34,35,36 and a rate-dependent homogenized elasto-plasticity constitutive model respectively. Cyclic loads representing fatigue loading are applied on the outer surface of the homogenized exterior domain. The constitutive relations used in the concurrent multiscale model are summarized in the Supplementary Methods Section of this paper.

Fig. 7: Schematic of the concurrent computational domain.
figure 7

Microstructural RVEs are embedded in homogenized exterior domains.

Concurrent multiscale model with embedded microstructural RVE

The application of periodicity or other boundary conditions directly on the microstructural RVE domain can result in significant inaccuracies due to improper representation of the domain and consequent boundary effects32,33,54. The concept of concurrent multiscale domains, in which the RVE is embedded in a homogenized exterior domain, has been successfully introduced in43,55,56 to alleviate this discrepancy. Boundary conditions are applied on the exterior boundary of the concurrent domain, while the RVE shares a common interface with the homogenized domain. Similar concurrent approaches have also been proposed for coupled continuum-polycrystalline domains in29,57.

The embedded RVE in Fig. 7 is generated by directly mapping the 2D EBSD micrograph of a polycrystalline microstructure shown in Fig. 2, and extruding it in the orthogonal direction to a suitable depth. Each section of the RVE contains about 50–200 grains. Previous studies on Ni-based superalloys in34,54,58 have demonstrated the effectiveness of microstructural RVEs of this size. Some of the EBSD micrographs containing cracked regions are not fully supported by first and second nearest neighbor grains. While this may lead to some error in the simulation results due to limited grain representation in RVE models, the relatively lower strains levels may not significantly alter the response fields as noted in59. The homogenized exterior medium in the concurrent model may also mitigate the surrounding boundary effect due to influence from neighboring grains. Furthermore, the final posterior model is constructed using multiple experimental microstructures and testing. With increasing data, the Bayesian model improves at distinguishing between all the sources of model noise and the meaningful relation existing between crack nucleation sites and state variable fields.

For fatigue crack modeling, it is necessary that the microstructural RVEs be constructed directly from EBSD scans of the region containing cracks for one-to-one correspondence with experimental observations. Statistically equivalent RVE (SERVEs), with matching morphological and crystallographic distributions, are not ideal for this one-to-one correspondence. 3D microstructures may consequently be created by focused ion beam (FIB)–based serial-sectioning, or synchrotron-based CT methods. However, it is challenging to obtain interrupted 3D experimental microstructures in the context of fatigue experiments. The decision to create 3D microstructural RVEs by extruding 2D EBSD scans was made since such types of data were not available for this study. In60, the authors have inferred that sub-surface crystallography and grain morphology, are needed to accurately predict surface strain distributions. In the present study, it is assumed the error due to the extrusion will not be significant since the thickness of the RVE is smaller than the average grain size. Additionally, the use of multiple microstructural instantiations may have a mitigating effect on the crack nucleation model.

The concurrent multiscale model can mitigate other spurious boundary effects on the evolving RVE variables. A sample microstructural RVE, shown in Fig. 7, has dimensions of 168 μm × 69 μm × 7.5 μm, while the homogenized exterior rectangular domain has dimensions of 336 μm × 138 μm × 42 μm. This corresponds to twice the length and width of the RVE and a depth of 7.5 μm plus half of the RVE width. Five of the six faces of the embedded RVE share an interface with the exterior homogenized domain, enforcing displacement compatibility across these interfaces. The remaining RVE face represents the free surface of the fatigue experiment and is unconstrained during the cyclic simulations.

The concurrent FE model is discretized into TET4 elements, generated using an automated meshing pipeline built within the Simmetrix code61. This minimizes the computational and memory requirements and provides a robust discretization in terms of numerical accuracy and geometric conformity. Since the state variables near twin boundaries are key to the crack model, a graded mesh with highly refined elements near the grain/twin boundaries, and coarser elements away from these critical regions, is developed as seen in Module II of Fig. 1. The exterior domain mesh matches the RVE mesh at the interface, but coarsens away from the interface.

Creating a microstructural state variable database related to fatigue nucleation

The following steps are executed for creating the database.

Experimentally obtained training data sets for supervised learning

The Bayesian modeling approach to fatigue crack nucleation requires a data set for supervised learning, consisting of experimentally observed points for which the binary cracked/uncracked states are known. To generate this data set, the locations of crack nucleation sites are recorded for the experimental micrographs, e.g. in Fig. 2. Each TET4 element in the concurrent model, for which three nodes are on a cracked interface is marked as belonging to the set of cracked elements. All other elements in the ensemble are designated as uncracked. This labeling method classifies the data set into two distinct states for training the probabilistic model.

Extracting core microstructural state variables for the crack nucleation model

The concurrent multiscale model in Methods Section I is now developed for each of the 12 experimental microstructures and subsequently loaded under fatigue loading conditions. Analogous to the experiments in51, a stress-controlled fully reversed (R = −1) cyclic loading with a triangular wave-form and 1 Hz frequency is applied on the computational model. The peak stress in this loading cycle corresponds to a 1% total strain. All state variables are recorded at the peak tensile load of each cycle. This process generates a complete database of mechanical fields for all the 12 experimental microstructures, which is used for training the Bayesian crack nucleation model.

Time-histories of state variables are calculated and recorded as quantities of interest in the development of the fatigue crack nucleation model. From a comprehensive set of mechanistic variables, including stress and strain measures, dislocation evolution indicators, fatigue indicator parameters, etc., a set of ten core state variables are deduced. Most of the other indicators are functions of these core parameters. The core variables are tabulated in Table 2.

Table 2 List of core state variables used in the probabilistic crack nucleation analysis.

The set of variables in Table 2 is designated as X. Additionally, a regularized difference operator (similar to unsharp masking) acting on each state variable is appended to the set X. Correspondingly the dimension of the set is k = dim(X) = 20. The values of X are assembled into a state variable database for use in the training and testing of the probabilistic model.

Observations on parallel slip configuration from experiments

Observations from fatigue experiments in9 suggest that dislocation slip parallel to the coherent twin boundaries results in persistent slip band (PSB) and thin shear band formation due to strain localization. Adjacent to the twin boundary, plasticity occurs due to dislocation glide on slip planes parallel to the twin boundary. The repeated passing of these dislocations and migration of dislocations to this region, shears the \(\gamma ^{\prime}\) precipitates, promoting additional plasticity by reducing the barrier to subsequent dislocation motion. It has been observed that these highly localized regions of plasticity nucleate cracks, several nanometers offset from the actual twin boundary.

Some twin boundaries are identified by grains that are oriented such that their slip systems favor dislocation motion along their shared boundary plane. This special grain-pair configuration is called a parallel slip configuration. It is characterized by a grain pair separated by a coherent twin boundary, where the resolved shear stress on one of the shared slip systems is significantly higher than that of the other slip systems. This configuration results in primary slip along that single shared plane and no new dislocation migration occurs from other planes to the slipped region. The twin boundary lengths are also relatively long and straight compared to other grain boundaries. Consequently, this geometric effect and the lack of out-of-plane dislocation motion leaves the shared slip system free of obstacles, which allows dislocations to glide unimpeded for long distances. In cyclic deformation, uninterrupted localized slip promotes the early stages of PSB formation.

Indicator of dominant parallel slip

Since parallel slip configurations are identified as precursors to crack nucleation, it is important to develop an indicator for dominant parallel slip. This indicator should identify coherent twin boundaries and subsequently search for grain-pairs, oriented such that the resolved shear stress favors slip along the shared boundary plane. The following procedure is adopted to account for these conditions and lead up to a definition for the parallel slip indicator.

The criteria used to determine if two grains are in a Σ3 twinned configuration is based on two conditions. The first condition is that the misorientation angle must be 60 ± 2, and the second requires that the absolute value of the inner product between the shared slip plane normals is >0.985. For each grain pair, the shared slip plane normal is defined by the unit vector average of the two slip plane normals, whose inner product has the largest absolute value. Specifically, these vectors are determined by the following process.

For two adjacent grains g1 and g2 in the microstructure, the two respective slip planes that most closely share the same normal are selected according to the criterion:

$$({l}^{* },{k}^{* })=\arg \mathop{\max }\limits_{l,k}(| {{{{\bf{n}}}}}_{l}^{{g}_{1}}\cdot {{{{\bf{n}}}}}_{k}^{{g}_{2}}| )$$
(4)

where l and k enumerate the four possible slip planes in the 〈111〉 family. Correspondingly, the shared slip plane normal for the grains g1 and g2 are defined as:

$${{{{\bf{n}}}}}_{\,{{\mbox{p}}}\,}^{{g}_{1},{g}_{2}}=\frac{{{{{\bf{n}}}}}_{l* }^{{g}_{1}}+{{{{\bf{n}}}}}_{k* }^{{g}_{2}}}{| | {{{{\bf{n}}}}}_{l* }^{{g}_{1}}+{{{{\bf{n}}}}}_{k* }^{{g}_{2}}| | }$$
(5)

where \({{{{\bf{n}}}}}_{l* }^{{g}_{1}}\) is the normal to the l* slip plane of grain g1 and \({{{{\bf{n}}}}}_{l* }^{{g}_{2}}\) is the normal to the k* slip plane of grain g2.

The existence of a shared plane allows for a partition of the standard octahedral slip systems of FCC materials. Let \({{{\mathcal{A}}}}\) denote the set of all its slip systems for a grain that shares a slip system with an adjacent grain. Also, let \({{{{\mathcal{A}}}}}_{{{\mbox{p}}}}\) be the set of slip systems that contain the shared slip plane characterized by \({{{{\bf{n}}}}}_{\,{{\mbox{p}}}\,}^{{g}_{1},{g}_{2}}\) in Eq. 5. Finally, let \({{{{\mathcal{A}}}}}_{{{\mbox{o}}}}\) be the complement of \({{{{\mathcal{A}}}}}_{{{\mbox{p}}}}\) in \({{{\mathcal{A}}}}\), i.e. the complementary set of slip systems that do not contain the shared slip plane.

Two special resolved shear stresses (RSS) are defined to determine when parallel slip is dominant for a given twin boundary. The maximum RSS is given as \({\tau }_{* }=\mathop{\max }\limits_{\alpha \in {{{\mathcal{A}}}}}\left(| {{{\bf{M}}}}:{{{{\bf{s}}}}}_{0}^{\alpha }| \right)\), where \({{{{\bf{s}}}}}_{0}^{\alpha }={{{{\bf{m}}}}}_{0}^{\alpha }\otimes {{{{\bf{n}}}}}_{0}^{\alpha }\) is the Schmid tensor for the α slip system, \({{{{\bf{m}}}}}_{0}^{\alpha }\) is the slip direction, and \({{{{\bf{n}}}}}_{0}^{\alpha }\) is the slip plane normal in the reference configuration. Correspondingly, the parallel slip RSS and the out-of-plane RSS are respectively defined as:

$${\tau }_{{{\mbox{p}}}}=\mathop{\max }\limits_{\alpha \in {{{{\mathcal{A}}}}}_{{{\mbox{p}}}}}\left(| {{{\bf{M}}}}:{{{{\bf{s}}}}}_{0}^{\alpha }| \right)\quad \,{{\mbox{and}}}\,\quad {\tau }_{{{\mbox{o}}}}=\mathop{\max }\limits_{\alpha \in {{{{\mathcal{A}}}}}_{{{\mbox{o}}}}}\left(| {{{\bf{M}}}}:{{{{\bf{s}}}}}_{0}^{\alpha }| \right)$$
(6)

The parallel slip RSS is the maximum RSS on the shared slip plane, while the out-of-plane RSS is the maximum over the remaining three slip planes. When τp > τo, primary slip occurs along the twin boundary, and out-of-plane dislocations do not disrupt this plastic flow. Therefore, parallel slip dominates the deformation process. This condition is combined with the coherent twin and a plastic flow requirement to define a binary parallel slip indicator Ξ as \(\xi \,or\,\bar{\xi }\) for a given material volume. The parallel slip indicator is activated with Ξ = ξ, if all of the following three requirements are met simultaneously:

  1. 1.

    The material volume is located on a coherent twin boundary;

  2. 2.

    The parallel slip plane is dominant, i.e. τp > τo;

  3. 3.

    The maximum plastic slip-rate exceeds a minimum threshold, i.e. \(\dot{\gamma } > {\dot{\gamma }}_{{{\mbox{m}}}}\), where \({\dot{\gamma }}_{{{\mbox{m}}}}\) is a minimum threshold slip-rate.

The threshold slip-rate is selected as \({\dot{\gamma }}_{{{\mbox{m}}}}=1{0}^{-5}{s}^{-1}\) to indicate that the material is undergoing plasticity. If any one of the three conditions are not met, then \({{\Xi }}=\bar{\xi }\). The parallel slip indicator is designed to signal time-evolving, point-wise occurrence of parallel slip throughout the microstructure. It is analogous to the time-independent, grain-wise-homogeneous parallel slip configuration definition in1,8,9,15. The importance of Ξ is that it acts as a sufficiency condition for crack nucleation within the experimental data set of microstructures. This empirical data is exploited in the probabilistic model framework.

Formulation of the Bayesian model for classifying crack nucleation

Let a discrete random variable Θ represent the possible observation of a crack event in a fixed, reference material volume V*, and let a continuous k-dimensional random vector X denote the set of core state variables pertaining to that particular region. The entire microstructural field is classified by one of two distinct states, viz. (i) with a crack event Θ = C, and (ii) without a crack event \({{\Theta }}=\bar{C}\). Additionally, let the previously defined parallel slip indicator Ξ in the microstructure be a discrete random variable with two possible states, viz. (i) Ξ = ξ indicating a parallel slip configuration, and (ii) \({{\Xi }}=\bar{\xi }\) otherwise. In accordance with experimental observations, it is assumed that all cracked states C require that a material volume is in a parallel slip configuration Ξ = ξ. Therefore, the probability of observing a crack nucleation event C in the microstructure with given state variables x, homogeneous over some small volume V*, is given by Bayes’ theorem. For the case of Ξ = ξ, this is stated as:

$$P({{\Theta }}=C\,| \,{{{\bf{X}}}}={{{\bf{x}}}})={p}_{{{\Theta }}| {{{\bf{X}}}}}(C\,| \,{{{\bf{x}}}})={p}_{{{\Theta }}| {{\Xi }},{{{\bf{X}}}}}(C\,| \,\xi ,{{{\bf{x}}}})=\frac{{p}_{{{{\bf{X}}}}| {{\Theta }},{{\Xi }}}({{{\bf{x}}}}\,| \,C,\xi ){p}_{{{\Theta }}| {{\Xi }}}(C\,| \,\xi )}{{p}_{{{{\bf{X}}}}| {{\Xi }}}({{{\bf{x}}}}\,| \,\xi )}$$
(7)

where pΘX(Cx) and pΘ  Ξ,X(Cξ, x) are both posterior probability density functions (PDF’s), pX Θ,Ξ(xC, ξ) is a likelihood function, pΘ  Ξ(Cξ) is the prior belief of observing a crack nucleation event, and pX Ξ(xξ) is the marginal state variable PDF, which normalizes the expression. For the complementary event \({{\Xi }}=\bar{\xi }\), the posterior PDF pΘX(Cx) = 0. Hereafter, the subscripts on the PDF’s will be omitted to eliminate notational clutter, but may be used occasionally for providing beneficial clarity. Applying the law of total probability to Eq. 7 yields:

$$p(C\,| \,{{{\bf{x}}}})=\frac{p({{{\bf{x}}}}\,| \,C,\xi )p(C\,| \,\xi )}{p({{{\bf{x}}}}\,| \,C,\xi )p(C\,| \,\xi )+p({{{\bf{x}}}}\,| \,\bar{C},\xi )p(\bar{C}\,| \,\xi )}$$
(8)

The posterior PDF p(Cx) is the primary quantity of interest to be determined for the crack nucleation model. From Eq. 8, it is seen that four PDF’s must be calculated to fully characterize this posterior, viz. p(xC, ξ), \(p({{{\bf{x}}}}\,| \,\bar{C},\xi )\), p(Cξ), and \(p(\bar{C}\,| \,\xi )\). Each of these distributions is systematically constructed from the experimental data set in the following sections.

Determination of the prior distributions \(p({C}\,| \,\xi )\) and \(p(\bar{C}\,| \,\xi )\)

The prior belief that a crack event occurs in a material volume is constructed by calculating the ratio of cracked volume to total volume over all experimental microstructures. For the available experiments, let M = 12 denote the number of independent microstructures in the data set, and let \({V}_{C}^{(i)}\) and \({V}_{\,{{\mbox{tot}}}\,}^{(i)}\) be respectively the total cracked volume and total overall volume of the ith microstructure, where Ξ = ξ. The crack nucleation prior is correspondingly defined as:

$$p(C\,| \,\xi )=\frac{\mathop{\sum }\nolimits_{i}^{M}{V}_{C}^{(i)}}{\mathop{\sum }\nolimits_{i}^{M}{V}_{\,{{\mbox{tot}}}\,}^{(i)}}$$
(9)

The uncracked prior is constrained by the complementary event, and is therefore determined as \(P(\bar{C}\,| \,\xi )=1-P(C\,| \,\xi )\). Additionally, let the reference volume for the model be defined as:

$${V}_{* }=\frac{1}{M}\mathop{\sum }\limits_{i}^{M}{V}_{C}^{(i)}$$
(10)

which corresponds to the average crack volume over all microstructures.

Determination of the likelihood distributions \(p({{{\bf{x}}}}\,| \,{C},\xi )\) and \(p({{{\bf{x}}}}\,| \,\bar{C},\xi )\)

In general, the likelihood PDF’s, p(xC, ξ) and \(p({{{\bf{x}}}}\,| \,\bar{C},\xi )\), are high-dimensional, highly covariant joint distributions. The Nataf transformation method62,63,64 is used in this work to construct explicit functional representations of both p(xC, ξ) and \(p({{{\bf{x}}}}\,| \,\bar{C},\xi )\). It has been shown through Sklar’s theorem that a multivariate joint distribution can be decomposed into a composition of component-wise univariate marginal distributions and a multivariate copula, which retains the covariance structure of the random vector. For the Nataf transformation, the joint distribution is approximated using a Gaussian copula to permit a straightforward analytical computation. The following construction is formulated in the context of p(xC, ξ). However for \(p({{{\bf{x}}}}\,| \,\bar{C},\xi )\), analogous steps are performed by substituting \(\bar{C}\) for C.

The Nataf transformation is comprised of three sequential transformations. First, an observed data point x = [x1, x2, . . . ,xk] of the random set X of state variables is transformed component-wise into a vector of uniformly distributed random variables u = [u1, u2, . . , uk], where

$${u}_{i}={F}_{{X}_{i}| {{\Theta }},{{\Xi }}}({x}_{i}\,| \,C,\xi )$$
(11)

\({F}_{{X}_{i}| {{\Theta }},{{\Xi }}}\) is the marginal cumulative distribution function (CDF) for the ith component of X. Secondly, each component of u is standardized by applying the inverse standard normal CDF, specified by a mean of zero and a variance of one. Let the standardized vector z = [z1, z2, . . . , zk] be characterized by this transformation as:

$${z}_{i}={{{\Phi }}}^{-1}\left({u}_{i}\right)$$
(12)

Finally, a k-dimensional normal distribution is fit to the observations in the standardized space. This sequential application of the ingredients of the Nataf transformation is summarized as:

$$F({{{\bf{x}}}}\,| \,C,\xi )={{{\Phi }}}_{k}({{{\bf{z}}}})={{{\Phi }}}_{k}({{{\Phi }}}^{-1}({u}_{1}),{{{\Phi }}}^{-1}({u}_{2}),...\,,{{{\Phi }}}^{-1}({u}_{k}))$$
(13)

where the joint CDF F(xC, ξ) corresponds to the likelihood p(xC, ξ) and Φk is the k-dimensional multivariate normal CDF. Differentiating Eq. 13 yields:

$$p({{{\bf{x}}}}\,| \,C,\xi )=p({{{\bf{z}}}}\,| \,C,\xi )\mathop{\prod }\limits_{i}^{k}\frac{{p}_{{X}_{i}| {{\Theta }},{{\Xi }}}({x}_{i}\,| \,C,\xi )}{\phi ({z}_{i})}$$
(14)

where \({p}_{{X}_{i}| {{\Theta }}}({x}_{i}\,| \,C)\) is the marginal PDF for the ith component, ϕ is the standard normal distribution, and p(zC) is a multivariate normal distribution ϕk given as:

$$p({{{\mathbf{z}}}}\,| \,C,\xi )={\phi }_{k}({{{\mathbf{z}}}})=\frac{1}{\sqrt{{(2\pi )}^{k}\det {{{\mathbf{\Sigma }}}}}}\exp \left(\frac{1}{2}{{{\mathbf{z}}}}\cdot \left({{{{\mathbf{\Sigma }}}}}^{-1}{{{\mathbf{z}}}}\right)\right)$$
(15)

with the covariance matrix Σ. The computation of p(xC, ξ) for the state variable data set only requires the determination of \({F}_{{X}_{i}| {{\Theta }},{{\Xi }}}\), \({p}_{{X}_{i}| {{\Theta }},{{\Xi }}}\), and Σ. \({F}_{{X}_{i}| {{\Theta }},{{\Xi }}}\) and \({p}_{{X}_{i}| {{\Theta }},{{\Xi }}}\) are represented by a regularized, weighted empirical CDF and a weighted kernel density estimation from the observed x data, respectively. The covariance matrix Σ, defined as: \({{{\mathbf{\Sigma }}}}={\mathbb{E}}[{{{\mathbf{z}}}}\otimes {{{\mathbf{z}}}}]\), is approximated by a weighted unbiased estimator of the observed data as:

$${{{\mathbf{\Sigma }}}}=\frac{1}{1-\omega }\mathop{\sum }\limits_{i}^{{N}_{{{\mbox{e}}}}}{w}^{(i)}{{{{\mathbf{z}}}}}^{(i)}\otimes {{{{\mathbf{z}}}}}^{(i)}$$
(16)

where Ne is the total number of discretized volumes (finite elements) over all M microstructures in a state of Ξ = ξ, w(i) are observation weights corresponding to the normalized material volume of each point such that \(\mathop{\sum }\nolimits_{i}^{{N}_{{{\mbox{e}}}}}{w}^{(i)}=1\), and \(\omega =\mathop{\sum }\nolimits_{i}^{{N}_{{{\mbox{e}}}}}{w}^{(i)2}\) is the sum of squares of the observation weights.

An example of the action of the Nataf transformation for the uncracked likelihood PDF \(p({{{\bf{x}}}}\,| \,\bar{C},\xi )\) is shown in Fig. 8a for two representative state variables. The von Mises stress and maximum plastic slip-rate data points \({{{\bf{x}}}}=[{\sigma }_{{{\mbox{VM}}}},\dot{\gamma }]\), are depicted in the original X space, transformed into the uniform U space, and then transformed again into the standardized Z space. Additionally, grid lines are added and transformed accordingly to illustrate the component-wise warping of each observation space.

Fig. 8: Illustration of the Nataf transformation for two variables.
figure 8

a Example action of the Nataf transformation for the bivariate \({{{\bf{x}}}}=[{\sigma }_{{{\mbox{VM}}}},\dot{\gamma }]\). Approximately 500,000 data points from all 12 microstructures are transformed from the observed X space to the marginally uniform U space to the standardized Z space. b Iso-probabilistic contours of the bivariate normal model fit for the uncracked likelihood \(p({{{\bf{z}}}}\,| \,\bar{C})\) in the bivariate example case of \({{{\bf{x}}}}=[{\sigma }_{{{\mbox{VM}}}},\dot{\gamma }]\) and c Iso-probabilistic contours of the model fit for \(p({{{\bf{z}}}}\,| \,\bar{C})\) by transforming the standardized model back to the X space.

Figure 8b demonstrates the application of the model described in Eq. 15 to determine a bivariate normal distribution as an approximation for the standardized observations z. This approximation is achieved through the computation of Σ, using Eq. 16. It is seen that the preservation of the iso-probabilistic property of the transformation and of the covariance structure are determined by how well the standard bivariate normal fits the data in the standardized Z space, shown in Fig. 8b. The original joint PDF, \(p({{{\bf{x}}}}\,| \,\bar{C},\xi )\), is recovered in the X space by applying Eq. 14. The outcome, shown in Fig. 8c, is an excellent model for a joint distribution, which would typically be complicated to analytically represent and infer from the data.