Bayesian Modeling of Inconsistent Plastic Response due to Material Variability

Bayesian Modeling of Inconsistent Plastic Response due to Material Variability

F. Rizzi, M. Khalil, R.E. Jones, J.A. Templeton, J.T. Ostien, B.L. Boyce Sandia National Laboratories, P.O. Box 969, Livermore, CA 94551, USA.Corresponding Author Email:

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 well-known for some time [Hill [1963], Nemat-Nasser [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 lower-bound 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 physically-realistic 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 design-build-test 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 bio-based, 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 well-known 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 least-squares 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 overly-complex 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 data-fit 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 real-world relevance that a synthetic dataset would not; however, we apply some pre-processing 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 high-throughput, micro-tension 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 stress-strain curves from the array of nominally identical dogbone-shaped 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 stress-strain curve, as shown in Fig. 1(b), is qualitatively similar and behaves in a classically elastic-plastic 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 pre-load 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 zero-stress, mildly worked material resulting from the pre-load 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 state-of-the-art 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 non-ideal microstructure. We estimate the stress-strain 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 zero-stress/zero-strain 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 pre-processing 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


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).

Figure 1: (a) An array of nominally identical micro-tension “dogbone” specimens, (b) experimental data from Boyce et al. [Boyce et al. [2017]] color-coded by batch, and (c) the reduced data set used in this work.

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


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


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


For the inelastic response, we employ a von Mises () yield condition between an effective stress derived from and an associated flow stress, , as


The rate independent, associative flow rule is written in the current configuration as the Lie derivative of the elastic left Cauchy-Green tensor (cf. Simo and Hughes [1998])


The Lagrange multiplier enforces consistency of the plastic flow with the yield surface, obeys the usual Kuhn-Tucker conditions, and can be interpreted as the rate of plastic slip. Finally, we make the flow stress


a function of the equivalent plastic strain


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 post-yield behavior; and , superpose a more gradual transition in stress-strain 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 coarse-grained 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, cross-sectional 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 least-squares 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 non-dimensional cross-section 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 cross-sectional 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 cross-sectional area per sample rather than utilizing a more accurate spatially-varying 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


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 finite-element based forward model is relatively expensive to query (each tension simulation takes approximately 1 cpu-hour), the inverse problem of parameter estimation and model selection via direct evaluation becomes infeasible. Instead, we construct inexpensive-to-evaluate, 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


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 polynomial-based surrogate model, can be written as


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 non-intrusive 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 stress-strain curves are independent since each experiment is a self-contained 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 stress-strain curve from the -th batch which consists of a sequence of stress observations obtained at the strain locations . A widely-adopted approach is to express the discrepancy between an observation and surrogate model prediction using an additive noise model, as in


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:


where we recall that represents the stress observations collected from the -th stress-strain curve of the -th batch. With the assumption that each experimental curve is independent from the others, the combined likelihood takes the form


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 non-measurement (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 first-order Legendre-Uniform PCE:


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)


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






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 cross-section 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 data-based approach that selects the optimal model as the one that strikes the right balance between data-fit and model simplicity [Beck [2010], Kass and Raftery [1995]]. This balance is monitored using the so-called model evidence, or marginal likelihood, for each model , given by


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 trade-off between the data-fit 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]]:


where the expectation is with respect to the posterior PDF . The first term on the right hand side of Eq. (21) quantifies the data-fit and is known as goodness-of-fit. The second term is equal to the relative entropy (or Kullback-Leibler 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


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 Gauss-Hermite 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.


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


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:


can now be approximated using Monte Carlo sampling strategy as


where are Latin Hypercube samples [Le Maître and Knio [2010]] drawn from a standard multi-variate normal distribution Gauss-Hermite 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):


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 stress-strain 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 elastic-plastic stress response. We construct surrogate model for this five parameter model, as outlined in Sec. 4.2. In order to construct Legendre-Uniform 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 first-order PCE:


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 stress-strain curves obtained from the corresponding plasticity simulations over the 5032 points in parameter space. Fig. 3 shows the normalized (relative) root-mean-square 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 over-fitting 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 piece-wise surrogate was superior to a global representation due to the lower smoothness of response relative to the 5-parameter model. The piece-wise surrogate was constructed using plastic strain as a classifier to switch between surrogates for the elastic and plastic regimes.

Figure 2: Stress-strain curves samples used to build the surrogate for the five parameter model.
Figure 3: Normalized root-mean-square error in the 5-parameter PCE surrogate for the engineering stress as a function of engineering strain. Results are shown for PCE surrogates of orders up to nine.

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 stress-strain 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 cross-section 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 5-parameter plasticity model with errors embedded in Young’s modulus and yield strength whereas is a 2-parameter plasticity model with an error embedded in the cross-section 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 cross-sectional 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 second-most 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 cross-sectional 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 near-yield behavior, refer to Fig. 1c.

Model Bayes Factor with respect to
Other models
Table 1: Bayes factor as defined in Eq. (22) with respect to

5.3 Calibration

The previous analysis has resulted in an optimal model that strikes a balance between data-misfit 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 Metropolis-Hastings 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 maximum-a-posteriori (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 cross-sectional 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 post-processing that was carried out in converting the applied loads to engineering stresses relied upon over-estimates of the sample cross-sectional 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 cross-sectional 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 cross-sectional 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 ).

Figure 4: 1D and 2D marginals of the posterior parameter PDF for model .

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 cross-sectional 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 (non-degenerate) 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 cross-sectional 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)
Figure 5: Comparison of posterior PDFs obtained using the additive (red, solid curves) and embedded (blue, dashed curves) error models for the parameters (a) yield strength and (b) area correction factor .

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 cross-sectional area correction factor, . We also examine the predictions obtained using the 5-parameter model with a purely classical additive model error, . We obtain 150 posterior predictive realizations of the stress-strain 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 over-estimation 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)
Figure 6: Sample posterior predictive realizations obtained using the 5 parameter model for the additive (a,c) and embedded (b,d) approach. The top row shows the results obtained by sampling the joint parameter posterior density and propagating these samples through the model only, while the middle row shows the results with the contribution of the additive error. The experimental observations are shown in (e) for comparison.

6 Conclusions

We have presented a method which can model the variability of a material which is well-described 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 well-calibrated 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 data-misfit and model complexity. The constitutive model we developed is amenable to non-intrusive 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, macro-scale 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 post-necking failure behavior of similar materials.


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 DE-NA0003525.


  • Battaile et al. [2015] C. C. Battaile, J. M. Emery, L. N. Brewer, and B. L. Boyce. Crystal plasticity simulations of microstructure-induced 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
  • 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. Extreme-value statistics reveal rare failure-critical defects in additive manufacturing. Advanced Engineering Materials, 2017.
  • Debusschere et al. [2017] B. Debusschere, K. Sargsyan, C. Safta, and K. Chowdhary. UQ Toolkit., 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 reduced-order 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 2-step chemical kinetics model for n-dodecane 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
  • 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 bluff-body 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 gauss-hermite 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 structure-property 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(1-2):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.
  • Nemat-Nasser [1999] S. Nemat-Nasser. 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 force-field 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 component-based 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 atomistically-informed 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 Wiener-Askey polynomial chaos for stochastic differential equations. SIAM Journal on Scientific Computing, 24(2):619–644, 2002.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description