{'fraction_non_alphanumeric': 0.030700956676730565, 'fraction_numerical': 0.007002291659088429, 'mean_word_length': 4.847389131128363, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Plasticity models of material variability based on uncertainty quantification techniques

The advent of fabrication techniques like additive manufacturing has focused attention on the considerable variability of material response due to defects and other micro-structural 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.


Introduction
Variability of material response due to defects and other micro-structural aspects has been well-known for some time [1,2,3,4]. In many engineering applications inherent material variability has been ignorable, 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 micro-structural features [5,6,7].
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 [8]; 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. [9] clearly shows pronounced variability in the resultant yield and hardening.
Given this current state of the technology, the need to enhance design methodology to account for this variability in order to meet performance thresholds with high confidence is clear. In this work, we leverage tools from uncertainty quantification (UQ) [10,11,12] 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 technique to embed the modeled stochasticity in distributions on the physical parameters of the model itself was developed by Sargsyan, Najm, and Ghanem [13], and in this work we adapt it to model the inherent variability of an AM metal [9]. 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. [7] applied the stochastic reduced order model (SROM) technique [14] to weld failure. The SROM technique has many of the basic components of 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 Sec. 2 of this work, we describe the selected experimental dataset [9] 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 develop the methods necessary to perform Bayesian calibration of the material parameters: selection of prior distributions to represent the state of knowledge prior to calibration, design of the likelihood function that determines how close the model response is to the calibration data, and the Markov chain Monte Carlo sampling needed to evaluate the posterior distribution of the parameters that quantifies their means and uncertainties. In particular, we adapt both the traditional additive error [15] and the newer embedded error [13] UQ methods to the representation of the observed mechanical response; and we develop surrogate models of the full finite element simulation tailored to the elastic-plastic response of interest to facilitate efficient Monte Carlo sampling. In Sec. 5, we provide the results of the surrogate response building and calibration processes in order to compare the two methods in light of the selected data. We also employ sensitivities provided by the surrogate in order to discuss model selection, and make assessments about the importance of the various parameters. In Sec. 6, we discuss the results in light of a simple analytic version of the representation problem that serves to illustrate the flow of the calibration process and emphasizes the attributes that make the embedded noise model particularly suitable to representing inherent material variability. We also describe how the variability models can be used in an enhanced design process. In Sec. 7, we emphasize the innovations of the proposed approach to modeling the mechanical response to microstructural material variability.

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. [9], we have six experimental datasets, each consisting of 120 stress-strain 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 stress-strain curve 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 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 zero strain reference configuration 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. We assume the measurement noise to be 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 thresh-old 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 N b = 6 batches of stress data with N i = {64, 77, 91, 79, 64, 46} curves, respectively. To make the data suitable for the inverse problem of parameter calibration, we interpolate each curve and extract n ε = 151 points over the interval (0,3)% to finally arrive at where i enumerates the batches, k enumerates the N i stress curves within the i-th batch, and σ (i,k) j = σ (i,k) (ε j ) represents the stress measured at the j-th strain value, ε j = 0.03j/(n ε − 1), for the k-th curve of the i-th batch. The resulting dataset is shown in Fig. 1(c).
To expedite development of the appropriate analysis and modeling of materials with significant intrinsic variations, we assume all variability beyond the nearly negligible measurement noise stems from the underlying material response. This will lead to conservative estimates of material variability; however, given relevant data, variations in the as-built geometry could be included in the variability analysis or corrected for in pre-processing of the stress-strain data.

Plasticity Theory
To model the observed behavior which resembles standard von Mises plastic response, we adopt a standard finite deformation framework [16] with a multiplicative decomposition of the deformation gradient into elastic and plastic parts where F e is associated with lattice stretching and rotation, and F p is associated with plastic flow. Following Ref. [17], we assume an additive stored energy potential written in terms of the elastic deformation Here, the elastic volumetric deformation is given by J e = det(F e ) = det(F) since plastic flow is assumed to be isochoric, and the deviatoric elastic deformation is measured byb e = J −2/3 e F e F T e . We associate the elastic constants κ and µ with the bulk modulus and shear modulus, respectively, and relate them to Young's modulus, E, and Poisson's ratio, ν, via the linear elastic relations κ = and µ = E 2(1+ν) . The Kirchhoff stress resulting from the derivative of the stored energy potential, W , is For the inelastic response, we employ a J 2 (von Mises) yield condition between an effective stress derived from s = dev[τ ] and an associated flow stress,  Figure 1: (a) An array of nominally identical micro-tension "dogbone" specimens, (b) experimental data from Boyce et al. [9] color-coded by batch, and (c) the reduced data set used in this work.
Υ, 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. Ref. [17]) 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, Y ; linear hardening coefficient, H, and nonlinear exponential saturation modulus K and exponent B. In tension, the yield strength, Y , determines the onset of plasticity; the hardening coefficient H determines the linear trend of the post-yield behavior; and K, B superpose a more gradual transition in stress-strain from the trend determined by Young's modulus E in the elastic regime to H in the plastic regime. These material parameters form the basis of our analysis of material variability. To be clear, this standard J 2 plasticity model is a coarse-grained representation of the microstructural variations that engender the variability in the mechanical response, with the plastic strain representating 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 [18] 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, original cross-sectional area and gauge length.

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 [15,19,20,21,13,22]. 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 distribution of the parameters p(θ|D, M ) that best fits the available observations D given the model choice M . The width of the distribution depends on the consistency of the model with the data and the amount of data. By using this probabilistic framework and physical interpretations of the parameters, we aim to quantify the material variability.

Bayesian inference for parameter calibration
Consider our model M for stress σ = M (ε; θ) comprised of Eqs. (4)(5)(6)(7)(8), where ε is the independent variable and θ = {E, Y, H, K, B} are the parameters of interest. By setting to {H, K, B} or {K, B} 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 ν = 0.3; however, we allow the Young's modulus, E, to vary, so that the locus of yield points is not constrained to a line.
Bayes rule relates the data and prior assumptions on the parameters into the posterior density of the target parameters as where p(D|θ, M ) is the likelihood of observing the data D given the parameters θ and model M , p(θ|M ) is the prior density on the parameters reflecting our knowledge before incorporating the observations, and p(D|M ) is the evidence. 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. Here, we will employ relatively uninformative uniform prior densities based on plausible ranges for each of the parameters (more details will be given in Sec.5). Experimental data influences the resulting posterior probability only through the likelihood p(D|θ, M ), which is based on the difference between the data D and the model predictions M (ε; θ). The likelihood plays an analogous role to the cost/objective function in traditional fitting/optimization in the sense that it defines what model predictions are close to the 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 which leads questions of model discrepancy, comparison and selection which will be briefly discussed in Sec. 4.3 and Sec. 5. In general, given the complexities of the model M , the posterior density p(θ|D, M ) is not known in closed form and one has to resort to numerical methods to evaluate it. Markov Chain Monte Carlo (MCMC) methods [23,24] provide a suitable way to sample from the posterior density, and to estimate it using, e.g. , kernel density estimation.

Surrogate Model
MCMC sampling of the posterior density involves many sequential evaluations of the model M (ε; θ). Since the finite-element based forward model M (ε; θ) is relatively expensive to query (each tension simulation takes approximately 1 cpu-hour), the inverse problem of determining the parameters θ becomes infeasible. We overcome this by building an efficient, sufficiently accurate surrogate modelM (ε; θ) of the physical response M (ε, θ) over the region of interest with a polynomial chaos expansion (PCe, see A for a brief review of polynomial chaos expansions). The fidelity of the surrogate with the full model will be discussed in detail in Sec. 5.1.
Since rough bounds of each of the parameters can be estimated from the data and knowledge of similar materials, we assume that the model parameters θ follow a uniform prior density p(θ). We construct a corresponding PC expansion of the random vector θ via where ξ is a vector of standard uniform random variables and P θ defines the number of terms in the expansion. A corresponding expansion of the model response can be written asM and serves as a suitable surrogate model which, given {σ I }, can be evaluated by drawing samples from the distribution of ξ and then evaluating the polynomial expansion, Eq. (11). Methods to obtain the coefficients {σ I } are discussed in A. A by-product of this expansion is the direct access to parametric sensitivities over the range of interest. The sensitivities of the response to the input parameters can be used to identify the most influential parameters. This analysis is particularly useful for parameter elimination in problems with large dimensionality, as described in more detail in Sec. 5.

Likelihood Formulation
As mentioned before, the likelihood is the term in Eq. (9) that accounts for the data. To formulate the likelihood, one needs to reason about what data is available and its relationship with the model predictions. From the discussion in Sec. 2, one can argue that the batches are independent. Furthermore, within a given batch, we assume all the N i 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.

Additive error formulation
Consider the k-th stress-strain curve from the i-th batch which consists of a sequence of stress observations {σ . A widely-adopted approach is to express the discrepancy between an observation and surrogate model prediction using an additive noise model as where {η capturing the discrepancy between observations and model predictions at a given ε j . This formulation is predicated on the assumption that the modelM (ε; θ) 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 to assume the errors to be independent and identically distributed (i.i.d.) Gaussian random variables with zero mean, i.e. η j ∼ N (0, ς 2 ), where ς 2 is the variance. This yields the following likelihood where we recall that D (k) i represents the stress observations collected from the k-th stress-strain curve of the i-th batch. By assuming that each curve is independent from another, we can write for the full dataset of the i-th batch. The standard deviation ς can either be fixed in advance, if some knowledge about the experimental process is available, or it can be inferred along with the target parameters θ. Moreover, it can be assumed to be either constant or varying with the data points. The basic additive error formulation can be enriched by adding a term capturing the discrepancy between the model prediction and truth represented by the physical data, leading to whereη represents the discrepancy between the model prediction and truth. A structure for the model error is more difficult to prescribe than that for the data error. In fact, the calibrated model now effectively becomesM (ε, θ) +η(ε), and not simply the originalM (ε, θ). Given that this additional term is not physically associated with the presumed sources of non-measurement variability its applicability outside the training regime is delicate. Lastly, this additive term can yield difficulties because it can lead to violations of physical laws [25,26].

Embedded model discrepancy
A more suitable approach to representing variability embedded in the physical model involves adding the model discrepancy error [27,28] to the parameters where η is an additive noise term akin to that in Eq. (13). In this case, θ +η is a random vector with density and moments to be estimated, whereas η (i,k) j is determined by a priori estimates of measurement noise. The random vector θ +η can be represented with a PCe. For instance, for a single parameter θ we can write θ +η = The problem is thus transformed into a density estimation problem, where our objective is now to estimate α = {α 0 , α 1 , . . .} that parametrize and define 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 strictly 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 aleatoric/irreducible rather than epistemic/reducible. In the conventional case, as more data is taken into account, the width of the posterior density shrinks, tending to a Dirac delta function at the true parameter value assuming informative data and negligible model discrepancy. On the other hand, in the present context of density estimation, the objects of inference are the parameters α, and the posterior density is thus on α. Thus, the more data is taken into account, provided the data is sufficiently informative, the distribution on α narrows while the width of the distribution of the parameters θ remains finite and conforms to the data. For this embedded technique, the model calibration problem thus involves finding the posterior distribution on α via Bayes theorem Eq. (9) where α has been substituted for θ and, again, p(α|D,M ) is the posterior, p(D|α,M ) is the likelihood, and p(α|M ) is the prior. Once the posterior for the parameters is characterized, it can be propagated through the model to obtain the posterior predictive distributions for quantities of interest (namely stress in this case). The key feature of these predictions is that their uncertainty is affected by both parameter and model uncertainties. For brevity, we leave the full mathematical details of this embedded approach, including the likelihood formulation, to B.

Results
In this section, we present the details of the construction of the particular surrogate models from the full finite element plasticity model, their calibration to the experimental data, and the physical interpretation of the resulting predictions and parameter estimates. Most of the numerical results presented below are obtained using the UQ Toolkit [29] package.

Surrogate Model
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, E; yield strength, Y ; hardening modulus, H; saturation modulus, K; and saturation exponent, B, which control the elastic-plastic stress response. We build Legendre-Uniform PC expansions of these parameters by assuming that they are uniformly distributed over a chosen range standard uniform random variables. 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. 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 e.g. Refs. [20,21,22]. The priors for α are constructed such that the target physical parameters have their target priors, Eq. (19).

Two parameter perfect plasticity model
For the two parameter model, the stress is expressed as a function of strain and Young's modulus, E, and yield strength, Y according to where we use the superscript "(2)" to identify this as the two parameter model. As mentioned before, the full tension simulation is expensive to evaluate, so we leverage the PC expansions of the inputs, Eq. (19), to create PCes of the stress at each strain {ε j } value To build this sequence of PCes, {M j , j = 0, n ε −1}, we employ regression using uniform random samples in the two-dimensional space (ξ 1 , ξ 2 ). In particular, we generate 1004 training (1000 samples in the inner domain and four additional ones for the corners) to build the surrogate model, and 250 validation samples to assess its accuracy and check for over-fitting. (We chose regression over the more computationally efficient stochastic collocation via sparse quadrature grids given the superior results in constructing the piece-wise smooth surrogates to be introduced forthwith.) Fig. 2 shows the training samples mapped back to the physical domain (E, Y ) using Eq. (19), and the corresponding plasticity simulations. Given the number of samples, the regression approach is suitable for constructing PC expansions of up to ninth order. Computing a higherorder expansion would lead to an under-determined problem. The ensemble of responses, shown in Fig. 2b, are clearly piece-wise linear, as expected, with a slight downward slope in the post-yield response due to finite deformation effects.
For the two parameter model, a key feature is that the stress-strain behavior has a discontinuity of the first derivative when the behavior switches from elastic to plastic as the stress, which depends on E, exceeds the particular yield Y value. This dependence on both parameters can be observed in Fig. 3, which shows the stress plotted as a function of E and Y at three different strain locations, ε = 0.12, 0.39, 1.95%. These points have been chosen such that the first, ε = 0.12%, is within a regime for which all points in the (E, Y ) space are fully elastic; the second value, ε = 0.39%, is in a mixed elastic-plastic regime, and the third point, ε = 1.95%, identifies a fully plastic regime. Fig. 3 shows that within the elastic and plastic regimes, the behavior is linear, and the separator is a line. We will exploit this observation, and use the plastic strain (which is zero in the elastic regime) as a classifier to sub-divide the training data.
To compare the accuracy of global PC expansions up to ninth order and a piecewise linear surrogate built over the elastic and plastic sub-domains, in   Fig. 4 we observe that a global linear PCe is accurate where the regime is either fully elastic or plastic, but inaccurate in the mixed region where the discontinuity in the response makes a global representation sub-optimal. Also, as we increase the order of the PCe from first to fifth order, the results do not change within the elastic and plastic regions, but improve in the mixed region. However, when the order of the expansion is at least fifth order, the errors do not decrease as rapidly which suggests over-fitting. Lastly, the low-order piece-wise linear surrogate has the lowest error in both norms across the strain range and, hence, it is more suitable than a high order global surrogate for this model response.

Three parameter linear hardening model
Augmenting the two parameter model with the post-yield phenomenology controlled by the hardening modulus (H) results in the three parameter model, σ(ε) = M (3) (ε; {E, Y, H}). In this case, we use a total of 2008 training samples (2000 in the inner domain and the 8 additional ones for the corners), and 500 validation points. Again, we build global polynomial surrogates of increasing order and a piecewise low order polynomial surrogate to represent the response surface. To capture the additional complexity in the post-yield response we used a quadratic PC in the plastic regime of the mixed, piece-wise surrogate which is linear in the elastic regime. The resulting 2 -norm and ∞ -norm errors (not shown for brevity) have the same trends as those for the two parameter model shown in Fig. 4 due to similarity in the slope discontinuity of the response and, likewise, the piece-wise surrogate is the best representation of the full simulation response for linear hardening.

Five parameter saturation hardening model
The five parameter model adds the saturation modulus (K) and saturation exponent (B), yielding σ(ε) = M (5) (ε; {E, Y, H, K, B}). To build and check this surrogate, we collected 5032 training and 500 validation samples. As shown in Fig. 5, the five parameter model adds a smoother elastic-plastic transition through the presence of the saturation modulus and exponent. Fig. 6 shows the relative error based on the 2 -norm and ∞ -norm for the various PCes as a function of the strain. As with the piece-wise surrogate for the three-parameter model, here we use a linear PCe for the elastic and a quadratic PCe for the plastic. Unlike the results for the two simpler models, the piece-wise surrogate does not outperform the global surrogates. In this case, a global PCe of order ≥6 gives the lowest errors. Since the errors are comparable for these

Sensitivity Analysis
As mentioned, one advantage of building a PCe surrogate is that one can obtain global Sobol sensitivities of a target quantity of interest with respect to the input parameters [30]. Here we compute the total sensitivities [30] of the stress using the surrogate built at each strain point. Fig. 7 shows the sensitivities obtained over the range of the surrogate for the: (a) two, (b) three, and (c) five parameter models. The sensitivities are influenced by the range chosen to build the surrogate model. The experimental stress-strain curves show little hardening and so the surrogate was constructed with over a narrow range of hardening parameters, Eq. (19), which appropriately minimizes their importance. Also, it is apparent that the relative importance of the parameters evolves with strain. As expected, Fig. 7 shows that the Young's modulus E is the dominant parameter within the elastic regime. At larger strains, in the transition between the elastic and plastic regimes, E becomes gradually less important and the yield Y starts to dominate. In the more complex models, the hardening H and the saturation parameters K and B play a relative minor role due to the fact that the data displays little hardening and the two parameter model is a good representation of the majority of the stress behavior. In both the three and five parameter models the sensitivity to H is slim and almost negligible, whereas the added K in the five parameter model is clearly not negligible apparently for its role in determining the "knee" in the stress-strain at the elastic-plastic transition. More discussion of this behavior will be given in Sec. 5.3.3.

Calibration
In this section, we discuss and contrast the results obtained from the inverse problem formulated the additive error and those obtained using the embedded  formulation. We use the three parameter, linear hardening model as a reference case, which we discuss in detail, and then show the main results for the other models.

Inversion with the additive error model
We assume the measurement noise to be constant, ς 2 , along the strain axis, i.e. the measurement error does not depend on the strain ε. Hence, the parameters to be inferred are θ = {E, Y, H, ς 2 }. As priors, we choose uniform densities with ranges coinciding with those chosen to build the surrogate model in Eq. (19). For the variance, we choose a uniform prior over the positive axis. This is typically appropriate because the surrogate might not be reliable outside the range where it was computed on. The prior plays a minor role if a substantial amount of data is available, making the problem likelihood-informed rather than prior-informed in the large data limit. We leverage this case to highlight some key features of Bayesian calibration applied to the present problem of representing material variability: first, the dependence on the type of surrogate model used, second, the batch-to-batch differences in the resulting parameters, third, the correlations between the target parameters, and finally, convergence of the results with the amount of data used in the inverse problem. Fig. 8 shows the joint posteriors between the physical parameters, Fig. 8(a,b,c), and the marginalized posterior for the standard deviation of the measurement error, Fig. 8(d), obtained using 10 stress-strain curves randomly chosen from D 3 . Results are shown for each of the surrogates we constructed. The convergence of the resulting PDFs with polynomial order gives us confidence that the higher order global surrogates and the piece-wise surrogate lead to sufficiently accurate posterior parameter distributions. Since the distributions are not skewed, it is apparent that the extracted elastic modulus E is not correlated with the yield Y , or the hardening H, whereas there is a weak negative correlation between yield and hardening. We conjecture that the pre-yield data informs E independent of the other parameters and, likewise, the yield point informs Y ; however, there is a trend in the experimental data, see Fig. 1c, for high yield points to lead to subsequent low post yield slopes and vice versa. The maximum a posteriori (MAP) value of the inferred standard deviation is around 0.027 GPa. This value is larger than the one estimated directly from the experimental data (0.020 GPa). This is expected because we are not accounting for model error and, therefore, the model discrepancy is lumped into the measurement error. Note that the range used for the plots are much smaller than those originally chosen for the surrogate construction, Eq. (19). Fig. 9 shows the joint posteriors among the physical parameter for the separate batches obtained using 40 randomly selected curves from each batch. Clearly, the mean parameters of the batches are quite variable and the distributions are, for the most part, distinct and well-separated. However, the correlation structures are similar, suggesting that the batches behave qualitatively in the same way. Fig. 10 illustrates the convergence trend in the posterior distributions as a function of the number N of curves used in the calibration. The panels of Fig. 10 show samples of the joint posteriors p(E, Y ), p(E, H) and p(Y, H) obtained from the third batch, as a function of the number N of curves used to run the problem. The densities shift and narrow as more data is taken into account, and, given the range of the graphs, it appears that the densities have more-or-less converged which suggests sufficient data has been obtained to characterize the parameters. Fig. 11 shows the predictions using the posterior distribution along with an ensemble of curves from the third batch used to run the inference. The error bars represent the posterior predictive uncertainty stemming from the data noise, while the black solid line represents the mean prediction and the gray band represents the ±2 standard deviations due to posterior uncertainty. In this case, the data set leads to a narrow posterior distribution around the mean curve, which is reflected in a tight gray band. Overall, we see that the model prediction and the superimposed error bars cover the experimental data and, hence, provide a good representation of it. However, the inferred value of the measurement noise is not representative of the actual noise that was estimated from the experimental data since the lack of model discrepancy term amplifies the inferred data noise. Hence, despite the predictive analysis showing that we are able to recover the variability of the data, our conclusion is this is an unsuitable inversion model since the full variability is represented as measurement error. Consequently, predictions made using this model would underestimate and inadequately account for the actual variability in the material.

Inversion with embedded error model
In this section, we discuss some key calibration results obtained using the embedded error model.  results.) For each plot, the blue bars represent the posterior uncertainty from the data noise, the black solid line represents the mean prediction, the dark gray band represents the ±2 standard deviations due to posterior uncertainty, and Figure 11: Posterior predictive results obtained for the additive inversion using the three parameter model along with 40 curves (green) from D 3 . The blue bars represent the posterior uncertainty from the data noise, while the black solid line represents the mean prediction and the gray band represents the ±2 standard deviations due to posterior uncertainty. the light gray band represents the ±2 standard deviations due to model error. We observe that in all cases the experimental data is captured and described well by the model predictions. Here, the variability of the data is mostly described by the prediction uncertainty due to parameter variability, while the data noise is small and comparable to the estimate of the measurement error obtained from the data itself. This is the key difference with respect to the additive results shown in Fig. 11. By accounting for model discrepancy through the embedded error terms, we are able to characterize the variability more properly because it is not artificially lumped in the data noise. The contribution stemming from the posterior uncertainty is again quite small, suggesting that we have accounted for sufficient amount of data. The results for three models show that overall they have similar predictive capabilities.

Comparison of additive and embedded error models
To illustrate the differences between the additive and embedded error-based inference results, Fig. 13 shows the posterior PDFs for each parameter using the data from the third batch obtained from the additive inference (left column) and embedded error approach (right column). We note in both cases, the contribution stemming from the posterior uncertainty of the parameter estimates is again quite small, as shown in Fig. 11 and Fig. 12, suggesting that we have accounted for a sufficiently large amount of data. The key difference in the parameter posterior distributions is the embedded term which explicitly enables the calibrated models to reflect the inherent variability in the material parameters. In contrast, calibration with the additive error term tends to attribute variability to that component of the model which cannot be interpreted as a specific level of variability in any particular parameter. Engineering Stress (GPa) Experimental Data Pred error due to model Pred error due to posterior 2 due to data noise (c) Figure 12: Posterior predictive results obtained with the embedded model-based inversion using the (a) 2, (b) 3, and (c) 5-parameter models using 91 curves (green) from D 3 . The blue bars represent the posterior uncertainty from data noise, the black solid line represents the mean prediction, the light gray band represents the ±2 standard deviations due to model inadequacy and the dark gray band is the posterior uncertainty. Posterior PDFs of the two parameters common to all the models, the Young's modulus E and yield strength Y , highlight this difference in uncertainty quantification strategies. Fig. 13 shows that the embedded approach attributes a significant amount of variability in E based on the data while the additive approach has a small uncertainty. In effect, the distribution in the additive approach converges to the best average value for E, since E is treated as a constant rather than a distribution. All variability arising from data variation, as opposed to measurement error and finite sample size, is attributed to the additive term where it is confused with all other uncertainties. With the embedded approach, the distributions of the elastic modulus E for the three models are broad, flat, and fairly consistent. Likewise, the difference in the two approaches is readily apparent in the distributions of Y shown in Fig. 13. The distributions of the yield strength Y for the two simpler plasticity models are essentially the same, as they treat yield as a well-defined point unlike the more complex 5-parameter, saturation hardening model. When interpreted in light of the model sensitivities shown in Fig. 7 and discussed in Sec. 5.2, Y is well informed by the data, resulting in the narrow distribution and small uncertainties in the additive approach and consistent broad distributions using the embedded approach. The low sensitivity to H in the 3 and 5-parameter models gives rise to broad Gaussian distributions with the additive formulation indicating more informative data is needed. This deficiency also contributes to the qualitative differences in the densities for H and Y with the additive formulation. In contrast to H, the 5 parameter model is sensitive to both the Y and K, as shown in Fig. 7, and yet the data is not sufficient to fully inform each independently. Instead, the calibration results in a broad joint PDF of the two parameters since they effect similar changes in the model's behavior near the elastic-plastic transition. This confounding effect also likely gives rise to the bimodal posterior distributions of these parameters. Lastly, given that the sensitivity to B is essentially negligible, it is not surprising that the prior, restricted by the range of the surrogate, exerts significant influence on the posterior distribution of this parameter in both formulations. Fig. 14 shows 100 posterior predictive realizations obtained using the 5 parameter model for the additive approach Fig. 14(a,c) and embedded approach Fig. 14(b,d) calibrated using D 3 . The top row shows the results obtained by sampling the joint posterior density of the parameters θ, and pushing these samples through the modelM (ε j ; θ) only; while the bottom row shows the results with the contribution of the measurement noise η (3,k) j . It is evident from comparing Fig. 14 to the data shown in Fig. 1b that the embedded error approach yields a suitable representation where the variability of the response is determined by the variability of the material parameters. On the contrary, the additive approach yields a tight envelope of predictions, and the full variability is only captured by the added (and over-estimated) measurement noise. Comparing individual realizations to the curves obtained experimentally, it is clear that the embedded approach with noise yields curves that match quite well the trends observed in the experiment, the classical method deviates considerably

Discussion
Following the discussion begun in Sec. 4.3.2, we will use a simplification of the model M (ε, θ) to summarize the key concepts in this work and help generalize the intuition needed to model variability. Let us consider only the elastic response so that the nominal model is: and limit our attention to data for a single batch in the elastic regime. The embedded model of the data is whereη and η j are mean zero random variables. The additive model omits αη which varies the slope of the stress-strain curve. In both cases, so, in the limit of infinite informative data, both formulations recover the cor- . This is illustrated in the comparison of Fig. 14a with the average of Fig. 14b in the elastic regime. The difference between the two formulations becomes clear when examining the variance at a given ε i and covariance across all samples ε i of the calibrated representations. Recall that the stress-strain data for a single specimen has virtually no measurement noise and yet the stress-strain curves are essentially lines with slopes that vary across a batch. Without the αη term, the simplest additive model, where the sequence of random variables {η j } are assumed independent and identically distributed, obtains a variance around the mean that is determined by the total variance of the dataset. Furthermore, the realizations are not smooth. They have a wide and constant variation around their mean trends and, consequently, have distinct offsets in stress at zero strain, as can been seen in Fig. 14c. In contrast to this model employing only uncorrelated noise η j , the embedded model accounts for most of the variation in the data with a distribution of slopes effected by the αη term. This is consistent with the variations in the dataset which is composed of highly correlated data for each test, i.e. each test gives essentially the same linear relationship between σ j and ε j at every j. The fan-like ensemble of realizations shown in Fig. 14d clearly represents the continuity and the particular type of variation seen in the data.
This basic illustration was motivated by our data where each individual experiment is well-described by the hypothesized model. When that is the case, the example shows that the embedded model better represents the intrinsic material variability and, by extension, the underlying physics. If, however, external measurement noise were the primary source of uncertainty, the additive model would provide a good representation. As in traditional constitutive modeling, a rational choice of how to formulate the representation of observed variability can only be assessed by examination of the data, and then confirming the validity of that choice by comparing synthetically replicated experiments to the observed behavior, as in Fig. 14. This choice can be guided by examining whether or not the apparent noise is correlated with the mean behavior and the model prediction, as was done in this work.
In general, there are three kinds of uncertainty that should be considered during calibration to experimental data: (A) external measurement error, (B) intrinsic variability, and (C) model form error. Given that measurement noise is typically uncorrelated with the underlying physical response it is typically modeled with white noise. Moreover, it is reducible by replicating the experiment and collecting more data in the sense that the posterior distributions of model parameters converge and narrow. In contrast, variability in the material properties cannot be reduced by increased data gathering, although more data will typically better inform the estimated joint distribution of material parameters.
The hallmark of inherent variability is individual experimental curves which are well explained by a model with appropriate physical parameters, but have systematic parametric discrepancies across the set of curves. The embedded error formulation is well-suited to represent this source of variability. Finally, model form error refers to relevant physics which are unincorporated in the model, and manifests itself through discrepancies between model predictions and the actual data. When present, a model of a single realization will display a systematic discrepancy from the data it is trying to emulate. Since this error can confound the determination of the other errors it is crucial to perform model selection, as was done in this work albeit for a dataset that was generally well-represented by all members of the model family. (It should be noted the embedded formulation [13] was originally developed to mitigate this type of error.)

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. [13] (which, as mentioned, was originally developed to address model discrepancy), we can mathematically represent the material variability as material parameters drawn from a well-calibrated joint distribution. This Bayesian approach is consistent with our understanding of microstructural variability and appropriate for UQ studies requiring forward propagation of this variability. In particular, we developed a constitutive model of variability that is amenable to non-intrusive sampling and adaptable to direct evaluation in simulation codes that handle fields of distributions. This enables a robust design methodology that can predict performance margins with high confidence.
Another important contribution of this work is to 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 understand how that variability would manifest in applications. We have demonstrated that the embedded method accurately characterizes the aleatoric uncertainty present in the experimental observations and enables "black-box" engineering UQ analysis. It gives insight into what aspects of a homogeneous, macroscale constitutive model are most strongly affected by microstructural variability and enables quantitative model selection. It is able to distinguish the variability of different batches and therefore assess their relative performance. It demonstrates convergence with increasing data and the convergence of common parameters in a nested hierarchy of models. 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.

A Polynomial chaos expansion
A polynomial chaos expansion (PCe) is a spectral representation of a random variable. Here we provide a brief description of the PCe construction and refer readers to Refs. [31,32,33] for more details. A PCe representation of any realvalued random variable λ with finite variance is an expansion of the form λ = ∞ |I|=0 λ I Ψ I (ξ 1 , ξ 2 , . . .), where ξ I are independent identically distributed (i.i.d.) standard random variables, λ I are the coefficients, I = (I 1 , I 2 , . . .) ∀I j ∈ N 0 is an infinite-dimensional multi-index, |I| = I 1 + I 2 + . . . is the 1 norm, and Ψ I are multivariate normalized orthogonal polynomials written as products of univariate orthonormal polynomials: The basis functions ψ Ij are polynomials of order I j in the independent variable ξ j orthonormal with respect to the probability density p(ξ j ). For instance, if the germ ξ is a standard Gaussian random variable, then the PCe is based on Hermite polynomials. Different choices of ξ j and ψ m are available via the generalized Askey family [32]. The PCe (25) converges to the true random variable λ in the mean-square sense [34]. For computational purposes, the infinite dimensional expansion (25) must be truncated: λ = I∈I λ I Ψ I (ξ 1 , ξ 2 , . . . , ξ ns ), (27) where I is some index set, and n s is some finite stochastic dimension that typically corresponds to the number of stochastic degrees of freedom in the system. For example, one popular choice for I is the total-order expansion of degree p, where I = {I : |I| ≤ p}, see e.g. Ref. [33]. Given the expansion for the input λ(ξ), the PCe for a target quantity of interest Q produced by the model evaluation Q = M (λ) can be wriiten in a similar form Q(ξ) = M (λ(ξ)) = I∈I q I Ψ i (ξ 1 , ξ 2 , . . . , ξ ns ).
Methods to compute PC coefficients are broadly divided into two groups, namely intrusive and non-intrusive [33]. The former involves substituting the expansions into the governing equations, and applying orthogonal projection to the resulting equations, resulting in a larger and modified system for the PCe coefficients. This approach is applicable when one has access to the full forward model and can readily modify the governing equations in the simulator. The other, non-intrusive approach is more generally applicable and involves finding an approximation in the subspace spanned by the basis functions by evaluating the original model many times.
One such non-intrusive method relies on orthogonal projection of the solution and is known as non-intrusive spectral projection (NISP). In general, this integral must be estimated numerically. An alternative method of non-intrusively obtaining PCe coefficients is regression, which involves solving the linear system:    Ψ I 1 (ξ (1) ) · · · Ψ I K (ξ (1) ) . . . . . .
where Ψ I n is the nth basis function, q I n is the coefficient corresponding to that basis, and ξ (m) is the mth regression point. In the regression matrix A each column corresponds to a basis element and each row corresponds to a regression point from the training set.

B Embedded discrepancy
As discussed in Ref. [27], the embedded discrepancy likelihood often involves highly nonlinear and near-degenerate features, thus forcing one to find an alternative way to approximate it in a computationally feasible manner. Sargsyan et al. [27] suggest several options based on the assumption of conditional independence between the data points. In this work, we rely on the Gaussian approximation to the marginalized likelihood, which for the i-th batch D i , can be written as where µ j (α) ≡ µ(ε j ; α) = E ξ [M (ε j ; (θ +η)(ξ)] and ς 2 j (α) ≡ ς 2 (ε j ; α) = V ξ [M (ε j ; (θ +η)(ξ)] are the mean and variance of the model at fixed α and strain point. These moments are computed by constructing a PCe for the outputs by propagating the PCe of the input argument in Eq. (17):