# Parametric Analysis of a Phenomenological Constitutive Model for Thermally
Induced Phase Transformation in Ni-Ti Shape Memory Alloys

P. Honarmandia,∗, A. Solomoub, R. Arroyavea, D. Lagoudasb a Department of
Materials Science and Engineering, Texas A&M University, College Station,
Texas, USA
b Department of Aerospace Engineering, Texas A&M University, College Station,
Texas, USA<EMAIL_ADDRESS>
###### Abstract

In this work, a thermo-mechanical model that predicts the actuation response
of shape memory alloys is probabilistically calibrated against three
experimental data sets simultaneously. Before calibration, a design of
experiments (DOE) has been performed in order to identify the parameters most
influential on the actuation response of the system and thus reduce the
dimensionality of the problem. Subsequently, uncertainty quantification (UQ)
of the influential parameters was carried out through Bayesian Markov Chain
Monte Carlo (MCMC). The assessed uncertainties in the model parameters were
then propagated to the transformation strain-temperature hysteresis curves
(the model output) using first an approximate approach based on the variance-
covariance matrix of the MCMC-calibrated model parameters and then an explicit
propagation of uncertainty through MCMC-based sampling. Results show good
agreement between model and experimental hysteresis loops after probabilistic
MCMC calibration such that the experimental data are situated within 95%
Bayesian confidence intervals. The application of the MCMC-based UQ/UP
approach in decision making for experimental design has also been shown by
comparing the information that can be gained from running replicas around a
single new experimental condition versus running experiments in different
regions of the experimental space.

  * November 2017

Keywords: Ni-Ti Shape Memory Alloys, Design of Experiment, Bayes’ Theorem,
Markov Chain Monte Carlo-Metropolis Hastings Algorithm, 95% Bayesian
Confidence Interval

##  1 Introduction

Nickle-Titanium ($Ni-Ti$) alloys are one of the most common shape memory
alloys (SMAs), which are widely used in different engineering applications,
such as biomedical devices and implants, microelectromechanical systems,
sensors and actuators, seismic protection tools, and aerospace products and
structures [1, 2]. Such applications are enabled by either the shape memory
effect or super-elasticity, which in turn are the result of a reversible
thermoelastic martensitic transformation that is triggered by applying thermal
and/or mechanical loads [3, 4].

The macroscopuc thermal actuation response of SMAs enables actuators with high
specific weight ratio compared to the conventional alternatives [5, 6].
Thermal actuation can occur during thermal cycles under isobaric condition
where inelastic strain is recovered through the forward and reverse
martensitic phase transformation induced during cooling and heating [7].
However, the forward and reverse transformations take place at different
temperatures due to energy dissipation during this cyclic thermal process [8],
which results in a hysteresis loop in strain-temperature space.

Under the framework of Integrated Computational Materials Engineering (ICME),
completing the loop along the process-structure-property-performance chain
through multiscale modeling is a necessary (albeit not sufficient) condition
for the efficient and robust design of materials for specific applications. In
the case of $Ni-Ti$ SMAs, different thermodynamics and kinetics modeling are
usually applied to connect process and structure, such as CALculation of PHAse
Diagram (CALPHAD) [9, 10], precipitation of secondary phases [11, 12], and
phase field models [13, 14, 15]. In regard to the structure-property linkages,
a number of constitutive models have been proposed over the recent decades to
predict the thermo-mechanical response of $Ni-Ti$ SMAs. These models can be
categorized into phenomenological [16, 17, 18, 19, 20, 21, 22, 23, 24, 25] or
micromechanics-based models [7, 26, 27, 28, 29] which have been reviewed
thoroughly by Patoor et al. [8] and Lagoudas et al. [30]. Generally, micro-
mechanical models are high-fidelity and expensive models that use
microstructural information to predict the macroscopic behavior of SMAs while
phenomenological models are cheap and fast since the free energy associated
with a non-microscopic homogenized material volume is taken into account for
the prediction of the macroscopic responses [21].

Despite a substantial number of constitutive models for the response of SMAs,
there are a very few studies associated with analysis, calibration, and
uncertainty quantification (UQ) of the model parameters and subsequent
uncertainty propagation (UP) from these parameters to the model responses.
These uncertainties can result from different sources, including natural
uncertainty (NU) from the random nature of the system, model parameter
uncertainty (MPU) due to lack of knowledge and data about the model
parameters, and Model structure uncertainty (MSU) as a result of incomplete or
simplified physics in the model [31]. The importance of UQ and UP, which tend
to be ignored in deterministic approaches towards model parameterization, is
noticeably highlighted in the case of robust design where it is sought that
the system response be relatively independent of the uncertainties of its
variables/parameters under the operating conditions. These inverse and forward
uncertainty analyses are discussed in more detail as follows.

Model parameter estimation based on available (experimental) data for the
system outputs is an inverse problem, which can be solved through
deterministic or probabilistic approaches [32]. Deterministic calibration such
as least squares approach yields a single best estimate based on the
minimization of the squared difference between the model results and data,
while probabilistic calibration provides probability distributions and
uncertainty bounds for the model parameters. In different engineering fields,
deterministic approaches have commonly been used to calibrate the model
parameters, even though it is already widely recognized that probabilistic
approaches can provide better predictions with the uncertainty/confidence
bounds [33].

In this regard, it is worth noting that many combinations of the model
parameters usually lead to similar model outputs, but, by definition,
deterministic calibration only arrives at one best estimate. Probabilistic
calibration, on the other hand, provides an ensemble of possible combinations
of model parameter values in the form of probability distributions [34].Among
probabilistic approaches, Bayesian-based Markov Chain Monte Carlo algorithms
have become the most commonly used tools due to the fast development in
computing capabilities [35]. These sampling techniques are more robust and
simple compared to the traditional analytic and numerical approaches to solve
high-dimensional intractable integrals in the probabilistic calibration
process [36, 37]. In addition, MCMC methods can be used to sample from complex
multivariate distributions that are usually very hard to sample using other
sampling techniques such as inversion or rejection sampling [37].

MCMC sampling, however, is an expensive method, particularly in high-
dimensional cases where the parameter convergence in the calibration process
can be very time-consuming. A reduction in the dimensionality of the parameter
space through sensitivity analysis (SA), however, can significantly improve
the MCMC calibration efficiency. Sensitivity analysis (SA) is often applied to
determine the impact of the parameters’ variability/uncertainty on the total
variability/uncertainty of the model results [38] in order to find the
insensitive model parameters that can be eliminated from the calibration
process. In this regard, design of experiments (DOE) is the most common
approach, which requires a description of experiments (such as complete
factorial design (CFD) offered by Fisher [39] or fractional factorial design
(FFD) to consider possible combinations of pre-defined parameters’ levels)
besides a subsequent analysis of the experiments (such as ANalysis Of
VAriance, ANOVA).

In robust design, the final results of the model and their corresponding
uncertainties are the main interest. Therefore, a UP technique should be
applied by considering the trade-off between cost and information loss to
propagate uncertainties from the model parameters to the model outputs [40].
This forward problem can approximately be solved through different methods
with different cost and precision, e.g., forward model analysis of the
parameters’ optimal combinations, statistical second moment methods,
polynomial chaos expansion, etc.

In the present work, the parameter sensitivity analysis and MCMC calibration
has been performed for the phenomenological constitutive model proposed by
Lagoudas et al. [21]. This model has been selected due to the consideration of
three important features in SMAs which are usually ignored in this type of
modeling: 1) the smooth transitions in thermo-mechanical responses of SMAs
during forward and reverse phase transformations that take place in a range of
temperatures and/or mechanical loadings, 2) the dependency of maximum
transformation strain after full transformation on the magnitude of the
applied stress although it does not hold at sufficiently high level of applied
stresses where the material system reaches a saturated value for the maximum
transformation strain resulting from the formation of all possible favored
martensitic variants, and 3) the definition of a critical driving force in
terms of the magnitude of the applied stress and the transformation direction
[21]. In section 2, this thermo-mechanical model is described thoroughly. In
section 3, parameter sensitivity analysis using the combination of CFD and
ANOVA are explained in details and subsequently utilized to perform the
sensitivity assessment of the model parameters. In section 4, the applied MCMC
sampling algorithm is described and used to probabilistically calibrate the
sensitive model parameters, and then the parameter uncertainties are
propagated to the model outcomes using different UP approaches. This is
followed by the summary and conclusion in section 5.

##  2 Description of the Constitutive Thermo-Mechanical Model

In this section, the phenomenological constitutive model proposed by Lagoudas
et al. [21] is described in detail. In this model, a total Gibbs free energy
is defined for the SMA material during the phase transformation in terms of
the thermos-elastic contributions of austenite ($G^{A}$), martensite ($G^{M}$)
and their interaction ($G^{mix}$) as follows:

$G(\sigma,T,\varepsilon^{t},\xi,g^{t})=(1-\xi)G^{A}(\sigma,T)+\xi
G^{M}(\sigma,T)+G^{mix}(\sigma,\varepsilon^{t},g^{t})$ (1)

where the external state variables $\sigma$ and $T$ are the applied stress
tensor and the absolute temperature, and the internal state variables
$\epsilon^{t}$, $\xi$, and $g^{t}$ account for the inelastic strain produced
during forward/reverse martensitic transformation, the total volume fraction
of martensite, and the transformation hardening energy, respectively.
Moreover, $G^{A}$, $G^{M}$, and $G^{mix}$ can be written as the following
equations:

$\displaystyle G^{A/M}$
$\displaystyle=-\frac{1}{2\rho}\sigma:S^{A/M}\sigma-\frac{1}{\rho}\sigma:\alpha(T-T_{0})$
(2)
$\displaystyle+c^{A/M}\Bigg{[}(T-T_{0})-T\ln{\Big{(}\frac{T}{T_{0}}\Big{)}}\Bigg{]}-s_{0}^{A/M}T+u_{0}^{A/M}$

$G^{mix}(\sigma,\varepsilon^{t},g^{t})=-\frac{1}{\rho}\sigma:\varepsilon^{t}+\frac{1}{\rho}g^{t}$
(3)

The constant parameters $\rho$ and $\alpha$ are the alloy density and the
thermal extension tensor, and the phase-dependent parameters $S$, $c$,
$s_{0}$, and $u_{0}$ corresponds to the compliance tensor, specific heat,
specific entropy, and specific internal energy, respectively. It should be
noted that the notation $(.:.)$ shows the inner product of two second-order
tensors.

Before applying the thermodynamics laws, an evolution relationship is defined
for the transformation strain in terms of the volume fraction of martensite
during forward and reverse transformations [21]:

$\dot{\varepsilon}^{t}=\Lambda^{t}\dot{\xi}\ \ ,\ \
\Lambda^{t}=\Bigg{\\{}\begin{array}[]{l}\Lambda_{fwd}^{t}\ \ ,\ \
\dot{\xi}>0\\\ \Lambda_{rev}^{t}\ \ ,\ \ \dot{\xi}<0\end{array}$ (4)

where $\Lambda^{t}$ represents the transformation direction tensor, which is
expressed for both directions as:

$\Lambda_{fwd}^{t}=\frac{3}{2}H^{cur}\Big{(}\frac{\sigma^{{}^{\prime}}}{\bar{\sigma}}\Big{)}\
\ ,\ \ \Lambda_{rev}^{t}=\frac{\epsilon^{t}}{\xi}$ (5)

$\bar{\sigma}=\sqrt{(3/2)\sigma^{{}^{\prime}}:\sigma^{{}^{\prime}}}$ is the
effective stress where $\sigma^{{}^{\prime}}$ is the deviatoric stress, and
$H^{cur}$ is the maximum transformation strain after full matensitic
transformation as a function of the effective stress:

$H^{cur}(\bar{\sigma})=H_{sat}(1-e^{-k\bar{\sigma}})$ (6)

where $H_{sat}$ is a saturated value for the maximum transformation strain at
sufficiently high value of the effective stress, and the parameter $k$
determines the exponential rate of $H^{cur}$ variation from 0 to $H_{sat}$ [7,
21].

Another evolution model can also be written between the variation rate of the
transformation hardening energy and the volume fraction of martensite during
forward and reverse transformation, as follows:

$\dot{g}^{t}=f^{t}\dot{\xi}\ \ ,\ \
f^{t}=\Bigg{\\{}\begin{array}[]{l}f_{fwd}^{t}\ \ ,\ \ \dot{\xi}>0\\\
f_{rev}^{t}\ \ ,\ \ \dot{\xi}<0\end{array}$ (7)

$f^{t}$ is a direction-dependent hardening function, which has been proposed
in this model such that it can capture the smooth transitions from the elastic
to transformation regime and vice versa:

$\Bigg{\\{}\begin{array}[]{l}f_{fwd}^{t}(\xi)=\frac{1}{2}a_{1}(1+\xi^{n_{1}}-(1-\xi)^{n_{2}})+a_{3}\\\
f_{rev}^{t}(\xi)=\frac{1}{2}a_{2}(1+\xi^{n_{3}}-(1-\xi)^{n_{4}})-a_{3}\end{array}$
(8)

where $n_{i}$ can change in the range of (0, 1]. Closer value of $n_{i}$ to
zero corresponds to more smooth transition at the beginning or the end of the
forward phase transformation (in the case of $i=1$ or $2$, i.e., $n_{1}$ or
$n_{2}$) or their counterparts for the reverse phase transformation (in the
case of $i=3$ or $4$, i.e., $n_{3}$ or $n_{4}$).

Now, the first and second law of thermodynamics can be applied. At a local
point in the material, the conservation of energy or first law of
thermodynamics can be expressed as:

$\rho\dot{u}=\sigma:\dot{\varepsilon}-\textnormal{div}(q)+\rho r$ (9)

where $\dot{u}$, $q$, and $r$ denote the internal energy rate, the heat flux
vector, and the internal heat generation rate, respectively. The second law of
thermodynamics at this local point can also be written as the Clausius-Planck
inequality [41]:

$\rho\dot{s}+\frac{1}{T}\textnormal{div}(q)-\frac{\rho r}{T}\geq 0$ (10)

It can turn into the following inequality through multiplying two sides of the
inequality by T and substituting div$(q)$ by its equivalent obtained from
equation 9:

$\rho\dot{s}T+\sigma:\dot{\varepsilon}-\rho\dot{u}\geq 0$ (11)

It should be noted that Gibbs free energy can be defined in terms of internal
energy using the Legendre transformation, as follows:

$G=u-\frac{1}{\rho}\sigma:\varepsilon-sT$ (12)

If $\dot{u}$ in equation 11 is substituted by the first derivative of obtained
u in Equation 12 with respect to time, the following thermodynamic constraint
is resulted for the Gibbs free energy rate:

$-\rho\dot{G}-\dot{\sigma}:\varepsilon-\rho s\dot{T}\geq 0$ (13)

Using the chain rule, the second law of thermodynamics can turn into:

$-\rho\Big{(}\partial_{\sigma}{G}:\dot{\sigma}+\partial_{T}{G}:\dot{T}+\partial_{\varepsilon^{t}}{G}:\dot{\varepsilon}^{t}+\partial_{\xi}{G}:\dot{\xi}+\partial_{g^{t}}{G}:\dot{g}^{t}\Big{)}-\dot{\sigma}:\varepsilon-\rho
s\dot{T}\geq 0$ (14)

Here, $-\rho\partial_{\xi}{G}=p$, $-\rho\partial_{\varepsilon^{t}}{G}=\sigma$,
and $-\rho\partial_{g^{t}}{G}=-1$ are three generalized thermodynamic forces.
In addition, the total infinitesimal strain and entropy can be expressed as
follows based on Coleman and Noll approach [42]:

$\varepsilon=-\rho\partial_{\sigma}{G}=S\sigma+\alpha(T-T_{0})+\varepsilon^{t}$
(15)

$s=-\partial_{T}{G}=\frac{1}{\rho}\alpha:\sigma+c\ln{\Big{(}\frac{T}{T_{0}}\Big{)}+s_{0}}$
(16)

Therefore, the final inequality for the second law of thermodynamics can be
expressed as:

$p\dot{\xi}+\sigma:\dot{\varepsilon}^{t}-\dot{g}^{t}\geq 0$ (17)

where the generalized thermodynamic force p can be obtained by taking
derivative of total $G$ in equation 1 with respect to $\xi$:

$\displaystyle p$
$\displaystyle=-\rho\partial_{\xi}{G}=G^{M}(\sigma,T)-G^{A}(\sigma,T)$ (18)
$\displaystyle=\frac{1}{2}\sigma:\Delta
S\sigma+\sigma:\Delta\alpha(T-T_{0})-\rho\Delta
c\Bigg{[}(T-T_{0})-T\ln{\Big{(}\frac{T}{T_{0}}\Big{)}}\Bigg{]}$
$\displaystyle+\rho\Delta s_{0}T-\rho\Delta u_{0}$

As observed in this equation, $p$ is proportional to $ΔG=G^{M}-G^{A}$ that is
the difference between the Gibbs free energy of pure martensite and pure
austenite phase.

Inserting equations 4 and 7 into equation 17 results in:

$(\sigma:\Lambda^{t}+p-f^{t})\dot{\xi}=\pi^{t}\dot{\xi}\geq 0$ (19)

where $\pi^{t}$ is the total thermodynamic force which must be positive during
the forward transformation (when $\dot{\xi}>0$) and negative during the
reverse transformation to hold the second law of thermodynamics (when
$\dot{\xi}<0$). In this model, this thermodynamic driving force must reach a
stress-dependent critical thermodynamic driving force ($Y^{t}$) in order to
start and continue the phase transformation. Different values of $Y^{t}$ for
forward and reverse phase transformation lead to a hyper-surface for each
transformation direction, as follows:

$\Phi^{t}(\sigma,T,\xi)=0\ \ ,\ \
\Phi^{t}(\sigma,T,\xi)=\Bigg{\\{}\begin{array}[]{l}\Phi_{fwd}^{t}=\pi_{fwd}^{t}-Y_{fwd}^{t}\
\ ,\ \ \dot{\xi}>0\\\ \Phi_{rev}^{t}=-\pi_{rev}^{t}-Y_{rev}^{t}\ \ ,\ \
\dot{\xi}<0\end{array}$ (20)

where,

$Y^{t}(\sigma)=\Bigg{\\{}\begin{array}[]{l}Y_{fwd}^{t}=Y_{0}^{t}+D\sigma:\Lambda_{fwd}^{t}\
\ ,\ \ \dot{\xi}>0\\\ Y_{rev}^{t}=Y_{0}^{t}+D\sigma:\Lambda_{rev}^{t}\ \ ,\ \
\dot{\xi}<0\end{array}$ (21)

$Y_{0}^{t}$ is a constant, and the model parameter $D$ indicates the
dependency of the critical value for driving force on the applied stress.

In the case of iso-baric condition ($\sigma$=cons), the value of $\xi$ can be
obtained at any specific $T$ during forward/reverse phase transformation using
equation 20. Since $\varepsilon^{t}$ is directly related to $\xi$, a
hysteresis loop can be obtained in transformation strain-temperature diagram
if $T$ alters from martensite start temperature ($M_{s}$) to martensite finish
temperature ($M_{f}$) and from austenite start temperature ($A_{s}$) to
austenite finish temperature ($A_{f}$) for forward and reverse phase
transformation, respectively.

It is worth noting that the model parameters can be correlated to the material
properties for the sake of the model calibration. In this regard, the
material’s properties $E^{A/M}$ and $\nu^{A/M}$ can represent the model
parameters $S^{A/M}$ and $\alpha$, and $M_{s}$, $M_{f}$, $A_{s}$, $A_{f}$,
$C^{A/M}$ can be associated with the model parameters $\rho\Delta s_{0}$,
$\rho\Delta u_{0}$, $a_{1}$, $a_{2}$, $a_{3}$, $Y_{0}^{t}$, and $D$, where
$E^{A/M}$ and $\nu^{A/M}$ are Young’s modulus and Poisson’s ratios of
austenite/martensite phase, and $C^{A/M}$ is the lines slope of
austenitic/martensitic transformation temperatures ($A_{s/f}/M_{s/f}$) in
stress-temperature phase diagram [7, 21].

##  3 Sensitivity Analysis of the Model Parameters

Before model calibration, SA is usually required to reduce the dimensionality
of the parameter space with minimum possible loss of information about the
model outcome in order to save calibration time and cost as a result of faster
convergence of the model parameters. In this work we carried a DOE based on a
complete factorial design (CFD) coupled to ANOVA in order to identify the
parameters most highly correlated to the output of the model.

###  3.1 Complete Factorial Design

CFD is an experimental design method that considers all possible combinations
of the levels ($R$) defined for factors/variables/parameters ($N$) in a given
system/model, which is equivalent to $R^{N}$ level-factor combinations. In the
case of parameter experimental design, model response should be obtained for
each level-parameter combination for the sake of sensitivity assessment.

Based on the expert’s intuition, 14 parameters have been considered as
candidates for model calibration. Two plausible values have been selected for
each parameter as lower and upper levels in the context of two-level CFD.
These levels are considered as the applied initial values $pm10\%$ of the
range of the parameters in order to hold the following temperature constraints
for their design combinations:

$M_{f}<M_{s}<A_{s}<A_{f}$ (22)

Therefore, $2^{14}=16384$ level-parameter combinations can be constructed
during CFD whose responses are determined through the difference between the
hysteresis loop obtained from running the model with each parameter
combination and a reference hysteresis loop obtained using the the average
value of lower and upper levels of the parameters. It is worth noting that any
reference curve can be applied in so far as it is consistent for all the
responses. Tschopp et al. [43] has evaluated various approaches to determine
the similarity/difference of two images by comparing two vectors which
represent the image features. In our case, these vectors were constructed
based on the coordinate of spatial points on the resulting hysteresis loops:

$d_{SE}=\displaystyle\sum_{i}(X_{i}-Y_{i})^{2}$ (23)

Now, the CFD results can be applied for SA of the given model parameters
through ANOVA.

###  3.2 Analysis of Variance

A 14-way ANOVA has been performed in this work using ”anovan” function in
Matlab Machine Learning toolbox, where N is the number of independent Model
Parameters–see Appendix A.

ANOVA is a statistical hypothesis testing technique based on the deviation
associated with each individual level of each factor/parameter/variable or
each level combination of each interaction between
factors/parameters/variables in the system from the overall mean of responses.
The resulting response for each one of $2^{14}$ factorial design combinations
has been used in the 14-way ANOVA to evaluate the sensitivity of 14 model
parameters introduced in section 2. The ANOVA results associated with each
parameter are shown in table 1. In this table, the model parameters have been
ranked based on their corresponding p-values in descending order, which
indicate the parameters from the highest to the lowest sensitivity. According
to the significance level 0.05, the first eight parameters have been selected
for further uncertainty analysis in order to reduce the computational cost for
the model calibration that will be discussed in section 4.1.

Table 1:  ANOVA table results showing the model parameters ranked based on
their p-values in descending order.

  * Source | Sum sq. | d.f. | Mean sq. | F | Prob$>$F
---|---|---|---|---|---
$H_{sat}$ | 0.65941 | 1 | 0.65941 | 800.62 | 5.3041e-172
$A_{f}(K)$ | 0.33278 | 1 | 0.33278 | 404.05 | 8.5195e-89
$M_{s}(K)$ | 0.08907 | 1 | 0.08907 | 108.14 | 3.0049e-25
$M_{f}(K)$ | 0.04402 | 1 | 0.04402 | 53.45 | 2.7735e-13
$C^{A}(MPa/K)$ | 0.03488 | 1 | 0.03488 | 42.35 | 7.8628e-11
$k(MPa^{-1})$ | 0.02462 | 1 | 0.02462 | 29.90 | 4.6223e-08
$E^{M}(GPa)$ | 0.01029 | 1 | 0.01029 | 12.50 | 4.0821e-04
$A_{s}(K)$ | 0.00323 | 1 | 0.00323 | 3.93 | 0.0476
$E^{A}(GPa)$ | 0.00018 | 1 | 0.00018 | 0.21 | 0.6443
$C^{M}(MPa/K)$ | 0.00008 | 1 | 0.00008 | 0.09 | 0.7612
$n_{1}$ | 0 | 1 | 0 | 0 | 1
$n_{2}$ | 0 | 1 | 0 | 0 | 1
$n_{3}$ | 0 | 1 | 0 | 0 | 1
$n_{4}$ | 0 | 1 | 0 | 0 | 1
Error | 13.4819 | 16369 | 0.00082 |  |
Total | 14.6805 | 16383 |  |  |

##  4 Bayesian Uncertainty Quantification and Propagation of the Model
Parameters

###  4.1 Parameter Estimation and Uncertainty Quantification using Markov
Chain Monte Carlo-Metropolis Hastings Algorithm

In order to sample the parameter space, one can use Bayesian-based MCMC
sampling approaches [37]. Among MCMC approaches, Metropolis-Hastings (M-H)
sampling is one of the most common applied tools to probabilistically solve
the inverse problem associated to probabilistic parameter calibration from
experimental data. In the case of model calibration, this technique updates
the existing prior knowledge about the parameters to their posterior
information based on the available data from the model responses.

In the present work, MCMC toolbox in Matlab has been used to probabilistically
calibrate the eight sensitive model parameters identified in section 3.2.
First, the prior parameter information, including the initial values
($\theta^{0}$), valid ranges, and probability density functions (PDFs) for the
parameter are considered. Then, a new parameter vector/candidate
$\theta^{cand}$ is randomly sampled from a proposal posterior PDF ($q$) which
is defined as a multivariate Gaussian distribution in the toolbox with a mean
value at $\theta^{0}$ and an arbitrary variance-covariance matrix. The
candidate is accepted/rejected using a probabilistic criterion defined based
on the M-H ratio that is expressed as follows:

$MH=\frac{p(\theta^{cand})p(D|\theta^{cand})}{p(\theta^{0})p(D|\theta^{0})}\frac{q(\theta^{0}|\theta^{cand})}{q(\theta^{cand}|\theta^{0})}$
(24)

In this equation, the first ratio is called Metropolis ratio which is the
ratio of the posterior probability of $\theta^{cand}$ to $\theta^{0}$ given
data ($D$). In the context of the Bayes’ theorem, the posterior probabilities
are proportional to the prior probability times the likelihood—probability of
the data given model results evaluated with the candidate model parameters—for
each case. We considered the likelihood to be a Gaussian distribution centered
at the data whose error is considered as the distribution variance. In the
case of multiple data sets, it should be noted that the product of the
likelihoods obtained from each data set can yield the total likelihood if data
sets are assumed to be conditionally independent.

The Metropolis ratio can be applied for the acceptance/rejection of the
candidate when the proposal distribution is symmetric; however, this is not
usually the case since the probability of moving from $\theta^{0}$ to
$\theta^{cand}$ is not necessarily equal to the probability of the reverse
move. For this reason, Hastings introduced the second ratio in equation 24 to
involve the above-mentioned non-symmetric property of the proposal
distribution [37, 44]. Now, $min(MH\times 100,100)$ is considered as the
acceptance probability of the candidate [44]. According to this criterion, If
the candidate is accepted, then $\theta^{1}=\theta^{cand}$; otherwise,
$\theta^{1}=\theta^{0}$. The sampling of the new candidate and its
acceptance/rejection process are sequentially repeated $n$ times by finding
the M-H ratio using the new sampled candidate and the previous accepted
parameter during each iteration. At the end of the procedure, N samples of the
parameter vector ($\big{\\{}\theta^{0},...,\theta^{N}\big{\\}}$) are generated
to represent a multivariate posterior distribution for the model parameters.

Before the convergence of the parameter samples, samples generated during the
”burn-in period”—in which the different samples are causally correlated—must
be removed before further analysis of the parameters. After the removal of the
burn-in period, the mean values and the square root of diagonal elements in
variance-covariance matrix of the remaining samples in the convergence region
can be introduced as the plausible optimal values and uncertainties of the
model parameters.

In our calibration work, initial values and their ranges have been determined
based on the expert’s beliefs, and an uniform (non-informative) probability
distribution has been selected for each model parameter due to the absence of
knowledge about the parameters’ PDF. In addition, three experimental datasets
for transformation strain-temperature hysteresis loop under different isobaric
conditions have simultaneously been utilized together to calibrate the model.
This has been performed through the calculation of the difference between
various isobaric hysteresis curves obtained from each parameter sample and the
corresponding experimental data using equation 23. For the sake of the
parameter calibration, these differences must simultaneously be compared with
zero in the likelihood function where $M(\theta)$ and $D$ are the mentioned
differences and zero, respectively.

In this work, the likelihood variance $\sigma^{2}$ (the squared error of the
distances between model and experimental hysteresis curves from zero) is set
as an unknown hyper-parameter in the calibration process and updated in
parallel with the other model parameters during MCMC sampling. The only
difference is that an inverse-gamma distribution is defined as the non-
informative prior PDF for this hyper-parameter, which yields the same form of
distribution for the posterior PDF after applying the Bayes’ theorem since the
inverse-gamma distribution is the corresponding conjugate prior density for
exponential likelihood functions. Therefore, $\sigma^{2}$ is sampled from an
inverse gamma posterior distribution during MCMC process [45]. In addition, it
should be noted that an adaptive proposal distribution has been employed in
this calibration work. In this scheme, an arbitrary positive covariance
($V_{0}$) is selected in the beginning of MCMC sampling, and then it is
adapted in each next iteration using the covariance obtained from the ensemble
of previous samples in the MCMC chain [46, 47].

In this work, a constrained MCMC approach has also been developed in order to
consider the constraints for the transformation temperatures mentioned in
equation 22. For this purpose, the likelihood is penalized by assigning an
infinitely high value for the distance between model and experimental
hysteresis curves in the cases that the constraints for the parameters are not
satisfied. In these cases, the induced penalty results in a likelihood close
to zero, which in turn results in an infinitesimally small acceptance
probability for new model parameter candidates that violate constraints.

In the calibration of the sensitive model parameters, high number of samples
(i.e., 200,000) are generated to ensure the convergence of the parameters and
the uniqueness of the optimal response in the parameter space. The multi-
variate posterior distribution after MCMC sampling can be assessed through
joint and marginal frequency distributions. In figure 1, some of the parameter
joint distributions have been shown in 2D color graphs as examples, where the
colors represents different density of samples in the parameter space,
increasing from blue to red spectra. These types of graphs can also
demonstrate the correlation between each pair of parameters qualitatively. For
instance, negative/positive linear correlations can be inferred from the
negative/positive slopes of the elliptical shapes in the graphs. The linearity
of the sample distribution in the parameter space can be determined
quantitatively through Pearson correlation coefficient ($-1<\rho<1$) that is
expressed as:

$\rho_{X,Y}=Corr(X,Y)=\frac{cov(X,Y)}{\sigma_{X}\sigma_{Y}}$ (25)

where $\sigma_{X}$ and $\sigma_{Y}$ are the standard deviations for parameters
$X$ and $Y$, respectively, which equal the square root of the corresponding
diagonal elements in $8\times 8$ variance-covariance matrix. $cov(X,Y)$ is
also the covariance of $X$ and $Y$ obtained from the corresponding non-
diagonal element in the variance-covariance matrix. Generally, a closer value
of $\rho$ to either -1 or 1 implies a higher linear correlation between $X$
and $Y$, whereas a closer value to zero suggests their lower linear
correlation. The negative or positive sign refers to the negative or positive
linear correlation. The linear coefficients for all pairs of parameters in the
thermo-mechanical model have been listed in table 2. As can be observed in
this table, most of the pair parameters are linearly uncorrelated although 6
out of 28 parameter pairs show some weak or moderate correlations with
$|\rho|$ greater than 0.2. In addition, the parameters $C^{A}$ and $A_{f}$
(1-d) and the parameters $E^{M}$ and $C^{A}$ (1-f) show the highest and the
lowest linear correlation, respectively, regardless of the correlation signs.

Figure 1: Joint frequency distributions for some of the pair model parameters
in the form of 2D color graphs with normalized color bars. Table 2:  Linear
correlation coefficient for each pair of model parameters.

$A_{s}(K)$ $A_{f}(K)$ $M_{s}(K)$ $M_{f}(K)$ $C^{A}(MPa/K)$ $E^{M}(GPa)$
$H_{sat}$ $k(MPa^{-1})$ $A_{s}(K)$ 1 -0.0721 0.2948 -0.1080 0.3813 -0.0167
-0.1330 -0.0378 $A_{f}(K)$ -0.0721 1 -0.0150 -0.0822 0.5836 -0.0553 0.0249
-0.0672 $M_{s}(K)$ 0.2948 -0.0150 1 -0.4628 0.1156 -0.0547 -0.0465 -0.0415
$M_{f}(K)$ -0.1080 -0.0822 -0.4628 1 -0.1280 -0.0289 -0.1355 0.0160
$C^{A}(MPa/K)$ 0.3813 0.5836 0.1156 -0.1280 1 -0.0015 0.0476 -0.0405
$E^{M}(GPa)$ -0.0167 -0.0553 -0.0547 -0.0289 -0.0015 1 0.2771 0.0091 $H_{sat}$
-0.1330 0.0249 -0.0465 -0.1355 0.0476 0.2771 1 -0.2103 $k(MPa^{-1})$ -0.0378
-0.0672 -0.0415 0.0160 -0.0405 0.0091 -0.2103 1

As mentioned above, the burn-in period at the beginning of MCMC sampling must
be removed before the calibration and UQ of the model parameters, which can be
found using the cumulative mean distributions for the parameters, as shown in
figure 2. In these figures, there are some noisy trends at the beginning of
MCMC sampling for the values of all the parameters, but they reach almost
plateaus after about 63,000 generations that show the convergence values (mean
values) for the parameters. Therefore, the first 63,000 samples are extracted
from the total MCMC sample collection before further analysis. The mean and
the square root of the diagonal elements in the variance-covariance matrix of
the remaining samples in the convergence region can be introduced as the
plausible optimal values and the uncertainty of the model parameters, which
are listed in table 3.

Table 3:  Plausible optimal values and uncertainties (standard deviations) of
the model parameters after MCMC calibration besides their initial values.

$A_{s}(K)$ $A_{f}(K)$ $M_{s}(K)$ $M_{f}(K)$ $C^{A}(MPa/K)$ $E^{M}(GPa)$
$H_{sat}$ $k(MPa^{-1})$ Initial Values 307.0 318.0 300.0 270.0 9.0 40.0 0.034
0.02 After MCMC Calibration 296.6$\pm$8.5 322.6$\pm$11.3 280.4$\pm$6.6
259.9$\pm$8.5 11.8$\pm$4.1 35.6$\pm$8.6 0.0517$\pm$0.0044 0.0595$\pm$0.0237

Figure 2: Cumulative mean distribution plot for each model parameter after
MCMC sampling.

After removal of the burn-in period, the marginal posterior frequency
distribution of the parameters can be plotted using the remaining MCMC samples
as observed in figure 3. These plots suggest skewed Gaussian distributions for
the marginal posterior frequency distribution (although it seems to be a
Gaussian distribution in the case of the parameter $H_{sat}$), which are
truncated from both sides by the lower and upper bounds of the parameters.
However, the parameters $E^{M}$ and $k$ show flatter distributions due to
their low sensitivities among the selected parameters, as shown in table 1. In
general, model parameters with lower sensitivities can suggest broader sample
distributions since their variations make less changes in model outcomes.

Figure 3: marginal posterior frequency distribution plot for each model
parameter after the removal of the burn-in period from MCMC samples.

###  4.2 Uncertainty Propagation from the Model Parameters to the Model
Outcomes

As mentioned earlier, only the final results of the model and the
corresponding uncertainties are considered important in the context of the
robust design. The reason is that such information can help the designers to
reduce the influenceof the system performance by the existing uncertainties in
the outcomes, regardless of how many variables and models have been involved.
Therefore, an appropriate UP technique is essential to propagate uncertainties
from the model parameters/input variables to the outputs with a relatively
high efficiency and a minimum or at least quantifiable loss of information
[40]. Among different UP methods, first/second order second moment (FOSM/SOSM)
approaches are probably the most common approximation procedures in
engineering [48]. In this work, FOSM approach has been used to perform non-
linear propagation of the variance-covariance matrix of the model parameters.

Let $\theta$ and $\bar{\theta}$ vectors denote the optimal MCMC samples for
the parameters (after removal of the burn-in period) and their mean values,
respectively, and $M(\theta)$ denote the model objective (output) function. If
the objective function is estimated by a first order Taylor expansion at
$\bar{\theta}$:

$M(\theta)\approx M(\bar{\theta})+\sum_{i=1}^{N}{\frac{\partial
M(\theta)}{\partial{\theta_{i}}}(\theta_{i}-\bar{\theta}_{i})}$ (26)

The expected values associated with first and second moment of the above first
order Taylor expansion approximates the mean value ($\bar{M}$) and variance
($\sigma_{M}^{2}$) of the objective function $M$ [48] as follows:

$\bar{M}=E[M(\theta)]\approx M(\bar{\theta})$ (27)

$\sigma_{M}^{2}=E[(M(\theta)-\bar{M})^{2}]\approx\sum_{i=1}^{N}{{|\frac{\partial
M(\theta)}{\partial{\theta_{i}}}|}^{2}\sigma_{i}^{2}}+\sum_{i=1}^{N}{\sum_{j\neq
i}^{N}\frac{\partial M(\theta)}{\partial{\theta_{i}}}\frac{\partial
M(\theta)}{\partial{\theta_{j}}}\sigma_{ij}}$ (28)

where $\sigma_{i}^{2}$ and $\sigma_{ij}$ are the variance of the parameter
$\theta_{i}$ and the covariance between the parameters $\theta_{i}$ and
$\theta_{j}$, respectively. The approximations in equation 27 and 28 has been
derived thoroughly by Kriegesmann [49]. It is worth noting that equation 28
can also be expressed in terms of the variance-covariance matrix of the model
parameters ($V$) [50] as follows:

$\sigma_{M}^{2}\approx g^{T}Vg$ (29)

where g is a column vector which consists of all partial derivative elements,
i.e., $g=(\frac{\partial M(\theta)}{\partial{\theta_{1}}},...,\frac{\partial
M(\theta)}{\partial{\theta_{N}}})^{T}$.

The uncertainties of parameters reported in table 3 have been propagated to
the uncertainty of the transformation strain along the hysteresis curves by
applying the variance-covariance matrix obtained from the optimal MCMC samples
(without burn-in period) in the above-mentioned FOSM approach. The results
have been shown in figure 4 for different isobaric conditions. In these
figures, the blue/red solid and dashed lines correspond to the model results
during cooling/heating process (forward/reverse martensitic transformation)
and their experimental counterparts, respectively. In addition, the blue/red
shaded regions indicate 95% Bayesian confidence interval (95% BCI) for model
results during cooling/heating. It should be noted that the interval is
equivalent to $\varepsilon^{t}(T)\pm 2\sigma_{\varepsilon^{t}}$, where
$\varepsilon^{t}(T)$ and $\sigma_{\varepsilon^{t}}$ can be calculated through
equations 27 and 28, respectively.

As can be observed, there are generally good agreements between the model and
experimental hysteresis loops, or at least the experimental data are located
inside the 95% BCIs. Although the given thermo-mechanical model cannot exactly
predict the slope of the curves during forward and reverse transformation,
there are no discrepancies between model results and experimental data for the
maximum transformation strain ($H^{cur}$) in each isobaric condition. The
difference in the curves’ slopes can be attributed to the slight effects of
other parameters that are not considered in the calibration, the missing
physics in the model, and the uncertainties in the experimental data all
together.

Another important feature in theses graphs are the unrealistic humps around
the transformation temperatures, which can be related to the deficiency of the
partial derivatives in the applied UP approach (28) at where the change
suddenly. For this reason, a direct UP has been performed using the model
forward analysis of optimal MCMC samples to find more rational regions for 95%
BCIs. After running the model for all the optimal samples, 2.5% of the
resulting hysteresis curves have been eliminated from each one of the top and
down margins to obtain the mentioned intervals. As shown in figure 5, this
approach yields smoother and more precise boundaries for 95% BCIs with no
humps. Although the direct UP approach is more expensive than the FOSM method,
it is required to be used in this case to provide more realistic and precise
results for 95% BCIs which are very important in robust design.

Figure 4: 95% BCIs for the model hysteresis curves obtained from the FOSM
approach for different isobaric conditions, a) 100, b) 150, and c) 200 MPa,
besides their corresponding experimental counterparts.

Figure 5: 95% BCIs for the model hysteresis curves obtained from the direct UP
approach for different isobaric conditions, a) 100, b) 150, and c) 200 MPa,
besides their corresponding experimental counterparts.

##  5 MCMC application in decision making for experimental design

In this section, The goal is to respond the experimentalists’ question about
which experimental design can gain more knowledge about the system?—performing
experiments at the same conditions or at different conditions in design input
space. For this purpose, it is assumed that the probabilistically calibrated
model in section 4.1 can yield the ground truth and its natural uncertainty
bounds at any design condition (i.e., any isobaric condition). According to
this assumption, two different synthetic experimental sets are designed
through sampling the hysteresis curves from the model uncertainty bounds
obtained at any isobaric condition of interest. The first set contains three
random samples (replicas) from the same experimental condition (i.e., $\sigma$
= 150 MPa), while the second set considers three random samples from three
different conditions (i.e., $\sigma$ = 175, 250, and 300 MPa). After the
generation of synthetic data sets, the calibrated results obtained in Table 3
for the influential parameters are considered as the parameter priors which
are separately updated against each set of synthetic data using the MCMC
technique in a sequential way of data training. In this way of training, the
parameter posterior distribution obtained after each training is considered as
the prior distribution for the next training. The relative entropy
(Kullback–Leibler (K-L) divergence) is calculated to find how diverged the
probability distributions are relatively, i.e., what is the distance between
the prior and final posterior probability distributions of the parameters
after the three sequential MCMC calibrations in each case. The calculated K-L
divergence can be introduced as a comparison measure for the amount of
information gained from each synthetic experimental sets. Generally, K-L
divergence for continuous probability distributions P and Q is defined as:

$D_{KL}(P||Q)=\int_{-\infty}^{+\infty}p(x)log(\frac{p(x)}{q(x)})dx$ (30)

where $p(x)$ and $q(x)$ are the densities of $P$ and $Q$ at any given random
point $x$, respectively. In the case of multivariate Gaussian distribution,
the integration in Equation 30 can be obtained as follows:

$D_{KL}(\mathcal{N}_{1}||\mathcal{N}_{2})=\frac{1}{2}\Big{[}ln(\frac{|\Sigma_{2}|}{|\Sigma_{1}|})-d+tr(\Sigma_{2}^{-1}\Sigma_{1})+(\mu_{2}-\mu_{1})^{T}\Sigma_{2}^{-1}(\mu_{2}-\mu_{1})\Big{]}$
(31)

The equivalent normal distributions obtained for parameter prior and posterior
distributions are compared for each set of experimental design using Equation
31. The K-L divergence values for experimental sets 1 and 2 are calculated
around 3.8 and 4.4, respectively. The higher value of the K-L divergence for
experimental set 2 indicates a greater difference between prior and final
posterior distributions in this case, which implies that the consideration of
different conditions in experimental design can gain more information about
the system rather than experimental replicas in the same condition.

##  6 Summary and Conclusion

In this work, the calibration and UQ of the sensitive parameters in a thermo-
mechanical model have been performed against three different isobaric
experimental data simultaneously using a constrained MCMC-MH algorithm in the
context of Bayesian statistics. Eight sensitive parameters have been found
before the calibration using a DOE approach which includes CFD and N-way
ANOVA. After MCMC sampling of these eight parameters, their convergence have
been checked by the joint frequency distributions and the cumulative mean
plots, which have also been utilized to identify the qualitative correlation
of each pair parameters and the burn-in period, respectively. In addition, the
linear correlation between each two parameters has been calculated through
Pearson coefficient. These coefficients suggest that most pair parameters are
linearly uncorrelated, although a few of them show weak or moderate linear
correlations.

After removal of the samples belonging to the burn-in period, skewed Gaussian
distributions have been obtained for the marginal posterior frequency
distributions of the model parameters, which are truncated in the range
defined for each parameter. Among the selected parameters, the parameters with
the lowest sensitivities show flatter distributions due to their less
influences on the model outcome. Moreover, it should be noted that the mean
and square root of the diagonal elements in the variance-covariance matrix of
the remaining optimal samples after the elimination of the burn-in period
provide the plausible optimal values and uncertainties of the selected
parameters, respectively.

At the end, the parameters’ uncertainties have been propagated to the
uncertainties of the transformation strain along the hysteresis curves using
two various UP techniques in order to find 95% BCIs during cooling/heating
process. In this regard, FOSM approach results in some unrealistic humps
around the transformation temperatures. For this reason, a direct UP technique
has been performed through the forward analysis of optimal MCMC samples in
order to achieve more rational and precise uncertainty bands that are very
important in the context of robust design. The hysteresis curves obtained from
the mean values of the parameters are in good agreements with their
experimental counterparts for different given isobaric conditions, or at least
it can be stated that the experimental data are situated in 95% BCIs. Although
the slopes of the forward and reverse transformation curves are not predicted
exactly, the maximum transformation strain associated with each isobaric
condition is very close to its corresponding data.

In this work, It has been also shown that the Markov Chain Monte Carlo
approach can be applied for decision making in experimental design. Assuming
the calibrated model can predict the ground truth and the natural
uncertainties in different experimental conditions, two synthetic experimental
sets have been sampled from the model response uncertainty bounds obtained at
the same and different isobaric conditions, respectively. After sequential
MCMC updates of parameters’ posterior distributions with each set, their
information gains are compared together through the calculation of K-L
divergence for each case which indicate the difference between prior and final
posterior distributions of the model parameters. A higher information gain has
been obtained from the experimental set 2 that are sampled from different
experimental conditions (different isobaric conditions).

Generally, this work introduces a complete framework for UQ of model
parameters and subsequent UP from model parameters to outputs as a guideline
for material design.

## References

## References

  * [1] Kang G and Kan Q 2017 Thermomechanical cyclic deformation of shape-memory alloys Cyclic Plasticity of Engineering Materials (John Wiley & Sons, Ltd) pp 405–530
  * [2] Kan Q and Kang G 2010 Int. J. Plast. 26 441–465
  * [3] Yu C, Kang G, Kan Q and Xu X 2017 Mech. Mater. 112 1–17
  * [4] Yu C, Kang G and Kan Q 2017 Acta Mech. Sin. 33 619–634
  * [5] Bhaumik S K, Ramaiah K V and Saikrishna C N 2014 Nickel-titanium shape memory alloy wires for thermal actuators Micro and Smart Devices and Systems Springer Tracts in Mechanical Engineering (Springer, New Delhi) pp 181–198
  * [6] Dilibal S, Sahin H, Dursun E and Engeberg E D 2017 Electr. Eng. 99 923–930
  * [7] Baxevanis T, Cox A and Lagoudas D C 2014 Acta Mech. 225 1167–1185
  * [8] Patoor E, Lagoudas D C, Entchev P B, Brinson L C and Gao X 2006 Mech. Mater. 38 391–429
  * [9] Povoden-Karadeniz E, Cirstea D C, Lang P, Wojcik T and Kozeschnik E 2013 Calphad 41 128–139
  * [10] Frost M, Sedlák P, Sippola M and Šittner P 2010 Smart Mater. Struct. 19 094010
  * [11] Johnson L and Arróyave R 2016 Mater. Des. 107 7–17
  * [12] Tapia G, Johnson L, Franco B, Karayagiz K, Ma J, Arroyave R, Karaman I and Elwany A 2017 J. Manuf. Sci. Eng 139 071002–071002–13
  * [13] Steinbach I and Shchyglo O 2011 Curr. Opin. Solid State Mater. Sci. 15 87–92
  * [14] Ke C b, Cao S s, Ma X and Zhang X p 2012 Trans. Nonferrous Met. Soc. China 22 2578–2585
  * [15] Ke C B, Cao S S and Zhang X P 2015 Modelling Simul. Mater. Sci. Eng. 23 055008
  * [16] Yu C, Kang G and Kan Q 2017 Acta Mech. Sin. 33 619–634
  * [17] Yu C, Kang G, Kan Q and Xu X 2017 Mech. Mater 112 1–17
  * [18] Haller L, Nedjar B, Moumni Z, Vedinaş I and Trană E 2016 Continuum Mech. Thermodyn. 28 957–975
  * [19] Kan Q and Kang G 2010 Int. J. Plast. 26 441–465
  * [20] Jiang D, Kyriakides S and Landis C M 2017 Extreme Mech. Lett. 15 113–121
  * [21] Lagoudas D, Hartl D, Chemisky Y, Machado L and Popov P 2012 Int. J. Plast. 32-33 155–183
  * [22] Saleeb A F, Padula S A and Kumar A 2011 Int. J. Plast. 27 655–687
  * [23] Bodaghi M, Damanpack A R, Aghdam M M and Shakeri M 2014 Int. J. Eng. Sci. 82 1–21
  * [24] Song Z and Dai H H 2015 Int. J. Eng. Sci. 88 40–63
  * [25] Zhu P, Feng P, Sun Q P, Wang J and Dai H H 2016 Int. J. Plast. 85 52–76
  * [26] Mirzaeifar R, DesRoches R, Yavari A and Gall K 2013 Acta Mater. 61 4542–4558
  * [27] Yu C, Kang G, Kan Q and Song D 2013 Int. J. Plast. 44 161–191
  * [28] Yu C, Kang G and Kan Q 2015 J. Mech. Phys. Solids 82 97–136
  * [29] Heinen R and Miro S 2012 Computer Methods in Applied Mechanics and Engineering 229-232 44–55
  * [30] Lagoudas D C, Entchev P B, Popov P, Patoor E, Brinson L C and Gao X 2006 Mech. Mater. 38 430–462
  * [31] Choi H J, Mcdowell D L, Allen J K and Mistree F 2008 Eng. Optimiz. 40 287–307
  * [32] Hill M C 2000 In Building Partnerships 1–10
  * [33] Kim Y J, Kim K C, Park C S and Kim I H 2014 In The 2nd Asia conference of International Building Performance Simulation Association 28–29
  * [34] Larssen T, Huseby R B, Cosby B J, Høst G, Høgåsen T and Aldrin M 2006 Environ. Sci. Technol. 40 7841–7847
  * [35] Jackman S 2000 Am. J. Political Sci. 44 375–404
  * [36] Olbricht W, Chatterjee N D and Miller K 1994 Phys. Chem. Minerals 21 36–49
  * [37] Lynch S M 2007 Introduction to Applied Bayesian Statistics and Estimation for Social Scientists (Springer Science & Business Media)
  * [38] Saltelli A, Annoni P, Azzini I, Campolongo F, Ratto M and Tarantola S 2010 Comput. Phys. Commun. 181 259–270
  * [39] Fisher Ronald Aylmer S 1966 The design of experiments 8th ed (Edinburgh ; London : Oliver & Boyd)
  * [40] Ghanem R and Red-Horse J 1999 Physica D: Nonlinear Phenomena 133 137–144
  * [41] Paglietti A 1977 Ann. Inst. Henri Poincaré, Section (A): Physique Theorique 27 207–219
  * [42] Coleman B D and Noll W 1963 Arch. Rational Mech. Anal. 13 167–178
  * [43] Tschopp M A and Hernandez-Rivera E 2017 Quantifying similarity and distance measures for vector-based datasets: Histograms, signals, and probability distribution functions ARL-TN-0810 US Army Research Laboratory Aberdeen Proving Ground United States, US Army Research Laboratory Aberdeen Proving Ground United States
  * [44] Billera L J and Diaconis P 2001 Statistical Science 16 335–339
  * [45] Gelman A, Carlin J B, Stern H S, Dunson D B, Vehtari A and Rubin D B 2014 Bayesian Data Analysis vol 2 (Boca Raton, FL: CRC press)
  * [46] Haario H, Laine M, Mira A and Saksman E 2006 Stat. Comput. 16 339–354
  * [47] Haario H, Saksman E and Tamminen J 2001 Bernoulli 7 223–242
  * [48] Putko M M, Taylor Arthur C I, Newman P A and Green L L 2001 J. Fluids Eng. 124 60–69
  * [49] Kriegesmann B 2012 Probabilistic Design of Thin-walled Fiber Composite Structures (ISD)
  * [50] Tellinghuisen J 2001 J. Phys. Chem. A 105 3917–3921
  * [51] Montgomery D C 2017 Design and Analysis of Experiments (John Wiley & Sons)

## Appendix A: Two-way ANOVA

In the case of two-way ANOVA where there are $L$ observations and two multi-
level factors/parameters/variables $A$ and $B$, system response can be defined
as a Means model as follows [51]:

$y_{ijl}=\mu_{ij}+\epsilon_{ijl}=\mu+\tau_{i}+\beta_{j}+(\tau\beta)_{ij}+\epsilon_{ijl}\
\ ,\ \ \Bigg{\\{}\begin{array}[]{l}i=1,...,I\\\ j=1,...,J\\\
l=1,...,L\end{array}$ (32)

where $I$ and $J$ are the number of levels associated with
factors/parameters/variables A and B, and $\tau_{i}$, $\beta_{j}$, and
$(\tau\beta)_{ij}$ denote the effects of the level $i$ of factor $A$, level
$j$ of factor $B$, and their interaction, respectively. $\epsilon$ is the
residual or the random error which is considered as an independent normally
distributed function with a fixed variance, $\mathcal{N}(0,\sigma^{2})$. In
this context, the following hypotheses are statistically evaluated to find the
most sensitive factors/parameters/variables and interactions in the system
[51]:

$\displaystyle H_{0}:\ \ \tau_{1}=\tau_{2}=...=\tau_{I}=0$ $\displaystyle
H_{1}:\ \ \textnormal{at least one }\tau_{i}\neq 0$ (33) $\displaystyle
H_{0}:\ \ \beta_{1}=\beta_{2}=...=\beta_{J}=0$ $\displaystyle H_{1}:\ \
\textnormal{at least one }\beta_{j}\neq 0$ (34) $\displaystyle H_{0}:\ \
(\tau\beta)_{ij}=0\ \ \textnormal{for all $i$ and $j$}$ $\displaystyle H_{1}:\
\ \textnormal{at least one }(\tau\beta)_{ij}\neq 0$ (35)

In each case, $H_{0}$ and $H_{1}$ are null and alternative hypothesis,
respectively. It should be noted that a component (either a
factor/parameter/variable or an interaction) is identified sensitive when the
corresponding null hypothesis is rejected after hypothesis testing. In this
regard, a variance partitioning into the components of the system is required;
however, this partitioning can be performed on the total sum of squares
instead due to the proportionality between total variance and total sum of
squares.

Let $y_{i..}$, $y_{.j.}$, $y_{ij.}$, and $y_{...}$ denote the sum of all the
factorial design responses which includes the effects of the level $i$ of the
factor/parameter/variable $A$, the level $j$ of the factor/parameter/variable
$B$, the combination of the level $i$ for $A$ and the level $j$ for $B$, and
all the level-factor/parameter/variable combinations, respectively. Moreover,
let $\bar{y}_{i..}$, $\bar{y}_{.j.}$, $\bar{y}_{ij.}$, and $\bar{y}_{...}$
denote the corresponding averages as follows [51]:

$\displaystyle
y_{i..}=\displaystyle\sum_{j=1}^{J}\displaystyle\sum_{l=1}^{L}y_{ijl}\ \ ,\ \
\bar{y}_{i..}=\frac{y_{i..}}{JL}$ (36) $\displaystyle
y_{.j.}=\displaystyle\sum_{i=1}^{I}\displaystyle\sum_{l=1}^{L}y_{ijl}\ \ ,\ \
\bar{y}_{.j.}=\frac{y_{.j.}}{IL}$ (37) $\displaystyle
y_{ij.}=\displaystyle\sum_{l=1}^{L}y_{ijl}\ \ ,\ \
\bar{y}_{ij.}=\frac{y_{ij.}}{L}$ (38) $\displaystyle
y_{...}=\displaystyle\sum_{i=1}^{I}\displaystyle\sum_{j=1}^{J}\displaystyle\sum_{l=1}^{L}y_{ijl}\
\ ,\ \ \bar{y}_{...}=\frac{y_{...}}{IJL}$ (39)

Total sum of squares ($SS_{T}$) can be decomposed to the sum of squares of the
system components, as follows:

$\displaystyle
SS_{T}=\displaystyle\sum_{i=1}^{I}\displaystyle\sum_{j=1}^{J}\displaystyle\sum_{l=1}^{L}(y_{ijl}-\bar{y}_{...})^{2}$
$\displaystyle=\displaystyle\sum_{i=1}^{I}\displaystyle\sum_{j=1}^{J}\displaystyle\sum_{l=1}^{L}\Big{[}(\bar{y}_{i..}-\bar{y}_{...})+(\bar{y}_{.j.}-\bar{y}_{...})$
(40) $\displaystyle\ \
+(\bar{y}_{ij.}-\bar{y}_{i..}-\bar{y}_{.j.}+\bar{y}_{...})+({y}_{ijk}-\bar{y}_{ij.})\Big{]}^{2}$
$\displaystyle=JL\displaystyle\sum_{i=1}^{I}(\bar{y}_{i..}-\bar{y}_{...})^{2}+IL\displaystyle\sum_{j=1}^{J}(\bar{y}_{.j.}-\bar{y}_{...})^{2}$
$\displaystyle\ \
+L\displaystyle\sum_{i=1}^{I}\displaystyle\sum_{j=1}^{J}(\bar{y}_{ij.}-\bar{y}_{i..}-\bar{y}_{.j.}+\bar{y}_{...})^{2}$
$\displaystyle\ \
+\displaystyle\sum_{i=1}^{I}\displaystyle\sum_{j=1}^{J}\displaystyle\sum_{l=1}^{L}({y}_{ijl}-\bar{y}_{ij.})^{2}$
$\displaystyle=SS_{A}+SS_{B}+SS_{AB}+SS_{E}$

where $SS_{A}$, $SS_{B}$, $SS_{AB}$, $SS_{E}$ are the sum of squares
associated with $A$, $B$, $A$-$B$ interaction, and the random error,
respectively. It is worth noting that $SS_{E}$ can be obtained by the
subtraction of $SS_{A}$, $SS_{B}$, and $SS_{AB}$ from $SS_{T}$. In equation
40, the first equality results from the addition and subtraction of the same
terms on the right side, and the second equality is a consequence of the
polynomial expansion of the function where the six cross products becomes
zero. If these sum of squares are divided by their corresponding degrees of
freedom, the mean squares ($MS$) are obtained. The expected values of these
mean squares when the levels are random variables can provide relationships to
estimate the total deviations of all the levels associated with $A$, $B$, and
$A$-$B$ interaction, as follows:

$\displaystyle
E(MS_{A})=E\Big{(}\frac{SS_{A}}{I-1}\Big{)}=\sigma^{2}+\frac{JL\displaystyle\sum_{i=1}^{I}\tau_{i}^{2}}{I-1}$
(41) $\displaystyle
E(MS_{B})=E\Big{(}\frac{SS_{B}}{J-1}\Big{)}=\sigma^{2}+\frac{IL\displaystyle\sum_{j=1}^{J}\beta_{j}^{2}}{J-1}$
(42) $\displaystyle
E(MS_{AB})=E\Big{(}\frac{SS_{AB}}{(I-1)(J-1)}\Big{)}=\sigma^{2}+\frac{L\displaystyle\sum_{i=1}^{I}\displaystyle\sum_{j=1}^{J}(\tau\beta)_{ij}^{2}}{(I-1)(J-1)}$
(43) $\displaystyle
E(MS_{E})=E\Big{(}\frac{SS_{E}}{IJ(L-1)}\Big{)}=\sigma^{2}$ (44)

It should be noted that the mean square of each factor/parameter/variable or
interaction can be used in the above equations in the problems that the levels
are fixed. In the case that the null hypothesis about the insensitivity of a
factor/parameter/variable or an interaction is true, the total corresponding
effects of the component becomes zero; therefore, the corresponding expected
value in equation 41, 42 or/and 43 equals the expected value for the mean
squares of the random error, i.e., its variance ($\sigma^{2}$). On the other
hand, if there are any effects from the system components, the corresponding
expected values are greater than $\sigma^{2}$. It is clear that bigger ratio
of the expected value of a component to $\sigma^{2}$ is equivalent to more
sensitivity of the component. This is also true for the ratio of the
corresponding mean squares, which is called $F$-ratio. However, a certain
criterion is required to fully reject the null hypotheses for the
identification of the sensitive components in the system. For this reason, a
$F$-distribution is created from anyone of the null hypotheses based on the
numerator and denominator degrees of freedom of $F$-ratio to calculate the
corresponding p-value. Generally, the area underneath the F-distribution
confined between the F-ratio value and infinity is introduced as the p-value.
Different significance levels can be set as the criteria to reject the null
hypotheses, i.e., 0.01, 0.05, or 0.1. For any system component, a p-value less
than the significance level rejects the corresponding null hypothesis, and
demonstrates the sensitivity of the given component.

In the case that there is just one observation/response for each level-
factor/parameter/variable combinations (generally for the models and
simulations where their response is deterministic such as our thermo-
mechanical model), the variance of the random error ($\sigma^{2}$) can not be
determined by the expected mean squares of the random error in equation 44
since $L$ equals $1$. Under this condition, no clear differentiation can be
established between the effect of the interaction and the random error in
equation 43. Therefore, the interaction effect in equation 32 must be
considered zero in order to be able to estimate $\sigma^{2}$ using equation
43; otherwise, the effect of each individual factor/parameter/variable can not
be estimated through equations 41 and 42. ANOVA table for this case can be
observed in table 4 [51].

Table 4:  Two-way ANOVA table for the case of one observation/response per
level-factor/parameter/variable combination [51].

  * Source of | Sum of | Degrees of | Mean | Expected Mean | F-
---|---|---|---|---|---
Variation | Squares | Freedoms | Square | Square | Ratio
$A$ | $J\displaystyle\sum_{i=1}^{I}(\bar{y}_{i.}-\bar{y}_{..})^{2}$ | $I-1$ | $MS_{A}$ | $\sigma^{2}+\frac{J\sum\tau_{i}^{2}}{I-1}$ | $\frac{MS_{A}}{MS_{E}}$
$B$ | $I\displaystyle\sum_{j=1}^{J}(\bar{y}_{.j}-\bar{y}_{..})^{2}$ | $J-1$ | $MS_{B}$ | $\sigma^{2}+\frac{I\sum\beta_{j}^{2}}{J-1}$ | $\frac{MS_{B}}{MS_{E}}$
Residual or $AB$ | Subtraction | $(I-1)(J-1)$ | $MS_{E}$ | $\sigma^{2}+\frac{\sum\sum(\tau\beta)_{ij}^{2}}{(I-1)(J-1)}$ |
Total | $\displaystyle\sum_{i=1}^{I}\displaystyle\sum_{j=1}^{J}(y_{ij}-\bar{y}_{..})^{2}$ | $IJ-1$ |  |  |