Bayesian Modeling of Inconsistent Plastic Response due to Material Variability
Abstract
The advent of fabrication techniques such as additive manufacturing has focused attention on the considerable variability of material response due to defects and other microstructural aspects. This variability motivates the development of an enhanced design methodology that incorporates inherent material variability to provide robust predictions of performance. In this work, we develop plasticity models capable of representing the distribution of mechanical responses observed in experiments using traditional plasticity models of the mean response and recently developed uncertainty quantification (UQ) techniques. We demonstrate that the new method provides predictive realizations that are superior to more traditional ones, and how these UQ techniques can be used in model selection and assessing the quality of calibrated physical parameters.
1 Introduction
Variability of material response due to defects and other microstructural aspects has been wellknown for some time [Hill [1963], NematNasser [1999], McDowell et al. [2011], Mandadapu et al. [2012]]. In many engineering applications inherent material variability has been insignificant, and traditionally the design process is based on the mean or lowerbound response of the chosen materials. Material failure is a notable exception since it is particularly sensitive to outliers in the distributions of microstructural features [Dingreville et al. [2010], Battaile et al. [2015], Emery et al. [2015]]. Even in cases of material failure, traditional engineering design approaches have been able to successfully ignore material variability through the use of empirical safety factors. As the engineering community moves towards more physicallyrealistic and efficient designs while maintaining needed safety, it is necessary to replace these empirical safety factors with confidence bounds that account for actual material variability.
Currently, additive manufacturing (AM) is of particular technological interest and provides strong motivation to not only model the mean response of materials but also their intrinsic variability. Additive manufacturing has the distinct advantages of being able to fabricate complex geometries and accelerate the designbuildtest cycle through rapid prototyping [Frazier [2014]]; however, currently, fabrication with this technique suffers from variability in mechanical response due to various sources, including defects imbued by the process, the formation of residual stresses, and geometric variation in the printed parts. As an example, high throughput tensile data from Boyce et al. [Boyce et al. [2017]] clearly shows pronounced variability in the resultant yield and hardening.
The current state of AM technology, and other manufacturing methods with intrinsic variability, e.g. nano and biobased, would clearly benefit from an enhanced design methodology accounting for this variability in order to meet performance thresholds with high confidence. In this work, we leverage tools from uncertainty quantification (UQ) [Le Maître and Knio [2010], Xiu [2010], Smith [2013]] to provide material variability models, realizations, and, ultimately, robust performance predictions.
It is wellknown that any model is an approximation of the physical response of a real system. Typically, models are characterized by many parameters, and thus appropriately tuning them becomes a key step toward reliable predictions. The most common approach to model calibration is leastsquares regression which yields a deterministic result appropriate for design to the mean. Bayesian inference methods provide a more general framework for model calibration and parameter estimation by providing a robust framework for handling multiple sources of calibration information as well as a full joint probability density on the target parameters. Traditionally, Bayesian techniques have been applied in conjunction with additive noise models that are appropriate for modeling external, uncorrelated influences on observed responses. Recently, a novel technique has been developed to embed the modeled stochasticity in distributions on the physical parameters of the model itself [Sargsyan et al. [2015b]], and in this work we adapt it to model the inherent variability of an AM metal [Boyce et al. [2017]]. This is not the only method available in this emerging field of probabilistic modeling of physical processes for engineering applications. There are commonalities between many of the methods. Notably, the work of Emery et al. [Emery et al. [2015]] applied the stochastic reduced order model (SROM) technique [Field et al. [2015]] to weld failure. The SROM technique has many of the basic components of the embedded noise model: a surrogate model of the response to physical parameters, a means of propagating distributions of parameters with Monte Carlo (MC) sampling and computing realistic realizations of the predicted response.
In practice, many plausible models exist to represent the trends in the noisy observational data. Typically, these models aim to capture different aspects of the relevant physical phenomenon or use different paradigms to represent similar aspects. Some of the available models may be overlycomplex representations of the system response in relation to the available observations. This issue is amplified by the use of models to capture the modeling error in addition to the mean trends. In this context, the optimal model can be obtained using Bayesian model selection, where optimality is measured in terms of datafit and model simplicity [Beck [2010], Kass and Raftery [1995], Berger and Pericchi [1996], Verdinelli and Wasserman [1996]].
In Sec. 2 of this work, we describe the selected experimental dataset [Boyce et al. [2017]] that motivates this effort and provides calibration data. This deep dataset provides realworld relevance that a synthetic dataset would not; however, we apply some preprocessing and simplifying assumptions to facilitate the task of developing the methodology. In Sec. 3, we review the basic plasticity theory that provides the basis for the material variability models developed in Sec. 4. In Sec. 4, we present the Bayesian framework for model calibration and selection, including the formulation of the likelihood functions that are required in this context. In particular, we adapt both the traditional additive error [Kennedy and O’Hagan [2001]] and a more novel embedded error [Sargsyan et al. [2015b]] model. In Sec. 5, we provide the results relating to surrogate model construction, parameter estimation and model selection with the available experimental data. In Sec. 6, we emphasize the innovations of the proposed approach to modeling the mechanical response to microstructural material variability.
2 Experimental Data
We focus this work on the analysis of highthroughput, microtension experimental measurements of additively manufactured stainless steel. From the experiments of Boyce et al. [Boyce et al. [2017]], we have six experimental data sets, each consisting of 120 stressstrain curves from the array of nominally identical dogboneshaped specimens shown in Fig. 1(a). (The data from distinct builds of the array are referred to as batches throughout the remainder of the manuscript.) Each stressstrain curve, as shown in Fig. 1(b), is qualitatively similar and behaves in a classically elasticplastic fashion; however, the material displays a range of yield strengths, hardening and failure strengths and some variability in its apparent elastic properties.
To simplify the data and remove some of the uncertainties associated more with the loading apparatus than the material, we omit the preload cycle to approximately 0.2% strain. The remainder of the mechanical response is monotonic tensile loading at a constant strain rate, see Fig. 1(c). We associate the zero strain reference configuration of each sample with the zerostress, mildly worked material resulting from the preload cycle. The resulting stress, , and strain, , values are derived from the customary engineering stress and strain formulas. It is important to note that the both stress and strain are dependent on measurements of the specimen geometry. The gauge section is nominally 1mm1mm4mm with variations of 2–10%, depending on batch. These variations are large relative to traditional manufacturing tolerances and reflect the current stateoftheart in additive manufacturing process and quality control. Physically, the variations in dimensions and effective properties result from incomplete sintering and compaction of the particulate feed material and the subsequent porosity, surface roughness, residual stress and nonideal microstructure. We estimate the stressstrain measurement noise to be approximately Gaussian with 0.009% standard deviation in the strain measurement and 20.0 MPa in the stress measurement based on the analysis of the random variations for individual curves and the noise in their zerostress/zerostrain intercepts. Since we do not try to model failure in this effort, we discard tests that do not reach at least 3% strain. This threshold was chosen to be sufficiently large that each sample curve is well within the plastic regime (and near peak stress), and yet retain sufficient data to enable calibration. This preprocessing yielded batches of stress data with curves, respectively. To make the data suitable for the inverse problem of parameter calibration, we interpolate each curve and extract points over the interval (0,3)% to finally arrive at
(1) 
where enumerates the batches, enumerates the stress curves within the th batch, and represents the stress measured at the th strain value, , for the th curve of the th batch. The resulting data set is shown in Fig. 1(c).
(a) 
(b) 
(c) 
3 Plasticity Theory
To model the observed behavior which resembles classical von Mises plastic response, we adopt a standard finite deformation framework [Simo [1988]] with a multiplicative decomposition of the deformation gradient, , into elastic and plastic parts
(2) 
where is associated with lattice stretching and rotation, and is associated with plastic flow. Following Simo and Hughes [1998], we assume an additive stored energy potential written in terms of the elastic deformation
(3) 
Here, the elastic volumetric deformation is given by since plastic flow is assumed to be isochoric, and the deviatoric elastic deformation is measured by . We associate the elastic constants and with the bulk modulus and shear modulus, respectively, and relate them to Young’s modulus, , and Poisson’s ratio, , via the linear elastic relations and . The Kirchhoff stress resulting from the derivative of the stored energy potential, , is
(4) 
For the inelastic response, we employ a von Mises () yield condition between an effective stress derived from and an associated flow stress, , as
(5) 
The rate independent, associative flow rule is written in the current configuration as the Lie derivative of the elastic left CauchyGreen tensor (cf. Simo and Hughes [1998])
(6) 
The Lagrange multiplier enforces consistency of the plastic flow with the yield surface, obeys the usual KuhnTucker conditions, and can be interpreted as the rate of plastic slip. Finally, we make the flow stress
(7) 
a function of the equivalent plastic strain
(8) 
and the following parameters: initial yield, ; linear hardening coefficient, , and nonlinear exponential saturation modulus and exponent . In tension, the yield strength, , determines the onset of plasticity; the hardening coefficient determines the linear trend of the postyield behavior; and , superpose a more gradual transition in stressstrain from the trend determined by Young’s modulus in the elastic regime to in the plastic regime. These material parameters form the basis of our analysis of material variability. To be clear, this standard J plasticity model is a coarsegrained representation of the microstructural variations that engender the variability in the mechanical response, with the plastic strain covering a wide variety of underlying inelastic mechanisms and the physical definitions of the material parameters shaping our interpretation of the underlying causes of the variable response.
We approximate the tensile test with a boundary value problem on a rectangular parallelepiped of the nominal gauge section with prescribed displacements on two opposing faces and traction free conditions on the remaining faces to effect pure tension. Finite element simulations are performed in Albany [Salinger et al. [2016]] using the constitutive model described in this section. The engineering stress and strain corresponding to that measured in the experiments are recovered from the reaction forces, prescribed displacements, crosssectional area and gauge length.
4 Calibration formulation
In general, a calibration problem involves searching for the parameters of a given model that minimize the difference between model predictions and observed data. In this work, we adopt a Bayesian approach to the calibration problem [Kennedy and O’Hagan [2001], Sivia [1996], Rizzi et al. [2012, 2013], Sargsyan et al. [2015b], Marzouk et al. [2007a]]. In contrast to leastsquares fitting resulting in a single set of parameter values, in a Bayesian perspective the parameters are considered random variables with associated probability density functions (PDFs) that incorporate both prior knowledge and measured data. The choice of Bayesian methods is well motivated by the data, which agree with the chosen model to a high degree but uncertainty is present in the model parameters both within and across all batches. Bayesian calibration results in a joint posterior probability density of the parameters that best fits the available observations given the model choice . The parametric uncertainty reflected in the posterior PDF depends on the consistency of the model with the data and the amount of data. We aim to quantify the material variability using this probabilistic framework and physical interpretations of the parameters.
4.1 Bayesian inference for parameter calibration
Consider our model for the engineering stress , being a full finite element plasticity model with an underlying plasticity model described in Section 3, where is the independent variable and is the vector containing physical parameters as well as auxiliary and nuisance parameters as will be defined later. By setting or to zero we can form a nested sequence of models with 2, 3, or 5 parameters with perfect plastic, linear hardening, or saturation hardening phenomenology, respectively. Given that we only have one dimensional tension data, we fix the Poisson’s ratio ; however, we allow the Young’s modulus, , to vary, so that the locus of yield points is not constrained to a line. We also allow for geometric variability through a nondimensional crosssection correction factor, , that influences the model output and include in . This geometric correction is motivated by the fact that the observed engineering stress data was computed using an average crosssectional areas based on the outer dimensions of each sample, and can be interpreted as the ratio of the effective load bearing area of the sample to its measured average area. The correction factor aims to mitigate the effect of utilizing one nominal value for crosssectional area per sample rather than utilizing a more accurate spatiallyvarying area profile for each sample and the fact that the outer dimensions lead to an overestimate of actual load bearing area of the AM tensile specimens due to the physical imperfections discussed in Sec. 2. Also, since the gauge length, and hence the strain, are relatively error free, it is not included in the calibration parameters.
In a Bayesian setting for model calibration and selection, Bayes’ rule is used to relate the information contained in the data and prior assumptions to the parameters in the form of a posterior probability density function as
(9) 
Here is the likelihood of observing the data given the parameters and model , is the prior density on the parameters reflecting our knowledge before incorporating the observations, and is the model evidence (which we will compute for model selection purposes). It is important to note that the denominator is typically ignored when sampling from the posterior since it is a normalizing factor, independent of , that ensures the posterior PDF to integrate to unity; however, this term, known as the model evidence, plays a central role in model selection, as will be described later. In this context, we will employ relatively uninformative uniform prior densities due to lack of prior knowledge of the model parameters in the present context of the response of AM tensile specimens. Experimental data influences the resulting posterior probability only through the likelihood , which is based on some normalized measure of the distance between the data and the model predictions . The likelihood plays an analogous role to the cost/objective function in traditional fitting/optimization in the sense that it describes the misfit between model predictions and observational data. Specific forms of the likelihood will be discussed in Sec. 4.3. As Eq. (9) suggests, the outcome is conditioned on the model chosen, leading to questions regarding model comparison and selection which will be discussed in Sec. 4.4. In general, given the complexities of the model , the posterior density is not known in closed form and one has to resort to numerical methods to evaluate it. Markov chain Monte Carlo (MCMC) methods [Gamerman and Lopes [2006], Berg and Billoire [2008]] provide a suitable way to sample from the posterior density, while kernel density estimation, for example, can be used to provide subsequent estimates of the posterior PDF.
4.2 Surrogate Modeling
The sampling of the parameters’ posterior probability density and the evaluation of the Bayesian model evidence involves many evaluations of the computational model. Since the finiteelement based forward model is relatively expensive to query (each tension simulation takes approximately 1 cpuhour), the inverse problem of parameter estimation and model selection via direct evaluation becomes infeasible. Instead, we construct inexpensivetoevaluate, accurate surrogates for the response of interest using polynomial chaos expansion (PCE) [Wiener [1938], Ghanem and Spanos [1991], Xiu and Karniadakis [2002]]. Marzouk et al. [Marzouk et al. [2007b]] have shown that such surrogates can be effectively constructed using UQ techniques with a presumed uniform density on the parameters’ range of interest.
Since rough bounds of each of the parameters can be estimated from the data and knowledge of similar materials, we represent the unknown parameters using a spectral PCE in terms of a set of independent and identically distributed standard uniform random variables , as in
(10) 
where represents the dimensionality of , are the orthogonal PC basis elements (Legendre polynomials in this case), and defines the number of terms in the expansion. Eq. (10) is, essentially, a linear transformation that maps standard uniform random variables to the unknown parameters over their range of interest. A corresponding expansion of the model response, acting as a polynomialbased surrogate model, can be written as
(11) 
and is constructed as functions of the physical parameters, at each value of engineering strain at which data is available. The PC coefficients for the inputs, , and outputs, , of the model can be obtained, using one of two approaches: Galerkin projection or stochastic collocation. We utilize a nonintrusive stochastic collocation method with regression [Xiu [2010]] to estimate the unknown PC coefficients as it does not require any modification to the existing computational models/simulators. Details on this procedure are given in [Le Maître and Knio [2010]], with an application of PCE surrogate modeling with computationally intensive numerical models in [Khalil et al. [2015]].
4.3 Likelihood Formulation
As mentioned, the likelihood is the term in Bayes’ rule, Eq. (9), that accounts for the data. Physical reasoning based on data is available and its relationship with the model prediction helps to formulate the likelihood. From the discussion in Sec. 2, one can argue that the batches are independent. Furthermore, within a given batch, we assume all the stressstrain curves are independent since each experiment is a selfcontained test, performed on separate specimens, i.e. the variability of each specimen is the result of its specific microstructure.
We consider two different formulations of the likelihood, which lead to different formulations of the inverse problem, and hence models of the material variability. The formulations differ by how they account for measurement noise and other variability, and how they are affected by (systematic) model discrepancies. Since each formulation leads to qualitatively different predictions, interpretations, and realizations, we are interested in how each is able to discriminate material variability from other sources of randomness. In this section and in the Results section we will discuss how, given that plastic strain is a coarse metric of the inelastic deformation in additively manufactured materials, discrepancies between the observed data and the model predictions can be interpreted physically. The results in Sec. 5 will illustrate how the posterior responds to the quantity of data and its variability.
4.3.1 Additive error formulation
Consider the th stressstrain curve from the th batch which consists of a sequence of stress observations obtained at the strain locations . A widelyadopted approach is to express the discrepancy between an observation and surrogate model prediction using an additive noise model, as in
(12) 
where is the sample from the set of random variables capturing the discrepancy between observations and model predictions at a given . This formulation is predicated on the assumption that the model accurately represents the true, physical process occurring with fixed, but unknown, parameters. This a strong assumption (and one of the main deficiencies of this approach) since models are, in general, only approximations of observed behavior. Nevertheless, this is a commonly used method due to its simplicity.
In lieu of a completely characterized measurement error model (which is rarely obtained in practice), it is reasonable and expedient to assume the errors are distributed in an independent and identically distributed (i.i.d.) manner, with an assumed Gaussian probability density of zero mean (unbiased) and a variance parameterized by . This yields the likelihood:
(13) 
where we recall that represents the stress observations collected from the th stressstrain curve of the th batch. With the assumption that each experimental curve is independent from the others, the combined likelihood takes the form
(14) 
for the full data set of the th batch.
The standard deviation can be either fixed in advance (based on prior knowledge) or inferred along with the other unknown parameters . Moreover, it can be assumed to be either constant or varying with the strain value. In this context, represents the joint error attributed to modeling error as well as observational noise.
One can enrich this additive error formulation by adding a term capturing the discrepancy between the model prediction and truth (i.e. model error) separately. A structure for the model error is more difficult to prescribe than that for the measurement error. Given that this additional term is not physically associated with the presumed sources of nonmeasurement (material and other physical) variability, its applicability outside the training regime is tenuous. Furthermore, the additive combination of the two error terms can lead to challenges in disambiguation in a parameter estimation context. Lastly, this additive term can yield difficulties because it can lead to violations of physical laws and constraints [Salloum and Templeton [2014a, b]].
4.3.2 Embedded model discrepancy
In order to avoid the difficulties associated with an additive error formulation, we embed the model error in key parameters, essentially converting the unknown parameters into random variables that introduce variability in model predictions due to their uncertain nature. This approach, as detailed in [Sargsyan et al. [2015b]], represents selected parameters using PCE. In our context, we will assume a uniform distribution for each parameter in Eq. (12), given by the firstorder LegendreUniform PCE:
(15) 
in which represents the mean term and dictates the level of variability in which contributes to the total model error. We will also consider the model selection problem as to whether or not to embed the error in , which amounts to keeping the term or setting it explicitly to zero. These coefficients (one or two per parameter) need to be inferred from the experimental observations. Such embedding of error, along with a classical additive error, may account for both model and measurement errors, respectively. A noticeable advantage of this approach, when compared to employing additive error alone, is that the model error is embedded in the model parameters and hence we can propagate the calibrated model error through numerical simulations to any output of interest. (Note that this use of a PCE is distinct from its use in the surrogate modeling, where the extent was selected to cover the feasible range of the parameters and not inferred, as it is in this context.)
The problem of calibrating such an embedded error model is equivalent to that of estimating the probability density of the parameters in which the error is embedded. Specifically, our objective is to estimate that parametrize the density of . This is in contrast to the conventional use of Bayesian inference for parameter estimation, i.e. additive error formulations, in which one infers the parameter and not its density. Also, the data for our present calibration problem motivates the embedded approach since it suggests the uncertainties are aleatory/irreducible rather than epistemic/reducible.
In this context, the model calibration problem thus involves finding the posterior distribution on via Bayes’ theorem Eq. (9)
(16) 
where has been substituted for , denotes the posterior PDF, is the likelihood PDF, and is the prior PDF. Note that reduces to the classical parameter vector when no error embedding is performed ( and for ). Among the different options detailed in [Sargsyan et al. [2015b]] for the likelihood construction, we employ the marginalized likelihood, which for the th batch , can be written as
(17) 
where
(18) 
and
(19) 
are the mean and variance of the model predicted stress at fixed and strain point . These moments are computed using the quadrature techniques that are commonly relied upon in uncertainty propagation, as in the input characterization, Eq. (15); see [Sargsyan et al. [2015b]] for detailed description of the methodology involved in likelihood evaluation.
4.4 Bayesian model selection
The expansion in Eq. (15) resulting from the embedded model discrepancy approach is general in that it does not impose limitations on which parameters should embed the error. Furthermore, as mentioned in Sec. 4.1, we can model the stress using 2, 3, or 5 physical parameters. We have also chosen to examine the inclusion of a crosssection correction factor, , which could also be chosen as another parameter in which to embed an error. Therefore, there are three physical models that are competing to fit the data, with up to 6 physical parameters, all of which are competing to be bestowed with an embedded error term. That amounts to a total of 88 plausible models that are competing to fit the given data. Since determining the optimal model and optimal parameters for model error embedding a priori is not possible, we perform model selection using Bayes factors [Berger and Pericchi [1996], Verdinelli and Wasserman [1996]]. This method of model selection is a databased approach that selects the optimal model as the one that strikes the right balance between datafit and model simplicity [Beck [2010], Kass and Raftery [1995]]. This balance is monitored using the socalled model evidence, or marginal likelihood, for each model , given by
(20) 
The computation of the model evidence is typically neglected in Bayesian calibration as it acts merely as a normalizing factor in Bayes’ rule, Eq. (9). The model evidence acts as a quantitative Ockham’s razor that, when maximized, performs an explicit tradeoff between the datafit and the model simplicity [Beck [2010]]. This is elucidated when examining the logarithm of the model evidence [Muto and Beck [2008], Sandhu et al. [2017]]:
(21) 
where the expectation is with respect to the posterior PDF . The first term on the right hand side of Eq. (21) quantifies the datafit and is known as goodnessoffit. The second term is equal to the relative entropy (or KullbackLeibler divergence) between the prior and posterior PDFs, also known as the information gain [Konishi and Kitagawa [2007]]. Given sufficient data , the information gain is normally higher for more complex models (models with more parameters or with parameters of greater prior PDF support reflecting poor prior knowledge as to their values). Hence, the information gain term quantifies model complexity by examining the relative difference between prior and posterior PDFs. This is in contrast to the frequentist approach where model complexity, in general, only depends on the number of model parameters rather than their relative prior supports [Konishi and Kitagawa [2007]].
Model selection using Bayes factors involves the direct comparison of two plausible models at a time, with Bayes factor defined as the ratio of their evidences given by
(22) 
with indicating that model is more likely while indicates that is the more likely model. One can supplement the evidence with prior belief on these models (called prior model probability) and in that case we would be interested in the ratios of posterior model probabilities. In our context, we assume that all models have equal prior probabilities ahead of analyzing the data and thus Bayes factor is a suitable way to compare the plausible models.
The key challenge is to compute the model evidence efficiently and accurately. In this investigation, we will modify a technique known as adaptive GaussHermite quadrature [Naylor and Smith [1982], Liu and Pierce [1994], Hakim et al. [2018], Khalil and Najm [2018]] by employing importance sampling for integration rather than Gauss quadrature. In an importance sampling framework, we utilize a multivariate normal distribution as a proposal/sampling distribution, with a mean vector and covariance matrix equal to the posterior mean and posterior covariance of , i.e.
(23) 
The model evidence , being an integral over the parameter space, can now be estimated using importance sampling with as the proposal distribution. To do this, we reformulate the evidence as
(24) 
and utilize the change of variable , with being the lower triangular matrix in the Cholesky decomposition of the covariance matrix, , to obtain the following formulation for the evidence:
(25) 
can now be approximated using Monte Carlo sampling strategy as
(26) 
where are Latin Hypercube samples [Le Maître and Knio [2010]] drawn from a standard multivariate normal distribution GaussHermite quadrature points where . For the estimator in Eq. (26) to be accurate, we need good estimates for the posterior mean vector and covariance vector, and ), which can be obtained adaptively using the following iterations ( denoting the iteration number):
(27)  
(28)  
(29) 
For each individual model, these iterations start with an initial guess for the posterior parameter mean vector and covariance matrix. We obtain those from an initial run of a Markov chain Monte Carlo sampler over the parameter posterior. The iterations stop with a termination criterion on . For this application, the termination criteria is met when the relative change in is less than .
5 Results
In this section, we present results relating to the construction of surrogate models from the full finite element plasticity model, their calibration to the experimental data, and the selection of the optimal predictive model using Bayesian model selection. Most of the results are obtained using the UQ Toolkit [Debusschere et al. [2017]].
5.1 Surrogate Model Construction
To describe the material stressstrain behavior we analyze, calibrate and compare three nested plasticity models of increasingly complex phenomenology, namely perfect plasticity, linear hardening and saturation hardening. As mentioned in Sec. 3, we focus on five parameters: Young’s modulus, ; yield strength, ; hardening modulus, ; saturation modulus, ; and saturation exponent, , which control the elasticplastic stress response. We construct surrogate model for this five parameter model, as outlined in Sec. 4.2. In order to construct LegendreUniform PCE surrogates for the model outputs (engineering stresses at a collection of engineering strain values), we assumed that the material parameters are uniformly distributed and characterized by the following firstorder PCE:
(30)  
where are independent identically distributed (i.i.d.) standard uniform random variables. The geometric parameter will be discussed in the next section.
We chose these parameters ranges to be large enough that the corresponding predictions can capture the variability of the experimental data shown in Fig. 1b. The ranges are also chosen so that the constructed surrogate can be used to emulate the two and three parameter model by setting the relevant parameters to zero. Also, we remark that the expansion with i.i.d. random variables is a common step to build the surrogate model. Any correlations between the parameters will then be discovered through the inverse problem, see [Marzouk et al. [2007a], Rizzi et al. [2012, 2013], Khalil et al. [2015], Safta et al. [2017]] for example. To construct and validate this surrogate, we collected training and validation samples.
Fig. 2 shows the stressstrain curves obtained from the corresponding plasticity simulations over the 5032 points in parameter space. Fig. 3 shows the normalized (relative) rootmeansquare error for the various PCE surrogates with increasing order as a function of the strain. In this case, a PCE surrogate of order 6 or greater exhibits the lowest errors. Since the errors do not drop significantly for PCE surrogates of order beyond 6, we are likely to be overfitting the full simulation data when using regression in fitting the coefficients of higher order PCE models. For the 2 (perfect plasticity) and 3 (linear hardening) parameter models, a piecewise surrogate was superior to a global representation due to the lower smoothness of response relative to the 5parameter model. The piecewise surrogate was constructed using plastic strain as a classifier to switch between surrogates for the elastic and plastic regimes.
5.2 Model selection using Bayes’ factor
For illustrative purposes, we will use the third batch of available data as described in Sec. 2. We will reserve the joint analysis of all available batches for possible future investigations due to the observed large and abrupt changes in stressstrain behaviors across batches (as opposed to a more smooth variability within batches) as illustrated in Fig. 1c. As was discussed in Sec. 4.4, the stress can be modeled using 2, 3, or 5 physical parameters, resulting in three physical models that are competing to fit the data. Furthermore, depending on the number of model being analyzed, there are up to 6 parameters (including the crosssection correction factor) to be bestowed with an embedded error term. That amounts to a total of 88 plausible models that are competing to fit the given data.
We apply the methodology outlined in Sec. 4.4 to compete the model evidence of each of these models. The competing models are labeled according to (a) the inclusion/exclusion of physical parameters and (b) whether or not an error is embedded in those physical parameters. For clarity, models will be distinguished by subscript labels. The appearance of the parameter in the label subscript indicates its inclusion in the model, while the appearance of the parameter with a subscript denotes the inclusion of the parameter, as well as an embedded error as in Eq. (15). For example, model is a 5parameter plasticity model with errors embedded in Young’s modulus and yield strength whereas is a 2parameter plasticity model with an error embedded in the crosssection correction factor. For all models, we assume the measurement noise intensity (or standard deviation), , to be an unknown constant along the strain axis, i.e. the measurement error does not depend on the observed strain nor stress values. We therefore infer along with the other parameters. To enforce positivity of in the inference step, we chose to infer the natural logarithm of with a uniform prior on [5,0] in log(). For the material property parameters, we choose uniform priors with ranges coinciding with those chosen to build the surrogate model in Eq. (5.1). As for the crosssectional area correction factor, , we utilized a uniform prior with support on [0,2].
Having computed the model evidence for all 88 plausible models, it was observed that the model evidence was maximized for model and thus we compute Bayes factor with respect to this model as reported in Table 1. We can see that model is the most likely model (based on model evidence), with the secondmost likely model, , being 40% less likely.
The following observations are noteworthy:

The 5 parameter model is the most likely physical model, in comparison to the 3 and 2 parameter models (as described in Sec. 3).

The variability in the mechanical responses as observed from experiments can be effectively and compactly modeled via an embedded error in the yield strength and the crosssectional area correction factor . Physically, this result is plausible given the dimensional variations and the fact that the AM process is likely to affect plastic response and unlikely to alter the effective elastic modulus.

The data provides relatively weak evidence to suggest that the variability can be modeled by embedding the errors in the hardening modulus, , saturation exponent , or Young’s modulus, . This too is plausible given the small variations in hardening and the trend in nearyield behavior, refer to Fig. 1c.
Model  Bayes Factor with respect to 

1  
0.60  
0.14  
0.05  
0.02  
0.01  
Other models 
5.3 Calibration
The previous analysis has resulted in an optimal model that strikes a balance between datamisfit and model simplicity (as indexed by the number of parameters and which are allowed to be distributions). As part of the previous analysis, we sampled from the posterior parameter PDF using adaptive MetropolisHastings Markov chain Monte Carlo sampling technique [Gamerman and Lopes [2006], Berg and Billoire [2008]]. Figure 4 shows the 1D and 2D marginal posterior parameter PDFs, in which and denote the PCE coefficients that account for the embedded errors in the respective parameters. Note the unimodality of the PDFs, implying the existence of a global maximumaposteriori (MAP) estimate for the parameters. There are additional insights that can be made based on further examination of the (marginal) PDFs in Fig. 4:

The crosssectional area correction factor has a mean value of 0.76 and varies with most probability between 0.7 and 0.8. This indicates that the postprocessing that was carried out in converting the applied loads to engineering stresses relied upon overestimates of the sample crosssectional areas, as expected given that the effective load bearing area is at most the area given by the outer dimensions of the gauge section.

There is strong negative correlation between Young’s modulus, , and the crosssectional area correction factor, . This is to be expected given the basic relationship between force, area and stress and, also, the fact that the stress at any strain value is bounded within the envelope of variability.

There is strong negative correlation between the intensity of the embedded errors in yield strength, , and the crosssectional area correction factor, . This is logical since the data suggests a fixed level of variability in the engineering stress at any strain value and both of these terms combine in contributing to that total variability.

Examining the marginal PDF of the intensity of the additive error, , we observe a relatively insignificant MAP estimate close to 0.02 GPa. This indicates that the embedded error terms are effective in capturing the variability in the mechanical response, with an embedded error in the yield strength having a MAP estimate close to 115 GPa (from the marginal posterior PDF of ).
We next examine the differences in the inferred posterior parameter PDFs between the additive and embedded error models. To reiterate, the key difference between the two approaches are the embedded errors which explicitly allow the calibrated models to reflect the inherent variability in the mechanical responses using parametric variability. In contrast, calibration with solely an additive error attributes variability in the observed quantities (i.e. engineering stress) to an additive error term which is not connected directly to the variability of any particular parameter. We focus on the two parameters that are selected to carry an embedded error in the optimal model according to the model selection study (see Sec. 5.2). Those parameters are the yield strength, , and crosssectional area correction factor, . Fig. 5 shows the posterior PDFs for the two parameters obtained using the additive and embedded error approach. We highlight the following key observations:

It is obvious that the embedded approach attributes a significant amount of variability in yield strength based on the data while the additive approach results in smaller uncertainty. In effect, the distribution in the additive approach converges to a delta distribution at the best average value for , since the parameter is treated as an unknown constant rather than a random variable. The width of the distribution in the additive approach characterizes the lack of data (not the inherent material variability), while the embedded approach converges to a (nondegenerate) distribution characterizing the range of that help explain the observed variability.

Likewise, the difference in the two approaches is readily apparent in the distributions of the crosssectional area correction factor shown in Fig. 5. Although the mean estimate for is more or less the same using the two approaches, the uncertainty (reflected in the support of the PDFs) is larger for the embedded error approach and is interpretable as the variation of a physical parameter responsible for the inconsistent mechanical response.
(a)  (b) 
5.4 Predictions under uncertainty
In this section, we will examine the predictions obtained using the optimal model (see previous model selection sections) with embedded model error in yield strength, , and the crosssectional area correction factor, . We also examine the predictions obtained using the 5parameter model with a purely classical additive model error, . We obtain 150 posterior predictive realizations of the stressstrain curves from both models as shown in Fig. 6 reflecting the uncertainty in the calibrated parameters and the embedded error terms (for model ). The top row shows the results obtained by sampling the joint posterior parameter PDF of the parameters , and pushing these samples through the model only; while the middle row shows the results with the contribution of the additive error term. These results highlight the differences between the additive and embedded error approach in characterizing material variability effects on mechanical responses:

It is evident from comparing Fig. 6a,c (additive error) and Fig. 6b,d (embedded error) to the data shown in Fig. 6e that the embedded error approach yields a suitable representation where the variability of the response is determined by the variability of the material parameters (with the aid of embedded parametric distributions).

The additive approach yields a tight envelope of predictions, and the full variability is only captured by the additive error term, which results in overestimation of the predictive uncertainty.

Comparing individual realizations to the curves obtained experimentally, it is clear that the embedded approach (with and without additive error contribution) yields curves that match quite well the trends observed in the experiment by appropriately characterizing material variability and measurement noise, while the classical method deviates considerably and is dominated by the high frequency, uncorrelated white noise.
(a)  (b) 
(c)  (d) 
(e) 
6 Conclusions
We have presented a method which can model the variability of a material which is welldescribed by existing plasticity models of mean response, but contains microstructural variability leading to different macroscopic material properties. By leveraging the embedded error method of Sargsyan et al. [Sargsyan et al. [2015a]] (which, as mentioned, was originally developed to address model discrepancy), we can mathematically represent the variability in mechanical response as material and geometric parameters drawn from a wellcalibrated joint distribution. This Bayesian approach was shown to be consistent with our understanding of microstructural variability and allows physical interpretation of the sources of variability through the parameters. These interpretations can be linked to processing through inference and intuition. The Bayesian methodology is appropriate for UQ studies requiring forward propagation of physical variability and also enables the utilization of Bayesian model selection to arrive at an optimal model that strikes a balance between datamisfit and model complexity. The constitutive model we developed is amenable to nonintrusive sampling (and adaptable to direct evaluation in simulation codes that handle fields of distributions), and thus enables a robust design methodology that can predict performance margins with high confidence.
Another important contribution of this work is the illustration of the contrast the proposed approach with commonly used uncertainty formulations. The standard, additive error formulation appropriately accounts for the uncertainty in the experiments arising from measurement error. Yet, in the case of inherent variability, it characterizes all the uncertainty as measurement error which results in unwarranted confidence in the material properties and an inability to correctly and physically understand how that variability would manifest in applications. We have demonstrated that the embedded method accurately characterizes the aleatory uncertainty present in the experimental observations and enables engineering UQ analysis of material and geometric sources of uncertainty. It gives insight into what aspects of a homogeneous, macroscale constitutive model are most strongly affected by microstructural variability and enables quantitative model selection. Moreover, it is capable of representing both significant external noise and inherent variability in a unified formulation.
In future work, we will extend the methodology to the postnecking failure behavior of similar materials.
Acknowledgments
This work was supported by the LDRD program at Sandia National Laboratories, and its support is gratefully acknowledged. B.L. Boyce would like to acknowledge the support of the Center for Integrated Nanotechnologies. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DENA0003525.
References
 Battaile et al. [2015] C. C. Battaile, J. M. Emery, L. N. Brewer, and B. L. Boyce. Crystal plasticity simulations of microstructureinduced uncertainty in strain concentration near voids in brass. Philosophical Magazine, 95(10):1069–1079, 2015.
 Beck [2010] J. Beck. Bayesian system identification based on probability logic. Structural Control and Health Monitoring, 17(7):825–847, Nov. 2010. ISSN 15452255. doi: 10.1002/stc.424. URL http://doi.wiley.com/10.1002/stc.424.
 Berg and Billoire [2008] B. A. Berg and A. Billoire. Markov chain Monte Carlo simulations. Wiley Online Library, 2008.
 Berger and Pericchi [1996] J. Berger and L. Pericchi. The intrinsic bayes factor for model selection and prediction. Journal of the American Statistical Association, 91:109–122, 1996.
 Boyce et al. [2017] B. L. Boyce, B. C. Salzbrenner, J. M. Rodelas, L. P. Swiler, J. D. Madison, B. H. Jared, and Y.L. Shen. Extremevalue statistics reveal rare failurecritical defects in additive manufacturing. Advanced Engineering Materials, 2017.
 Debusschere et al. [2017] B. Debusschere, K. Sargsyan, C. Safta, and K. Chowdhary. UQ Toolkit. http://www.sandia.gov/UQToolkit, 2017.
 Dingreville et al. [2010] R. Dingreville, C. C. Battaile, L. N. Brewer, E. A. Holm, and B. L. Boyce. The effect of microstructural representation on simulations of microplastic ratcheting. International Journal of Plasticity, 26(5):617–633, 2010.
 Emery et al. [2015] J. M. Emery, R. V. Field, J. W. Foulk, K. N. Karlson, and M. D. Grigoriu. Predicting laser weld reliability with stochastic reducedorder models. International Journal for Numerical Methods in Engineering, 103(12):914–936, 2015.
 Field et al. [2015] R. Field, M. Grigoriu, and J. Emery. On the efficacy of stochastic collocation, stochastic Galerkin, and stochastic reduced order models for solving stochastic problems. Probabilistic Engineering Mechanics, 41:60–72, 2015.
 Frazier [2014] W. E. Frazier. Metal additive manufacturing: a review. Journal of Materials Engineering and Performance, 23(6):1917–1928, 2014.
 Gamerman and Lopes [2006] D. Gamerman and H. F. Lopes. Markov chain Monte Carlo: stochastic simulation for Bayesian inference. CRC Press, 2006.
 Ghanem and Spanos [1991] R. Ghanem and P. Spanos. Stochastic Finite Elements: A Spectral Approach. Springer Verlag, New York, 1991.
 Hakim et al. [2018] L. Hakim, G. Lacaze, M. Khalil, K. Sargsyan, H. Najm, and J. Oefelein. Probabilistic parameter estimation in a 2step chemical kinetics model for ndodecane jet autoignition. Combustion Theory and Modelling, 22(3):446–466, 2018.
 Hill [1963] R. Hill. Elastic properties of reinforced solids: some theoretical principles. Journal of the Mechanics and Physics of Solids, 11(5):357–372, 1963.
 Kass and Raftery [1995] R. E. Kass and A. E. Raftery. Bayes Factors. Journal of the American Statistical Association, 90(430):pp. 773–795, 1995. ISSN 01621459. URL http://www.jstor.org/stable/2291091.
 Kennedy and O’Hagan [2001] M. C. Kennedy and A. O’Hagan. Bayesian calibration of computer models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 63(3):425–464, 2001.
 Khalil and Najm [2018] M. Khalil and H. N. Najm. Probabilistic inference of reaction rate parameters from summary statistics. Combustion Theory and Modelling, 22(4):635–665, 2018.
 Khalil et al. [2015] M. Khalil, G. Lacaze, J. C. Oefelein, and H. N. Najm. Uncertainty quantification in LES of a turbulent bluffbody stabilized flame. Proceedings of the Combustion Institute, 35(2):1147–1156, 2015.
 Konishi and Kitagawa [2007] S. Konishi and G. Kitagawa. Information Criteria and Statistical Modeling (Springer Series in Statistics). Springer, oct 2007.
 Le Maître and Knio [2010] O. Le Maître and O. Knio. Spectral Methods for Uncertainty Quantification. Springer, New York, NY, 2010.
 Le Maître and Knio [2010] O. Le Maître and O. M. Knio. Spectral methods for uncertainty quantification: with applications to computational fluid dynamics. Springer Science & Business Media, 2010.
 Liu and Pierce [1994] Q. Liu and D. A. Pierce. A note on gausshermite quadrature. Biometrika, 81(3):624–629, 1994.
 Mandadapu et al. [2012] K. K. Mandadapu, A. Sengupta, and P. Papadopoulos. A homogenization method for thermomechanical continua using extensive physical quantities. Proc. R. Soc. A, 468(2142):1696–1715, 2012.
 Marzouk et al. [2007a] Y. M. Marzouk, H. N. Najm, and L. A. Rahn. Stochastic spectral methods for efficient bayesian solution of inverse problems. Journal of Computational Physics, 224(2):560 – 586, 2007a.
 Marzouk et al. [2007b] Y. M. Marzouk, H. N. Najm, and L. A. Rahn. Stochastic spectral methods for efficient Bayesian solution of inverse problems. Journal of Computational Physics, 224(2):560–586, 2007b.
 McDowell et al. [2011] D. McDowell, S. Ghosh, and S. Kalidindi. Representation and computational structureproperty relations of random media. JOM Journal of the Minerals, Metals and Materials Society, 63(3):45–51, 2011.
 Muto and Beck [2008] M. Muto and J. Beck. Bayesian Updating and Model Class Selection for Hysteretic Structural Models Using Stochastic Simulation. Journal of Vibration and Control, 14(12):7–34, 2008.
 Naylor and Smith [1982] J. C. Naylor and A. F. Smith. Applications of a method for the efficient computation of posterior distributions. Applied Statistics, pages 214–225, 1982.
 NematNasser [1999] S. NematNasser. Averaging theorems in finite deformation plasticity. Mechanics of Materials, 31(8):493–523, 1999.
 Rizzi et al. [2012] F. Rizzi, O. Knio, H. Najm, B. Debusschere, K. Sargsyan, M. Salloum, and H. Adalsteinsson. Uncertainty Quantification in MD Simulations. Part II: Inference of forcefield parameters. SIAM J. Multiscale Model. Simul., 10(4):1460–1492, 2012.
 Rizzi et al. [2013] F. Rizzi, R. E. Jones, B. J. Debusschere, and O. M. Knio. Uncertainty quantification in md simulations of concentration driven ionic flow through a silica nanopore. ii. uncertain potential parameters. The Journal of Chemical Physics, 138(19):194105, 2013.
 Safta et al. [2017] C. Safta, M. Blaylock, J. Templeton, S. Domino, K. Sargsyan, and H. Najm. Uncertainty quantification in LES of channel flow. International Journal for Numerical Methods in Fluids, 83(4):376–401, 2017.
 Salinger et al. [2016] A. G. Salinger, R. A. Bartlett, A. M. Bradley, Q. Chen, I. P. Demeshko, X. Gao, G. A. Hansen, A. Mota, R. P. Muller, E. Nielsen, et al. Albany: Using componentbased design to develop a flexible, generic multiphysics analysis code. International Journal for Multiscale Computational Engineering, 14(4), 2016.
 Salloum and Templeton [2014a] M. Salloum and J. A. Templeton. Inference and uncertainty propagation of atomisticallyinformed continuum constitutive laws, part 1: Bayesian inference of fixed model forms. International Journal for Uncertainty Quantification, 4(2), 2014a.
 Salloum and Templeton [2014b] M. Salloum and J. A. Templeton. Inference and uncertainty propagation of atomistically informed continuum constitutive laws, part 2: Generalized continuum models based on gaussian processes. International Journal for Uncertainty Quantification, 4(2), 2014b.
 Sandhu et al. [2017] R. Sandhu, C. Pettit, M. Khalil, D. Poirel, and A. Sarkar. Bayesian model selection using automatic relevance determination for nonlinear dynamical systems. Computer Methods in Applied Mechanics and Engineering, 320:237–260, 2017.
 Sargsyan et al. [2015a] K. Sargsyan, H. Najm, and R. Ghanem. On the statistical calibration of physical models. International Journal of Chemical Kinetics, 47(4):246–276, 2015a.
 Sargsyan et al. [2015b] K. Sargsyan, H. N. Najm, and R. Ghanem. On the Statistical Calibration of Physical Models. International Journal for Chemical Kinetics, 47(4):246–276, 2015b.
 Simo and Hughes [1998] J. Simo and T. Hughes. Computational Inelasticity. Springer New York, New York, NY, 1998.
 Simo [1988] J. C. Simo. A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition: Part i. continuum formulation. Computer methods in applied mechanics and engineering, 66(2):199–219, 1988.
 Sivia [1996] D. Sivia. Data Analysis: A Bayesian Tutorial. Oxford Science, 1996.
 Smith [2013] R. C. Smith. Uncertainty quantification: theory, implementation, and applications, volume 12. Siam, 2013.
 Verdinelli and Wasserman [1996] I. Verdinelli and L. Wasserman. Bayes factors, nuisance parameters, and imprecise tests. In J. M. Bernardo, J. O. Berger, A. P. Dawid, and A. F. M. Smith, editors, Bayesian Statistics 5, pages 765–771. Oxford University Press, London, 1996.
 Wiener [1938] N. Wiener. The Homogeneous Chaos. Am. J. Math., 60:897–936, 1938.
 Xiu [2010] D. Xiu. Numerical methods for stochastic computations: a spectral method approach. Princeton university press, 2010.
 Xiu and Karniadakis [2002] D. Xiu and G. Karniadakis. The WienerAskey polynomial chaos for stochastic differential equations. SIAM Journal on Scientific Computing, 24(2):619–644, 2002.