Bayesian methods for characterizing unknown parameters of material models

ARTICLE IN PRESS

JID: APM

[m3Gsc;February 16, 2016;21:33]

Applied Mathematical Modelling 000 (2016) 1–17

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Bayesian methods for characterizing unknown parameters of material models J.M. Emery a,∗, M.D. Grigoriu b, R.V. Field Jr. a a b

Sandia National Laboratories, Albuquerque, NM, United States Cornell University, School of Civil & Environmental Engineering, 220 Hollister Hall, Ithaca, NY, United State

a r t i c l e

i n f o

Article history: Received 19 May 2015 Revised 7 January 2016 Accepted 20 January 2016 Available online xxx Keywords: Stochastic reduced order models Bayesian method Uncertain material parameters

a b s t r a c t A Bayesian framework is developed for characterizing the unknown parameters of probabilistic models for material properties. In this framework, the unknown parameters are viewed as random and described by their posterior distributions obtained from prior information and measurements of quantities of interest that are observable and depend on the unknown parameters. The proposed Bayesian method is applied to characterize an unknown spatial correlation of the conductivity field in the definition of a stochastic transport equation and to solve this equation by Monte Carlo simulation and stochastic reduced order models (SROMs). The Bayesian method is also employed to characterize unknown parameters of material properties for laser welds from measurements of peak forces sustained by these welds. © 2016 Elsevier Inc. All rights reserved.

1. Introduction Material properties at the fine scale vary randomly in space, e.g., microstructural conductivity and mechanical properties of weld fusion zones. Mechanical systems are often subjected to actions that vary randomly in space and time, e.g., pressure on the skin of aircrafts and features of most biological tissues. The determination of the response of these materials and systems involves solutions of equations with random coefficients, input, and end conditions, referred to as stochastic equations. Monte Carlo simulation is the only general method for solving stochastic equations irrespective of their complexity. However, its use in applications encounters two obstacles. First, the method is computationally unfeasible if the time for calculating a single solution sample is excessive, especially when many samples are needed to construct reliable solution statistics. Second, the implementation of the method requires full probabilistic characterization of the random entries of these equations, which may not be available in many applications. This paper addresses the first concern by application of stochastic

∗ Corresponding author at: Sandia National Laboratories, P.O. Box 5800 MS 0346, Albuquerque, New Mexico, United States. Tel.: +1 603 488 1555; fax: +1 505 844 4616. E-mail addresses: [email protected] (J.M. Emery), [email protected] (M.D. Grigoriu), rvfi[email protected] (R.V. Field Jr.).

http://dx.doi.org/10.1016/j.apm.2016.01.046 S0307-904X(16)30042-7/© 2016 Elsevier Inc. All rights reserved.

Please cite this article as: J.M. Emery et al., Bayesian methods for characterizing unknown parameters of material models, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.01.046

JID: APM 2

ARTICLE IN PRESS

[m3Gsc;February 16, 2016;21:33]

J.M. Emery et al. / Applied Mathematical Modelling 000 (2016) 1–17

reduced-order models (SROMs) and it addresses the second through development of a Bayesian method for characterizing unknown parameters of material models. It is assumed that the probability laws of the random entries of a stochastic equation have known functional forms but some of their parameters are unknown and unobservable, e.g., parameters of spatial correlations of random field models for material microstructures. We develop a Bayesian framework in which the unknown parameters are viewed as random variables. The initial information on the unknown parameters is captured by prior distributions. Measurements of quantities of interest that are generated numerically and that depend on these parameters are used to update the prior distributions. Updated distributions, referred to as posterior distributions, summarize all available information on the unknown parameters. Stochastic equations depending on unknown parameters characterized by posterior distributions are hierarchical. They are classical stochastic equations that can be solved by existing methods if conditioned on these parameters. We are not alone in employing a Bayesian framework in a mechanics application to inform a model with uncertain parameters. Applications span from heat conduction, to structural health monitoring, to structural dynamics, and include damage of composite materials such as concrete and fiber reinforced structures. Wang and Zabaras used a Bayesian approach for inverse heat conduction problems to ascertain uncertain material properties and boundary conditions [1]. Beck and Yuen use the approach to choose a class of models based on system response data for structural reliability applications [2]. Similarly, Grigoriu and Field use a Bayesian approach for model selection in a frame validation round-robin challenge problem [3]. Bogdanor et al. use a Bayesian approach to introduce uncertainty in initial material defects for multiscale modeling of failure in composites [4]. Vanik et al. use modal data to update model stiffness parameters for structural health monitoring [5]. Simeon et al. describe a Bayesian inference approach to model updating for material parameters in a reinforced concrete beam that used modal predictions and measurements to account for model and measurement uncertainty [6]. This is only a very brief survey of the literature and there are a multitude of other articles discussing the application of Bayes’ theorem to mechanics problems that are naturally awash with uncertainty. In our work, we introduce the following novelties. First, the unknown model parameters are not observable and they need to be inferred from material properties at different scales, e.g., the correlation parameter λ in the first example is inferred from measurement of apparent conductivity. Moreover, when our method is applied to an engineering problem with two variables simultaneously, we discover an unexpected non-uniqueness and we develop an explanation. Second, surrogate models based on SROMs [7] are used to solve stochastic differential equations (SDE) rather than brute force Monte Carlo simulation or approximate methods that cannot capture sample properties. The essentials of the proposed method are outlined in the following section. The method is applied in subsequent sections to solve two stochastic problems. The first is a one-dimensional transport equation with random conductivity whose spatial correlation depends on an unknown parameter. Measurements of apparent conductivity are used to construct posterior distributions for the unknown correlation parameter. The second is a laser weld problem whose performance depends on various mechanical properties, some of which are assumed to be unknown. Measurements of peak weld loads are used to construct posterior distributions for these parameters.

2. Bayesian method Suppose the probability laws of the random entries of a stochastic equation are defined up to a set of unknown parameters that are collected in a vector λ. Conditional on λ, the equation is a classical stochastic equation and can be solved by existing methods, e.g., Monte Carlo simulation. It is assumed that λ cannot be observed directly. Quantities of interest that are observable and depend on λ are used to characterized this vector. We develop a Bayesian framework for characterizing λ. In this framework λ is viewed as a random vector denoted by . Information on the distribution of  is used to construct a prior density f (λ). If only the range of  is known, f (λ) is assumed to be uniform in this range. Since λ is unobservable, quantities of interest that can be measured and depend on λ are used to update f (λ) from:

f  (λ ) ∝ f  (λ ) (λ | data ),

(1)

where f (λ) is the posterior density of , (λ|data) denotes the likelihood function of  corresponding to measurements of the selected quantity of interest, and the symbol ∝ means the ratio of the left- and right-hand sides of Eq. (1) is constant. In this framework, the original stochastic equation becomes a hierarchical stochastic equation that can be solved in two steps. First, independent samples {λi } of  are generated from its posterior density f (λ). Second, solutions of the stochastic equation conditional on { = λi } can be obtained by existing methods for solving stochastic equations [7–13]. They can be subsequently used to find unconditional solution statistics. In Section 3, we illustrate the full implementation of this method for solving stochastic equations with uncertain parameters and develop statistics of quantiles of interest using Monte Carlo simulation and stochastic reduced-order models (SROM). In Section 4, we apply the Bayesian framework to a problem of practical relevance.

Please cite this article as: J.M. Emery et al., Bayesian methods for characterizing unknown parameters of material models, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.01.046

ARTICLE IN PRESS

JID: APM

J.M. Emery et al. / Applied Mathematical Modelling 000 (2016) 1–17

[m3Gsc;February 16, 2016;21:33] 3

3. Random conductivity Let U(x) be the solution of the stochastic differential equation:



d d A (x ) U (x ) dx dx



= 0,

0 < x < l,

(2)

with the boundary conditions U (0 ) = α and U (l ) = β , where A(x) is a real-valued homogeneous random field defined on a probability space (, F, Pr ) that takes values in a strictly-positive, bounded interval on the real line. The condition A(x) > 0 a.s. guarantees that the SDE in Eq. (2) remains elliptical for almost all samples of A(x). The solution of Eq. (2) exists and is unique [14] and has the form:

x

U (x ) = α + (β − α ) 0l

dy/A(y )

0 dy/A (y )

x

= α + (β − α ) 0l

B(y ) dy

0 B (y ) dy

,

0 ≤ x ≤ l,

(3)

where B(x ) = 1/A(x ) is bounded since almost all samples A(x) are strictly positive. The apparent conductivity of a specimen is a random variable defined by:



Aapp (l ) = (1/l )



l 0

B(y ) dy

−1

,

(4)

and depends on properties of B(x), specimen length l, and boundary conditions. It represents the conductivity of a specimen with the same length and boundary conditions as that in Eq. (2) but with space invariant conductivity Aapp (l) that lets through the same current. The solution U(x) in Eq. (3) for this specimen is U (x ) = α + (β − α ) x/l. Since the current through l the specimens with conductivities A(x) and Aapp (l) are A(x ) U  (x ) = (β − α )/ 0 B(y ) dy and Aapp (l ) U  (x ) = Aapp (l ) (β − α )/l, l we have (β − α )/ 0 B(y ) dy = Aapp (l ) (β − α )/l, which yields Eq. (4). 3.1. Parameter characterization Suppose B(x) is a homogeneous Beta translation field defined by:





−1 B(x ) = a + (b − a ) FBeta ◦  G (x ) ,

0 ≤ x ≤ l,

(5)

where 0 < a < b < ∞, FBeta is the distribution of the standard Beta variable with shape parameters (p, q),  denotes the distribution of the standard Gaussian random variable N(0, 1), and  G(x) is a homogeneous Gaussian random field with mean 0, variance 1, and correlation function E[G(x + ξ ) G(x )] = exp − λ |ξ | , where λ > 0 is the correlation decay parameter and ξ is the spatial lag term. The marginal distribution of B(x) is Beta with range [a, b] and shape parameters (p, q) [15, 16, Section 3.1]. The Beta model for B(x) is chosen because it obeys certain physical constraints, unlike other possible choices, since its samples are strictly positive provided the range of the Beta distribution is in (0, ∞). Fig. 1 shows with dash and solid lines the covariance function of G(x) and the covariance function of B(x) scaled by its variance, i.e., c˜(ξ ) = E[B˜(x + ξ ) B˜(x )]/Var[B(x )], where B˜(x ) = B(x ) − E[B(x )]. The covariance functions of these fields nearly coincide for (λ, p, q ) = (1, 2, 3 ) (left panel). For (λ, p, q ) = (5, 1/2, 3 ) (right panel), the covariance function of B(x) is approximately exponential but decays faster than that of G(x). The covariance function of the random field B(x) can be obtained from that of G(x) by using the relationship between these fields given by Eq. (5) and illustrated by Fig. 1. To simplify calculations, we assume c˜(ξ ) = exp − λ |ξ | for all values of λ.

Fig. 1. Covariance functions of G(x) (dash lines) and B(x) scaled by its variance (solid lines) for l = 10, (λ, p, q ) = (1, 2, 3 ) (left panel) and (λ, p, q ) = (5, 1/2, 3 ) (right panel).

Please cite this article as: J.M. Emery et al., Bayesian methods for characterizing unknown parameters of material models, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.01.046

ARTICLE IN PRESS

JID: APM 4

[m3Gsc;February 16, 2016;21:33]

J.M. Emery et al. / Applied Mathematical Modelling 000 (2016) 1–17

Fig. 2. Histograms of B¯ (l ) for l = 10, (λ, p, q ) = (1, 2, 3 ) (left panel) and (λ, p, q ) = (5, 1/2, 3 ) (right panel). Solid lines are densities of standard Gaussian variable.

The probability law of the random field B(x) in Eq. (5) depends on the range [a, b], the shape parameters (p, q), and the correlation decay parameter λ. It is assumed that λ is the only uncertain parameter – a scalar and thus not denoted with bold font. It is an unobservable parameter in the sense that it cannot be measured. However, λ can be inferred from the specimen apparent conductivity Aapp (l), a quantity of interest that can be measured and depends on λ. As previously stated, we view λ as a random variable and construct its distribution from prior information and measurements of apparent conductivity. The apparent conductivity defined by Eq. (4) is the inverse of the spatial average of B(x), or equivalently,

B¯ (l ) =

1 l



l 0

B(x ) dx =

1 , Aapp (l )

(6)

i.e., the spatial average of B(x) over [0, l] coincides with the inverse of the apparent conductivity Aapp (l). The mean and standard deviation of random variable B¯ (l ) are:

μ¯ (l ) = E[B¯ (l )] = E[B(x )] = μ and   2 λ l + exp(−λ l ) − 1 )

= σ g(l, λ ), σ¯ (l ) = Var[B¯ (l )] = σ λl

(7)

where σ denotes the standard deviation of B(x) [3]. These equations show explicitly the dependence of the moments of B¯ (l ) on λ. We note that the variance of B¯ can be made as small as desired by increasing the specimen length for a fixed, unknown λ and that the number of samples needed to estimate accurately λ depends on the specimen length. These properties indicate that robust estimates of λ are always available for any ratio l/λ. Generally, the distribution of the functional B¯ (l ) of the random field B(x) cannot be obtained analytically but can be estimated from samples of B(x). Fig. 2 shows histograms of B¯ (l ) centered and scaled by their means and standard deviations given by Eq. (7) for l = 10, (λ, p, q ) = (1, 2, 3 ) (left panel), and (λ, p, q ) = (5, 1/2, 3 ) (right panel). The solid lines are densities of the standard Gaussian variable. The plots show that B¯ (l ) is approximately Gaussian for a broad range of values of λ and (p, q). To simplify calculation, we assume that B˜(l ) is a Gaussian variable, i.e., B˜(l ) ∼ N (μ ¯ (l ), σ¯ (l )2 ). The implementation of the Bayesian method requires one to specify a prior density f (λ) for , calculate the likelihood function based on the available measurements, and aggregate them to find the corresponding posterior density f (λ) for . We assume that the prior information is limited to the range (λl , λh ), 0 < λl < λh < ∞, of so that f (λ) is constant in this range. The likelihood function is constructed under the assumption that B¯ (l ) is a Gaussian variable. Let {ξ j, k }, k = 1, . . . , n j , denote nj noiseless, independent measurements of B¯ (l j ) for specimens with length lj , j = 1, . . . , J. The likelihood function corresponding to these measurements has the expression:

(λ | {ξ j,k } ) ∝

J 

σ g( l j , λ )

−n j

j=1



exp −

nj 1 2 k=1

 ξ − μ 2  j,k , σ g( l j , λ )

(8)

with the previous notation. The posterior density f (λ) of has the expression:

f  (λ ) ∝ 1(λl , λh ) (λ | {ξ j,k } ).

(9)

Please cite this article as: J.M. Emery et al., Bayesian methods for characterizing unknown parameters of material models, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.01.046

ARTICLE IN PRESS

JID: APM

J.M. Emery et al. / Applied Mathematical Modelling 000 (2016) 1–17

[m3Gsc;February 16, 2016;21:33] 5

under the postulated prior, where 1(C ) = 1 if C is true and zero otherwise. The posterior density f (λ) captures all available information on and can be used to quantify our uncertainty in the unknown parameter λ. It is used in the following section to characterize probabilistically the solution U(x) of Eq. (2). 3.2. Stochastic equations with uncertain parameters ˆ of λ, i.e., It is common to characterize U(x) in Eq. (2) by the solution of this equation corresponding to an estimate λ ˆ ). The approach is computationally efficient but provides no information on the the conditional random field U (x ) | ( = λ sensitivity of the solution of Eq. (2) to the uncertainty in λ. This is a severe limitation when the information on λ is limited so that the range of values of this parameter is large. Our procedure accounts explicitly for the uncertainty in λ and develops credible intervals on solution statistics. For example, the probability Pr(U (x ) ≤ u | ) for fixed but arbitrary x and u is a function of so that it is random variable. The 5%- and 95%-quantiles of the distribution of this random variable define the 90%-credible interval for Pr(U (x ) ≤ u | ). Since x and u are arbitrary, similar calculations yield credible intervals for any other values of these arguments. This section outlines two algorithms for finding statistics of U(x) in Eq. (2) and constructing credible intervals on statistics of U(x) by Monte Carlo simulation and SROM methods. 3.2.1. Monte Carlo simulation We use a hierarchical Monte Carlo algorithm to capture and quantify the uncertainty in λ and its effects on the solution U(x). The algorithm has the following three steps. – Step 1. Uncertainty in : Generate nλ independent samples {λi }, i = 1, . . . , nλ , from the posterior density f (λ) of that characterizes the unknown parameter λ in the Bayesian framework considered in this study. – Step 2. Conditional solutions: Calculate nu independent samples {ui, j (x)}, j = 1, . . . , nu , of U(x) for each sample λi , i = 1, . . . , nλ , of , i.e., samples of the conditional solutions {U (x ) | ( = λi )}. – Step 3. Solution statistics: Use the samples of and U(x) to obtain statistics of quantities of interest. For example, the distribution of U(x) at x = l/2 can be obtained from:





F (u ) = Pr(U (l/2 ) ≤ u ) = E F (u | ) =

 λl λh

F (u | λ ) f  (λ ) dλ,

(10)

where F (u | λ ) = Pr(U (l/2 ) ≤ u | = λ ) denotes the marginal distribution of the conditional random field U (x ) | ( = λ ). The Monte Carlo estimate of F(u) corresponding to the samples generated in the first two steps is: nλ 1 Fˆ (u ) = nλ i=1



nu 1 1(ui, j (l/2 ) ≤ u ) nu



j=1

=

nλ 1 Fˆ (u | λi ), nλ

(11)

i=1

where Fˆ (u | λi ) is an estimate of F(u|λi ). Similar algorithms can be developed for other quantities of interest, e.g., moments of U(x) and apparent conductivity Aapp (l) or other functionals of U(x). Alternatives to Steps 2 and 3 can be used to characterize the solution of Eq. (2). We may generate a single sample ui (x) of U (x ) | ( = λi ) for each sample λi of and use the resulting solutions {ui (x)} to characterize U(x). For example, the distribution in Eq. (10) can be given in the form:

 

F (u ) = E 1 U (l/2 ) ≤ u =



=

Pr( ∈ Ir )

r

Pr( ∈ Ir )F¯ (u | Ir ),



1 Pr( ∈ Ir )





Ir

F (u | λ ) f  (λ ) dλ

(12)

r

by properties of conditional expectation [17, Section 2.17], where {Ir } is a partition of the range of . The scaled integral F¯ (u | Ir ) in the above square bracket is the local average of the conditional distribution F(u|λ) over Ir and can be estimated by:

Fˆ¯ (u | Ir ) =

nλ

 1 1 λi ∈ Ir F (u | λi ), nλ,r /nλ

(13)

i=1

where nλ, r denotes the number of samples {λi } in Ir and Pr( ∈ Ir ) nλ,r /nλ . Both Monte Carlo algorithms can be used to construct credible intervals on statistics of U(x). Suppose our objective is to construct (1 − cu )%-credible intervals for the distribution F(u) in Eq. (10), i.e., which constitute the counterpart of classical confidence intervals in the Bayesian framework [18]. The estimates Fˆ (u | λi ) in Eq. (11) are samples of random variable F(u| ) so that the cu /2-smallest and the cu /2-largest samples {Fˆ (u | λi )} are the limits of the (1 − cu )%-credible interval on F(u). For example, if nλ = 100 and (1 − cu )% = 90%, the credible intervals are given by the fifth smallest and the fifth largest value of {Fˆ (u | λi )}. The Monte Carlo algorithm in Eqs. (12) and (13) can also be applied to develop credible intervals on F(u) by following similar calculations with Fˆ¯ (u | Ir ) in place of Fˆ (u | λi ). The credible intervals are given by the cu /2-smallest and the c /2-largest samples {Fˆ¯ (u | I )}. u

r

Please cite this article as: J.M. Emery et al., Bayesian methods for characterizing unknown parameters of material models, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.01.046

ARTICLE IN PRESS

JID: APM 6

[m3Gsc;February 16, 2016;21:33]

J.M. Emery et al. / Applied Mathematical Modelling 000 (2016) 1–17

3.2.2. SROM method The SROM and Monte Carlo methods are conceptually similar. Differences between these methods relate to the manner in which solution samples are constructed. Following is the algorithm used for calculating solution statistics by the SROM method. – Step 1. Uncertainty in : Use the posterior density f (λ) of or independent samples of this random variable to construct ˜ } and Voronoi cells { k }, k = 1, . . . , m, i.e., non-overlapping intervals centered on {λ ˜ } ˜ for with samples {λ a SROM k k that partition the range  of [7]. It is intended that m be much smaller than the number of samples nλ used in the Monte Carlo simulation described in Section 3.2.1, i.e., m <
k

k

{u˜k, j (x )} and {v˜ k, j (x )}, j = 1, . . . , nu , by following, e.g., developments in [7]. ˜ , viewed as a function of λ and denoted – Step 3. Surrogate model in : The first order Taylor expansion of U(x) about λ k by U(x, λ) to mark its dependence on λ, has the form:





˜ k + (λ − λ ˜ k ) U (x, λ ˜ k ) + V (x, λ ˜ k ) (λ − λ ˜ k ), U (x ) = U (x, λ ) = U x, λ for fixed but arbitrary x, which yields the surrogate model [19]:

U˜ (x ) = U˜ (x, ) =

m





˜ k) , 1( ∈ k ) U˜k (x ) + V˜k (x ) ( − λ

x ∈ (0, l ),

(14)

k=1

in for the solution U(x) of Eq. (2). The model differs from current surrogate models that are deterministic functions for samples of random elements, i.e., samples of in our case. In contrast, U˜ (x ) in Eq. (14) is a random field for samples of . Samples of U˜ (x ) are produced in the following manner. First, we generate nu independent samples, {u˜k, j (x )} and ˜ . Second, samples of U˜ (x ) {v˜ k, j (x )}, j = 1, . . . , nu , of U˜k (x ) and V˜k (x ) for each sample λ˜ k , k = 1, . . . , m, of the SROM

corresponding to nλ independent samples {λi }, i = 1, . . . , nλ , of are calculated from Eq. (14). If λi is in  l , the corre˜ )}, j = 1, . . . , nu . Denote the sample of U˜ (x ) corresponding sponding set of nu samples of U˜ (x ) is {u˜l, j (x ) + v˜ l, j (x ) (λi − λ l to = λi and sample j of the random fields {U˜k (x ), V˜k (x )} by u˜i, j (x ). The algorithm produces nu nλ samples {u˜i, j (x )} of U˜ (x ). – Step 4. Solution statistics: Use the samples of U˜ (x ) to obtain statistics of quantities of interest. For example, the distribution F(u) in Eq. (10) can be estimated as in Eq. (11). Credible intervals for statistics of quantities of interest can be constructed by the approach in Eqs. (10) and (11) since for each sample λi of there are nu samples of U˜ (x ).

In summary, the implementation of the surrogate model in Eq. (14) involves nu deterministic solutions of Eq. (2) and nu ˜ of ˜ . These deterministic solutions of a version of Eq. (2) obtained by differentiation with respect to λ for each sample λ k solution are samples {u˜k, j (x )} and {v˜ k, j (x )} of U˜k (x ) and V˜k (x ), k = 1, . . . , m. Once the surrogate model U˜ (x ) is implemented, elementary calculations yield samples of U˜ (x ). The surrogate model in Eq. (14) is constructed with respect to the random variable . The uncertainty in the random ˜ is represented by samples of {U˜ (x )} and {V˜ (x )} to minimize calculations. field A(x) for = λ k k k The surrogate model in Eq. (14) is for the solution U(x) of Eq. (2). Similar surrogate  models can be developed for various quantities of interest. For example, the random variable F (u | ) = Pr U (l/2, ) ≤ u can be characterized by the surrogate model:

F˜ (u, ) =

m





˜ k ) + Fλ (u | λ ˜ k ) ( − λ ˜ k) , 1( ∈ k ) F (u | λ

k=1

˜ ) = ∂ F (u | λ )/∂λ | ˜ where the sensitivity factor Fλ (u | λ k λ=λ





k

can be approximated by

(15)







˜ + λ ) ≤ u − Pr U (l/2, λ k

˜ ) ≤ u / λ, λ > 0. The probabilities in this finite difference approximation of F (u | λ ˜ ) can be estimated Pr U (l/2, λ k λ k ˜ ˜ from samples of solutions U˜k (x ) = U (x ) | ( = λk ) and U (x ) | ( = λk + λ ). We conclude with the observation that the surrogate model in Eq. (14) is not unique. Alternative surrogate models can be developed by viewing the solution of Eq. (2) as a parametric random field, i.e., a deterministic function of x ∈ (0, l) that depends on and a random vector  whose distribution has -dependent parameters, e.g., we may have  ∼ N( , 1). This surrogate model has to be constructed over a higher dimensional space containing the range of random vector ( , ). 3.2.3. Numerical results To illustrate these methods, let the true value of the correlation decay parameter be λ = 5, but assume it is unknown. The known parameters in the definition of B(x) are a = 1, b = 20, p = 1/2, q = 3 and l = 10 and the boundary conditions are U (0 ) = 0 and U (l ) = 1. The prior information on λ is limited to its range (λl , λh ) = (0.01, 8.0 ). Measurements of apparent conductivity Aapp (l) are available on two specimens with lengths l = 10 and l = 20. These two measurements are samples of B¯ (l ) in Eq. (6) corresponding to λ = 5, i.e., the true value of this parameter. Please cite this article as: J.M. Emery et al., Bayesian methods for characterizing unknown parameters of material models, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.01.046

JID: APM

ARTICLE IN PRESS J.M. Emery et al. / Applied Mathematical Modelling 000 (2016) 1–17

[m3Gsc;February 16, 2016;21:33] 7

Fig. 3. Posterior density of for λ = 5, a = 1, b = 20, p = 1/2, q = 3, and l = 10 based on generated measurements of Aapp (l) on samples with size l = 10 and l = 20.

Fig. 4. Moments of posterior density (left panel) and posterior distributions (right panel) of by Monte Carlo and SROM methods for λ = 5, a = 1, b = 20, p = 1/2, q = 3, and l = 10 based on generated measurements of Aapp (l) on samples with size l = 10 and l = 20.

Fig. 3 illustrates a histogram of the posterior density f (λ) for the prior density f  (λ ) = (λh − λl )−1 1(λl ≤ λ ≤ λh ) and measurements of apparent conductivity on two specimens of length l = 10 and l = 20, i.e., J = 2 and n j = 1 in Eq. (8). The ˜ for . The solid line in the figure is a fit to this histogram that is used to generate samples of and construct a SROM ˜ . The dash and solid lines in the right left panel in Fig. 4 shows with stars and circles the first five moments of and ˜ . The SROM ˜ in the figure has m = 10 samples. panel are the posterior distribution of and of Fig. 5 shows Monte Carlo estimates of the mean (red) and standard deviation (green) of U(x) computed several ways. The heavy solid lines correspond to the case with random with the density f (λ) generated with the algorithm in Section 3.2.1. The heavy dashed lines in the left panel of these figures corresponds to estimates with deterministic λ = 5. The thin solid lines in the right panel of these figures are corresponding SROM estimates generated following the algorithm in Section 3.2.2. The heavy red lines at 45 degrees in Fig. 5 are the Monte Carlo estimate for E[U(x)] with random. The other estimates of E[U(x)] are completely covered by these because they are indistinguishable from the aforementioned at the figure’s scale. Fig. 6 shows Monte Carlo estimates of the skewness and kurtosis of U(x) for random with the density f (λ) and for λ = 5. The heavy solid lines in Fig. 6 are Monte Carlo estimates of skewness (blue) and kurtosis (green) coefficients for random with density f (λ) computed with the algorithm outlined in Section 3.2.1. The heavy dash lines in the left panel are Monte Carlo estimates of skewness (red) and kurtosis (cyan) for λ = 5. The thin lines in the right panel are SROM estimates of skewness (red) and kurtosis (cyan) based on the algorithm in Section 3.2.2. We note that the first two moments of U(x) are not sensitive to the uncertainty in . However, the estimates of the skewness and kurtosis of U(x) for random and λ = 5 differ, which indicates that higher order properties of U(x) are sensitive to the uncertainty in . Fig. 7 shows a Monte Carlo estimate of the distribution Pr(U (l/2, λ ) ≤ u ) for λ = 5 with heavy dashed lines. When is uncertain, Pr(U (l/2, ) ≤ u ) is a random variable as a function of . There are other approaches that could be used. For example, one could formulate a statistical inverse problem and use the maximum likelihood method. In this case, we prefer Please cite this article as: J.M. Emery et al., Bayesian methods for characterizing unknown parameters of material models, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.01.046

JID: APM 8

ARTICLE IN PRESS

[m3Gsc;February 16, 2016;21:33]

J.M. Emery et al. / Applied Mathematical Modelling 000 (2016) 1–17

Fig. 5. Monte Carlo estimates of mean (red) and standard deviation (green) of U(x) computed three ways: with random (heavy solid lines) by the algorithm in Section 3.2.1; with λ = 5 (heavy dashed line, left panel); and by the algorithm in Section 3.2.2 (thin solid lines, right panel). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 6. Monte Carlo estimates of skewness (blue and red) and kurtosis (green and cyan) U(x) computed three ways: with random (heavy solid lines) by the algorithm in Section 3.2.1; with λ = 5 (heavy dashed line, left panel); and by the algorithm in Section 3.2.2 (thin solid lines, right panel). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 7. 90%-credible intervals on the distribution F(u| ) (thin solid lines), E[F(u| )] (heavy solid lines), and F(u| ) (heavy dash lines) by Monte Carlo (left panel) and SROM (right panel).

Please cite this article as: J.M. Emery et al., Bayesian methods for characterizing unknown parameters of material models, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.01.046

JID: APM

ARTICLE IN PRESS

[m3Gsc;February 16, 2016;21:33]

J.M. Emery et al. / Applied Mathematical Modelling 000 (2016) 1–17

9

Fig. 8. The observed force-displacement response of a partial-penetration, laser butt weld under tension (data and inset photo from [20]).

the Bayesian approach because we want to construct credible intervals with the SROM method. Here, we estimate and plot the expected value E[F(u| )] of the distribution Pr(U (l/2, ) ≤ u ) using the algorithms proposed in Sections 3.2.1 and 3.2.2 because we see this as consistent with our Bayesian framework. The heavy solid lines in Fig. 7 are estimates with Monte Carlo (left panel) and SROM (right panel). The thin solid lines are 90%-credible intervals on the distribution of U(l/2) that capture effects of the uncertainty in λ. These credible intervals are relatively large for some thresholds u. 4. Random laser welds Laser welding is popular in engineering applications because it does not require physical contact and, due to the coherence of the energy source, it is well suited for joints near heat sensitive materials. Laser welds are formed by generating enough heat with the focused energy of the laser beam to melt the base metals, with no additional material added during the process. The result is highly variable material properties that contribute to significant variability in the structural response of the welded component. A demonstration of this variability is shown in Fig. 8 where the force-displacement curves are plotted for forty tensile specimens cut from butt welded plates. The butt welding process results in a partial penetration weld at the specimen midplanes and we are careful to note that the curves in Fig. 8 do not represent stressstrain relations. The specimens were nominally identical, but the peak load varies by 13% and the displacement at peak load varies by 59%. In this section, we apply the Bayesian framework to estimate properties of a finite element material model used to predict the strength of laser welds. We begin with a description of the finite element and constitutive models. The strength of laser welds depends on various mechanical properties. We only consider effects of uncertainty in the material initial yield stress y and the hardening constant h on the peak force P in the force-displacement relationship providing a global characterization of weld performance. The peak load is a quantity of interest that is observable and depends on y and h. Two cases are considered. The first assumes that mean value of the hardening parameter is unknown. The second assumes that the mean values of both hardening and initial yield stress are unknown. Our objective is to construct posterior densities for the unknown parameters, i.e., posterior density for mean hardening and initial yield stress, based on measurements of peak force and prior information on the unknown parameters. To achieve this objective, we follow the approach in Eqs. (8) and (9). 4.1. Mechanical and finite element models We develop a three-dimensional finite element model that predicts the load bearing capacity of the weld coupon given information about the uncertainty of parameters in the weld constitutive model. The geometry of the finite element model replicates the nominal dimensions of the free-span tensile specimens, 38.1 × 6.35 × 1.6 mm, accounting for planes of symmetry, and has the nominal 0.76 mm weld penetration as shown in Fig. 9. This finite element mesh contains 6440 linear hexahedral elements. Appropriate kinematic boundary conditions are used for the planes of symmetry and displacements are applied such that tensile loads are generated, opening the weld (refer to Fig. 8 inset). The finite element calculations are completed with Adagio, Sandia National Laboratories’ SIERRA Solid Mechanics quasi-static finite element code [21]. An elastic-plastic material model is used for both the weld region and the bulk material away from the weld. However, in this study, only the model parameters in the weld region were considered uncertain. The plastic flow is idealized to be governed by a yield surface having the form:

σy = y + κ ,

(16)

Please cite this article as: J.M. Emery et al., Bayesian methods for characterizing unknown parameters of material models, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.01.046

ARTICLE IN PRESS

JID: APM 10

[m3Gsc;February 16, 2016;21:33]

J.M. Emery et al. / Applied Mathematical Modelling 000 (2016) 1–17

Fig. 9. (a) Three sides of the 3-D geometry and mesh used for model calibration showing the weld region in light gray and (b) a detail of the mesh and notch at the weld.

where y is the initial yield strength and κ is the isotropic hardening. The evolution of κ is expressed as:

h r

κ ( p ) = [1 − exp (−r p )],

(17)

with h the linear hardening term, r the recovery term and  p the equivalent plastic strain. Only the initial yield strength y and the hardening term h of the weld region are considered to be uncertain with their probability models known up to their mean values. The plasticity parameters for the base material were estimated through prior work [20] and held constant. We assume the elastic properties are fixed with modulus of elasticity = 180 GPa and Poisson’s ratio = 0.27. Details about numerical convergence and treatment of the weld geometry can be found in [22]. 4.2. Unknown hardening Suppose the probabilistic model for material hardening is known to follow a Beta distribution; however, the mean hardening λh is unknown and must be inferred from measurements. We view the unknown parameter λh as random and denote it by := h in our Bayesian framework. It is assumed that (1) takes values in a known interval (ah , bh ), 0 < ah < bh < ∞, and (2) n measurements {pi }, i = 1, . . . , n, of the peak force are available. The measurements are independent samples of the peak force P corresponding to perfectly known probabilistic models for material properties. Since the information on is limited to its range, the prior density f (λh ) of this variable is uniform in (ah , bh ). The measurements of the peak load are used to update f (λh ) to a posterior density f (λh ) via the likelihood function (λ|{pi }) by following the approach in Eqs. (8) and (9). The construction of (λ|{pi }) requires the densities of the conditional random variable P| . For the conductivity problem discussed in Section 3, we bypassed the construction of the density of the conditional variable Aapp (l)| by noting that this variable is approximately Gaussian (see Fig. 2). The Gaussian approximation is not valid for the weld problem under consideration, so that the density of P| has to be constructed numerically. The following three step algorithm has been used to construct the density of P| . – Step 1. Samples of peak force: A mesh with points {λh, i }, i = 1, . . . , nh , has been selected in the postulated range (ah , bh ) ¯ of P generated under the assumption that the mean hardening is λh, i have of . Independent samples {pi, s }, s = 1, . . . , s, been used to construct histograms for the conditional random variables P | ( = λh,i ). Please cite this article as: J.M. Emery et al., Bayesian methods for characterizing unknown parameters of material models, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.01.046

ARTICLE IN PRESS

JID: APM

J.M. Emery et al. / Applied Mathematical Modelling 000 (2016) 1–17

[m3Gsc;February 16, 2016;21:33] 11

Table 1 Parameters for Eq. (18).

λh, i (MPa)

α

β

124864.9 129761.6 134658.3 139554.9 144451.6 149348.3 154244.9 159141.6 164038.2 168934.9

1.78 1.95 2.02 2.24 2.46 2.66 2.76 2.82 2.69 3.00

3.00 3.00 3.00 3.00 3.00 3.00 3.00 2.74 2.52 2.61

Fig. 10. Posterior density f (λh ) with n = 10, 20, 50 and 100.

– Step 2. Model for f(p|λh ): Beta densities,

f ( p | λh ) =

( p − a )α−1 (b − q )β −1 , (b − a )α+β −1 B(α , β )

a ≤ p ≤ b,

(18)

with range (a, b) = (280, 480 ) are fitted to these histograms, where B(α , β ) denotes the Beta function. It is possible to have the same range (a, b) for all λh ∈ (a, b) since the histograms have similar range. Table 1 gives discrete values λh, i for nh = 10 of λh used for calculations and the corresponding estimates of the shape parameters (α , β ). – Step 3. Interpolate: The variation of the shape parameters for an arbitrary λ ∈ (ah , bh ) is obtained from their estimates in Table 1 by interpolation. The resulting shape parameters are used to define the conditional density f(p|λh ) for all λh in (ah , bh ). This conditional density is needed to calculate the likelihood function. We note that f(p|λh ) is the analogue of the densities in Fig. 2. Let f(p|λh ) be the conditional density of P|λh , λ ∈ (ah , bh ), obtained by the above procedure. It has the form in Eq. (18) with shape parameters obtained by interpolating between values in Table 1. The likelihood function corresponding to the  measurements {pi }, i = 1, . . . , n, of the peak load is (λh | { pi } ) = ni=1 f ( pi | λh ) so that the posterior density of has the   form f (λh )∝f (λh ) (λh |{pi }). Numerical results are for nh = 10, s¯ = 100, y = 6.4710e + 04 MPa, r = 9.2394e − 01, and various samples n = 10, 20, 50 and 100. Fig. 10 shows the posterior density f (λh ) for increasing sample size n. These distributions are approximately centered on the true mean hardening 146,825 MPa as determined with the experimental data and become narrower as n increases. Table 2 tabulates estimates of the mean for each n using the mean of the posterior density. 4.3. Unknown hardening and yield stress Suppose the probabilistic models for material hardening and yield stress are known to follow Beta distributions; however, the mean hardening λh and mean yield stress λy are unknown and must be inferred from measurements. We view λh and λy as random variables and denote them by h and y , respectively, so that the uncertainty in the parameters of the probabilistic models of the weld material is captured by the bivariate random vector  = ( h , y ). It is assumed that (1)  takes values in a known two-dimensional rectangle (ah , bh ) × (ay , by ), 0 < ah < bh < ∞, 0 < ay < by < ∞, and (2) n Please cite this article as: J.M. Emery et al., Bayesian methods for characterizing unknown parameters of material models, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.01.046

JID: APM 12

ARTICLE IN PRESS

[m3Gsc;February 16, 2016;21:33]

J.M. Emery et al. / Applied Mathematical Modelling 000 (2016) 1–17 Table 2 Estimated mean found using the posterior density f (λh ) and percent difference from the true mean. The true mean value of hardening is 146,825 MPa.

λh (MPa) n

mean

%-difference

10 20 50 100

153,905 147,931 145,728 146,942

4.8 0.8 −0.7 0.08

Fig. 11. Range of samples of P |  = (λh,i , λy, j ).

measurements {pi }, i = 1, . . . , n, of the peak force are available. The measurements are generated and constitute independent samples of the peak force P corresponding to perfectly known probabilistic models for material properties. Since the postulated prior information on  is limited to the range of this vector, the prior density f (λh , λy ) of  is taken uniform in (ah , bh ) × (ay , by ). Measurements {pi } of the peak force P are used to construct the likelihood function (λh , λy |{pi }) of the unknown parameters and derive the posterior density of  from f (λh , λy )∝f (λh , λy ) (λh , λy |{pi }). As for the previous case, the likelihood function (λh , λy |{pi }) has to be constructed numerically since the distribution of the conditional random variable P| is not available analytically. We follow the algorithm in the previous section to find ¯ j = 1, . . . , j, ¯ is set in the range (ah , bh ) × (ay , by ) of . For (λh , λy |{pi }). First, a mesh with nodes (λh, i , λy, j ), i = 1, . . . , i, ¯ ¯ the calculations here we use i = j = 10. Histograms are constructed for the conditional random variables P |  = (λh,i , λy, j )

¯ of P generated under the assumption  = (λh,i , λy, j ). In our example, we use from independent samples {pij, s }, s = 1, . . . , s, s¯ = 100, resulting in a total of 10,000 calculations for peak load. Second, probability density functions f(p|λh, i , λy, j ) are fitted to the histograms of P |  = (λh,i , λy, j ). The output of this step consists of estimates of the parameters of f(p|λh, i , λy, j ) at the nodes of the mesh in (ah , bh ) × (ay , by ). Third, the estimates of the parameters of f(p|λh, i , λy, j ) are extended to all points in the range of  by interpolation so that f(p|λh, i , λy, j ) becomes available at all points (λh , λy ) in (ah , bh ) × (ay , by ). The construction of meaningful likelihood functions requires that the densities of the conditional random variables P| have the same support for all values of  in (ah , bh ) × (ay , by ). Beta models of the type in the previous section with range covering the samples of P |  = (λh,i , λy, j ) corresponding to all (λh, i , λy, j ) are unsatisfactory since the range of these samples varies significantly in (ah , bh ) × (ay , by ), as illustrated in Fig. 11. Alternative models need to be considered for the density of P|. The left and right panels in Figs. 12 and 13 show estimates of the mean and standard deviation and of skewness and kurtosis coefficients of P |  = (λh,i , λy, j ). All estimates of the skewness coefficient are positive taking values in the range (0.2600, 1.3538) with an average value of 0.7311. The range of kurtosis estimates is (2.1821, 6.6830) with an average value of 3.6272. It is clear that the accuracy of the estimates for these moments is a function of the sample size s.¯ In our example, s¯ = 100 is a practical maximum given the computational demands associated. Due to physical constraints of the system, we propose to use the Gumbel distribution for the densities of the conditional random variables P |  = (λh,i , λy, j ). The density has the form:

f ( p | λh , λy ) = ρ exp[−ρ ( p − ψ ) − e−ρ ( p−ψ ) ],

q ∈ R,

(19)

where ρ and ψ are shape and location parameters [23]. We note that the support of the Gumbel model in Eq. (19) is the entire real line, so that it can be used to construct a meaningful likelihood function, and that the skewness and kurtosis coefficients, 1.14 and 5.4, of this model are in the ranges (0.2600, 1.3538) and (2.1821, 6.6830) of the estimates of these coefficients. The left and right panels in Fig. 14 show estimates of ρ and ψ obtained from sets of 100 independent samples Please cite this article as: J.M. Emery et al., Bayesian methods for characterizing unknown parameters of material models, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.01.046

JID: APM

ARTICLE IN PRESS

[m3Gsc;February 16, 2016;21:33]

J.M. Emery et al. / Applied Mathematical Modelling 000 (2016) 1–17

13

Fig. 12. Estimates of the mean and standard deviation of P |  = (λh,i , λy, j ).

Fig. 13. Estimates of the skewness and kurtosis of P |  = (λh,i , λy, j ).

Fig. 14. Shape and location parameters, ρ and ψ , of the Gumbel model proposed for P |  = (λh,i , λy, j ).

of P |  = (λh,i , λy, j ) corresponding to all pairs (λh, i , λy, j ). The estimates are tabulated in Tables 3 and 4. These estimates, available at the nodes of the mesh in (ah , bh ) × (ay , by ), are extended to all points of this rectangle by interpolation so that the density of P| becomes available at all values of . All numerical results are for nh = 10, r¯ = 100, r = 9.2267e − 01, and various samples n. The second and third columns of Table 5 give the mean of the posterior density f (λh , λy ) for samples of size n = 10, 20, 50, and 100. While the means of f (λh , λy ) for different n are in the vicinity of the (λh , λy ) = (146, 825 MPa, 64, 717 MPa ), they do not seem to lock on the true value of (λh , λy ). This observation is consistent with the contour lines of the posterior density f (λh , λy ) shown in Fig. 15. The evolution of f (λh , λy ) with n illustrated in this figure is rather different from that of f (λh ) in Fig. 10. The posterior density f (λh ) is centered on λh = 146, 825 MPa and becomes more and more concentrated on this value as n increases suggesting that it converges to a Dirac delta function centered on λh = 146, 825 MPa as n → ∞. In contrast, the posterior density f (λh , λy ) seems to converge to a set of equally likely values of (λh , λy ), rather than a point. This prevents Please cite this article as: J.M. Emery et al., Bayesian methods for characterizing unknown parameters of material models, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.01.046

ARTICLE IN PRESS

JID: APM 14

[m3Gsc;February 16, 2016;21:33]

J.M. Emery et al. / Applied Mathematical Modelling 000 (2016) 1–17 Table 3 ρ -parameter for Eq. (19).

λy, j λh, i

55009.2

57166.4

59323.6

61480.8

63638.1

65795.3

67952.5

70109.7

72266.9

74424.2

124800.9 129695.0 134589.2 139483.3 144377.5 149271.6 154165.8 159059.9 163954.1 168848.2

0.0703 0.0569 0.0558 0.0570 0.0551 0.0567 0.0622 0.0548 0.0510 0.0596

0.0715 0.0674 0.0650 0.0621 0.0569 0.0600 0.0666 0.0644 0.0562 0.0627

0.0650 0.0706 0.0627 0.0546 0.0633 0.0544 0.0659 0.0681 0.0557 0.0574

0.0745 0.0761 0.0629 0.0683 0.0718 0.0610 0.0580 0.0592 0.0556 0.0637

0.0601 0.0764 0.0782 0.0679 0.0659 0.0611 0.0633 0.0616 0.0586 0.0600

0.0689 0.0662 0.0758 0.0659 0.0742 0.0665 0.0615 0.0587 0.0644 0.0550

0.0837 0.0712 0.0623 0.0707 0.0705 0.0627 0.0778 0.0618 0.0687 0.0647

0.0810 0.0753 0.0678 0.0702 0.0715 0.0681 0.0646 0.0663 0.0694 0.0701

0.0744 0.0685 0.0651 0.0744 0.0823 0.0566 0.0704 0.0674 0.0670 0.0660

0.0910 0.0757 0.0754 0.0820 0.0719 0.0637 0.0758 0.0534 0.0739 0.0593

Table 4 ψ -parameter for Eq. (19).

λy, j λh, i

55009.2

57166.4

59323.6

61480.8

63638.1

65795.3

67952.5

70109.7

72266.9

74424.2

124800.9 129695.0 134589.2 139483.3 144377.5 149271.6 154165.8 159059.9 163954.1 168848.2

313.5812 318.4115 321.4394 324.5826 330.7019 336.2830 338.0300 342.6851 346.6685 352.7137

319.6717 324.8193 326.4027 331.9098 334.5646 338.8289 343.1455 354.3391 354.4961 357.1680

329.2164 330.3423 338.2309 340.8321 342.1645 345.5706 350.9496 354.1666 362.8863 368.5353

334.9084 336.0162 340.5122 343.3445 346.9049 354.7713 360.5588 360.8818 364.4881 371.0857

340.0535 347.3560 348.8428 353.5720 358.4167 359.2523 365.2896 369.0974 374.1151 377.0911

348.5087 351.8462 352.4214 360.6215 363.9662 370.5042 373.8323 377.3527 377.9014 386.9674

353.1120 355.9409 361.1091 365.4469 370.7983 372.8261 378.7587 384.0073 385.1016 392.5455

360.1857 366.9317 368.5024 373.4443 376.2257 381.5918 386.2165 390.3631 395.4078 393.1355

367.6991 372.7716 377.1450 379.2682 380.5185 390.1921 391.9800 395.2313 399.7547 405.3954

372.7674 377.0613 380.8547 386.3024 390.7058 393.2711 397.6234 403.8614 405.7864 410.5863

Table 5 Estimated means found using the posterior density f (λh , λy ) and root-mean-square difference from the true values. Note the true mean values are λy = 64, 717 MPa and λh = 146, 825 MPa. Mean

Root-mean-square diff.

n

λh (MPa)

λy (MPa)

RMS (MPa)

10 20 50 100

150166 151178 146945 143985

65486 63490 64414 65102

2424.2 3198.0 230.4 2026.6

accurate estimates of the mode of f (λh , λy ). The fourth column in Table 4 gives the root-mean-square difference of the posterior means from the actual means. We attempted to find an explanation for the lack of uniqueness of our procedure. Let P be a real-valued quantity of interest whose distribution depends on a vector λ of uncertain parameters, that we view as a random vector  in the Bayesian framework. Suppose  takes values in a finite, known set {λ1 , . . . , λm }, the members of {λ1 , . . . , λm } are equally likely, and n independent observations {pi }, i = 1, . . . , n, of P are available. The density of the conditional random variable  P| is f ( p | λ ) = m k=1 f k ( p) δ (λ − λk ), where fk (p) is the density of P |  = λk and δ (λ ) = δ (λ1 ), . . . ,δ (λd ). The likelihood functions corresponding to the data {pi } has the expression:

 ( λ | { pi } ) =

n i=1



f ( pi

| λ) =

n m

i=1

k=1



f k ( pi ) δ ( λ − λ k )



=

m n

k=1



f k ( pi )

δ ( λ − λk ),

(20)

i=1

and the posterior density f (λ) of  is proportional with (λ|{pi }) since the prior density f (λ) is uniform in {λ1 , . . . , λm } by assumption. Suppose (λ | { pi } ) = a1 for k = 1, . . . , m , m < m, and (λ | { pi } ) = a2 for k > m and that a1 > a2 . This as    sumption holds if ni=1 fk ( pi ) ni=1 fl ( pi ) for all k, l ∈ {1, . . . , m } and ni=1 fk ( pi ) ni=1 fl ( pi ) for all k, l ∈ {m + 1, . . . , m}.  Under this assumption, f (λ) takes two distinct values, one in {λ1 , . . . , λm } and one in {λm +1 , . . . , λm }. Since f (λ) is maximized by any λ ∈ {λ1 , . . . , λm }, there is no unique solution. We can only state that the actual value of the uncertain parameter λ is a member of the set {λ1 , . . . , λm }. The lack of uniqueness can be eliminated by considering non-uniform prior Please cite this article as: J.M. Emery et al., Bayesian methods for characterizing unknown parameters of material models, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.01.046

JID: APM

ARTICLE IN PRESS J.M. Emery et al. / Applied Mathematical Modelling 000 (2016) 1–17

a

b

c

d

[m3Gsc;February 16, 2016;21:33] 15

Fig. 15. Contour plots for the posterior density f (λh , λy ) with increasing n = 10, 20, 50, and 100. The open circle indicates the mode of the posterior density.

a

b

Fig. 16. (a) The log-likelihood contour plot with two points marked with ‘x’ and ‘o’. (b) The Gumbel distribution plotted for the points indicated with ‘x’ and ‘o’.

Please cite this article as: J.M. Emery et al., Bayesian methods for characterizing unknown parameters of material models, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.01.046

JID: APM 16

ARTICLE IN PRESS

[m3Gsc;February 16, 2016;21:33]

J.M. Emery et al. / Applied Mathematical Modelling 000 (2016) 1–17

Fig. 17. Force-displacement curves for the points marked by ‘x’ and ‘o’.

distributions, provided there is information justifying such an assignment, and/or by considering additional quantities of interest in the analysis to break the tie. The above considerations are consistent with the posterior density f (λh , λy ) in Fig. 15 that, for large sample sizes n, takes large nearly constant values in a subset of its support. Fig. 16(a) shows the natural logarithm of the likelihood function ln((λh , λy |{pi })) and is marked with an ‘x’ and an ‘o’, (67,780 MPa, 133,385 MPa) and (59,608 MPa, 161,645 MPa) respectively, for two points that give nearly the same peak value (mode). Fig. 16(b) plots the conditional densities f(p|λh, i , λy, j ) for these two points. It is evident that the conditional densities are nearly identical for these points. Finite element calculations were run with these (yield, hardening) pairs and their force-displacement are plotted in Fig. 17. This shows that the two samples follow different paths to nearly the same peak loads, 3268 and 3261 N. This suggests differentiation between them could be achieved using a second metric, e.g., the displacement associated with peak force or the energy dissipated up to peak force. 5. Conclusions We developed a Bayesian framework to characterize unknown parameters of probabilistic models for material properties. In this framework, the unknown parameters are viewed as random and described by their posterior distributions obtained from prior information and measurements of quantities of interest that are observable and depend on the unknown parameters. The original stochastic equation becomes a hierarchical stochastic equation that can be solved in two steps. First, independent samples of the unknown parameters are generated from its posterior density. Second, solutions of the stochastic equation conditional on those samples can be obtained by existing methods for solving stochastic equations. In the numerical examples, we first demonstrated the Bayesian framework using a one-dimensional transport equation with random conductivity whose spatial correlation depends on an unknown parameter. Measurements of apparent conductivity that were generated numerically were used to construct posterior distributions for the unknown correlation parameter. Subsequently, the stochastic transport equation was solved by two methods and statistics of the apparent conductivity were computed. This demonstrated the efficacy of stochastic reduced-order models for solving stochastic equations and outlined the hierarchical nature of the Bayesian framework. The second numerical example focused on using the Bayesian framework for a practical engineering problem of interest. Here, the framework was used to develop estimates for uncertain parameters in the constitutive model of a laser weld. It was assumed that information about the uncertain parameters, the hardening modulus and yield stress in a plasticity model, was known up to but excluding their mean values. When one uncertain parameter was considered the framework provided an accurate estimate of the mean value. However, surprisingly, when two values were considered uncertain, the framework appeared to highlight a non-uniqueness. We developed an explanation for the non-uniqueness and suggest it may be circumvented with better prior information or by considering additional quantities of interest. In summary, the following summarizes the novel contributions from our work: •

•

•

The unknown model parameters are not observable and they need to be inferred from material properties at different scales, e.g., the correlation parameter λ in the first example is inferred from measurement of apparent conductivity. Surrogate models based on SROMs are used to solve SDEs rather than Monte Carlo simulation or approximate methods that cannot capture sample properties. When our method was applied to ascertain two means of material properties, we discovered a non-unique solution. This was rather unexpected and we developed an explanation in the discussion around Eq. (20).

Please cite this article as: J.M. Emery et al., Bayesian methods for characterizing unknown parameters of material models, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.01.046

JID: APM

ARTICLE IN PRESS J.M. Emery et al. / Applied Mathematical Modelling 000 (2016) 1–17

[m3Gsc;February 16, 2016;21:33] 17

Acknowledgments Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Company, for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-AC04-94AL85000. References [1] J. Wang, N. Zabaras, Hierarchical Bayesian models for inverse problems in heat conduction, Inverse Probl. 21 (1) (2005) 183–206. [2] J. Beck, K. Yuen, Model selection using response measurements: Bayesian probabilistic approach, J. Eng. Mech. 130 (2) (2004) 192–203. [3] M. Grigoriu, R.V. Field Jr., A solution of the static frame validation challenge problem using Bayesian model selection, Comput. Methods Appl. Mech. Eng. 197 (29–32) (2008) 2540–2549. [4] M.J. Bogdanor, C. Oskay, S.B. Clay, Multiscale modeling of failure in composites under model parameter uncertainty, Comput. Mech. 56 (2015) 389–404. [5] M. Vanik, J. Beck, S. Au, Bayesian probabilistic approach to structural health monitoring, J. Eng. Mech. 126 (7) (2000) 738–745. [6] E. Simoen, G. De Roeck, G. Lombaert, Dealing with uncertainty in model updating for damage assessment: A review, Mech. Syst. Signal Process. 56-57 (2015) 123–149. [7] M. Grigoriu, A method for solving stochastic equations by reduced order models and local approximations, J. Comput. Phys. 231 (19) (2012) 6495–6513. [8] I.M. Babuška, F. Nobile, R. Tempone, A stochastic collocation method for elliptic partial differential equations with random input data, SIAM J. Numer. Anal. 45 (3) (2007) 1005–1034. [9] I.M. Babuška, R. Tempone, E. Zouraris, Galerkin finite element approximations of stochastic elliptic partial differential equations, SIAM J. Numer. Anal. 42 (2) (2004) 800–825. [10] J. Foo, X. Wan, E. Karniadakis, The multi-element probabilistic collocation method (ME-PCM): Error analysis and applications, J. Comput. Phys. 227 (2008) 9572–9595. [11] P. Frauenfelder, C. Schwab, R.A. Todor, Finite element for elliptic problems with stochastic coefficients, Comput. Methods Appl. Mech. Eng. 194 (2005) 205–228. [12] C. Schwab, R.A. Todor, Karhunen–Loéve approximation of random fields by generalized fast multipole methods, J. Comput. Phys. 217 (2006) 100–122. [13] D. Xiu, Numerical Methods for Stochastic Computations: A Spectral Method Approach, Princeton University Press, Princeton, 2010. [14] I.M. Babuška, R. Tempone, E. Zouraris, Solving elliptic boundary value problems with uncertain coefficients by the finite element method: The stochastic formulation, Comput. Methods Appl. Mech. Eng. 194 (12–14) (2005) 1251–1294. [15] R.V. Field Jr., M. Grigoriu, Model selection for random functions with bounded range, applications in science and engineering, in: A. d‘Onofrio (Ed.), Bounded Noises in Physics, Biology and Engineering, Springer, Birkhäuser Basel, 2013, doi:10.1007/978-1-4614-7385-5. [16] M. Grigoriu, Applied Non-Gaussian Processes: Examples, Theory, Simulation, Linear Random Vibration, and MATLAB Solutions, Prentice Hall, Englewoods Cliffs, NJ, 1995. [17] M. Grigoriu, Stochastic Calculus. Applications in Science and Engineering, Birkhäuser, Boston, 2002. [18] A. Zellner, An Introduction to Bayesian Inference in Econometrics, John Wiley & Sons, Inc., New York, 1971. [19] R.V. Field Jr., M. Grigoriu, J.M. Emery, On the efficacy of stochastic collocation, stochastic Galerkin, and stochastic reduced order models for solving stochastic problems, Probab. Eng. Mech. 41 (2015) 60–72. [20] B.L. Boyce, P.L. Reu, C.V. Robino, The constitutive behavior of laser welds in 304l stainless steel determined by digital image correlation., Metall. Mater. Trans. A 37A (2006) 2481–2492. [21] SIERRA Solid Mechanics Team, Adagio 4.18 User’s Guide, Technical report sand2010-6313, Sandia National Laboratories, Albuquerque, NM, 2010. [22] J.M. Emery, R.V. Field Jr., J.W. Foulk III, K.N. Karlson, M.D. Grigoriu, Predicting laser weld reliability with stochastic reduced-order models„ Int. J. Numer. Methods Eng. 103 (2015) 914–936. [23] N.L. Johnson, S. Kotz, N. Balakrishnan, Continuous Univariate Distributions, vol. 2, John Wiley & Sons, Inc., 1995.

Please cite this article as: J.M. Emery et al., Bayesian methods for characterizing unknown parameters of material models, Applied Mathematical Modelling (2016), http://dx.doi.org/10.1016/j.apm.2016.01.046