A Bayesian Framework for the Estimation of the Single Crystal Elastic Parameters from Spherical Indentation Stress-Strain Measurements

This paper presents a two-step Bayesian framework for the estimation of the intrinsic single crystal elastic stiffness parameters from the measurements of spherical indentation stress-strain responses in multiple individual grains of a polycrystalline sample, whose crystal lattice orientations have been measured using electron back-scattered diffraction technique. The first step requires the establishment of the functional dependence of the indentation elastic modulus given the lattice orientation and the intrinsic single crystal elastic stiffness parameters. Previous efforts for this step required a large database of computationally expensive finite element (FE) simulations in order to establish this function with adequate accuracy. In this paper, it is shown that the introduction of a Bayesian framework can greatly reduce the number of simulations necessary to establish this function, while introducing practically useful measures of uncertainty which can guide the selection of specific additional simulations that are expected to best improve the predictive accuracy of the function. The second step involves a Markov Chain Monte-Carlo (MCMC) sampling of the distribution of possible values for the single crystal elastic stiffness parameters based on a given set of experimentally measured elastic indentation moduli in individual grains of different lattice orientations. This second step is accomplished by calibrating the available experimental data to the function established in the first step. This novel framework is presented and demonstrated in this paper for an as-cast cubic polycrystalline Fe-3% Si sample and a hexagonal polycrystalline commercially pure (CP-Ti) titanium sample.


Introduction
Continued development and application of physics-based multiscale materials models is largely hampered by the lack of protocols for reliably estimating the intrinsic material properties at the microscale (e.g., the grain-scale properties in modeling of polycrystalline materials). In recent years, instrumented indentation techniques have been demonstrated to be capable of providing consistent and reliable measurements at the lower length scales (up to submicron length scales) [1][2][3][4][5]. Although small-scale mechanical measurements are now quite reliable, it has not been a straightforward process to extract the intrinsic material properties from such measurements. As specific examples, one would hope to estimate the values of the single crystal elastic constants and the critical resolved shear strengths from the instrumented nanoindentation measurements.
Currently employed strategies for extracting intrinsic material properties from indentation tests have generally involved the calibration of physics-based finite element (FE) models of these tests to the corresponding set of experimental measurements [6,[9][10][11]. In this regard, it has been pointed out in recent work [7] that these protocols are much more robust when the calibration is attempted in the form of the normalized indentation stress-strain curves as opposed to directly matching the load-displacement curves. This is mainly because the initial elastic response and the elastic-plastic transition occur over a very short early portion of the load-displacement curve that is not easily identified and isolated, resulting in a very high sensitivity of the extracted values of the intrinsic material properties to small changes in the calibration procedures.
The calibration of the FE simulated indentation stress-strain curves to the experimentally measured indentation stress-strain curves for any selected material system essentially involves solving an inverse problem. In other words, the guessed values of the intrinsic material properties of interest become inputs to the FE simulations. Typically, one has to search over a large multidimensional space to find the best-fits between the FE predictions and the measurements. The main challenge comes from the high computational expense of FE simulations of the indentation experiments. It should be noted that establishing each data point on the FE predicted indentation stress-strain curve needs the simulation of a suitable unloading segment [7], and this drives up the cost of the simulation significantly. Given all of the complexity described, the only logical path forward is to establish a reduced-order model for the FE simulations of the indentation test, and to use the reduced-order model in solving the inverse problem described above. In recent work [9], we have formalized this approach as a two-step process: (1) establishing a reduced-order model calibrated to FE simulations of indentations that takes the relevant intrinsic material properties as inputs and predicts indentation properties (defined suitably on an indentation stress-strain curve), and (2) the extraction of the intrinsic material properties from the available measurements (typically performed on grains of different orientations in a polycrystalline sample) through calibration with the reduced-order model established in step (1). The second step described above typically involves the solution to an optimization problem (i.e., minimizing the difference between the measurements and the predictions from the reduced-order model). The viability of this two-step protocol for extracting the values of the single crystal elastic constants and the critical resolved shear strengths in Fe-3%-Si has been demonstrated in recent work [7,9].
The main difficulty with the two-step protocol described above lies in building the reduced-order model (i.e., step (1)). Because of the need to cover a large space (for example for extracting single crystal elastic constants, the input space of interest is the product space spanning all combinations of the single crystal elastic constants, 11 , 12 , 44 , and all possible grain orientations), one needs to generate a large amount of the FE simulation data in order to establish a high-fidelity reducedorder model. The difficulty of this task is amplified significantly in dealing with hcp crystals, where the numbers of the intrinsic properties is significantly larger (for example, modeling the elastic deformation in hcp crystals requires specification of five independent single crystal elastic constants). In prior work [9], the reduced-order models were built using standard regression approaches. Although these regression approaches produced excellent results, they do not scale well to problems with larger numbers of the intrinsic properties (because of the need to generate a large amount of data spanning the entire input domain).
The primary goal of this paper is to demonstrate the utility of Bayesian strategies for (i) optimizing the reduced-order model building effort involved in step (1), and (ii) providing estimates of the desired intrinsic material parameters (single elastic constants specifically) with uncertainty measures from available experimental data (spherical indentation measurements). Towards these goals, we will develop and present a Bayesian inference framework for both steps of the two-step protocol described above. Bayesian inference has been instrumental in model-building tasks with limited amount of data [12][13][14][15]. The adoption of a Bayesian inference framework for the extraction of the intrinsic material properties from indentation measurements offers the following main advantages: (i) it is expected to dramatically reduce the number of FE simulations needed to produce the reduced-order model generated in step (1), and (ii) it provides a much more rigorous quantification of the uncertainty in the estimates of the intrinsic material properties obtained in step (2), while accounting for the uncertainty in the measurements as well as other sources. In this paper, we first develop the framework, and subsequently demonstrate its application to the extraction of single crystal elastic properties in selected cubic and hexagonal metals.

New Bayesian Inference Framework for the Estimation of Intrinsic Material Properties from Indentation Measurements
Let denote the set of intrinsic material properties to be established. For cubic crystals, this represents the set of three elastic constants, i.e., = { 11 , 12 , 44 }. Let denote an available set of observations of the indentation properties corresponding to the set of crystal orientations . This set of observations could come from either FE simulations or the physical experiments. We shall note the source of the data using subscripts sim and exp on these variables. Furthermore, in the notation employed in this paper, a set of values for a variable is denoted by an upper-case symbol, while an individual element of the set is denoted by its non-bold counterpart. As an example, a single value of the indentation property will be denoted by . Furthermore, a collection of variables is also denoted by bold symbols. As an example, a single crystal orientation would be denoted by , as it denotes a set of three Bunge-Euler angles [16]. However, a collection of grain orientations would be represented by . Likewise, a single set of intrinsic material parameters is denoted by , whereas a collection of the sets of intrinsic material parameters is denoted . Employing this notation, the central tasks in the two-step protocol developed in this work are the following: the sample. It is noted that the indentation properties are measured by the spherical indentation protocols mentioned earlier, while the orientations are measured using electron back-scattered diffraction (EBSD) techniques [17].
Prior experimental work [37] in single-phase polycrystalline metals has focused on exploring the dependence of indentation modulus on the lattice orientation of the indented grains (i.e., individual crystals). These findings were verified by suitable FE simulations [9]. Recently, a reduced-order model which captures the dependence of indentation modulus on both orientation and an arbitrary set of intrinsic material parameters has been established from FE simulations. The mathematical form of the reduced-order model for the present application is adopted from this prior work [9] as where K ( ) denote the symmetrized Surface Spherical Harmonics basis over the relevant orientation space of interest, and P ( ) denote a multivariate Legendre polynomial product basis.
In other words, one can express P (̅) = P 1 ( ̅ 1 )P 2 ( ̅ 2 ) … P ( ̅ ), where = ( 1 , 2 … ) forms a multi-index array, each element of which is a nonnegative integer allowed to vary from 0 to the selected maximum degree, , i.e., where ( ) enumerates the spherical harmonics that implicitly reflect the crystal symmetries of interest [16,18]. The integers and denote the truncation levels adopted in the use of Eq. (1). It is emphasized here that the model form used in Eq. (1) denotes a Fourier representation using an orthonormal basis that has been previously shown to produce compact representations for mechanical responses of crystalline solids [7,9,[19][20][21][22]. One of the central features of a Fourier representation is that the Fourier coefficients are completely independent of each other. The goal of the reduced-order modeling task here is to estimate the values of , expressed in a vector notation as , from the sparse amount of available data, as it is being generated from the expensive FE simulations. Even more importantly, our goal is to drive the model building in an optimal way by identifying the specific set of inputs for the next FE simulation such that it maximizes the improvement to the reduced-order model being built.

Building the Reduced-Order Model
The reduced-order model (see Eq. (1)) needs to be built such that it makes good predictions for the indentation modulus over a large domain of input parameters ( , ). Given the large domain of the input parameters (e.g., covering the range of values for the three independent parameters defining cubic elasticity and the two independent parameters defining the indentation direction in the crystal reference frame) and the high cost of executing a FE simulation for generating each data point, it is highly desirable to explore Bayesian regression approaches for estimating the unknown Fourier coefficients in Eq. (1). Let the corresponding sets of intrinsic parameters, , used as inputs to simulations be denoted as . The data generated from FE simulations will be denoted { , , } following the notation introduced earlier.
Bayesian approaches treat model parameters (e.g., Fourier coefficients in Eq. (1)) as stochastic variables exhibiting a distribution of values. Most importantly, Bayes' theorem allows one to update the distributions for the model parameters given new data (i.e., observations) and is commonly expressed as It is expedient to treat the distributions associated with all the stochastic variables in Eq. (3) as normal (i.e., Gaussian) distributions. As a specific example, the i th observed value of the indentation modulus is modeled as being generated from a deterministic model, with added stochastic noise, as where (0, −1 ) denotes a normal distribution with a zero mean and a variance of −1 . Note that the stochastic noise is assumed to be independent of location in the parameter space, i.e., homoscedastic. The likelihood for a set of independently observed indentation moduli can be established using the product rule as As noted earlier, the model parameters are also treated as stochastic variables. The prior belief on these variables is assumed to be specified by a normal distribution with a zero mean and a large variance of −1 as The application of Bayes' rule (Eq. (3)) to the problem at hand results in where ( | , , , , ) denotes the posterior (updated) distribution on the model parameters. The denominator in Eq. (7) reflects the probability of the observed outcomes irrespective of the model parameters chosen, and can be described by the marginalization of the likelihood with respect to the model parameters as In a fully Bayesian approach, the precision parameters, , , may also be treated as stochastic variables [24]. This allows for a separate application of Bayes' theorem expressed as Alternately, one can use point estimates from the maximization of the likelihood in Eq. (9), denoted as ̂,̂. This is equivalently interpreted as the maximization of the evidence of the observed data in Eq. (8) [25]. With this approach, the posterior distributions of model coefficients in Eq. (7) can be solved analytically (while assuming normal distributions for the various variables involved) [25][26][27]. The updated posterior distribution computed using the approach described above is generally expected to be sharper (i.e., lower variance) compared to the prior belief.
Obviously, the available observations may not produce a posterior distribution that is sharp enough (i.e., the uncertainty associated with the posterior is still too high for a given application). In such cases, one needs to examine carefully where one should produce additional data points (i.e., new observations) in order to maximize the sharpening of the posterior distributions. The general approach to solving this problem (i.e., identifying the new data points exhibiting the maximum potential for improving the model accuracy and reliability) involves making predictions for new inputs, and identifying the specific inputs that exhibited the highest variance (i.e., uncertainty) in their predictions as the locations where new observations should be generated [28,29]. This kind of a rational approach for deciding where to generate new data points is critical for situations where where ( , ) denote the new inputs. Therefore, the specific set of inputs which exhibit the highest variance for the prediction can be readily identified. Once the set of inputs are identified, and corresponding FE simulation performed, the next step is updating the distribution of model coefficients with the newly acquired observation. The update step to the distribution of the model coefficients is natural using a Bayesian framework in the sense that any knowledge acquired previously can be incorporated through the prior.
The posterior distribution of the parameters can continually be updated as incoming data is sequentially added by setting the prior as the previously inferred posterior distribution of model coefficients as shown in Eq. (11). Updates to the posterior distribution of model coefficients are performed until sufficient model convergence and prediction performance is attained. Model convergence is determined through the change in values of the model coefficients and parameters as data is added. Model performance is evaluated through various error metrics such as the leaveone-out-cross-validation (LOOCV) error [25,26]. Building the reduced-order model and critically evaluating its reliability and robustness completes the first step of the two-step protocol. It should be noted that this is intended to be performed only once for a given class of materials.

Estimating Intrinsic Material Properties from Indentation Measurements
For the second step of the protocol, our goal is to employ the reduced-order model built in the first step together with indentation measurements obtained from a given sample to estimate its intrinsic This likelihood is expressed as The evaluation of the likelihood described in Eq. (13) is performed using the reduced-order model, In practice, a class of algorithms have been developed in order to define these transitions and are referred as Metropolis-Hastings algorithms [24]. In this work, Single Component Metropolis Hastings (SCMH) is applied, which considers component wise transitions [33]. In the algorithm below for a given step , partial updates are performed for the sample for each component until all components are updated.
The basic steps for the implementation of the SCMH algorithm are as follows: 1. Initialize a starting point, 0 , using the best available information While the probability of a proposed transition is described by the proposal distribution ( * ), the probability of accepting the transition is given by ( * ). By assuming a flat prior for ( ), the acceptance probability of a proposed transition is completely specified by the posterior probability of the states evaluated within a normalizing constant using Eq. (13) [24,34]. The variances of the proposal distributions 2 are tuned during the "burn-in" period in order to meet an acceptance rate around ~0.23. Ensuring the acceptance rate lies around 0.23 has been shown to provide efficient convergence of the Markov chain for gaussian posteriors [35]. All of the computations described above were realized using functions readily available in MATLAB [36].

Problem Statement
For our first case study, we revisit the extraction of the single crystal elastic constants { 11 , 12 , 44 } of the bcc metal Fe 3%-Si, which was previously attempted using standard regression techniques. In the previous study [13], a total of 2286 simulations were needed to establish a high-fidelity reduced-order model in the first step of the two-step protocol. The simulated database consisted of the indentation modulus corresponding to 300 distinct sets of cubic stiffness constants within the domain 50 GPa ≤ 11 ≤ 250 GPa, 40 GPa ≤ 12 ≤ 150 GPa and 15 GPa ≤ 44 ≤ 120 GPa across 9 orientations selected within the fundamental zone of the relevant orientation space [9]. We note, these ranges encompass a vast number of cubic metals [40]. It is anticipated that the proposed Bayesian framework will need significantly less number of FE simulations to adequately capture the FE predicted indentation modulus within the same parameter space in a robust reduced-order model.

Model Building Process
The Bayesian model building process enables sequential design strategies through the identification of high value simulations which will best improve the predictive capability of the model. Since a database of simulations is already available, simulations are treated as "unseen" and are sampled based on the determined utility of performing the simulation. Before beginning the sequential design process, an initial set of simulations must be performed to establish an initial model.

Figure 1. Cross validation error of the reduced-order model built in the first step of the two-step protocol for different truncation levels in Eq. (1). The dashed line indicates 300 training points.
For the present study, a set of 123 FE simulations were selected from the previously performed 2286 simulations as this initial set. This initial set was selected to correspond to the boundaries of the intrinsic material parameter space. Following initialization, the reduced-order model in Eq. (1) was considered with different truncation levels of L = 8, 10, 12 for the symmetrized Surface spherical harmonics (differently shaped symbols in Figure 1) and Q = 1, 2, 3 for the maximal degree of the respective Legendre Polynomials [16] (different colors of symbols in Figure 1). The truncation levels of the reduced-order model can be treated as hyperparameters, and must be selected so that we produce the most robust and accurate reduced-order models. Leave-one-outcross-validation (LOOCV) was performed at various times during the update process and plotted in Figure 1 for the different truncation levels considered.  Figure 3 as a parity plot. The resulting mean absolute prediction error of the reduced-order model was found to be 2.16 GPa (see Figure 3) while the LOOCV error was found to be 2.22 GPa (see Figure 1). This is comparable to previous efforts based on standard regression techniques and utilizing the full database of 2286 FE simulations, where the LOOCV error was reported to be in the range of 2-2.5 GPa [9].

Extracting Intrinsic Material Parameters
At any point during the model building process, the Bayesian framework presented in Section 2.1 can be used to sample the posterior distribution on the material parameters via the MCMC approach. In order to accomplish this second step of the proposed framework, one needs to evaluate the likelihood function (see Eq. (13)); this requires the use of the reduced-order model obtained in step (1) as well as the relevant experimental indentation data. The reduced-order model with truncations Q=2, L=10 obtained after using 300 training points (described in Section 3.2) was selected for this example case study. Experimental data, including the mean and associated variance of measured indentation moduli, were previously reported for 11 different grains in a polycrystalline sample of Fe 3%-Si [37]. Using the MCMC procedure described in Section 2.2, 50000 samples were drawn. The resulting multivariate distribution is shown in Figure 4 as three univariate distributions. To recap, in Step (1) of the protocol used a minimal number of finite element simulations to establish a high fidelity reduced-order model. Using experimental data previously reported, [37] the established reduced-order model was used to sample the distribution of elastic constants in Step (2)  techniques. This small difference could be attributed to the experimental measurement errors (in both the indentation protocols employed here as well as the more conventional measurement protocols employed in literature). We further note that it should be possible to further refine the methodology presented here (i.e., Step (2) of the protocol) to identify specific additional grain orientations for indentation measurements that might improve specifically the estimates of 44 by reducing its variance. Such refinements will be pursued in future work.

Problem Statement
In order to demonstrate the versatility of the proposed framework, attention is now turned to the extraction of the elastic constants, = { 11 , 12 , 44 , 33 , 13 }, for the hcp metal CP-Ti (commercially pure titanium) [40]. Unlike the previous case study, a database of previously performed FE simulations was not readily available for this case study. Therefore, FE simulations were designed and performed specifically for this study as demanded by the Bayesian inference framework in the Step (1) of the protocol.
The FE model used for this study is the previously validated Finite Element model [9] developed using the commercial software ABAQUS [41]. The sample mesh consisted of 12,610 C3D8 continuum 3-D elements and is shown in Figure 5. The simulated indents were performed using an analytically defined rigid indenter with a tip radius of 16 m, consistent with the size used in the experiments on single crystal CP alpha-Ti grains reported in literature [1]. The dimensions of the sample mesh were taken as 9.6 m X 9.6 m X 4.8 m. The these were chosen to encompass a large number of hcp metals of future interest to our research [6]. The transverse elastic isotropy of the hcp symmetry implies that the elastic indentation response is dependent solely on the declination angle (Φ) between the indenter axis and c-axis of the hcp crystal. Therefore, one only needs to explore the orientation space defined by 0 ≤ Φ ≤ 2 radians. Our goal will be to employ the sequential design strategy once again to efficiently explore the multi-dimensional parameter space identified above in establishing a reliable and robust reduced-order model for the FE indentation simulations over the entire parameter space of interest.

Model Building Process
As with the previous case study, the truncation parameters (Q, L) are important hyper-parameters in the model building process. Since, these are not known a priori, we need to build reduced-order models with different values of these hyper-parameters and make suitable selections. The basic strategy employed here as follows: (1)    It is important to recognize that the parameter space was purposefully chosen to be applicable to many hcp metals of future interest to our research [6]. Predictions are very good over the chosen parameter space as shown in Figure 8. Therefore, within the defined parameter space, future extraction efforts would no longer necessitate the generation of a new model. Furthermore, there is little value in performing additional simulations within the defined parameter space to attempt to significantly improve the reduced-order model. The convergence of the associated model parameters in Figure 7 provides evidence that the reduced-order model is unlikely change drastically with the introduction of new simulations.
It should be noted that the significantly larger training set needed for this case study compared to the previous case study can be attributed to the following reasons: (i) the present case study involved a six-dimensional input space whereas the previous one involved a five-dimensional input space, (ii) the range of values for each input in this case study were selected to be significantly larger than the previous one, and (iii) the degree of elastic anisotropy and contrast captured in this case study is significantly larger compared to the previous case study. The degree of single crystal elastic anisotropy, can be quantified by the universal elastic anisotropic index, , [44,45] defined as where and are the bulk and shear moduli provided by Voigt and Reuss estimates (indicated by subscript and respectively) of a macroscopically homogenous polycrystalline material with uniform texture [46]. A maximum universal elastic anisotropic index of 7.2 was noted for the earlier cubic case study discussed in this paper, compared to 66.2 encountered in the current hcp case study. It is therefore quite reasonable that the number of training data points needed is significantly higher.

Extracting Intrinsic Material Parameters
The focus is now turned to the sampling of the posterior distribution of the elastic constants, = { 11 , 12 , 44 , 33 , 13 } , via MCMC. Similar to the previous case study, in order to sample from the posterior distribution of the intrinsic material parameters, the likelihood function in Eq. (13) must be computed using the available experimental data and the reduced-order model established in Step (1) (corresponding to truncation levels Q=2, L=4 using a training set of 2200 FE simulation data points). The experimental data for this case study was obtained from a prior openly shared dataset [1]. This data set included indentation moduli for 50 different crystal orientations on a CP-Ti sample. Following the procedure described in Section 2.2, 50000 samples were drawn using the MCMC approach. The resulting posterior distributions are shown in Figure 9 for   This indicates the relative low sensitivity of the indentation modulus to changes in 13 , when compared to the other elastic stiffness constants. As noted in the previous case study, it should be possible to extend the framework presented here to focus exclusively on improving the estimation of 13 [13]. However, such an effort could only be justified after the uncertainty in the literature reported values is rigorously quantified. The variance in the predictions of the surrogate model at selected orientations is compared with the experimental data in Figure 10. relatively much more data is available. Furthermore, this observation suggests that there is much more value in conducting additional tests at the lower declination angles, specifically in the range of 0-0.2 radians, compared to conducting them at the higher declination angles. This could be highly valuable input to the experimentalists for their future studies.

Conclusions
A statistical framework has been presented for the robust extraction of the intrinsic material Correspondence: Surya R. Kalidindi, surya.kalidindi@me.gatech.edu