Bayesian Calibration and Uncertainty Quantification for a Physics-Based Precipitation Model of Nickel–Titanium Shape-Memory Alloys

Uncertainty quantification (UQ) is increasingly receiving attention due to advances in computational infrastructure and computer simulation models in multiple applications. These models are used to help understand the behavior of complex systems through representation of the phenomena that govern these systems using mathematical and physical formulations, and subsequently predicting specific quantities of interest (QoIs). It is widely accepted that models are imperfect since they involve modeling very complex phenomena that typically occur across different length and time scales [1]. Even if the correct values of the model inputs are used in the simulation, the output of the model will differ from the true value of the physical phenomenon being simulated due to incomplete understanding of the model parameters, incomplete knowledge of the physics, and/or intrinsic stochastic behavior of the system being simulated [2]. Consequently, the challenging task of characterizing and quantifying the uncertainties associated with these models is crucial toward leveraging their benefits and streamlining their applicability. Applications of UQ span a wide range across several science, engineering, and social disciplines such as climate and weather [3], forestry [4], nuclear engineering [5], computational fluid dynamics [6], biological phenomena [7], medicine [8], and econometrics [9].

One field in which the application of UQ is highly desired is computational materials science in general [10] and the emerging field of Integrated Computational Materials Engineering (ICME) in particular [11–13]. ICME is a new approach within the broad field of materials science and engineering that aims to integrate computational materials models to enable the optimization of the materials, manufacturing processes, and component design long before components are fabricated [14]. It has evolved due to the tremendous advances in computational materials models over the past three decades. In a recent roadmapping study by The Minerals, Metals & Materials Society (TMS) [15], UQ of these computational models was identified as a fundamental issue that needs to be addressed within the ICME framework. As stated in many recent reports [10,12,13], all computational materials models have uncertainties due to the stochastic nature of materials structures in addition to the levels of uncertainties associated with simulation methods [14]. To date, the focus in ICME has been on the development of physics-based computational models, but UQ has not been a central question despite its crucial importance [10]. In this work, we conduct one among the early efforts in applying formal UQ to an important class of computational materials models known as precipitation evolution models.

Precipitation evolution models are concerned with understanding and predicting the formation, growth and evolution of secondary phase particles (precipitates) that form from an initially supersaturated matrix phase during prescribed thermal treatments. Precipitates can play an important role in enhancing and modifying material properties (e.g., strength and hardness), and hence part performance [16,17]. Accurate modeling of the evolution of precipitate populations requires considerable understanding not only of the thermodynamics (i.e., driving forces) and kinetics (i.e., diffusion, growth rates) of the materials of interest, but also of the coupling of sophisticated large-scale simulations with difficult to acquire experimental data. Consequently, uncertainties are typically associated with these interacting components which results in hampering our understanding of the underlying process and our efforts to develop meaningful predictions [1]. These uncertainties originate from many different sources, some of which are related to the choice of model itself, while others relate to model inputs and parameters [6]. For example, the model may include oversimplified assumptions about some complex and intricate physics. Another possible source of uncertainty might originate from the lack of knowledge about the true value of intrinsic physical properties of the material, such as density or diffusion coefficient, that serve as inputs to the model. Kennedy and O'Hagan [18] provide a classification of the sources of uncertainty in computer simulation models.

The literature on formally characterizing uncertainties associated with ICME computational materials models is sparse. Traditionally, research efforts have utilized Monte Carlo methods to achieve this task [19–23]. Although Monte Carlo methods represent a powerful and well-studied approach, they require a large number of simulation runs, which becomes impractical and often infeasible with ICME models which are typically computationally expensive. In contrast, we present a UQ approach based on surrogate modeling to conduct Bayesian calibration of a physics-based computational model used to predict the evolution of precipitates in advanced shape-memory alloys subject to heat treatment. The calibration problem (sometimes known as the inverse UQ problem) refers to the adjustment of a set of model parameters such that the model's predictions are in agreement with experimental data [24].

Shape-memory alloys (SMAs) are a class of advanced metal alloys attractive to the industrial and academic communities due to their ability to remember their original shape. In other words, after deformation they have the ability to retain their original shape upon being heated [17]. The application scope of SMAs is broad across various sectors such as the automotive, aerospace, manufacturing, and energy exploration sectors [25]. Recently, it has been shown that control of the distribution of secondary strengthening precipitate phases in some SMAs can enable further control strategies to tune their transformation behavior to specific engineering applications [16,26].

The paper is organized as follows: Sec. 2 presents an overview of the literature. The current work represents a unique intersection of the fields of UQ, Bayesian calibration, physics-based precipitation models, and shape-memory alloys. Next, we present the precipitation model for nickel–titanium shape-memory alloys (NiTi SMAs) and its attributes in Sec. 3. The Bayesian calibration model is described in Sec. 4, and subsequently, Sec. 5 reports the application of the calibration model into the precipitation model and the obtained results. We finalize with conclusions, insights, and directions for future work in Sec. 6.

Modeling and simulation of precipitate evolution during material processing is of great interest to material scientists and engineers since it can have considerable impact on mechanical properties and part performance. Research on precipitation modeling typically starts with the development of thermodynamic models—using, for example, the CALPHAD framework—in order to identify the thermodynamic conditions that yield specific phase constitution states [27–29]. One of the earliest numerical precipitation models upon which many other studies have relied is the model provided by Wagner et al. [30]. The existing literature on precipitation modeling focuses on microstructural features and mechanical properties. For instance, Robson et al. [31] develop a model to account for the effect of nucleation processes on dislocations, and predict precipitation behavior in aluminum [32]. A different study conducted by Deschamps et al. [33,34] presents research on microstructure evolution and modeling of precipitation kinetics by integrating nucleation growth and coarsening under different heat treatments schemes, and then maps them to mechanical properties. Other works focus on phase field models to assess the morphology of precipitates at different sizes [35] and to analytically predict the evolution of the precipitate volume fraction using calorimetry techniques [36].

While there exists a plethora of research efforts on the topic of computational materials models, there is a corresponding lack of efforts on uncertainty quantification of these models. Chernatynskiy et al. [10] present a review of the few existing works on UQ in computational materials model with emphasis on multiscale simulations. One example includes Cai and Mahadevan [37] who conducted uncertainty propagation in a manufacturing process through leveraging modeling using Gaussian processes. Another example is given in Crews and Smith [38] where the authors use UQ to predict heat transfer parameters of SMA bending actuators.

One possible explanation for the lack of research on UQ in the context of computational materials models is the fact that advances in these models are recent and have focused on aspects relevant to the underlying physics responsible for the phenomena of interest. As these models continue to develop and mature, the need for UQ receives more attention as noted in multiple reports [14,15]. An important note, however, is that despite the lack of this integration, UQ as an independent field is quite mature due to the broad scope of applications that it addresses. Numerous UQ frameworks and methods have been developed to solve different problems such as regression, subset selection, hypothesis testing, and basis expansions [39,40]. In this paper, we focus on Gaussian process (GP) models to conduct UQ in ICME computational materials models.

Gaussian processes have a long history of application in mining, agriculture, and forestry since 1920 [41]. Due to their versatility, they have experienced significant growth over the past three decades, especially with the advances in low-cost high-speed computing [42]. GP models are well suited for studying and quantifying the uncertainty of functions and for making predictions given observed data. They have a very wide range of applications that include spatial modeling [43], computer model emulation [44], image analysis [45], and supervised classification and prediction in machine learning [46]. Rasmussen and Williams [46] present the theory and applications of GPs, while Cressie [47] and Gelfand et al. [42] focus on GP applied to spatial statistics (often called geostatistics), its implementations and generalization to different problems. Kennedy and O'Hagan [18,48] set foot for GP in the UQ field with their seminal article on calibration of computer models, a study that served as base for many subsequent works, including the current work, such as Refs. [49–51] for surrogate modeling and [8,52,53] for calibration.

The current work is first to conduct systematic UQ of a precipitation evolution computational model in shape-memory alloys. By conducting this task, we improve the predictive capability of the model which is very useful in the design and processing of materials. We focus on NiTi SMAs as a model material due to their unique features highlighted in Sec. 1. Two excellent review articles on SMAs are presented by Lagoudas [54] and Ma et al. [25] summarizing modeling, engineering, and high temperature applications. The nickel–titanium (NiTi) alloy considered in this paper (commercially known as nitinol) is a popular class of SMAs that has been broadly investigated for medical applications due to its biocompatibility [55,56], superelasticity, and strain control [57,58].

The focus of the current work is on conducting uncertainty quantification of a physics-based model, rather than the development of the model itself. In other words, the physics-based model is treated as a black box and the focus is on the development of the statistical UQ framework. In this section, we present a high-level overview of the physics-based model of interest but omit technical details since they are out of the paper's scope.

Material modeling tools are generally used to determine the evolution of a material over time based on the environmental conditions imposed upon it. One approach, precipitation modeling, is used to predict the nucleation and growth of a secondary phase within the matrix of a primary phase. These small precipitates (secondary-phase particles dispersed throughout the material) are compositionally and/or structurally different than the surrounding matrix phase, allowing their process-dependent morphologies to be tracked as they evolve in the presence of various environmental conditions. The sizes and shapes of these particles can drastically alter the properties of the overall material [59,60]. In the context of this research, we are interested in the nickel composition of the matrix which is calculated at any step in the process through stoichiometric considerations and a mass-balance equation involving nickel in the precipitates and the initial nickel content of the material before thermal processing. In turn, the composition of nickel in the matrix, the volume fraction of the precipitates, and their size distribution control the final shape-memory behavior in NiTi alloys [16].

MatCalc, the microstructural modeling tool used in this research, utilizes the numerical Kampmann Wagner (NKW) approach, coupled with proprietary thermodynamic and kinetic databases, to model the phase transformation of a material [61–63]. The model is unique in that it calculates the evolution of a discrete size distribution of particles in a material system when given a set of times and temperatures representing a thermal processing schedule. This model strikes a balance between the fast but overly simplified Johnson–Mel–Avrami–Kolmogorov (JMAK) approach [64] and the accurate but computationally intensive phase field approach [65]. The size distribution provides a more complete physical representation of the system than JMAK and the mean-field approximation for each size class leads to much faster calculations overall when compared to more detailed microstructural evolution techniques, such as phase field modeling which explicitly simulates the evolution of the entire microstructure, including secondary phases.

The evolution of the precipitate size distribution is tackled by applying the nucleation, growth, and coarsening equations to each bin at each time step, then updating the size distribution bin sizes to reflect the evolution that occurred during said time step. The nucleation of each bin in the distribution is calculated by governing equations based on classical nucleation theory. These nucleated bins are then individually subjected to growth-based governing equations rooted in the thermodynamic extremal principle. Finally, the coarsening process is approximated by establishing a minimum precipitate radius through an analysis of the system as a whole. Precipitates smaller than the minimum size value are then dissolved and redistributed to the next nearest bins larger than the minimum radius. These three evolution calculations are repeated at every time step, leading to the overall evolution of the size distribution of particles. Each of these sections has a different set of governing equations which are not included in this paper for the sake of brevity, but can easily be found at the reference provided [63]. Although each set of governing equations has been derived to explain a different portion of the material's overall evolution, they predominantly rely on the same parameters and variables meaning that all three sets of equations can be calibrated holistically.

Three terms (interfacial energy, diffusion correction prefactor, and diffusion correction exponential factor) significantly influence the kinetic behavior of the overall model and as such can be used to calibrate the model to experimental data. The interfacial energy term influences the behavior of the nucleation portion of the model and the diffusion correction factors equally influence all three stages of precipitate evolution listed above by modifying the substitutional diffusivity of the matrix phase. Values for interfacial energy are known with much less accuracy than they should be, especially when considering how sensitive the entire model is to this term. In fact, MatCalc attempts to calculate an interfacial energy value using the generalized broken bond (GBB) method if no calibrated interfacial energy is given [63,66]. This GBB approach works well for many systems but failed to accurately reproduce experimental results for $Ni4Ti3$ precipitation. Previous parametric calibration attempts lead to an interfacial energy value of $0.0505 J/m2$⁠. The diffusion correction terms can be seen as unknown values because they are essentially fitting parameters for the diffusivity of the matrix. In the remainder of the paper, we focus on conducting statistical calibration of the three parameters described above.

We start by introducing some definitions and notations. We distinguish between the model presented in this section to conduct statistical calibration and the physics-based precipitation model (Sec. 3) that is to be calibrated. We will refer to the former as the statistical model and the latter as the computer model to avoid confusion. It is common practice in the literature to refer to computational models as computer models since they involve running a computer simulation code.

with mean vector $μ$ and covariance matrix C defined by a mean function $μ(·)$ and a covariance function $C(·,·)$ respectively. We will denote this as $f(·)∼GP(μ,C)$⁠, with $μ(x)=E[f(x)]$ and $C(xi,xj)=cov(f(xi),f(xj))$⁠. It is evident that the choice of mean and covariance functions will influence the distribution in Eq. (1); therefore, they need to be carefully selected. Useful guidelines for their selection are given in Refs. [46] and [67] for the interested reader.

In the calibration problem, we distinguish between two groups of inputs to the computer model:

• Calibration parameters denoted by the vector $θ$ represent physical parameters that can be specified as inputs to the computer model, and whose values are unknown or not measurable [18,49]. Examples include properties of materials (such as conductivity or interfacial energy, among many) which are not easily determined for materials models, or saturation constants in ecosystem models [68]. The goal of the calibration problem is to estimate these parameters such that the computer model simulations agree with the experimental observation of the real process being simulated.

• Control inputs (or design variables) denoted by the vector x are variables that are set to known values by the user [18]. Examples of these inputs include temperature, pressure, force, or any other quantities that are known and controlled by the user.

where ρ is a scaling factor of the computer model, $δ(·)$ is a discrepancy function or model inadequacy function whose role is to account for the missing physics in the computer model, and $ε(·)$ is the observation or measurement error.

The calibration problem involves a synthesis of computer simulations and experimental observations. More specifically, to solve the calibration problem, one needs to collect experimental observations, $Z(x)$⁠, and run a sufficiently large number of computer model simulations $Y(x,θ)$⁠. The main task is to use these to infer the values of the calibration parameters $θ$ such that the computer model output most closely matches the real process represented by experimental observations. An important challenge within this context is the computational burden associated with generating the large number of computer model simulations needed to effectively conduct calibration. For example, some computer models might require the numerical solution to systems of coupled nonlinear partial differential equations. In this case, even generating a few simulations can be a challenging task.

When the above holds true, it is common practice to construct a surrogate model (sometimes known as an emulator or a meta-model). This represents a computationally efficient approximation that encapsulates the behavior of the original computer model and can thus be used to generate a large number of simulations. Of course, when the computer model is not computationally expensive, it can then be directly used to conduct calibration, without the need to construct this computationally efficient surrogate description. We will refer to this as direct calibration. It is important to note that although our main focus is on calibrating the precipitation model using a surrogate model, we will also present the direct calibration approach for benchmarking and comparison purposes.

We model the discrepancy function defined in Eq. (2) as a Gaussian process, $δ(·)∼GP(μ,C)$⁠. The mean function $μ(x)=H(x)β$⁠, where $H(x)=[1,x⊤]$ is a row vector of basis functions and $β=[β0,β1,…,βqx]⊤$ is a vector of regressors. In other words, $μ(x)=β0+β1x1+⋯+βqxxqx$⁠. The covariance function C will depend on an additional set of parameters $ψ$ (commonly known as hyperparameters). These will be presented in detail in Sec. 5. Finally, the measurement error terms are modeled as independent and identically distributed (iid) normal random variables, $ε(x)∼N(0,τ2)$ for all $x∈X$⁠. Hence, the statistical model is fully characterized by the set of parameters $ϕ={ρ,θ,β,ψ,τ2}$⁠.

where $p(Z|ϕ)$ is the likelihood function defined as the probability of observing Z conditioned on the model parameters (⁠$ϕ$ in this case) [69]. The term $p(ϕ)$ is the prior distribution of $ϕ$ that captures our prior belief or knowledge, and $p(ϕ|Z)$ is the posterior distribution of $ϕ$ conditioned on experimental observations.

Selection of these prior distributions plays an important role in the parameter estimation process. This selection is driven by our prior knowledge or belief regarding the parameters. It is common practice to assign uninformative prior distributions to parameters for which little knowledge is available (for example, a uniform distribution or normal distribution with large variances). Finally, utilization of conjugate priors such that the prior and posterior distributions of the parameter belong to the same family is usually recommended and offers practical and computational advantages [41].

Several methods can be utilized to obtain this posterior distribution. We rely on Markov chain Monte Carlo (MCMC) methods, specifically Gibbs sampler and Metropolis–Hastings (MH) algorithm [70,71].

The BLUP predictor will be subsequently used in Sec. 5 to assess the prediction accuracy of the calibrated model using a test set or cross-validation procedures.

for $x∈X$ and $θ∈T$⁠.

where the star notation for the control inputs x in the set D₂ indicates that they might be of different values than the control inputs used in D₁.

where $H1(x,θ)=[1,x⊤]$ and $H2(x)=[1,x⊤]$⁠. Similar to the direct calibration case, the covariance functions $C1((·,·),(·,·))$ and $C2(·,·)$ are fully characterized by the sets of hyperparameters $ψ1,ψ2$⁠, respectively, and henceforth, all of the parameters from this statistical model can be grouped into two sets: $β=[β1⊤,β2⊤]⊤$ and $ϕ={ρ,ψ1,ψ2,τ2}$⁠, which along with calibration parameters $θ$ represent all the statistical parameters that will be estimated using the Bayesian procedure. Notice that the measurement error is still following a normal distribution $N(0,τ2)$⁠.

thus, analogous to Eqs. (9) and (10), respectively.

Bayesian estimation of the parameters is conducted in two stages. The first stage is to estimate the hyperparameters $ψ1$ by maximizing the posterior probability distribution $p(ψ1|Y)$⁠. Subsequently, in the second stage, the remainder of the parameters ${ρ,ψ2,τ2}$ is estimated by maximizing the conditional posterior distribution $p(ρ,ψ2,τ2|d,ψ1)$⁠. A detailed procedure to carry out these optimizations is presented in the seminal work by Kennedy and O'Hagan [48].

Following these two stages, the remaining step is to obtain the posterior distribution of the calibration parameters $θ$⁠. We employ MCMC methods in a similar fashion to the direct calibration case.

In this section, we implement both calibration approaches presented in Sec. 4 to calibrate the computer model of interest. As detailed in Sec. 3, the computer model in our case predicts the evolution of precipitates in NiTi shape-memory alloys subject to heat treatment. More specifically, the model predicts the atomic percent (at. %) of nickel in the matrix following heat treatment for a specified period of time at a preset temperature.

The control inputs and calibration parameters are given as follows:

• Three control inputs x:

  1. —

     Initial nickel content in the alloy (at. %)

  2. —

     Aging (heat treatment) temperature (°C)

  3. —

     Aging (heat treatment) time (s)

• Three calibration parameters $θ$⁠:

  1. —

     Interfacial energy (⁠$J/m2$⁠)

  2. —

     Diffusion correction prefactor (dimensionless)

  3. —

     Diffusion correction activation (⁠$J/mol$⁠)

The output of the process $Z(x)$ and of the computer model $Y(x,θ)$ correspond to the final nickel content of the surrounding matrix (at. %). The posterior distribution of the unknown calibration parameters will be derived using a synthesis of experimental observations and computer model simulations.

The execution time of the computer model is highly dependent on its inputs. For example, although simulation under a heat treatment setting with longer aging time is generally expected to take longer than with shorter aging time, some certain values of the calibration parameters can result in a different outcome. From running multiple simulations in the initial model testing phase, execution times ranged between 3 s and 5 min. It is important to point out that although this might not represent significant computational burden from the perspective of obtaining model prediction for a particular input combination, in conducting statistical calibration a large number of computer simulations is necessary to perform MCMC simulation (10,000 or more). For example, the median execution time is 2.5 min, thus generating a sufficient number of simulations to conduct MCMC calibration can take up to 2 weeks. This makes the task of conducting subsequent sensitivity analysis, uncertainty analysis, and what-if scenarios that are typically needed by the user quite challenging. Our focus in this work is thus on conducting surrogate model calibration. We will implement the direct MCMC calibration for benchmarking purposes.

The experimental observations consisted of 36 data points. Each data point represents the final nickel content of the matrix (at. %) corresponding to a specific heat treatment regime (aging temperature and aging time) for a particular NiTi alloy with initial composition. It is obvious that these data are consistent with the computer model output and control inputs, respectively.

To avoid oxygen contamination, the samples were individually sealed into quartz tubes under a high-purity argon atmosphere. Heat treatments were performed in a muffle furnace followed by subsequent water quenching. A technique called differential scanning calorimetry (DSC) was first used to directly measure the resulting martensitic transformation temperature of each heat treated alloy. The martensitic transformation temperature is a property of NiTi alloys that is highly sensitive to the amount of nickel in the surrounding matrix. The formation of nickel-rich precipitates causes a change in the composition of the alloy, which in turn reflects into a change in the martensitic transformation temperatures. In other words, the martensitic transformation temperature determined through DSC can be used to indirectly measure the matrix composition.

The relationship between NiTi composition and the martensitic transformation temperatures has been extensively studied in literature, such as Refs. [73–75]. In the current work, the DSC measurements of the martensitic transformation temperatures were converted to nickel composition (at. %) using curve fitting to data presented in Frenzel et al. [76]. This yielded the 36 experimental observations data set used henceforth. From this set, 31 will be used to calibrate (or train) the precipitation model as detailed next and five will be used for test and validation to assess the predictive performance of the calibrated model. Figure 1 and Table 3 in Appendix presents the whole experimental data set.

Following the notation of Sec. 4, we have that q_x = 3 and control inputs $x=[x1,x2,x3]⊤$⁠, $qθ=3$ and calibration parameters $θ=[θ1,θ2,θ3]⊤$⁠. We define region $X$ to include all values from the experimental set; that is the hypercube with bounds $Xmin={50.5,350,35000}$ and $Xmax={52,600,380000}$⁠. Similarly, region $T$ is also defined as a hypercube; however, its bounds are specified by the modelers based on their experience and prior knowledge of the process to a region within which they have reason to believe that the true values of $θ$ might lie. These bounds were defined to be $Tmin={0,0,0}$ (not included) and $Tmax={0.1,100,1000}$⁠.

is a distance measure between $xi$ and $xj$ known as the Mahalanobis distance [42], with length scale parameters ω_k, which control the relative influence of each input dimension to the process. It can be seen that the covariance function is fully defined by hyperparameters $ψ={σ2,ω1,ω2,ω3}$⁠. The smoothness parameter was kept constant at $ν=1/2$⁠, based on preliminary work by the authors and common practice in the literature [41,42].

The reasoning behind this selection is straightforward. With the exception of the variance parameters $σ2$ and $τ2$⁠, uninformative uniform distributions were selected for all the parameters due to the lack of prior information regarding the shape of the distributions. The scaling factor ρ is defined to take values between 0 and 1, while β_j can take on any real value and ω_j is limited to positive values due to the intrinsics of Mahalanobis distance. The calibration parameters $θ1,θ2$⁠, and θ₃ take values within the previously defined bounds of the study region $T$⁠. Notice that an indicator function $1{A}(x)$ has been used for the prior distributions of $θ1,θ2$⁠, and θ₃. This function is equal to 1 if $x∈A$ and 0 otherwise. It was used to reject values of $θ1,θ2$⁠, and θ₃ for which the final (output) and initial (input) content of nickel are equal, which indicates that precipitation did not occur. This represents an anomaly that needs to be discarded. Finally, inverse gamma prior distributions were selected for the variance parameters $σ2,τ2$⁠. These represent conjugate distributions with respect to the normal likelihood distribution.

Next, we conduct Bayesian calibration of the model by obtaining the posterior distributions of the parameters above using MCMC simulations. We employ the Metropolis–Hastings algorithm, which requires the evaluation of the posterior distribution given in Eq. (6) at each iteration. This implies that Y needs to be evaluated, and therefore the computer model $Y(x,θ)$ is executed $n×qθ=31×3=93$ times at each of the MCMC iterations. We utilize parallel programming in the 862-node Intel x86-64 Linux cluster supercomputer at Texas A&M high performance research computing facilities. MCMC simulations were implemented over ten cores in the supercomputer and run for 15,000 iterations, with a 25% burn-in period and thinned every fifth iteration. Results of the posterior distributions of the calibration parameters are shown in Fig. 2, where each plot displays the histogram with 20 bins and their corresponding kernel density estimate.

In the graphs, we can see that θ₁ and θ₂ have defined informative and unimodal posterior distributions, with θ₁ having a unique peak and θ₂ being skewed to the right. The posterior distribution of θ₃ shows some uniformity, which implies that all the values across its domain are equally likely to characterize the real process. Table 1 presents mean, mode (most frequent value), and standard deviation of the posterior distribution for the calibration parameters.

On its own, the distribution of θ₁ can be attributed to the fact that it appears as a cubic term within an exponential function of the governing equations for precipitate evolution. The narrow nature of its distribution, in comparison to the other two parameters, can be explained by its appearance in the governing equation for precipitate nucleation, which is the materials science equivalent of initial conditions drastically influencing the solution to a differential equation. The interfacial energy is so dominant in precipitation models due to the fact that relatively small values imply almost no barrier to nucleation, leading to the precipitation of all the precipitates that can form in the material (ultimately limited by thermodynamics) right at the beginning of the precipitation process. On the other hand, slightly high values for this property may result in a barrier to nucleation that effectively results in no precipitation of secondary phases at all during reasonable heat treatment times. Thus, only a narrow range of interfacial energy values can yield physically meaningful precipitation predictions. In fact, the estimated mean and mode of θ₁ are well within the range of values reported in the literature and very close to the value calculated through parametric sweep-experimental data coupling.

It is more difficult to ascertain direct significance from specific values of θ₂ and θ₃ distributions since they represent corrections to the diffusivity values taken from the MatCalc NiTi kinetic database. These corrections are necessary as one does not really know the state of the microstructure before the precipitation starts. Of particular effect is the presence, for example, of grain boundaries, dislocations, and other defects that can catalyze (facilitate) the nucleation of the second phase precipitate particles. The ranges for each parameter were selected such that only physically realistic values would be calculated for these two parameters. Specifically, θ₂ and θ₃ correspond to corrections to the diffusion coefficient and activation energy values, respectively. Although they represent correction factors and as such have no direct physical significance themselves, it can be inferred that these two diffusion correction values are less influential (i.e., have broader distributions) than θ₁ for two main reasons. First, the equations that contain them are not as drastic as a cubic in an exponential, and second, it is because they are implemented in a later stage of the precipitate evolution which reduces their influence with respect to θ₁.

We now calibrate the computer model using the surrogate calibration to account for the computational burden associated with direct calibration. Calibration will be conducted based on the redefined statistical model in Eq. (10).

The first step toward surrogate calibration is to construct the computationally efficient surrogate model that can be used in lieu of the original computer model. We start by generating a data set of computer model simulations over different combinations of $(x,θ)$⁠. A data set D₁ with size N = 3025 data points was generated. Latin Hypercube Sampling was used to select values of $(x,θ)$ uniformly over the $X×T$ space at which the computer model is evaluated. This data set was again implemented using parallelization in the supercomputer to reduce the execution time. The same data set of experimental observations used for direct calibration was utilized for surrogate model calibration. We denote this as D₂ with size n = 31.

are now squared Mahalanobis distances. Similar to the previous case, the covariance function C₁ is fully defined by hyperparameters $ψ1={σ12,ω1,…,ω6}$⁠, and C₂ by $ψ2={σ22,ω7,ω8,ω9}$⁠. Notice that these covariance functions represent special cases of the Matérn function when $ν→∞$ [42].

Next, we estimate the hyperparameters of the GP models by maximizing their marginal posterior distributions in stages 1 and 2 as explained in Sec. 4.2. Subsequently, estimation of the calibration parameters was performed using MCMC with the same specified simulation procedure; that is, 15,000 iterations with 25% burn-in period and thinning set to every fifth iteration. The posterior distributions for the calibration parameters are shown in Fig. 3.

By inspecting the obtained posterior distributions, we see some differences from those obtained using direct calibration. The distribution of θ₁ is an informative bimodal behavior with two distinctive peaks. The distribution of θ₂ is skewed to the left, in contrast to the direct calibration case, and the distribution of θ₃ shows similar uninformative behavior to the direct calibration case. Posterior values are presented in Table 2 in a similar manner as with previous case.

As with the previous case, θ₁ is the parameter with direct physical significance while θ₂ and θ₃ represent correction factors for the kinetic equations of the model. The bimodal nature of θ₁ is most likely a result of the inability of the surrogate model to directly compare initial and final nickel content of the matrix, which is used as a sanity check in the direct model. To elaborate on this, note that while the surrogate model is to replace the computer model, it does not capture any physics, and thus does not recognize that cases where the initial and final nickel content are identical should be discarded.

To make this more explicitly clear: in the direct calibration model, an unreasonably small value for the interfacial energy would result in instantaneous saturation of the precipitate phase to volume fractions predicted through thermodynamics (lever rule). On the other hand, unreasonably large values for this quantity would result in a practically infinite barrier for nucleation, leading to no precipitate phases regardless of the heat treatment time. The former case would result in complete nickel depletion at times close to t = 0, while the latter case would correspond to no nickel depletion at all. In the direct model, these two problematic cases could be detected and discarded from the sampled distribution. That being said, both peaks are within the acceptable range of values seen in the literature.

The explanation of θ₂ and θ₃ distributions remains essentially the same for the surrogate model as it was for the direct case. These two parameters are correction factors and as such have no direct physical connection. All that can be said is that the mean and mode values seen in Fig. 3 will result in physically realistic diffusion coefficient and activation energy values. Again, these two distributions are broader because they influence later stages of the precipitate evolution which have less of an impact on overall precipitate morphology than early stage parameters like θ₁.

In terms of hyperparameters, the estimated variances are as follows: $σ12=0.1624, σ22=0.3324, τ2=0.0092$⁠. These results are consistent with what we would expect. For instance, the marginal variance of the surrogate model $σ12$ is smaller than that of the discrepancy function $σ22$⁠. This is justified since the discrepancy function accounts for the missing physics between the computer model and the experimental observations; therefore, it handles more noisy data, which translates to larger uncertainty for the discrepancy function. Additionally, the surrogate model replaces a deterministic computer model; therefore, it is expected that the surrogate will have smaller uncertainty given that the simulation dataset D₁ does not involve any noise besides numerical errors.

Regarding the length scale parameters, a large value translates to low influence of that specific input dimension on the covariance function and vice versa. The estimated values are $ω1:9={1.6759,0.2967,12.9291,0.1317,4.7732,12.6946,29.9991,29.9424,5.9599}$ which show several interesting results. Within the surrogate model (Eq. (17)), the most important input dimension corresponds to the interfacial energy (related to ω₄), which makes sense as it highly influences the precipitation process of nickel. On the opposite side, Aging time (ω₃) has been estimated to be less impactful on the surrogate model. In Sec. 5.5, we present a variance-based sensitivity analysis of the model, which further explains the dependences of the output of the process on each of the inputs.

In Secs. 5.2 and 5.3, 31 experimental observations were used to calibrate the precipitation model using both the direct and surrogate modeling approaches. We now use five additional experimental observations from the test set to assess their performance.

where $npost$ is the total number of posterior distribution samples calculated during MCMC (after burn-in period and thinning) and $ϕ(i)$ is the ith posterior sample. A similar procedure is utilized for the predictive distribution parameters in the remaining Eqs. (8), (14), and (15).

First, we run a formal ten-fold cross-validation procedure solely focused on assessing the performance of the surrogate model as an approximation for the computer model and how accurate of an approximation it represents. Figure 4 shows the results of this cross-validation, where the abscissa in the plot is the simulated values (the computer model output), while the ordinate is the mean of the predictions obtained by the surrogate model; therefore, a point in the plot compares a simulation $Y(x,θ)$ with a prediction $E[Y(x)|Y,θ]$ for a specific input $(x,θ)$⁠. In this case, the straight line is a reference indicating ideal predictions in the case when simulations and predictions are identical.

The surrogate model appears to adequately capture the behavior of the computer model. This is demonstrated by the points being reasonably scattered close to the ideal red line. It is worth noting that we are comparing 3025 points from dataset D₁ and that the most dense regions in the plot are those close the red line. Numerically, the average deviation from the ideal line is 0.0661 (at. %), which further confirms the results.

Next, we use the predictive distributions from Eq. (20) to assess the performance for both approaches, where we employ the set of data points not used for training the models and compare their predictive distribution to the actual observed values (which might be noisy). Figure 5 displays model predictions for the test set consisting of five points using the direct and surrogate calibration approaches, respectively. The abscissa in the plot is the experimental observation while the ordinate shows the mean of the predicted value, $E[Z(x0)|Z]$⁠, at the squared marker. The confidence interval bars represent ±1 standard deviation, $σ0=var[Z(x0)|Z]$⁠.

where $Zobs,i$ is the ith experimental observation from the test set, $Zpred,i$ is the ith predicted value (mean of the predictive distribution in this case), and $npred$ is the number of data points in the test set (five in our case). The standard error of prediction was equal to 0.0656 and 0.0675 for direct and surrogate calibrations, respectively. These values are approximately equal to 3% of the full range of experimental observations, which is acceptable given all the experimental uncertainties in these types of studies.

where $xi(1),xi(2)∼iidU(Xmin,Xmax)$⁠. We denote a hybrid point $xu(1):x−u(2)$ where the jth element equals $xj(1)$ if $j∈u$⁠, and $xj(2)$ if $j∉u$⁠.

where a large value of $κ¯̂u2$ represents a higher influence of variables subset u on function f. For reporting purposes, the indices are frequently normalized by the variance of the function evaluated at dataset 1, $var(f(X(1)))$⁠.

We apply this procedure into our application and first calculate the indices for the computer model itself, in order to analyze which of the six inputs (three control inputs and three calibration parameters) impact it the most. First pane of Fig. 6 shows the calculated total sensitivity indices for the computer model, where the subsets u were selected to be single inputs (i.e., leftmost bar represents input x₁, hence its subset is $u={1}$⁠). As we can see, initial nickel content (x₁) dominates the variance of the computer model followed by interfacial energy (θ₁), which is consistent with expert knowledge and results from the hyperparameters estimation in Sec. 4.2.

Similarly, we calculated Sobol indices for the calibrated surrogate model. In this case, we use the posterior distributions of the calibration parameters reported in Fig. 3 to calculate the total sensitivity indices for each of the inputs as well. The results can be visualized in the second pane of Fig. 6 with consistent results once again, where input x₁ is computed as the most controlling input to the surrogate model variance.

Uncertainty quantification (UQ) in materials simulation models has not been a central question in the UQ literature despite their great importance. In this paper, we have presented an inverse UQ problem (commonly known as the calibration problem) to a materials simulation model that predicts the evolution of precipitates in nickel–titanium shape memory alloys (NiTi SMAs) subject to aging heat treatment. We used a Gaussian process (GP) based surrogate modeling approach to conduct Bayesian calibration. The surrogate modeling approach was used in order to account for the computational burden associated with the precipitation model. The predictive performance of this surrogate modeling approach was benchmarked against the classical Markov chain Monte Carlo (MCMC) approach using a case study that involves real-world experimental observations.

The Bayesian approach used is invaluable both for solving the calibration problem and for quantifying the uncertainties in model predictions. Results show that the calibrated model using both approaches has good prediction accuracy, making it suitable for ICME-based design of shape memory alloys.

There are multiple follow-up studies currently undertaken by the co-authors. Namely, employing adaptive sampling techniques to accelerate and improve the predictive performance of the calibrated model, and coupling the calibrated precipitation model with other physics-based microstructure evolution models to quantify the uncertainty throughout the entire material processing chain, with particular emphasis on laser-based additive manufacturing of shape memory alloys.

This work was supported by an Early Stage Innovations grant from NASA's Space Technology Research Grants Program, Grant No. NNX15AD71G. Portions of this research were conducted with High Performance Research Computing resources provided by Texas A&M University.² L. J. and R. A. also acknowledge the support of NSF through the NSF Research Traineeship (NRT) program under Grant No. NSF-DGE-1545403, “NRT-DESE: Data-Enabled Discovery and Design of Energy Materials (D³EM).” RA and IK also acknowledge the partial support from NSF through Grant No. NSF-CMMI-1534534.

 
• BLUP =

  best linear unbiased predictor

 
• GP =

  Gaussian process

 
• ICME =

  Integrated Computational Materials Engineering

 
• MCMC =

  Markov chain Monte Carlo

 
• QoI =

  quantities of interest

 
• SMA =

  shape-memory alloy

 
• UQ =

  uncertainty quantification