Microstructure-sensitive uncertainty quantification for crystal plasticity finite element constitutive models using stochastic collocation methods

Uncertainty quantification (UQ) plays a major role in verification and validation of computational engineering models and simulations, and establishes trust in the predictive capability of computational models. In the materials science and engineering context, where the process-structure-property-performance linkage is well known to be the only road mapping from manufacturing to engineering performance, numerous integrated computational materials engineering (ICME) models have been developed across a wide spectrum of length-scales and time-scales to relieve the burden of resource-intensive experiments. Within the structure-property linkage, crystal plasticity finite element method (CPFEM) models have been widely used since they are one of a few ICME toolboxes that allow numerical predictions, providing the bridge from microstructure to properties and performances. Several constitutive models have been proposed to capture the mechanics and plasticity behavior of materials. While some UQ studies have been performed, the robustness and uncertainty of these constitutive models have not been rigorously established. In this work, we apply a stochastic collocation (SC) method to quantify the uncertainty of the three most commonly used constitutive models in CPFEM, namely phenomenological models (with and without twinning), and dislocation-density-based constitutive models, for three different types of crystal structures, namely face-centered cubic (fcc) copper (Cu), body-centered cubic (bcc) tungsten (W), and hexagonal close packing (hcp) magnesium (Mg). Our numerical results not only quantify the uncertainty of these constitutive models in stress-strain curves, but also analyze the global sensitivity of the underlying constitutive parameters with respect to the initial yield behavior, which may be helpful for robust constitutive model calibration works in the future.


INTRODUCTION
Uncertainty quantification (UQ) has been a cornerstone in applied mathematics to verify and validate forward models, with application ranging from subsurface flow, climate change, integrated computational materials engineering (ICME) developments, advanced manufacturing, and many more.In the context of materials design, process-structure-property-performance relationship plays a critical role in establishing the linkage between manufacturing and desired properties, to which materials can be tailored for specific applications.Across the spectrum of length-and time-scales, from quantum to macro-scale, multiple ICME models have been developed in the quest of accurate prediction of materials properties and performance.Notable ICME models for deformation of metals includes, but are not limited to, concurrent or hierarchical couplings of density functional theory, molecular dynamics, kinetic Monte Carlo, dislocation dynamics, crystal plasticity, phase-field, finite element, and any sort of multi-physics and hybrid approaches.Despite its success, there are still much room and open questions for UQ research with respect to microstructure-induced mechanical response in metal alloys.
The concept of ICME models is particularly important because with the predictive capability, engineers can numerically approximate material properties and performance under different operating conditions, without actually performing physical experiments.With the rise of high-performance parallel computing, this step is often hailed as the third paradigm of science Agrawal and Choudhary (2016), supporting prior empirical and experimental, as well as theoretical research.However, as in any forward computational model that relies on deterministic calculations, there is a necessary need to quantify the uncertainty that is associated with the numerical predictions, in order to sufficiently enhance its fidelity and robustness.For the case of crystal plasticity finite element method (CPFEM) as an ICME model that bridges microstructure and materials properties and performance, the question of uncertainty is even more relevant because microstructures are well known to be stochastic, high-dimensional in image or volume representations, and may be anisotropic and heterogeneous.
Numerous UQ studies in computational solid mechanics have been conducted over the last decade.Given the critical importance of optimization and UQ for a wide variety of problems in materials science, several frameworks have been developed, see e.g., Panchal et al. (2013); McDowell (2007); Kalidindi et al. (2016), to provide robust predictions under uncertainty.A comprehensive review of UQ applications in ICME-based simulations can be found in Honarmandi and Arróyave Honarmandi and Arróyave (2020).For example, Zhao et al.Zhao et al. (2022) incorporated measurement and parametric uncertainty to quantify the uncertainty of critical resolved share stress for hcp Ti alloys from nano-indentation.Lim et al. Lim et al. (2019) investigated the mesh sensitivity and polycrystalline representative volume element (RVE), where initial textures, hardening models, and boundary conditions are uncertain, and showed that an adequate polycrystalline RVE is obtained by capturing 1000 or more grains.Tran and Wildey Tran and Wildey (2020) applied data-consistent inversion method to infer a distribution of microstructure features from a distribution of yield stress, where the push-forward density map via a heteroscedastic Gaussian process approximation is consistent with a target yield stress density.Kotha et al. Kotha et al. (2019a,b, 2020b,a) developed uncertainty-quantified parametrically homogenized constitutive models to capture uncertainty in microstructure-dependent stress-strain curve, as well as stochastic yield surface, which has been broadly applied for modeling multi-scale fatigue crack nucleation in Ti alloys Ozturk et al. (2019b,a) and for single-crystal Ni-based superalloys with support vector regression as an underlying machine learning model Weber et al. (2020).Sedighiani et al. Sedighiani et al. (2020, 2022) applied genetic algorithm and polynomial approximation to various constitutive models, including phenomenological and dislocation-density-based models.Tran et al. Tran et al. (2019b) applied stochastic collocation (SC) method to quantify uncertainty for dendrite morphology and growth via phase-field model.Acar et al. Acar et al. (2017) proposed a linear programming approach to maximize a mean of materials properties under the assumption of Gaussian distribution for both inputs and outputs.Fernadez et al. Fernandez-Zelaia et al. (2018) utilized Bayesian inference to quantify the uncertainty in stressstrain curves, where model parameters are treated as random variables.Tallman et al.Tallman et al. (2019Tallman et al. ( , 2020) ) applied Gaussian process regression and the Materials Knowledge System framework to predict a set of homogenized materials properties with uncertainty from a distribution function for crystallography orientations and textures.Inductive design exploration method (IDEM) Ellis and McDowell (2017);McDowell et al. (2009); Choi et al. (2008) has been introduced as a materials design methodology to identify feasible and robust design for microstructure features, which has been broadly applied to many practical problems.From a methodological perspective, numerous applied mathematical UQ techniques have been developed over a few decades.While both intrusive and non-intrusive options are available, non-intrusive polynomial chaos expansion (PCE) Xiu and Karniadakis (2002); Najm (2009), which is also known as non-intrusive spectral projection PCE, is one of the most widely used UQ methods to propagate uncertainty in physical models and computational simulations.Global sensitivity analysis methods, which is arguably constructed on top of the high-dimensional model representation, have also well studied.And apparently, numerous studies connecting PCE and global sensitivity analysis, for example, Sudret Sudret (2008), Crestaux et al. Crestaux et al. (2009), andSaltelli et al. Saltelli et al. (2010), have been conducted.Despite the fact that numerous UQ application works have been carried out in the materials literature and more specifically, in the structure-property relationship, to the knowledge of the authors, none has been applied to quantify the uncertainty related to the underlying constitutive models used in CPFEM.
The rest of the paper is organized as follows.Section 2 introduces PCE and SC as the UQ backbone methodology used in this paper.Section 3 describes the integrated UQ workflow and how they are implemented in practice.Section 4 shows the first case study of 5d UQ for fcc Cu using phenomenological constitutive model with slipping.Section 5 shows the second case study of 16d UQ for hcp Mg using phenomenological constitutive model with slipping and twinning.Section 6 shows the last case study of 7d UQ for bcc W using dislocation-density-based constitutive model.Section 7 discusses and Section 8 concludes the paper, respectively.

UNCERTAINTY QUANTIFICATION BACKGROUND
In this section, we describe UQ background for CPFEM.In Section 2.1, we summarize the theoretical foundation for generalized polynomial chaos expansion as a non-intrusive spectral projection method.In Section 2.2, we provide the mathematical background for Smolyak sparse grid construction for high-dimensional interpolation and integration, as well as some comparison to full tensor grid highlighting the computational advantage in SC method.We refer interested readers to Babuška et al. (2007); Nobile et al. (2008); Xiu (2009) for a more rigorous mathematical characterization of the SC method.

Generalized polynomial chaos expansion
The generalized Wiener-Askey PCE Xiu and Karniadakis (2002) represents the second-order random process f pθq as where I n pξ i 1 , ¨¨¨, ξ in q denotes the Wiener-Askey polynomial chaos of order n in terms of the random vector ξ " pξ i 1 , ξ i 2 , . . ., ξ in q, and c's are polynomial chaos expansion coefficients to be determined.Without loss of generality, Equation 1 can be rewritten as where there is a one-to-one correspondence between the function I n pξ i 1 , ¨¨¨, ξ in q and Φ j pξq.Φ j pξpθqq are orthogonal polynomials in terms of ξ :" tξ i pθqu d i"1 , i.e., where δ ij is the Kronecker delta and x¨, ¨y denotes the weighted-average, which is defined as the inner product in the Hilbert space of the variable ξ with respect to the weighting function W pξq, described in Table 1, as Here, p f j are the coefficients to be determined.In practice, the number of terms in (2) are truncated after a finite term P , where P `1 " pp `nq! p!n! , where p is the order of PCE, and n is the dimensionality of the problem, resulting in an approximation for finite PCE, as The PCE coefficients p f j is determined by non-intrusive spectral projection of (5) onto the orthogonal polynomial basis tΦ j u as Table 1 describes the relationship between the types of Wiener-Askey polynomial chaos and their corresponding underlying random variables.For uniformly distributed variables ξ used in this paper, the Wiener-Askey scheme Xiu and Karniadakis (2002) requires Legendre polynomials as the polynomial basis tΦ j u.

Stochastic collocation
Sparse grid methods Novak and Ritter (1996, 1997, 1999); Barthelmann et al. (2000) are a cornerstone in high-dimensional interpolation and integration that have been used in a variety of disciplines.In concert with the generalized polynomial chaos expansion Xiu and Karniadakis (2002) as a non-intrusive spectral projection approach, SC methods Babuška et al. (2007); Nobile et al. (2008); Xiu (2009) are developed to improve the efficiency of the generalized polynomial chaos expansion on high-dimensional problems using Smolyak sparse grids for integration.In the nutshell, the polynomial chaos expansion coefficients are computed based on the sparse grid framework that significantly reduces the effect of the curse-of-dimensionality.
Following Nobile et al. (2008), let i ě 1 and ξ i 1 , . . ., ξ i m i ( Ă r´1, 1s be a sequence of abscissas, we begin by introducing the one-dimensional Lagrange interpolation operator as where L i j pξq are the Lagrange polynomials of degree m i ´1, i.e.L i j pξq " . The full tensor product formula is perhaps the most straightforward, as which requires m i functions evaluations.Although simple, a major drawback of full tensor product is that the total number of points grows very fast in high dimensions.Numerous choices of collocation points are possible, such as Gauss-Legendre, Clenshaw-Curtis, Leja, and Gauss-Patterson Nobile et al. (2016).In this work, the weakly-nested Gaussian abscissas (see Section 3.6.2 of Dalbey et al. (2021) and Eldred and Burkardt (2009)), which are zeros of orthogonal polynomials, are utilized for quadrature.
Let U 0 " 0 and for i ě 1, define the isotropic Smolyak quadrature formula Smolyak (1963) is given by or equivalently as, where q ě n is an integer denoting the level of the construction Wasilkowski and Woźniakowski (1995).To compute the operator Apq, nq, one needs to evaluate f on the set of points Hpq, nq " where ξ i " tξ i 1 , . . ., ξ i m i u Ă r´1, 1s is the collection of abscissas used by the univariate interpolating operator U i .This set is a much smaller subset of those required by the full tensor product rule.If the sets are nested, i.e. ξ i Ă ξ i`1 , then Hpq, nq Ă Hpq `1, nq.
To illustrate the benefits in using the Smolyak sparse grid, compared to the full tensor grid, Table 2 compares the number of simulations required to achieve the same level of accuracy.The dimensionalities are chosen according to the case studies in this paper.The equivalent number points on full tensor grid point is computed as p2 p `1q ´1q n , where is the corresponding level of sparse grid.

Variance-based global sensitivity analysis
Following Sudret (2008), Crestaux et al. (2009), andSaltelli et al. (2010), we summarize the variance-based global sensitivity analysis based on Sobol' decomposition as follows.In the spirit of generalized polynomial chaos expansion (i.e.Equation 1 after finite truncation), the Sobol' decomposition of f pξq into the summands of increasing dimensions as Given a model of the form y " f pξ 1 , ξ 2 , . . ., ξ n q, with y as a scalar, a variance-based first order effect for a generic factor ξ i can be written as , where ξ "i is the vector ξ without the i-th element, i.e. ξ "i " pξ 1 , . . ., ξ i´1 , ξ i`1 , . . ., ξ n q.The main effect sensitivity index (first-order sensitivity coefficient) is written as It is relatively well-known that and therefore, the total effect sensitivity index can be obtained as In global SA, the importance of parameter ξ i is measured by comparing its variance of the conditional expectation V ξ i " E ξ "i ry|ξ i s ‰ against the total variance Vrys.S i measures the effect of ξ i by evaluating the variance contribution of the basis function p f i that depends strictly on the set of

UQ for CPFEM constitutive models
variables in ξ i , while T i measures the total effect of ξ i by evaluating the variance contribution of all basis function whose dependencies include ξ i .For mathematical and implementation details, interested readers are referred to Tang et al. (2010), Weirs et al. (2012), andCrestaux et al. (2009), where most of the computations are based on Monte Carlo sampling ξ.In the context of this manuscript, we can understand ξ as the set of parameters for the underlying constitutive model, whether it is phenomenological or dislocation-density-based, and y as the quantity of interest from the CPFEM model.It should be emphasized that the method used in this study is the mainstream global sensitivity analysis that is widely used in structural reliability studies.

UNCERTAINTY QUANTIFICATION WORKFLOW FOR CRYSTAL PLASTICITY
In this paper, we limit the scope of the UQ studies to cases with a unique RVE.As discussed in Section 7, relaxing this restriction will be the subject of future work.DREAM.3DGroeber and Jackson (2014)  Figure 1 describes the integrated framework coupling DAKOTA uncertainty quantification code and DREAM.3D+ DAMASK workflow.Based on the sparse grid construction specifications, such as anisotropic / isotropic, sparse grid level, as well as other sensitivity analysis options, DAKOTA sets up a list of input parameters to be determined and evaluated by the coupled DREAM.3D+ DAMASK workflow.The sets of simulations are then deployed on high-performance computing systems and evaluated in parallel to accelerate the process.Typically, for a fixed input parameter vector, an ensemble of microstructural RVEs are used; however, to reduce the computational cost in this study, we limit the scope of our investigation to one microstructure RVE.It should be noted that, if the initial microstructure is fixed, then DAKOTA would interact directly with DAMASK, and the role of DREAM.3D can be conveniently ignored.In this section, Section 3.1 briefly describes the fundamentals of CPFEM model, whereas Section 3.2 describes the UQ workflow for CPFEM based on DAKOTA UQ package as the wrapper and DREAM.3D and DAMASK as a forward ICME model.

Crystal plasticity model
Consider each point X in the reference configuration being mapped to the current configuration x by a linear transformation with the deformation gradient tensor F, where F " Bx BX .The Lagrangian strain tensor is defined as The total deformation gradient F can be multiplicatively decomposed into an elastic and plastic parts, The velocity gradient, which measures the deformation rate, is defined as where L p is the plastic velocity gradient evaluated in the intermediate configuration.The second Piola-Kirchoff stress measure S is defined as where E e is the elastic Green-Lagrange strain tensor, and C is the fourth-order stiffness tensor.The plasticity velocity gradient L p , driven by the second Piola-Kirchoff S, controls the evolution of the plastic deformation gradient as 9 Constitutive equations representing the flow stress, such as the phenomenological slip-based hardening model and the dislocation density-based hardening model, differ on L p are calculated based on a specific microstructure and a set of internal state variables.The grain size d in the RVE is a random variable characterized by a log-normal distribution, i.e., where µ D and σ D are materials-dependent parameters.

Forward uncertainty quantification problem
In this work, we consider a forward UQ problem for a fixed set (or ensemble) of microstructures using the SC method, where the set of internal state variables for constitutive model parameters are considered stochastic with some inherent uncertainty.While the non-intrusive PCE approach allows an arbitrary probability density of parameters, in this work, uniform distributions are imposed on the constitutive model parameters due to lack of prior knowledge.With the choice of uniform distributions on bounded intervals, according to the Wiener-Askey scheme shown in Table 1, Legendre polynomials are used to approximate the quantities of interest.
In the initial yield regime, we focus on the the estimated yield strain ε Y and yield stress σ Y .From the homogenized stress-strain curve ε vM ´σvM obtained from a CPFEM simulation, the monotonic cubic interpolation via PCHIP method Fritsch and Butland (1984) is utilized to interpolate ε vM ´σvM curve.From the approximated stress-strain curve, an estimation of modulus of elasticity is obtained by simple linear regression.An offset at ε " 0.002 with the estimated modulus of elasticity is drawn, where the coordinates of the intersection are pε Y , σ Y q.Statistics of these two quantities of interest are obtained and returned to DAKOTA package.
To set up the UQ study, a pre-processing compilation of constitutive model parameters are obtained from DAKOTA using a numerical toy model.The sets of constitutive model parameters are then appropriately parsed into DAMASK using Python scripts, along with the output of DREAM.3D for setting up the geometric file.With the correct setup of DAMASK simulations, the set of DAMSASK simulations are then performed in a massively parallel manner on a high-performance This is a provisional file, not the final typeset article computing cluster.The post-processing results are then collected from DAMASK, and parsed back to DAKOTA package using a Python interface.DAKOTA then performs the UQ and sensitivity analysis, concluding the UQ workflow for CPFEM based on DREAM.3D and DAMASK.

Constitutive law
We adopt the summary and tabulated parameters from Sedighiani et al. (2020Sedighiani et al. ( , 2022) ) (Tables 1 and  2).In the phenomenological constitutive model, the shear on each slip system α is modeled as where τ 0 is the slip resistance, 9 γ 0 is the reference shear rate, and n determines the strain rate sensitivity of slip.The influence of other slip system α 1 on the hardening behavior of the slip system α is modeled as where h αα 1 is the hardening matrix, which captures the micromechanical interaction among different slip systems h 0 , a, and τ 8 are slip hardening parameters for all 12 slip systems in fcc materials.q αα 1 is a measure for latent hardening with value of 1.0 for coplanar slip system α and α 1 and 1.4 otherwise.

Design of numerical experiments
We conducted our analysis with one RVE shown in Figure 3a, where the crystallographic texture is shown in Figure 3b.Average grain size of 60.9467µm is used, where an RVE of size 256µm 3 is generated.A finite element mesh of 16 3 is created to approximate the microstructure RVE.In DREAM.3D, the texture crystallography of Copper-type with Euler angles pφ 1 , θ, φ 2 q " p90 ˝, 35 ˝, 45 ˝q is used, the grain size parameters are set as µ D " 4.09, σ D " 0.2, which results   in a RVE with 182 grains, shown in Figure 3a.For DAKOTA, we set the sparse grid level " 3, dimensionality n " 5, which results in 351 inputs.For each set of input parameters, a CPFEM simulation is performed, followed by the post-process.The results are analyzed in the following section.A uniaxial loading condition is applied in the [100] direction with 9 ε " 10 ´3s ´1.
To compare with the default parameter, a single CPFEM simulation is performed with constitutive parameters described in Table 3.The stress-strain equivalent curve is shown in Figure 4.As a reference to experimental data 1 , modulus of elasticity for polycrystalline Cu is reported at 110 GPa, whereas its yield strength is reported as 33.3 MPa.Compared to the experiment, the computed 1 https://www.matweb.com/search/DataSheet.aspx?MatGUID=9aebe83845c04c1db5126fada6f76f7e This is a provisional file, not the final typeset article

UQ for CPFEM constitutive models
modulus and yield strength σ Y are on the same scale, but off roughly by a factor of 1.5ˆ, possibly due to different processing conditions resulting in slightly different alloys.

Numerical results
Figure 5 shows a compilation of stress-strain curve for 351 simulations, where each corresponds to a unique set of constitutive parameters for fcc Cu.As shown in this figure, the constitutive model effects not only the initial yield behavior, but also the modulus of elasticity and the hardening behavior.Figures 6a and 6b show the probability density function for ε Y and σ Y , respectively.The mode for ε Y is approximately 0.00207, whereas the mode for σ Y is approximately 11.65 MPa.It should be noted that with the current sparse grid level " 3, the approximation for ε Y may be imprecise.One of the possible reasons is that the elastic regime of copper is very small (as shown in Figure 5), and therefore, a more accurate approximation may be required to accurately capture the yield strain.
Figure 7 and Figure 8, respectively, show the Sobol' indices for ε Y and σ Y .Ranking from the most influential parameters to the least influential parameters for ε Y from the Sobol indices for main effects, T τ 0 " 0.7858, T h 0 " 0.7035, T n " 0.1922, T a " ´0.01038, and T τ8 " ´0.001983.Ranking from the most influential parameters to the least influential parameters for σ Y from the Sobol indices for main effects, T τ 0 " 0.8258, T h 0 " 0.4194, T n " 0.07659, T a " ´0.01091, and T τ8 " ´0.0009990.The order of influential parameters for fcc Cu, regarding the initial yield behavior, is τ 0 ą h 0 ą n ą a ą τ 8 .Since the main scope of this paper is about the initial yield behavior, it is not surprising that Figure 8

Constitutive law
Firstly introduced by Hutchinson Hutchinson (1976) and extended for twinning by Kalidindi Kalidindi (1998), the resistance on α " 1, . . ., N s slip systems evolve from ξ 0 to a system-dependent saturation value and depend on shear on slip and twin systems according to ) where h denotes the components of the slip-slip and slip-twin interaction matrices, h s-s 0 , h int , c 1 , c 2 are model-specific fitting parameters and ξ 8 represents the saturated resistance evolution.
The resistances on the β " 1, . . ., N tw twin systems evolve in a similar way, where h tw-s 0 , h tw-tw 0 , c 3 , and c 4 are model-specific fitting parameters.Shear on each slip system evolves at a rate of where slip due to mechanical twinning accounting for the unidirectional character of twin formation is computed slightly differently, 9 γ " p1 ´f tot tw q 9 γ 0 ˇˇˇτ ξ ˇˇˇn Hpτ q, (29)

Frontiers
where H is the Heaviside step function.The total twin volume is calculated as where γ char is the characteristic shear due to mechanical twinning and depends on the twin system.Interested readers are referred to Section 6.2.2 from Roters et al Roters et al. (2019).

Design of numerical experiments
Similar to the previous section we restrict the scope to one RVE shown in Figure 9a, where the crystallographic texture is shown in Figure 9b.Average grain size of 204.037µm is used, where an RVE of size 2048µm 3 is generated.A mesh of 64 3 is created to approximate the microstructure RVE.In DREAM.3D, the texture crystallography with Euler angles pφ 1 , θ, φ 2 q " p90 ˝, 0 ˝, 0 ˝q is used, the grain size parameters are set as µ D " 5.2983, σ D " 0.2, which results in a RVE with 1706 grains, shown in Figure 9a.For DAKOTA, we set the sparse grid level " 2, dimensionality n " 16, which results in 577 inputs.For each set of input parameters, a CPFEM simulation is performed, followed by the post-processing steps.The results are analyzed in the following section.A uniaxial loading condition is applied in the [100] direction with 9 ε " 10 ´3s ´1.
To compare with the default parameter, a single CPFEM simulation is performed with constitutive parameters described in Table 4.The stress-strain equivalent curve is shown in Figure 10.As a   reference to experimental data2 , the modulus of elasticity for polycrystalline Mg is reported at 44 GPa, whereas its yield strength is reported as 90-105 MPa.Compared to the experiment, the computed modulus and yield strength σ Y are well calibrated, as the computational results agree very well with the experimental data.

Numerical results
Figure 11 shows a compilation of stress-strain curves for 577 simulations, where each corresponds to a unique set of constitutive parameters for hcp Mg.As shown in this figure, the constitutive model has a minor effect on the effective modulus of elasticity, and more profound effect on the yield stress σ Y .Figure 13 and Figure 14, respectively, show the Sobol' indices for ε Y and σ Y .Ranking from the most influential parameters to the least influential parameters for ε Y from the Sobol indices for main effects, T τ 0,basal " 0.5668, T τ 0,C2 " 0.4772, T ntw " 0.2439, T h s´s 0 " 0.1021, T τ 8,basal " 0.07249, T τ 0,pyrxay " 0.6131, T ns " 0.02082, T τ 8,pyrxay " 0.01091.Ranking from the most influential parameters to the least influential parameters for σ Y from the Sobol indices for main effects, T τ 0,C2 " 0.3729 T ntw " 0.3684, T τ 0,basal " 0.3566, T τ 0,pyrxay " 0.1181, T τ 8,basal " 0.1064, T h s´s 0 " 0.1061, T τ 0,pris " 0.03861.Compared to Sedighiani et al. Sedighiani et al. (2020, 2022), our analysis shows some agreements, but mostly differ in the set of sensitive parameters.Possible explanations are due to (1) different quantities of interest and (2) methodological approach: Sedighiani et al.Sedighiani et al. (2020Sedighiani et al. ( , 2022) ) studies are conducted based on ANOVA, whereas our approach relies on global sensitivity analysis with Sobol' indices.

Constitutive law
For the sake of completeness, we adopt the dislocation-density-based constitutive law description from Cereceda et al.Cereceda et al. (2016Cereceda et al. ( , 2015)), Stukowski et al.Stukowski et al. (2015), and summarize it here.Interested readers are further referred to Cereceda et al. Cereceda et al. (2016, Figure 12a.SC probability density function of ε Y for hcp Mg.where α is a slip system, m α and n α are unit vectors in the normalized slip direction and the normal to the slip plane of the system α, respectively, P α S " m α b n α is a Schmid geometric projection tensor.The resolved shear stress of slip system α include both Schmid and non-Schmid factors as , where P α T/AT " a 1 m α b n α 1 is a non-Schmid tensor representing twinning and anti-twinning asymmetry and the effects due to non-glide stress components, P α ng " a 2 pn α b m α q b α à3 pn α 1 b m α q b n α 1 .n α 1 forms an angle of ´60 ˝with the reference slip plane defined by n α and changes sign with the direction of slip on each glide plane Koester et al. (2012).a 1 , a 2 , a 3 are calibrated material-dependent constants.The shear rate 9 γ α is given by the Orowan equation where b " a 0 ?3{2 is the magnitude of the Burgers vector, a 0 is the lattice parameter, T is the absolute temperature, ρ α is the density of mobile screw dislocations in slip system α, and v s pτ α , T q is the screw dislocation velocity, which captures the thermally activated character of dislocation motion.
Under the assumption that kink relaxation is significantly faster than kink-pair nucleation, the total time t t required for a kink pair to form and sweep a rectilinear screw dislocation segment of length λ α is t t " t n `tk " Jpτ α , T q ´1 `λα ´w Tran et al.

UQ for CPFEM constitutive models
where t n is the mean time to nucleate a kink pair, t k is the time needed for a kink to sweep half a segment length, J is the kink-pair nucleation rate, w is the kink-pair separation, v k is the kink velocity.
The kink-pair nucleation rate is modeled by an Arrhenius formulation as where v 0 is an attempt frequency, ∆H kp is the activation enthalpy of a kink pair stress τ α , k is Boltzmann's constant.The kink velocity is modeled as where B is simplified to a constant.The dislocation velocity can be modeled as where h " a 0 ?6{3 is the distance between two consecutive Peierls valleys.When t k !t n , v 0 Bpλ α ´wq 2 exp ˆ´∆H kp kT ˙Ñ 0, and the common diffusive velocity expression is recovered It is further elaborated in Sedighiani et al. (2020) that where p and q determine the shape of the short-range activation energy.
Following the Kocks-Mecking family of dislocation density evolution models Mecking and Kocks (1981), the mobile dislocation density on slip system α evolves in time is modeled as with the initial dislocation density ρ α pt " 0q " ρ α 0 .Dislocation multiplication is proportional to the inverse mean free path of the dislocation λ α and the plastic strain rate, as λ α is defined as where d g is the grain size, c is a hardening constant, ρ α f is the forest dislocation density and calculated as ρ α f " Dislocation annihilation occurs spontaneously when dipoles approach within a spacing of d edge The resolved shear stress τ α is corrected by replacing with τ α 1 and accounting for the latent and self-hardening as where τ h is the hardening stress, ξ αα 1 are the coefficients of the interaction matrix between slip system α and α 1 , respectively, as six possible independent interactions: self, coplanar, collinear, orthogonal, glissile, and sessile.

Design of numerical experiments
For dislocation-density-based CPFEM simulations, we used an RVE shown in Figure 15a, with the crystallographic texture shown in Figure 15b and Euler angles pφ 1 , θ, φ 2 q " p0, 0, 45q.Average grain size of 3.0339µm is used, where an RVE of size 32µm 3 is generated.A mesh of 32 3 is created to approximate the microstructure RVE.In DREAM.3D, the texture crystallography with Euler angles pφ 1 , θ, φ 2 q " p0 ˝, 0 ˝, 45 ˝q is used, the grain size parameters are set as µ D " 1.0986, σ D " 0.15, which results in a RVE with 2089 grains, shown in Figure 15a.For DAKOTA, we set the sparse grid level " 3, dimensionality n " 17, which results in 799 inputs.For each set of input parameters, a CPFEM simulation is performed, followed by the post-process.The results are analyzed in the following section.A uniaxial loading condition is applied in the [100] direction with 9 ε " 10 ´3s ´1.
To compare with the default parameter, a single CPFEM simulation is performed with constitutive parameters described in Table 5.The stress-strain equivalent curve is shown in Figure 16.As a   reference to experimental data3 , the modulus of elasticity for polycrystalline W is reported at 400 GPa, whereas its yield strength is reported as approximately 750 MPa.Compared to the experiment, the computed modulus agrees very well, but the yield strength σ Y differs, possibly due to the constitutive model is calibrated for single crystal W, whereas the experimental data is reported for polycrystalline W.

Numerical results
Figure 17 shows a compilation of stress-strain curves for 799 simulations, where each corresponds to a unique set of constitutive parameters for bcc W. As shown in this figure, the constitutive model has a minor effect on the effective modulus of elasticity, and more profound effect on the yield stress σ Y .Figure 19 and Figure 20, respectively, show the Sobol' indices for ε Y and σ Y .Ranking from the most influential parameters to the least influential parameters for ε Y from the Sobol indices for main effects, T p " 0.7356, T ∆H 0 " 0.212, T C λ " 0.2089, T q " 0.1942, T ρ α 0 " 0.1556, T τ Peierls " 0.05712, T ν 0 " ´0.03444.Ranking from the most influential parameters to the least influential parameters for σ Y from the Sobol indices for main effects, T p " 0.5668, T ∆H 0 " 0.2095, T C λ " 0.2079, T q " 0.1795, T ρ α 0 " 0.1560, T τ Peierls " 0.06172, T ν 0 " ´0.03663.Since the ranking does not change when the quantity of interest changes from ε Y to σ Y , we conclude that the importance of parameters is p ą ∆H 0 ą C λ ą q ą ρ α 0 ą τ Peierls ą ν 0 .

DISCUSSION
In this paper, we conducted several UQ studies for constitutive models in CPFEM with a single microstructure RVE for each case study.Three case studies are performed with different crystal structures, namely fcc, hcp, and bcc, for Cu, Mg, and W, respectively.In this paper, three materials systems with different crystal structures (fcc, bcc, and hcp) are studied.Depending on the crystal structure, there may be different slipping and twinning systems in terms of slipping and twinning directions in plastic deformation, leading to interesting materials behaviors and mechanisms.The    UQ such as those described in this manuscript, play an important role in constitutive model calibration for unknown material system in the future.Since there are only a limited number of physical constitutive models, it is important to conduct a UQ study to observe the range of quantities of interest, and to numerically rank the influence of constitutive model parameters.Based on the stress-strain compilation curve conducted for various constitutive model parameters, the material behaviors can be rigorously quantified.The obtained UQ results provide a foundational step for further constitutive model calibration for future works, mostly conducted via digital image correlation techniques Turner et al. (2015); Reu et al. (2018Reu et al. ( , 2021)).
Compared to polynomial approximation with full tensor grid, sparse grid approaches have a significant computational advantages, where this advantage grows with increasing dimensionality thanks to a slower growth rate Nobile et al. (2008).In the context of constitutive models, the computational reduction is mostly profound in the case of hcp system (as opposed to bcc and fcc), such as Mg and Ti, and in the case of dislocation-density-based constitutive model (as opposed to phenomenological model), where many parameters require careful calibration to obtain a sufficient agreement with experiments.For simple system with a relatively simple phenomenological constitutive model, the computational reduction is less severe.It is noteworthy that the level of the Smolyak sparse grid in this study has a little effect on the resulting probability density function of QoIs.This implies that the underlying function is perhaps mostly low-order.This observation can also be confirmed by the Sobol' indices, where the first-order Sobol' indices are much more dominant, compared to higher-order Sobol' indices.
To construct the response surface model, stochastic collocation provides a significant advantage for reducing the curse of dimensionality.However, when it comes to accuracy, Gaussian process regression, which is also the underlying surrogate model for Bayesian optimization, is arguably one of the best approaches in shallow machine learning.The direction of Bayesian optimization, e.g.Tran et al. (2019aTran et al. ( , 2020Tran et al. ( , 2021)); Tran (2019); Tran et al. (2022), for robust constitutive model calibration remains open for future research.
The scope of this manuscript is to quantify the microstructure-sensitive uncertainty.Obviously, it can be expanded to account for the entire stress-strain curve.However, due to the number of parameters involved in each constitutive model, there are hundreds to thousands of runs needed for a single microstructure RVE.Such computationally expensive numerical experiments require careful planning and execution, and therefore, remain a potential topic for future studies.It is important to point out that by restricting to one RVE per case study, this work does not address microstructure-sensitive uncertainty that either is related or induced by the underlying stochastic nature of microstructures.The direction of investigating a microstructure ensemble with many RVEs remain open for future work.

CONCLUSION
In this paper, we applied SC to quantify uncertainty associated with the initial yield behavior, mainly the estimated yield strain ε Y and the estimated yield stress σ Y for fcc Cu, hcp Mg, and bcc W. A variety of constitutive models are used, resulting with different parameterization and dimensionalities for the constitutive models considered.To mitigate the curse of dimensionality, Smolyak sparse grid is employed for high-dimensional integration to evaluate the PCE coefficients.Variance-based global sensitivity analysis is used to study the sensitivity analysis of the constitutive model parameters.
In light of the computational results presented in previous sections, there are several influential parameters that may have a significant effect on the initial yield behavior.For the phenomenological constitutive model, the slip resistance τ 0 , the slip hardening parameter h 0 , and the strain rate sensitivity parameter n are the most influential parameters, ranking in descending order.For the dislocation-density-based constitutive model, the p-exponent in glide velocity is the most influential parameter, followed by the activation energy for dislocation glide ∆H 0 , the dislocation mean free path parameter C λ , the q-exponent in glide velocity, and the initial dislocation density ρ α 0 .We conclude that in both constitutive models considered in this study, i.e. phenomenological (with and without twinning) and dislocation-density-based constitutive models, regarding the initial yield behavior, some parameters may have a profound effect on the QoI, while some others may not have a significant effect.The observation could potentially pave way for dimensionality reduction in constitutive model calibration in the future.
Figure 3a.Representative volume element for Cu.

Figure 4 .
Figure 4. Stress-strain equivalent curve for Cu with default parameters and determination of yield point.Modulus of elasticity is estimated as 188.5919GPa, while pε Y , σ Y q are estimated as (0.002052, 9.8527 MPa)
Figure7and Figure8, respectively, show the Sobol' indices for ε Y and σ Y .Ranking from the most influential parameters to the least influential parameters for ε Y from the Sobol indices for main effects, T τ 0 " 0.7858, T h 0 " 0.7035, T n " 0.1922, T a " ´0.01038, and T τ8 " ´0.001983.Ranking from the most influential parameters to the least influential parameters for σ Y from the Sobol indices for main effects, T τ 0 " 0.8258, T h 0 " 0.4194, T n " 0.07659, T a " ´0.01091, and T τ8 " ´0.0009990.The order of influential parameters for fcc Cu, regarding the initial yield behavior, is τ 0 ą h 0 ą n ą a ą τ 8 .Since the main scope of this paper is about the initial yield behavior, it is not surprising that Figure8agrees with Figure7in terms of Sobol' indices.

Figure 6a .
Figure 6a.SC probability density function of ε Y for fcc Cu.Figure 6b.SC probability density function of σ Y for fcc Cu.
Figure 6a.SC probability density function of ε Y for fcc Cu.Figure 6b.SC probability density function of σ Y for fcc Cu.
Figure 9a.Representative volume element for Mg.

Figure 10 .
Figure 10.Stress-strain equivalent curve for Mg with default parameters and determination of yield point.Modulus of elasticity is estimated as 45.1172 GPa, while pε Y , σ Y q are estimated as (0.004375, 107.2 MPa)

Figure
Figure 12a and Figure 12b, respectively, show the probability density function for ε Y and σ Y .The mode for ε Y is approximately 0.0054, whereas the mode for σ Y is approximately 99 MPa.The uncertainty explained in σ Y reasonably agree with experimental data.
Figure 12a.SC probability density function of ε Y for hcp Mg.Figure 12b.SC probability density function of σ Y for hcp Mg.
Figure 15a.Representative volume element for W.

Figure 16 .
Figure16.Stress-strain equivalent curve for W with default parameters and determination of yield point.Modulus of elasticity is estimated as 411.3528GPa, while pε Y , σ Y q are estimated as (0.01098, 3695 MPa).

Figures
Figures 18a andFigure18b shows the probability density function for ε Y and σ Y , respectively.The mode for ε Y is approximately 0.01010, whereas the mode for σ Y is approximately 3350 MPa.The uncertainty explained in σ Y reasonably agree with experimental data.For experimental data with σ Y « 750 MPa, fairly a few set of parameters can reproduce and calibrate accordingly, as shown in Figure 17.

Figure 18a .
Figure 18a.SC probability density function ε Y for bcc W.

Figure 18b .
Figure 18b.SC probability density function of σ Y for bcc W.

Table 1 .
Relationship between the types of Wiener-Askey polynomial chaos and their underlying

Table 2 .
The number of collocation points used by sparse grid and full tensor grid.
Adams et al. (2009) polycrystalline RVE with a specific crystallographic texture, depending on the material considered.DAKOTAAdams et al. (2009)and Python scripts are used to generate inputs for constitutive models, whereDAMASK Roters et al. (2019)is employed as the CPFEM forward model.Results are collected and post-processed in DAKOTA.
‚ microstructure generation ‚ crystallographic texture ‚ crystal plasticity finite element Figure 1.Integrating DAKOTA uncertainty quantification workflow to DREAM.3D and DAMASK crystal plasticity finite element simulations.In this framework, DAKOTA queries input parameters to the DREAM.3D+ DAMASK automatic workflow, and receives output(s)/quantity(quantities) of interest from DREAM.3D + DAMASK.