Stochastic uncertainty analysis of damage evolution computed through microstructure-property relations

Jump to: navigation, search

Journal Probabilistic Engineering Mechanics 25 (2010) 198-205
Authors Erdem Acar, Kiran N. Solanki, Masoud Rais-Rohani, Mark F. Horstemeyerd
Paper PDF File:Acar et al 2009.pdf


Uncertainties in material microstructure features can lead to uncertainty in damage predictions based on multiscale microstructure-property models. This paper presents an analytical approach for stochastic uncertainty analysis by using a univariate dimension reduction technique. This approach is used to analyze the effects of uncertainties pertaining to the structure-property relations of an internal state variable plasticity-damage model that predicts failure. The results indicate that the higher the strain the greater the scatter in damage, even when the uncertainties in the material plasticity and microstructure parameters are kept constant. In addition, the mathematical sensitivity analysis results related to damage uncertainty are consistent with the physical nature of damage progression. At the beginning, the initial porosity and void nucleation are shown to drive the damage evolution. Then, void coalescence becomes the dominant mechanism. And finally, when approaching closer to failure, fracture toughness is found to dominate the damage evolution process.



Accurate predictions of damage progression and failure in ductile materials require the capturing of history effects and the modeling of correlations among the various physical scales present in the material, ranging from the atomic-level interactions and the microstructure features to the macroscale behavior at the continuum level. With the help of mathematical models that can capture the multiscale microstructure-property relations, it would be possible to more accurately relate the structural responses, such as stress, strain, and damage, to key material parameters such as grain size, particle size, interfacial strength, and porosity.

Previous works on multiscale microstructure-property modeling includes those of Bammann et al., [1] [2] Horstemeyer [3], Ganapathysubramanian and Zabaras [4] and Shilkrot et al. [5] . For a survey of recent progress in multiscale microstructure-property relationship modeling and simulations, the reader is referred to [6] .

Although microstructure-property relations enable the modeling of history effects as well as the damage progression and failure, the presence of uncertainties in material microstructure features can lead to considerable variation in failure predictions. Recently, Horstemeyer et al. [7] and Solanki [8] used a first-order Taylor series (FOTS) uncertainty analysis to investigate the effects of stochastic uncertainties in the microstructure features and the boundary conditions that characterize the damage evolution in AA356-T6 cast aluminum alloy. In particular, void-nucleation, void-growth, and void-coalescence equations were evaluated and quantified in terms of the sensitivity and stochastic uncertainty of various parameters in the constitutive equations. However, the accuracy of the Taylor expansion method largely depends on the scale of uncertainty in the random parameters and the nonlinearity in the corresponding random response.

This paper extends the work of Horstemeyer et al. [7] by performing a more accurate stochastic uncertainty analysis using the univariate dimension reduction (UDR) technique [9] . UDR is an additive decomposition technique that evaluates the multidimensional integral of a random function by solving a series of one dimensional integrals. As such, UDR offers an efficient approach for the evaluation of statistical moments such as the mean, variance, skewness, and kurtosis of a random response. To fully characterize the uncertainty in a random response such as damage, the estimated moments can be used to find its probability distribution using a suitable probability distribution fitting technique.

Commonly used distribution fitting techniques include the Pearson and Johnson families of distributions, saddlepoint approximations, and generalized lambda distributions. In this work, we use the extended generalized lambda distribution (EGLD), which combines the generalized lambda distribution (GLD) with the generalized beta distribution (GBD). The efficiency and accuracy of distribution fitting via the GLD and the GBD can be found in [10] [11] [12] [13] . After the first four statistical moments of damage are calculated, the parameters of the EGLD are estimated by minimizing the differences between the moments of the EGLD and those obtained through UDR.

The overall procedure of UDR+EGLD uncertainty analysis is as follows. First, damage calculations are performed at some selected points in the random variable space. Next, the first four statistical moments of damage are estimated using UDR. Then, the probability distribution of damage is approximated with the EGLD. The accuracy of this analysis depends on the accuracies of the estimated moments and distribution fitting. Acar et al. [14] have also successfully used the UDR+EGLD approach in the field of structural reliability.

The main focus of this paper is to analyze the effects of stochastic uncertainties in microstructure features of the material on the uncertainty in damage, which is calculated using the microstructure -property relations in a finite element analysis (FEA) of the selected component. Moreover, we aim to quantify the influence of uncertainty of the individual parameters in the constitutive equations on the uncertainty in damage. In metallurgical studies, it is difficult to independently quantify the effects of microstructural parameters when complex interactions are inherent. The quantitative predictions of the numerical study are presented in light of metallurgical findings to assure that the numerical results provide insights that are consistent with physical observations.

This paper is organized as follows. A brief description of the microstructure -property relationship model of Horstemeyer et al. [3] is given in the next section. Section 3 discusses material microstructure characterizations and their uncertainties. Section 4 presents the calculation of statistical moments of damage using univariate dimension reduction. Section 5 gives a brief description of probability distribution fitting by using the extended generalized lambda distribution. Section 6 presents the results of the uncertainty analysis for damage, where sensitivities of various factors are also compared. Finally, concluding remarks appear in Section 7.

Microstructure-property relationships

An effective method to capture the microstructure-property relationships is by use of internal state variable (ISV) evolution equations, which are formulated at the macroscale level. The ISVs reflect lower length scale microstructural rearrangements so that history effects can be modeled. With the help of such a material model, it would be possible to relate structural responses of interest, such as stress, strain, and toughness, to key material parameters such as particle size, interfacial strength, and spacing at the lower length scale.

The microstructure-property relationship modeling framework used here is that developed by Bammann et al. [1] [2] and extended by Horstemeyer [3] to account for stress-state-dependent damage evolution. The pertinent equations in this model are denoted by the rate of change of the observable and the internal state variables. For the sake of clarity and completeness, a listing of these equations and their relation to the material microstructure is briefly given here with additional details provided in [3]. The first equation is the modified Hooke's law that includes damage, and is given as

\underline{\overset{\circ}{\sigma}} =


where \mu' and \lambda' are the temperature dependent shear modulus and the Lamé constant given as follows:

\mu'=\mu \left[ 1 - \frac{T}{T_{melt}}\exp \left(aa \left(1-\frac{T}{T_{melt}} \right) \right)\right],

\lambda'=K-bb\frac{T}{T_{melt}} - \frac{2}{3}\mu'


where aa is the shear modulus temperature dependent coefficient, bb is the Lamé constant temperature dependent coefficient, T_{melt} is the melting temperature in Kelvin, T is the current temperature in Kelvin, \mu and K are shear and bulk moduli of base material, \underline{\sigma} and \underline{\overset{\circ}{\sigma}} are the Cauchy stress and the co-rotational rate of the Cauchy stress, respectively, \phi is an ISV that represents the damage fraction or state of material within a continuum element in the context of FEA with \dot{\phi} representing its material time derivative, \underline{D^e} is the elastic deformation tensor, and \underline{I} is the second-order identity tensor. The underscore symbol indicates a second-rank tensor. Recognizing that \underline{D^e}=\underline{D}-\underline{D^p}, the ISV representing the plastic deformation tensor or inelastic flow rule, \underline{D^p}, is given by the relationship

\frac{\underline{\sigma'}-\underline{\alpha}}{\left\Vert \underline{\sigma'}-\underline{\alpha}\right\Vert}


where \underline{\sigma'} is the deviatoric part of stress tensor, T is the temperature in Kelvin, \underline{\alpha} is the kinematic hardening (an ISV reflecting the effect of geometrically necessary dislocation) and R is the isotropic hardening (an ISV reflecting the effect of statistically stored dislocation density). The function V(T) determines the magnitude of rate-dependence on yielding, Y(T) is the rate-independent yield stress, and f(T) determines when the rate-dependence affects initial yielding. The functions f(T), V(T), and Y(T), which are related to yielding with Arrhenius-type temperature dependence, are given as

V(T)=C_1 e^{(-C_2/T)}


Y(T)=C_3 e^{(-C_4/T)}


f(T)=C_5 e^{(-C_6/T)}


where C_1 through C_6 are the yield-stress-related material parameters obtained through different monotonic stress-state tests (tension, compression, and torsion) at different temperatures and strain rates. The evaluation of \underline{D^p} in Eq. (2) also requires the corotational rate of the kinematic hardening, \overset{\circ}{\underline{\alpha}} , and the material time derivative of isotropic hardening, \dot{R}, expressed in a hardening-recovery format as

\underline{\overset{\circ}{\alpha}} = \left\{ 


\dot{R} = \left\{


where DCS0, DCS and z capture the microstructure effect of grain size. In Eqs. (4) and (5), r_d(T) and R_d(T) are scalar functions that describe dynamic recovery whereas r_s(T) and R_s(T) are scalar functions that describe thermal (static) recovery with h(T) and H(T) representing the anisotropic and isotropic hardening modulus, respectively. Hence, the two main types of recovery that are exhibited by populations of dislocations within crystallographic materials are captured in the ISVs. The temperature-dependent functions in Eqs. (4) and (5) are found as





\right] e^{(-C_8/T)}


\right] e^{(-C_{14}/T)}


\right] e^{(-C_{8}/T)}-C_{10}T


\right] e^{(-C_{8}/T)}-C_{16}T


where J'_2=\frac{1}{2}(\underline{\sigma'}-\underline{\alpha})^2,  J'_3 = \frac{1}{3}(\underline{\sigma'}-\underline{\alpha})^3 C_7 through C_{12} are the material plasticity parameters related to kinematic hardening and recovery terms, C13 through C18 are the material plasticity parameters related to isotropic hardening and recovery terms, and Ca and Cb are the material plasticity parameters related to dynamic recovery and anisotropic hardening terms. Constants C1 through C18 are found from tension, compression and shear tests at different temperatures and strain rates.

The mechanical properties of a material depend upon the microdefects within its structure that can change as a result of deformation. When the number of microdefects accumulates, damage is said to have grown. The three components of damage progression mechanism are void nucleation, growth and coalescence from second-phase particles and pores. In this regard, the material time derivative of damage, \dot{\phi}, is expressed as

\dot{\phi}=\left(\dot{\phi}_{particles}+\dot{\phi}_{pores}\right)C +


\dot{\phi_{particles}} = \dot{\eta}\upsilon + \eta\dot{\upsilon}


e^{\left(-C_{{\eta T}/T}\right)}




\right] \times
\frac{2(2^{V(T)}/Y(T)^{-1})}{(2^{V(T)}/Y(T)^{+1})}\frac{\sigma_H}{\sigma_{\upsilon m}}


\dot{C}=[Cd_1 +Cd_2 (\eta\dot{\upsilon}+\dot{\eta}{\upsilon})]e^{(C_{CT}T)}(DCS_0/DCS)^z


where \upsilon is the void-growth, \eta is the void-nucleation, d is the particle size, f is the particle volume fraction, K_{IC} is the fracture toughness, \sigma_H and \sigma_{vm} are the hydrostatic and von Mises stresses, J1; J2; J3 and I1 are stress invariants, C_{coeff} is the void-nucleation coefficient parameter, C_{CT} and C_{\eta T} are the void-coalescence and the void-nucleation temperature-dependent parameters, Cd1 and Cd2 are related to the first and second normalized nearest-neighbor distance parameters, respectively, and constants a, b, and c are the stress-state-based void-nucleation constants.

In Eqs. (12) and (13), the damage progression is divided into void-nucleation and void-growth from silicon particles and pores. Eqs. (14) and (15) describe void nucleation evolution and the void growth related to silicon particles, respectively. For the porosity evolution, the void-growth rule given in Eq. (16) is used. Void coalescence is introduced to reflect pore-pore interactions and silicon particlepore interactions, as expressed in Eq. (17). In the microstructure-property relationship model, bulk and shear moduli, material plasticity parameters (C1 to C18) and microstructure damage parameters (z; DCS; DCS0; a; b; c; f ; d; KIC ; C_{coal}, C_{CT}, C_{\eta T}, m, initial void volume fraction, R_0) are treated as random variables. The time integral form of Eq. (12), which is the damage state, is used as a damage index to assess the failure. Additional descriptions related to material selection and the random variables are provided next.

Material microstructure features and associated uncertainties

Fig. 1. Optical micrograph of cast AA356-T6 showing second-phase particles [3]

The microstructure of a typical metallic material contains a large number of microdefects such as microcracks, dislocations, pores, and decohesions. Some of these defects are induced during the manufacturing process and are present before the material is subjected to mechanical loads and thermal fields. In general, these defects are small and distributed throughout most of the volume.

In this paper,we focus on cast AA356-T6. Cast aluminum-silicon (Al-Si) alloys contain a mixture of a eutectic phase (Si particles embedded in an Al-1% Si matrix) and a proeutectic Al-1% Si phase. Depending on the alloy type, trace amounts of other elements are also present to promote precipitation hardening or improve other casting properties. The microstructure of this material consists of the primary (Al--1.6 wt% Si) and eutectic (Al--12.6 wt% Si) phases. In the eutectic regions, large silicon particles and clusters form a dendritic substructure while the Si remains solutionized in the primary phase. The microstructural alterations have a strong influence on the monotonic mechanical properties of Al castings through changes in void nucleation, growth, and coalescence characteristics. The optical microscope image of the cast alloy sample, as shown in Fig. 1, suggests that the material is more isotropic in nature. The secondary constituent particles within the matrix were found to range in size from 3 to 6 \mu m, with the silicon particles ranging from 3 to 5 \mu m, with an average size of 4 \mu m and particle volume fraction of approximately 7% (5.25% to 8.75%), as shown in Fig. 2.

Table 1. Elastic-plastic model constants for cast AA356-T6.
Fig. 2. AA356-T6 microstructure probability distribution function plots for (a) particle size and (b) particle volume fraction.

A nonlinear regression algorithm that was developed in the Cast Light Metals Project [3] was used to model the uncertainty in the plasticity model constants listed in Table 1, for the A356-T6 aluminum alloy. The experimental uncertainty, UE , in each model constant is calculated based on random and systematic uncertainties in the measured items such as force, strain, and specimen size (Eqs. (17)-(19)) using the formula

U_E = \sqrt{U_r^2+U_s^2}


where Ur and Us represent the random and systematic uncertainty, respectively, in the measurements, and are calculated as

U_r = 2\sqrt{\frac{1}{M-1}\sum_{i=1}^M (r_i - r_{mean})^2}


U_s = \sqrt{U_L^2+U_{daq}^2}


where M is the number of experiments used to measure item r_i (e.g., force, strain), U_L is the percentage uncertainty in the load cell (1%), extensometer, or strain gauges, and U_{daq} is the percentage uncertainty in the data acquisition (0.25%) [8].

Fig. 3. Uncertainty quantification methodology for the ISV model.

As shown in Fig. 3, the experimental uncertainties and uniaxial experimental data are propagated through a model correlation routine to predict uncertainties in the plasticity material parameters. The calculated first-order coefficients of variation (COVs) of the plasticity material parameters related to the yield strength, the kinematic and the isotropic hardenings, and the shear and the bulk moduli along with their temperature and strain-rate dependences (Eqs. (1)-(10)) are as shown in Table 1. The COV of each parameter is calculated from

COV_s = \frac{\sqrt{U_s}}{\mu_s}


where \mu_s is the mean value of the parameter. The calculated COVs are assumed to be normally distributed and rounded off due to limited test data. The proposed COVs are significant indicators of scatter in the given parameters that enable the study of their influences on damage as well as the overall product design. Eqs. (12)-(17) include the damage progression, which is multiplicatively decomposed into void-nucleation, void-growth, and void-coalescence variables. The quantitative fractographic measurements were used to quantify the variability (COV values) in microstructure spatial clustering, as given in Table 2.

Table 2. Damage model constants for cast AA356-T6.

A probabilistic model for the 38 random variables listed in Tables 1 and 2 is better based on both expert judgment and experimental data by utilizing a Bayesian framework. However, for simplicity, we only considered uncertainties measured via experiments as indicated in the block diagram shown in Fig. 3. Furthermore, we assumed that all the random variables are independent and follow a normal distribution; therefore, the joint probability density function (PDF) of these variables can be obtained using only the marginal PDFs. It is worth noting that these assumptions by no means diminish the generality of the proposed approach for stochastic uncertainty analysis.

We also recognize that more extreme values than indicated in Fig. 2 and Tables 1 and 2 could exist in the material for all of the different parameters. However, the COVs in Fig. 2 and Tables 1 and 2 represent a statistically significant scatter in the selected parameters. The proposed COVs are sufficient to screen the relative influence of the parameters on the debonding and fracture characteristics of the silicon particles. Once these first-order defects are understood, the more dominant parameters can be studied over a wider range of values and distributions in future studies.

Calculation of statistical moments of damage using univariate dimension reduction

As noted earlier, uncertainty analysis requires the calculation of the statistical moments of the response function. Note that, in this case, the response function of interest is damage. The easiest way to estimate the first two statistical moments of a response function Y(X) is to use a first-order Taylor series (FOTS) approximation, which gives

\mu_Y = Y(\mu_X)


\operatorname{Var}_Y = \sum_{i=1}^{N} \left(\frac{\partial Y}{\partial X_i}\right)^2 \operatorname{VarX}_i


where \mu and Var are the mean and variance of Y(X), respectively, and N is the number of input random variables described by the vector X^T=(X_1, X_2,...,X_N). This formulation assumes a linear response function and normally distributed uncorrelated input variables. Hence, the estimations from Eqs. (21) and (22) are subject to error when the response function is nonlinear.

The accuracy of uncertainty prediction using a Taylor series can be improved by including the second-order terms. In this case, the mean and variance of the response function can be approximated as

\mu_Y = Y(\mu_X)+
\frac{1}{2}\sum_{i=1}^{N} \left(\frac{\partial^2 Y}{\partial X_i^2}\right)


\operatorname{Var}_Y = 
\left(\frac{\partial Y}{\partial X_i}\right)^2 \operatorname{VarX}_i
-\frac{1}{4}\sum_{i=1}^{N} \left(\frac{\partial Y}{\partial X_i}\right)\operatorname{Var}^2 X_i 
\right .

\left .
+ \frac{\partial Y}{\partial X_i} \frac{\partial^2 Y}{\partial X^2_i}
E(X_i - \mu_{x_i})^3 
+ \frac{1}{4}\left(\frac{\partial Y}{\partial X_i}\right)E(X_i - \mu_{x_i})^4

where E[ ] is the expectation operator. Notice that since the second-order derivatives of the response function are needed, this approximation is numerically expensive. Even though the secondorder derivatives are used in the formulation, uncertainty analysis by Taylor series expansion can still lead to erroneous predictions for nonlinear responses.

For uncertainty analysis, a computationally efficient alternative is the use of stochastic response surface methodology (SRSM) [15] [16]

Here, a functional form is assumed for the response in terms of the random input parameters (e.g., polynomial chaos expansion of Hermite polynomials). Upon determining the parameters of the functional approximation, the statistical moments of the responses can easily be computed. The main drawback of the SRSM is that if the response is highly nonlinear, higher-order Hermite polynomials must be used, and the number of terms in the polynomials grows rapidly as the degree of polynomial increases.

The statistical moments (about the origin) of a response function, Y(X) can be defined as

m_l \equiv E[Y^i(X)]=\int_{R^N} y^i(x)f_x(x)\, dx


where f_x(X) is the joint distribution function of all random variables in the vector X,E[ ] is the expectation operator, R^N is the N- dimensional random variable domain, and l represents the order of the statistical moment.

Depending on the functional form of the response function and the number of random variables involved, the exact calculation of the multi-dimensional integral in Eq. (25) can become very challenging, if not impossible. However, it is possible to estimate the value of this integral in various ways. The most accurate way is to use Monte Carlo simulations (MCSs). Since an MCS requires a large number of response samples, the computational cost can quickly escalate when each response evaluation requires the use of an expensive high-fidelity model. A more efficient way to estimate the multi-dimensional integral is to use the univariate dimension reduction (UDR) technique [9]. Many researchers [17] [18] [14] have successfully used UDR for the estimation of statistical moments.

in UDR, the response function Y(X) is first estimated as

\hat{Y} = \sum_{j=1}^N Y(X_j)-(N-1)Y_0



Y(X_j) = Y(\mu_1,\cdot,\mu_{j-1}, X_j, \mu_{j+1},\cdot, \mu_N)



Y_0 = Y(\mu_1,\cdot,\mu_N)


By using the binomial formula, Eq. (25) can be rewritten as

m_l \cong \sum_{i=0}^l\binom{l}{i} E \left\{Y(\mu_1,\cdot,\mu_{j-1}, X_j, \mu_{j+1},\cdot, \mu_N)\right\}^i

\times {-(N-1)Y(\mu_1,\cdot,\mu_N)}^{l-i}

Eq. (29) can be calculated recursively by computing all one-dimensional integrals simultaneously. The procedure for recursive calculation of single dimension integrals can be found in [14]. When the input random variables follow normal distributions, the one-dimensional integrals can be calculated effectively using the Gauss-Hermite quadrature rule. When the random variables follow non-normal distributions, a quadrature rule [9] can be used. Xi et al. [18], however, showed that, for some cases, the quadrature rule can cause instabilities. If so, then the function Y(X) in a specified dimension can be approximated using a metamodeling technique (e.g., moving least square regression, Kriging, radial basis functions) and the numerical integration can be performed using more effective techniques (e.g., adaptive Simpson's rule).

After the statistical moments of damage are calculated, the next step is to approximate the PDF of damage using the extended generalized lambda distribution (EGLD), which is briefly described in the next section

Estimating the probability distribution of damage using the EGLD

Approximating the distribution of a random variable using a few of its statistical moments has been of interest to many researchers. Usually, the first four moments are used for this purpose. Popular distribution fitting techniques include the Johnson distribution, Pearson distribution, saddlepoint approximations, and generalized lambda distributions. The main disadvantage of the Johnson distribution is that it is not very easy to determine the four parameters from the moments of the sample data. The major drawback of the Pearson family of distributions can lead to unstable results near the boundaries of families in the skewness-kurtosis plane [18], and for the saddlepoint approximation, the main shortcoming is that singularities can arise during the numerical procedure [17]

Another powerful tool to estimate distribution of random variables is the generalized lambda distribution (GLD), which is very flexible and can fit a wide variety of distributions [19] [20] However, [10] and [11] showed that the GLD may fail to provide a valid distribution in some regions of the moment space. They proposed using the EGLD, which is the extended version of the GLD, where the GBD is used when the accuracy of the GLD is not satisfactory. Motivated by their successful application of the EGLD in distribution fitting as well as the previous experience [14], we decided to use the EGLD for distribution fitting. The probability distribution fitting (that is, the estimation of the distribution parameters) is performed by matching the first four moments of the EGLD with the first four moments of the response function obtained through univariate dimension reduction. Since we use UDR for moment calculation and the EGLD in distribution fitting, we refer to this approach as UDR+EGLD [14].

Uncertainty analysis of damage evolution

Uncertainty characterization of damage at different strain values

Fig. 4. Damage-strain for cast AA356-T6.

An un-notched tension specimen made of castAA356-T6 is loaded at a strain rate of 10^{-3} 1/s until the point of failure. The plot of damage growth as a function of strain, ", is shown in Fig. 4. The plots of damage probability distributions at different points along the damage-strain curve (Fig. 4) are shown in Fig. 5. The scatter (COV) in damage (C_{\phi} and mean value of damage appear to increase as the strain value increases. The comparison of UDR+EGLDbased PDF curves with those found using MCS (of sample size 10,000) indicates a good agreement between the two distributions. In addition, normal distribution fits to damage are also shown in Fig. 5. Since the number of random variables is large (38 random variables), in view of the central limit theorem, the probability distribution of damage can also be represented reasonably well with a normal distribution.

Fig. 5. Probability distributions of damage at different points along the stress-strain curve.

Wealso compared the uncertainty analysis results of UDR+EGLD with those from FOTS uncertainty analysis. Table 3 shows that both the FOTS and UDR uncertainty predictions are close to the Monte Carlo simulations (with a sample size of 10,000). The observation that the FOTS predictions are in reasonably good agreement with the MCS results indicates the dominance of the first-order effects. A closer look at the values in Table 3, however, shows that the UDR predictions are slightly better than those obtained using FOTS. In addition, UDR+EGLD is computationally more efficient than FOTS. Even though FOTS requires calculation of sensitivity derivatives of the limit-state function, UDR+EGLD provides a sensitivity-free uncertainty analysis.

Table 3. Standard deviation of damage obtained through FOTS, UDR and MCS.

Carlo simulations (with a sample size of 10,000). The observation that the FOTS predictions are in reasonably good agreement with the MCS results indicates the dominance of the first-order effects. A closer look at the values in Table 3, however, shows that the UDR predictions are slightly better than those obtained using FOTS. In addition, UDR+EGLD is computationally more efficient than FOTS. Even though FOTS requires calculation of sensitivity derivatives of the limit-state function, UDR+EGLD provides a sensitivity-free uncertainty analysis.

It is also worth noting here that the results presented in Table 3 depend on the normality and independence assumptions. The use of non-normal probability distributions for the random variables and the inclusion of correlation may have a strong effect on the quantitative values of the results, while the qualitative indications still hold.

Sensitivity of damage uncertainty at different strain values

The sensitivity of damage with respect to an input random variable X_i is measured using the formula

S_i = \frac{\partial C_\phi}{\partial C_i}\frac{C_i}{C_\phi}


where C_i is the coefficient of variation of X_i. Notice that the second term is used to normalize the sensitivity factor. The sensitivity of damage uncertainty to the uncertainties of the microstructure- property parameters are depicted in Fig. 6. For instance, Fig. 6(a) shows that, at the very beginning of the damage evolution, the initial radius of a spherical void, R_0, (25th random variable; see Tables 1 and 2 for the random variable list) and nucleation coefficient (29th random variable) are the most influential parameters. Fig. 6(b)-(c), on the other hand, shows that as damage progresses, the initial temperature (23rd random variable), coalescence temperature dependence term (38th random variable), coalescence factor cd1 (33rd random variable), and initial void volume fraction (37th random variable) become more important. Finally, Fig. 6(d) displays the effect of parameters on the last stages of damage progression. The most influential terms can be ordered according to their importance as the fracture toughness (30th random variable), initial temperature (23rd random variable), and coalescence temperature dependence (38th random variable), respectively.

Fig. 6. Sensitivity of damage uncertainty to the uncertainties of the parameters in the microstructure-property relationship model.
Fig. 7. Eventual failure of a damaged material caused by nucleation, growth and coalescence [21].

The sensitivity plots shown in Fig. 6 are also consistent with the physics of the damage progression shown in Fig. 7. At the beginning, the void properties drive the damage evolution. Then, voids combine with each other and coalescence becomes the main driver. Near failure, macroscopic properties such as fracture toughness, KIC , determines the damage evolution process.

Concluding remarks

This paper has presented a new approach for stochastic uncertainty analysis using the combination of univariate dimension reduction (UDR) technique and the extended generalized lambda distribution (EGLD). This approach was used to perform an uncertainty analysis on an un-notched cast AA356-T6 specimen using an internal state variable plasticity-damage model. The effects of uncertainty in material microstructural features (i.e., voids, cracks, inclusions) on damage initiation and evolution (or accumulation) were investigated. The uncertainty analysis results based on UDR+EGLD were compared with those obtained from first-order Taylor series (FOTS) expansion and Monte Carlo simulation. The proposed approach is found to be computationally efficient and provides more accurate estimates of parametric uncertainty than the FOTS. The influence of uncertainty in individual parameters appearing in the constitutive equations on damage uncertainty was also studied.

From the results obtained in this study, we can draw the following conclusions: (i) the scatter in damage, as measured by the coefficient of variation of damage, can increase with strain even though the uncertainties in the input variables are kept fixed; (ii) the sensitivities of damage uncertainty to the uncertainties in the input random variables depend on the strain values. As the strain value changes (i.e., as damage evolves), the importance of the random variables changes; (iii) the sensitivities are found to be consistent with the physics of the damage progression. At the very beginning, the void properties (initial size and growth parameters) are found to drive the damage evolution. Then, void coalescence becomes the main driver, and finally, near the failure condition, macroscopic properties such as fracture toughness dominate the damage evolution process.


This material is based upon work supported by the Department of Energy under Award Number DE-FC26-06NT42755. This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed therein do not necessarily state or reflect those of the United States Government or any agency thereof.


  1. 1.0 1.1 Bammann DJ, Chiesa ML, Horstemeyer MF, Weingarten LI. Failure in ductile materials using finite element methods, structural crashworthiness and failure. In: Wierzbicki T, Jones N, editors. Elsevier applied science, The Universities Press (Belfast) Ltd.; 1993.
  2. 2.0 2.1 Bammann DJ, Chiesa ML, Johnson GC. Modeling large deformation and failure in manufacturing processes. In: Tatsumi T, Wannabe E, Kambe T, editors. Theoretical and applied mechanics, Elsevier Science; 1996. p. 359-76.
  3. 3.0 3.1 3.2 3.3 3.4 3.5 Horstemeyer MF. From atoms to autos a new design paradigm using microstructure-property modeling, Part 1: Monotonic loading conditions. Sandia National Laboratories, Sand2000-8662, March 2001.
  4. Ganapathysubramanian S, Zabaras N. Design across length scales: A reducedorder model of polycrystal plasticity for the control of microstructuresensitive material properties. Computer Methods in Applied Mechanics and Engineering 2004;193:5017-34.
  5. Shilkrot LE, Miller RE, Curtin WE. Multiscale plasticity modeling: Coupled atomistics and discrete dislocation mechanics. Journal of the Mechanics and Physics of Solids 2004;52(4):755-87.
  6. Graham-Brady LL, Arwade SR, Corr DJ, Gutierrez MA, Breysse D, Grigoriu M, Zabaras N. Probability and materials: From nano- to macro-scale: A summary. Probabilistic Engineering Mechanics 2006;21:193-9.
  7. 7.0 7.1 Horstemeyer MF, Solanki K, Steele WG. Uncertainty Methodologies to Characterize a Damage Evolution Model. In: Proceeding of international journal of plasticity conference. 2005.
  8. 8.0 8.1 Solanki KN. Physically motivated internal state variable form of a higher order damage model for engineering materials with uncertainty. Ph.D. dissertation. Mississippi State University (Eds.); 2008.
  9. 9.0 9.1 Rahman S, Xu H. A univariate dimension-reduction method for multidimensional integration in stochastic mechanics. Probabilistic Engineering Mechanics 2004;19(4):393-408.
  10. 10.0 10.1 Karian ZE, Dudewicz EJ, McDonald P. The extended generalized lambda distribution system for fitting distributions to data: History, completion of theory, tables, applications, the final word on moment fits. Communications in Statistics-Computation and Simulation 1996;25(3):611-42.
  11. 11.0 11.1 Karian Z, Dudewicz E. Fitting statistical distributions to data: The generalised lambda distribution and the generalised bootstrap methods. Boca Raton (Florida): CRC Press; 2000.
  12. Lampasi DA, Nicola FD, Podesta L. The generalized lambda distribution for the expression of measurement uncertainty. IEEE instrumentation and measurement conference. 2005.
  13. King RAR, MacGillivray HL. Fitting the generalized lambda distribution with location and scale-free shape functionals. In: Proceedings of the symposium on fitting statistical distributions to data. 2006.
  14. 14.0 14.1 14.2 14.3 King RAR, MacGillivray HL. Fitting the generalized lambda distribution with location and scale-free shape functionals. In: Proceedings of the symposium on fitting statistical distributions to data. 2006.
  15. Isukapalli SS, Roy A, Georgopoulos PG. Stochastic response surface methods (SRSMs) for uncertainty propagation: Application to environmental and biological systems. Risk Analysis 1998;18(3):351-63.
  16. Kim NH, Wang H, Queipo NV. Efficient shape optimization under uncertainty using polynomial chaos expansions and local sensitivities. AIAA Journal 2006; 44(5):1112-5.
  17. 17.0 17.1 Huang B, Du X. Uncertainty analysis by dimension reduction integration and saddlepoint approximations. Journal of Mechanical Design, ASME 2006; 126(1):26-33.
  18. 18.0 18.1 18.2 Xi Z, Youn BD, Gorsich DA. Reliability-based robust design optimization using the EDR method. SAE international conference, Paper no. 2007-01-0550. 2007. have successfully used UDR for the estimation of statistical moments.
  19. Freimer M, Mudholkar G, Kollia G, Lin C. A study of the generalized Tukey Lambda family. Communications in Statistics: Theory and Methods 1988; 17(10):3547-67.
  20. Lakhany A, Mausser H. Estimating the parameters of the generalized lambda distribution. Algo Research Quarterly 2000;3(3):47-58.
  21. Hammi Y, Horstemeyer MF. A physically motivated anisotropic tensorial representation of damage with separate functions for void nucleation, growth, and coalescence. International Journal of Plasticity 2007;23:1641-78.
Personal tools

Material Models