Uncertainty-quantified parametrically homogenized constitutive models (UQ-PHCMs) for dual-phase α/β titanium alloys

This paper develops an uncertainty-quantified parametrically homogenized constitutive model (UQ-PHCM) for dual-phase α/β titanium alloys such as Ti6242S. Their microstructures are characterized by primary α-grains consisting of hcp crystals and transformed β-grains consisting of alternating laths of α (hcp) and β (bcc) phases. The PHCMs bridge length-scales through explicit microstructural representation in structure-scale constitutive models. The forms of equations are chosen to reflect fundamental deformation characteristics such as anisotropy, length-scale dependent flow stresses, tension-compression asymmetry, strain-rate dependency, and cyclic hardening under reversed loading conditions. Constitutive coefficients are functions of representative aggregated microstructural parameters or RAMPs that represent distributions of crystallographic orientation and morphology. The functional forms are determined by machine learning tools operating on a data-set generated by crystal plasticity FE analysis. For the dual phase alloys, an equivalent PHCM is developed from a weighted averaging rule to obtain the equivalent material response from individual PHCM responses of primary α and transformed β phases. The PHCMs are readily incorporated in FE codes like ABAQUS through user-defined material modeling windows such as UMAT. Significantly reduced number of solution variables in the PHCM simulations compared to micromechanical models, make them several orders of magnitude more efficient, but with comparable accuracy. Bayesian inference along with a Taylor-expansion based uncertainty propagation method is employed to quantify and propagate different uncertainties in PHCM such as model reduction error, data sparsity error and microstructural uncertainty. Numerical examples demonstrate the accuracy of PHCM and the relative importance of different sources of uncertainty.


INTRODUCTION
Structural analysis of heterogeneous materials using phenomenological constitutive models, is often faced with inaccuracies stemming from the lack of connection with the material microstructure and underlying physics. Pure micromechanical analysis, on the other hand, is computationally prohibitive on account of the large degrees of freedom needed to represent the entire structure. To overcome these shortcomings, hierarchical multiscale models based on computational homogenization, have been proposed to determine the homogenized material response for heterogeneous materials that can be used in component-scale analysis. A detailed account of these approaches is given in ref. 1 . However, for nonlinear problems involving history-dependent constitutive relations, many multiscale methods incur prohibitive computational costs from solving the micromechanical problem for every macroscopic point in the computational domain. It is imperative to develop robust macro-scale constitutive models, incorporating characteristic microstructural features as well as underlying physical mechanisms of deformation.
The concept of parametric homogenization has been introduced in refs. [2][3][4][5] to: (i) overcome prohibitive computational costs of other homogenization methods, and (ii) bridge length-scales through explicit microstructural representation in structure-scale constitutive models. PHCMs for polycrystalline α-phase titanium alloys have been developed in refs. 1,6,7 . The physics-informed PHCMs differ from conventional phenomenological models in their unambiguous depiction of constitutive parameters and their dependencies. Coefficients in these thermodynamically consistent, reduced-order constitutive models are functions of representative aggregated microstructural parameters (RAMPs), corresponding to distributions of key morphological and crystallographic descriptors in the microstructure. The forms of equations in PHCMs are chosen to reflect fundamental deformation characteristics of the material, made in micromechanical observations. For the Ti alloys, characteristics like objectivity, anisotropy, tension-compression asymmetry, and history/path-dependence are represented in the micromechanical crystal plasticity FE (CPFE) models. Constitutive parameters in PHCMs are calibrated from CPFE analysis, and their functional forms in terms of RAMPs are determined using machine learning tools. Details of computational development and validation of the PHCMs for single phase near-α Ti alloys is given in refs. 1,6 . The present paper extends the PHCMs to two-phase polycrystalline α/β Titanium alloys like Ti-6Al-4V and Ti-6Al-2Sn-4Zr-2Mo that are widely used in aircraft engine components due to their superior mechanical, thermal and fatigue properties 8 .
Microstructures of dual phase Ti alloys are typically characterized by primary α and transformed β phases, the latter consisting of colonies of alternating laths of secondary α (hcp) and β (bcc) lamallae in a matrix of equiaxed primary α grains as shown in Fig. 1a. The anisotropic phases, colony structures and crystallographic orientation mismatch at grain and lath boundaries significantly influence their deformation, creep, and fatigue characteristics [9][10][11][12][13][14] . Under dwell loading conditions, the near-α and α/β alloys undergo time-dependent load shedding at boundaries between grains with moderate to large lattice misorientation, leading to facet formation signaling the nucleation of fatigue cracks [15][16][17] .
Nucleation and short crack growth in these alloys is strongly influenced by local microstructural features such as crystallography and size of micro-textured regions 9,16,18,19 . These observations call for material microstructure-based modeling to account for morphology and crystallography in the reliable prediction of their mechanical properties and fatigue behavior. CPFE models have been extensively used 9,13,20 to characterize deformation behavior associated with α and α + β alloys. However, component-scale analysis with these micromechanical models are computationally expensive for macroscopic analysis. The PHCMs in refs. 1,6,7 enable full-scale component analysis with explicit representation of microstructural features and characteristics. They enjoy high accuracy, while being many orders of magnitude more efficient than conventional multiscale analysis methods. In this paper, the PHCMs are extended to two-phase polycrystalline α/β titanium alloys like Ti6242S, accounting for the presence of the α and β phases in the granular structure. In addition, the PHCM constitutive relations incorporate a kinematic hardening law to account for the Bauschinger effect, an important aspect of the mechanical response under reversed cyclic loading.
In addition, uncertainty quantification and propagation comprise important tasks in the PHCM development process to account for different sources of uncertainties, as depicted in Fig. 2. Uncertainties in constitutive coefficients are attributed to variability in the material microstructure, microstructural characterization, calibration of CPFE model, calibration of PHCM, etc. For example, errors inherently occur in microstructural characterization of electron back-scatter diffraction (EBSD) images, used for the generation of microstructure-based statistically equivalent representative volume elements or M-SERVEs for polycrystalline microstructures. Generation methods of the M-SERVEs for different alloys have been developed in refs. [21][22][23] , where the statistical distributions of microstructural descriptors are matched with experimental data obtained from EBSD or scanning electron microscopy (SEM) images. This process also induces uncertainty with respect to matching the higher order statistics of distributions. Consideration of all uncertainties in the multi-scale framework is a very detailed process. Three sources of uncertainty in PHCM constitutive coefficients, considered in this paper (see Fig. 2), are: • Model reduction error due to the functional forms used to represent the microstructure-dependent constitutive parameters using machine learning tools; • Data-sparsity error due to the finite size of the calibration data obtained from a finite number of CPFE simulations used in PHCM calibration; • Uncertainty in the microstructural descriptors or RAMPs due to natural variability in the microstructure.
For predictive analysis of material response, the uncertainties in PHCMs must be adequately quantified and propagated with additional deformation. Bayesian inference, along with maximum likelihood estimation, has been used in ref. 14 to quantify (i) the uncertainty associated with constitutive parameters in PHCM functional forms obtained from machine learning, and (ii) the uncertainty due to data-sparsity associated with a limited number of microstructures and micromechanical simulations used in calibrating the PHCMs. Furthermore, with evolving deformation, the uncertainty in PHCM constitutive parameters is propagated to response variables like stresses and plastic strains. Monte-Carlo based uncertainty propagation is not suitable since it requires repeated evaluations of the PHCMs. Multiple evaluations are generally not desirable for macroscopic analysis under given loading conditions. A Taylor-series expansion-based uncertainty propagation method, originally developed in ref. 14 , is used to propagate the uncertainty to response variables. The resulting uncertainty-quantified PHCM has a unique advantage of being incorporated in any commercial software through user material windows without having to perform expensive Monte-Carlo simulations.
The paper is organized as follows. The accuracy of PHCM developed in this paper along with the effect of different sources of uncertainty on mechanical response of dual-phase titanium alloys is demonstrated through validation examples in "Results" section. The micromechanical crystal plasticity constitutive model and the framework to develop PHCMs for dual-phase titanium alloys is given in "Methods" section. Different sources of uncertainty present in the PHCM development are also discussed in this section.

Overview
The PHCM models, developed and calibrated for primary α and transformed β phases in "Methods" section, are used in the volume fraction based equivalent model in Eq. (15) to obtain the equivalent stress-strain response of M-SERVEs for a given volume fraction of transformed β grains. The PHCM constitutive equations in "Methods" section are implemented in the user subroutine UMAT of the commercial code Abaqus 24 . Numerical implementation and integration of the PHCMs in the UMAT have been described in details in ref. 6 . In this section, PHCM predictions are compared with micromechanical CPFE simulations under various loading conditions to demonstrate its accuracy and efficiency. This is followed by a numerical example that studies the effect of different uncertainties on the time-dependent stochastic mechanical response.
Validating PHCM with micromechanics-based CPFE simulations Experimental validation of the PHCMs for a near-α Ti6242S alloy containing~95% primary α phase has been successfully conducted in ref. 6   Therefore, the PHCMs for Ti6242S alloys containing primary α and transformed β phases are validated with micromechanical CPFE simulation results. The microstructure used in the validation process corresponds to the EBSD scan shown in Fig. 3a. It consists of 70% primary α and 30% transformed β grains by volume. The (0001) pole figure of the crystallographic texture for this microstructure is shown in Fig. 3b. The corresponding RAMPs derived from the figures are given in Table 1.
A 503 grain M-SERVE is constructed from the EBSD scan of Fig. 3a, for which the (0001) pole figure is depicted in Fig. 3c. This M-SERVE is meshed into 492733 linear tetrahedral elements. The PHCM predictions are compared with results of CPFE simulations for four different loading conditions. First, the CPFE model of the M-SERVE is subjected to different biaxial stresscontrolled loading conditions. The initial yield stresses along different directions are extracted following the procedure outlined in "Methods" section. The biaxial yield stresses are plotted in Fig.  4a, along with the PHCM yield surface projected on XZ plane. The latter is obtained using anisotropy coefficients calculated from functional forms in Box 1. Excellent agreement is found between the tensile and compressive yield stresses by CPFE simulations and the PHCM. Similarly, uniaxial yield stresses in tension and compression along different directions in each of the XY, XZ, and YZ plane are extracted from results of CPFE and PHCM simulation and compared in Fig. 4b. The difference between the PHCM and CPFE results for the uniaxial yield stresses and their tensioncompression asymmetry is quite small. These comparisons demonstrate the ability of the PHCM to accurately describe material anisotropy, as well as tension-compression asymmetry.
To assess the accuracy of post-yield predictions by PHCM, the M-SERVE is subjected to constant strain-rate loading along the X and Z directions, both in tension and compression. The homogenized stress-strain response from the CPFE and PHCM simulations is compared in Fig. 4c. The CPFE simulations of the M-SERVE take approximately 6 h with 24 CPUs for a true strain of 7%, while PHCM-based ABAQUS simulations of a single element, corresponding to the same microstructure, takes~1 s on single CPU for the same strain. It is observed that the anisotropy in flow stress is correctly predicted by PHCM both in tension and compression, with a small error in flow stress along the Zdirection under compression.
Finally, the PHCM performance is assessed for cyclic loading by subjecting the M-SERVE to strain-controlled, triangular wave and dwell loading as shown in Fig. 5b, d. The triangular cyclic loading profile has a peak strain of ϵ 0 = 1.2% with a ramp time of 1 s. The dwell loading profile corresponds to a ramp time of 1 s followed by a hold strain ϵ 0 = 1.2% for 120 s and unloading to zero strain in 1 s. These loadings are chosen to represent typical loading conditions used in fatigue testing of coupon specimens. CPFEMbased dwell simulation of the M-SERVE takes approximately 48 h with 24 CPUs, while the corresponding PHCM simulation takes about 3 s on a single CPU for the same number of dwell cycles. Similarly, the triangular-waveform cyclic simulation using CPFEM takes~35 h with 24 CPUs, while the same simulation using PHCM takes about 3 s on a single CPU. The maximum homogenized tensile and compressive stresses in each triangular cycle by CPFE and PHCM simulations are compared in Fig. 5a for 140 cycles. For dwell loading, the tensile stress at the end of the hold period and the compressive stress at the end of each dwell cycle, obtained by CPFE and PHCM simulations are plotted in Fig. 5c for 90 cycles. The PHCM underestimates the peak tensile stresses and overestimates the peak compressive stresses under both triangular cyclic and dwell loading for this M-SERVE. While the error in PHCM prediction is small under dwell loading, it is comparatively high along the Z-direction under triangular cyclic loading. These errors are observed as the applied peak strain lies in a region of the monotonic stress-strain response shown in Fig. 4c, where flow stresses in PHCM do not accurately match those from CPFE simulations. These errors may be attributed to a small mismatch in the elasto-plastic transition region during PHCM calibration. This can be improved through incorporating additional functions to better accommodate this transition in PHCM.
Effect of microstructural uncertainty on mechanical response of dual-phase Ti alloys The effect of various sources of uncertainty on the stochastic, time-dependent response of the dual-phase Ti-6242S alloy is studied in this section. The RAMPs for this alloy are given in Table 1.
To study their impact on the accuracy of predictions, the model reduction and data sparsity errors are first considered in stochastic PHCM simulations. These are introduced at the model development and calibration stages of the PHCM. A stochastic PHCM simulation is performed for a uniaxial tension test at a constant strain-rate of 1 × 10 −4 s −1 and the corresponding stress-strain response is depicted in Fig. 6a. The black solid line and the gray interval respectively correspond to the mean and the standard deviation of the Cauchy stress. The uncertainty due to model reduction and data sparsity error has a 1σ interval of~31 MPa corresponding to~3% of flow stress. This predicted modeling error interval is consistent with the deviation of the PHCM predictions from the CPFE model response, which is shown with markers in the figure. The predicted level of model uncertainty is also consistent with the set of PHCM validation tests in previous section.
Next, the effect of microstructural uncertainty on the mechanical response is analyzed. Microstructural variability, as shown in Fig. 1a, naturally arises during thermo-mechanical processing of the alloy. The microstructural uncertainty is represented by m Fig. 3 Details of validation microstructure. a Surface EBSD scan and b (0001) pole figure of the Ti6242S alloy containing 70% primary α and 30% transformed β phases by volume as given in ref. 13 ; c (0001) pole figure of the simulated 503 grain M-SERVE constructed from the EBSD scan showing material symmetry axes (v 1 : Red, v 2 : Blue, v 3 : Green).  statistical distributions of RAMPs, which are calculated from spatial variability of the microstructure in EBSD scans using a sampling method that has been discussed in ref. 7 . The resulting statistical moments of RAMPs from this sampling method are given in Table 2 and incorporated in the stochastic material evolution by the expansion-based UP method.
To understand the effect of the dual-phase α/β microstructure on a representative fatigue problem, the uncertainty-quantified PHCM is used to simulate a fully reversed cyclic loading with a strain amplitude of Δε = 2% and a cycle period of T = 40 s. Only microstructural uncertainty is considered in this analysis. Figure 6b shows the predicted evolution of the mean and 1σ interval of the stress-strain response, arising from microstructural uncertainty represented by the moments of RAMPs in Table 2. At the end of the first cycle, corresponding to the time marked with (*) in the figure, the uncertainty in the stress and plastic strain responses have standard deviation of 23.6 MPa and 1.62 × 10 −4 , respectively. This represents the variability of the mechanical fields in the material due to crystallographic and morphological variabilities in the microstructure. For example, the fatigue-driving stress at a critical point in a structural component like a bolt or a spot weld, may vary by this amount, depending on the particular microstructure of this critical volume of material that is subject to alleatoric uncertainty.
Multiple simulations are next performed, while individually turning on random variations in different features of the microstructure to identify the important microstructural features that control the uncertainty of material response. Table 3 shows the list of microstructural features and the corresponding stochastic RAMPs used in this analysis. For each class of microstructural uncertainty considered, the table shows the standard deviation of stress and plastic strain response (right) calculated at time (*) marked in Fig. 6b. The relative contribution of each source to the stochastic response is also given in the last column. The two constituents of the alloy, viz., the primary-α grains and the transformed-β colonies, contribute approximately  54% and 46% respectively to the total uncertainty in the response. This is shown in the last row of the Table 3 as "all sources" for each constituent. In the break-down of microstructural features for the primary-α grains, it is seen that the most important determinant of the stochastic response is crystallography, especially the crystallographic texture (~89%) for this alloy. This dominant effect mainly results from the micro-textured regions within the material, i.e., clusters of similarly oriented crystallographic regions. These are visualized in the EBSD scan of Fig. 1a and represented by the moments of OMA ij in Table 2. Due to anisotropy of the hcp Ti alloy, this variability results in large uncertainty in the stress response. The presence of a high degree of micro-texture is associated with poor fatigue performance in Ti alloys 16,25,26 , due to the production of coherent, continuous slip bands 27 that induce stressconcentrations at boundaries of the clusters 28 and propagates short-cracks 29 .
For the transformed-β colonies, Table 3 shows that two aspects of the microstructure account for most of the uncertainty in the mechanical response. These are (i) the crystallographic texture (58%) due to the presence of micro-texture, and (ii) α/β lath thickness (33%). As discussed before, the size-effect associated with the α/β interfaces is an important strengthening mechanism in this alloy. The natural variability in the interface thickness results in pockets of high and low strength material, resulting in a variance in the stress and plastic flow. The results indicate that uniformity of the microstructure through reduction of micro-texture and having α/β interfaces with similar thickness, can alleviate extreme values of stress and plastic strain in the microstructure. This can improve the fatigue performance of dual-phase Titanium alloys.

DISCUSSION
This paper develops an uncertainty-quantified parametrically homogenized constitutive model (UQ-PHCM) for dual-phase α/β Titanium alloys such as Ti6242S. It is an extension of the PHCMs for polycrystalline α-phase Titanium alloys developed in refs. 1,6,7 . In this paper a class of dual-phase α/β Ti alloys containing a single variant of crystallographic orientation corresponding to parallel α/ β lamellae in each colony is considered. The authors recognize that this is a limiting assumption, given that colony structures can be quite complex, e.g., in Widmanstätten or "basket-weave" microstructures, which are characterized by multiple sets of α lamellae with different variants co-existing within the β grains. This simplifying assumption avoids complex interaction of multiple BOR variants and may affect the polycrystalline deformation mechanisms and response. However, this assumption is deemed sufficient given that the focus of this paper is on the development of an overall multiscale framework.
The micromechanical deformation-informed PHCMs overcome prohibitive computational costs of other homogenization methods, and bridge length-scales through explicit microstructural representation in structure-scale constitutive models. The forms of equations in these thermodynamically consistent models are chosen to reflect fundamental deformation characteristics such as objectivity, anisotropy, tension-compression asymmetry, and history/path-dependence. Constitutive coefficients are functions of RAMP that represent distributions of key morphological and crystallographic descriptors in the microstructure. The functional forms are determined by machine learning tools operating on a data-set generated by micromechanical analysis. The PHCMs are readily incorporated in FE codes like ABAQUS through userdefined material modeling windows such as UMAT. Significantly reduced number of solution variables in the PHCM simulations compared to micromechanical models, make them several orders of magnitude more efficient, but with comparable accuracy.
An equivalent PHCM is developed to model dual-phase polycrystalline α/β Titanium alloys like Ti-6Al-2Sn-4Zr-2Mo. The equivalent PHCM assumes a Taylor model of uniform deformation gradient for the two phases, and employs a phase volume-fraction based weighted-averaging rule to obtain the equivalent material response from the individual PHCM responses of primary α and transformed β phases. Deterministic PHCMs are developed individually for the primary α and transformed β phases from Uncertainty-quantified PHCM-predicted stress-strain response of the dual-phase Ti-6242 alloy for various sources of uncertainty and comparison with the CPFE model response (shown with markers). a for uniaxial constant strain-rate loading of _ ε ¼ 1 10 À4 s À1 , showing the mean predicted stress and its 1σ interval corresponding to the combination of model reduction and data sparsity errors, b for strain-controlled cyclic loading with strain amplitude of Δε = 2% and period of 40 s with the 1σ interval corresponding to the microstructural uncertainty delineated in Table 2, and c zoomed-in views of boxes A and B showing details of the evolving 1σ interval calculated by the UP method.  13 and yields accurate results in comparison with CPFEM. A major advantage of this volume fraction-based weighting method over one that couples the effect of both phases into a single unified model is that it avoids the need to create a very high-dimensional micromechanical data-set for machine learning-based coefficient function evaluation. This requires a prohibitively large number of CPFE solutions for different M-SERVEs by simultaneously varying RAMPs belonging to both phases. The equivalent PHCMs for dual-phase polycrystalline microstructures capture important mechanical behavior observed in Ti alloys, viz., anisotropy, length-scale dependent flow stresses, tension-compression asymmetry in initial yield stress and hardening, grain size, lath size and strain-rate dependent plastic flow, and cyclic hardening under reversed loading conditions. The resulting constitutive equations consists of an anisotropic elastic model, an anisotropic yield function with tension-compression asymmetry and an anisotropic hardening model with a non-linear kinematic hardening rule. The machine learning toolkit 30 is used to express constitutive parameters as functions of RAMPs that characterize microstructural features such as distributions of crystallographic orientation and misorientation, grain size and lath size. The PHCM performance is assessed by comparing its predictions with CPFE simulations for a validation M-SERVE under different loading conditions. The PHCMs are found to yield satisfactory results for all loading conditions. While experimental validation of the PHCMs for a near-α Ti alloy has been successfully conducted in ref. 6 , there is a lack of experimental results for validation of the dual phase α/β Ti6242S alloys considered here, in the open literature. Such experimental datasets will comprise results of macroscopic coupon experiments such as tensile tests, along with relevant microstructure characterization data sets. This will be pursued in future studies.
An uncertainty quantification formulation is subsequently developed for the PHCMs to quantify and propagate uncertainties that exist at different stages of its development. Three sources of uncertainty, viz., model reduction error, data sparsity error, and microstructural uncertainty are considered. Consequently, the constitutive coefficients in PHCM are assumed to be random. The probability distributions of these constitutive coefficients are determined using Bayesian inference. The uncertainty in the constitutive coefficients is propagated to response variables such as stresses, plastic strains and state variables using a Taylor series expansion-based uncertainty propagation method. The relative contribution of different sources of uncertainty and the RAMPs towards total uncertainty in stresses and plastic strains is studied through numerical examples.
The uncertainty-quantified PHCMs present a unified framework in which microstructure-dependent material response with quantified uncertainties can be efficiently generated. The efficiency and accuracy of PHCMs, along with their easy integration in commercial software, make them highly promising tools for large scale structural analysis.

METHODS Overview
The systematic development of UQ-PHCM from micromechanical CPFE simulations is described below. Crystal plasticity constitutive model is summarized in the next section for dual-phase alloys following which a deterministic PHCM is developed. Different sources of uncertainty in PHCM development are discussed at the end and Bayesian inference is employed to quantify these uncertainties.

Crystal plasticity model for dual phase Titanium alloys
A dual-phase Titanium alloy Ti6242S is modeled in this paper. Samples of this alloy have been experimentally characterized in earlier work 9,13 and exhibits a bimodal (duplex) microstructure as shown in Fig. 1a. The αcolonies in Fig. 1a nucleate from prior-β grain boundaries during cooling of the alloy from the β phase. In general, colony structures can be quite complex such as Widmanstätten or "basket-weave" microstructures, which are characterized by multiple sets of α lamellae with different variants coexisting within the β grains [31][32][33] . These are typically achieved with higher cooling-rates or higher content of β-stabilizing elements 34 . The Widmanstätten microstructures involve interaction of a large number of slip systems. Relative crystallographic orientation of colonies belonging to different Burgers orientation relationship (BOR) variants within the transformed-beta phase can introduce complex mechanisms of dislocation transmission or pile-up at the colony interfaces. The boundaries between different oriented alpha variants provide strong barrier to dislocation activities, which can affect deformation mechanisms and thereby the mechanical response. Figure 1a shows different crystallographic orientation variants in the transformed beta phase region. However, there is preponderance of α lamellae in each colony that have identical crystallographic orientations, and thus belong to a single variant of BOR 32 . With this general observation, the present paper develops a model for a single variant of crystallographic orientation. This simplifying assumption that avoids complex interaction of multiple BOR variants and may have some effect on the polycrystalline deformation mechanisms and response. Since the focus of this paper is on the development of an overall multiscale framework, this simplifying assumption is deemed sufficient. The bimodal microstructure is characterized by equiaxed primary α grains consisting of hcp crystals and transformed β grains consisting of alternating laths of α (hcp) and β (bcc) phases, as shown in the schematic of Fig. 1a. Depending on the cooling-rate, the α colonies may be as large as the prior β grains 34 , and of comparable size to the primary-α grains in the duplex microstructure 9 , as shown in Fig. 1a. The volume fraction of the α and β lamellae in the transformed β colonies are assumed to be 88% and 12% respectively, from experimental observations discussed in refs. 13 During processing, the nucleation and growth of the α laths within the transformed β matrix follows one of 12 possible crystallographic orientations, given by BOR expressed as (101) β ||(0001) α , ½111 β jj½2110 α

36
. An example of this relationship is shown in Fig. 1b, which aligns a 1 ð½2110Þ slip direction of hcp crystal with the bcc b 1 (½111) slip direction. This constitutes a soft slip mode since plastic slip is easily transmitted at the lath boundary 9 . On the other hand, significant misalignment exists between α phase a 2 ð½1210Þ and β phase b 2 slip directions, and also between the a 3 ð½1120Þ and all 111 h i β directions in the β phase. These slip systems constitute hard slip modes as plastic slip is arrested at the lath boundary. As a result of the soft and hard slip modes, the ease of slip transmission for a 1 , a 2 , a 3 basal and prism slips varies significantly 13 . Depending on the hard and soft modes of slip transmission between the α and β laths, characteristic lengths are developed for Hall-Petch type, size-dependent slip system resistances due to dislocation barriers. These developments have been detailed in ref. 9 and summarized in the following sections. The colonies are segmented in DREAM.3D software as transformed β grains 37 and their BOR variants are identified for subsequent M-SERVE construction and CPFE analysis.
In the crystal plasticity model, the primary α phase is explicitly modeled using 30 slip systems. A computationally efficient equivalent homogenized model has been developed for transformed β grains 13 , based on the Taylor model assumption of uniform deformation gradient for all phases with known volume fractions. The model is discussed in the subsequent sections. The equivalent homogenized model has 78 aggregated slip systems corresponding to 30 secondary α slip systems and 48 β slip systems in the transformed β grain. While the α and β lath sizes and shapes can change across the transformed β grains 38,39 , they are assumed to be constant for all grains in the M-SERVE. In addition to the soft and hard slip modes discussed above, the bcc slip systems in {112}〈111〉 and {123}〈111〉 slip families are further divided into soft and hard slip directions. The {110} 〈111〉 slip systems exhibit symmetry with respect to the slip direction, while {112}〈111〉 and {123}〈111〉 slip systems are asymmetric with respect to the slip direction 40 . This leads to a slip direction-dependent shear strength for these slip systems, which are considered either soft or hard depending on the slip direction 13,40 . The crystal plasticity parameters are calibrated using simulations with an image-based CPFE model and matching the simulation response with those from single crystal and polycrystalline experiments. The general crystal plasticity formulation is the same for both hcp and bcc phases with the only difference introduced in the hardening laws.
Crystal plasticity constitutive relations for hcp and bcc phases. The crystal plasticity model introduces a multiplicative decomposition of the deformation gradient F into a plastic component F p in the stress-free intermediate configuration and an elastic component F e as: The second Piola-Kirchhoff stress S expressed in the intermediate configuration is related to the elastic Green-Lagrange strain E e as: Here C is a fourth order elastic stiffness tensor that is assumed to represent transversely isotropic and cubic symmetries for the hcp and bcc phases respectively. The plastic velocity gradient L p is related to the plastic slip-rate _ γ α and the Schmid tensors s α 0 ¼ m α 0 n α 0 for a slip system α with slip normal m α 0 and slip direction n α 0 , given as: Here nslip is the number of slip systems (30 for hcp phase and 48 for bcc phase). The plastic slip-rate is given by a power-law flow rule for both the hcp and bcc phases, as: where _ γ α is the reference slip rate, m is the strain rate sensitivity exponent, τ α is the resolved shear stress and χ α is the slip system back-stress. Evolution of χ α is given by the Armstrong-Frederick type non-linear kinematic hardening law 41 as, with c and d representing direct hardening and dynamic recovery coefficients respectively. τ α GP and τ α GF are the hardening contribution from geometrically necessary dislocations (GND) due to short and long range stresses, given by: c α 1 and c α 2 are material constants and G α , Q α and b α correspond to shear modulus, activation energy and Burgers vector respectively. ρ α GP and ρ α GF represent GND densities parallel and normal to slip plane and are obtained from the curl of plastic deformation gradient 17,28 . The slip system hardening resistance g α due to statistically stored dislocations is represented using the following hardening law with stress saturation.
where g α HP ¼ K α ffiffiffiffi D α p is the slip system resistance due to the Hall-Petch effect and D α is the characteristic length-scale governing the size-effect 9 . D α for different slip systems of hcp and bcc phases are based on the soft or hard slip for a given slip system and are given in Tables 4 and 5 respectively.
The hardening-rate on individual slip systems of the hcp components evolves as, Here h β 0 andg β s correspond to the initial hardening rate and saturation stress respectively. r and n are the hardening and saturation exponents respectively. For the bcc phase, the evolution of self-hardening rate is given by, where γ a ¼ R t 0 P nslip β¼1 j_ γ β jdt. The parameters h β 0 and h β s are the initial and asymptotic hardening rates, while τ β 0 and τ β s correspond to initial and asymptotic saturation stresses. γ a is the accumulated plastic slip on all the bcc slip systems at a given point. Finally, a yield-point phenomenon is attributed to the resistance due to precipitate shearing τ α p on the basal and prismatic slip systems in hcp crystals, and is represented using the following relationship 7 .
Hereτ α p andε α p are the yield-point phenomenon stress and plastic strain magnitudes respectively, and ϵ p is the effective plastic strain. For the bcc phase however, the yield-point phenomenon is not seen in experimental observations and is not considered.
It has been observed in experiments that Ti alloys exhibit tensioncompression asymmetry in their overall response. This has been attributed to mechanisms such as residual stresses in colonies during the growth processes and differential slip transmission mechanisms, depending on the loading direction 9 . The crystal plasticity parameters are separately calibrated in tension and compression and a stress-state dependent, weighted-averaging scheme is employed to determine different slip system parameters. The weight W t for tensile parameters is determined based on the lode angle θ, corresponding to the local stress S state in the π-plane. The lode angle is measured clockwise from the third positive principal axis 42 . The π-plane is divided into 6 regions by pure tensile and S. Kotha et al.
compressive stress axes. The weighting parameter W t for tensile parameters changes smoothly from 0 to 1 corresponding to pure compressive and tensile stress axes respectively. The weighted averaging ensures that the hardening parameters change smoothly for all the stress states. The weight for tension is given as: where S dev 1 ; S dev 2 and S dev 3 are the principal values of the deviatoric part of the second P-K stress tensor S and J 2 is its second invariant. A typical hardening parameter, e.g., g α 0 is determined as, where g α 0 ðTÞ and g α 0 ðCÞ correspond to initial slip system resistances calibrated from uniaxial tensile and compression experiments.
Equivalent homogenized model for transformed β beta colonies. In CPFE simulations, each grain in the M-SERVE is assumed to consist of either primary α or transformed β phase. In transformed β grains consisting of alternating laths of secondary α (hcp) and β (bcc) phases, an equivalent homogenized model has been developed in ref. 13 using assumptions of a uniform deformation gradient in the Taylor model for both hcp and bcc phases. The resulting model for colonies consists of an equivalent single crystal with 78 slip systems. The Cauchy stresses σ in each individual phase is: Using a volume-fraction based weighted-averaging rule, the homogenized stress in transformed beta colonies σ TB is obtained from the stresses in individual phases and their volume fractions as, where v hcp f and v bcc f are the volume fractions of the α and β laths in the colonies that are weights. In ref. 13 , these are taken to be 88% and 12% respectively.
Calibration of the crystal plasticity model parameters is an important step in the development of UQ-PHCMs. Deterministic calibration of these parameters from limited number of experimental data sets often results in non-unique parameters. Recently, this uncertainty has been accounted for in refs. 43,44 by employing Bayesian calibration, which treats the crystal plasticity parameters as random variables. However, random model parameters necessitates performing stochastic CPFE simulations, and this would significantly increase the computational burden of generating the PHCM calibration dataset. Therefore, deterministic calibration that uses genetic algorithm based optimization has been used to calibrate the crystal plasticity parameters for primary α Ti6242S from single crystal and polycrystalline tension/compression and creep tests in refs. 9,13,17,28,45,46 . The calibrated model, along with its validation for the primary α phase are given in ref. 1 . Similarly, tension/ compression tests for single α − β colonies have been used to calibrate crystal plasticity parameters for secondary α and β phases in refs. 9,13 . The CPFE model for the polycrystalline Ti6242S alloy (70% primary α and 30% transformed β phase) has been validated in ref. 9 for constant strain-rate and creep tests. The calibrated parameters for secondary α and β phases are summarized in Tables 4 and 5 respectively. The resulting calibrated crystal plasticity model is taken to be the reference model, and is used for all simulations in this paper.
Equivalent homogenized model for dual phase polycrystals with Taylor assumption. The idea of the equivalent homogenized model for transformed β colonies given in previous section is extended in this work Nature of slip (T) and (C) denote parameters for tension and compression tests respectively. D g is the grain size, D c is the colony size that is taken to be equal to D g , and l α is the α-lath thickness.
to model polycrystalline microstructural SERVEs or M-SERVEs that contain a combination of primary α and transformed β grains. A typical polycrystalline M-SERVEs is illustrated in Fig. 7a, where each grain is either a primary α (white) or a transformed β (black). The volume fraction of transformed β grains is 52% in this M-SERVE. The corresponding (0001) and (1120) pole figures of crystallographic texture for primary and secondary hcp phase is shown in Fig. 7b.
To examine the effectiveness of the equivalent homogenized model, the response of the polycrystalline M-SERVEs is simulated with three different volume fractions of the transformed β grains, viz., V TB = 22, 52 and 75%. The orientation of the bcc lath in a transformed β grain is determined from the adjacent secondary hcp lath orientation through the BOR. The M-SERVE is subjected to a constant, true strain-rate loading of 0.001 s −1 along the Xdirection. The resulting volume-averaged Cauchy stress-true strain plots are shown in Fig. 7c for the three different volume fractions, designated as (CPFE).
For the equivalent model with Taylor assumption, the M-SERVE is also simulated by successively assuming only primary α grains (V TB = 0) and transformed β grains (V TB = 1) in the M-SERVE. The volume-averaged stresses in the M-SERVE for these two limiting cases are denoted as σ PA ij and σ TB ij respectively. Using the volume-fraction based weighted-averaging rule with phase volume fractions on the volume-averaged stresses for each phase, the effective stress in the equivalent model is expressed as: The corresponding stresses using Eq. (15) are plotted in Fig. 7c and designated as (RoM). The results of the equivalent model match those of the two-phase CPFE model accurately for all volume fractions. This approach is consequently used to explicitly represent the effect of the primary α and transformed β phases in the PHCMs of the dual-phase Ti alloys discussed next.
Parametrically homogenized constitutive models for dual-phase titanium alloys Physics-based PHCMs have been developed for polycrystalline microstructures of single phase crystalline materials, e.g., the near-α Ti6242S alloy, in refs. 1,6 . The general forms of equations representing the evolution of state variables are chosen to be consistent with microscopic mechanisms of deformation, e.g., anisotropy, tension-compression asymmetry, hardening-softening behavior etc. The first law of thermodynamics is used to bridge length-scales and express constitutive coefficients as functions of RAMPs. Thermodynamic equivalence, according to the Hill-Mandel condition 47 defines a homogenized material as being energetically equivalent to a heterogeneous material with polycrystalline microstructures. This paper extends the previous work to dual phase Ti6242S alloys containing primary α and transformed β phases. Steps in the development of the corresponding PHCMs are discussed next.   D c is the colony size assumed to be equal to the grain size D g and l β is the β-lath thickness.

S. Kotha et al.
Sensitivity analysis and identification of RAMPs. The first step in PHCM development is a detailed sensitivity analysis to identify important microstructural descriptors and their distributions that govern the homogenized material response. These microstructural descriptors are characterized by a set of RAMPs, which relate the constitutive parameters in PHCM with the underlying microstructure. Functional dependencies of constitutive parameters in terms of the RAMPs is an important feature of the PHCMs. Different microstructural descriptors and RAMPs that influence the homogenized elasto-plastic response of primary α and transformed β M-SERVEs are given in Table 6. The associated RAMPs are defined below.
(i) Texture tensor (I tex ): The crystallographic c-axis orientation distribution of a M-SERVE is compactly represented by a texture tensor obtained from the weighted c-axes orientations of individual grains. It is defined as: whereĉ ðiÞ and V (i) correspond to the unit vector along the c-axis orientation and the volume of i th grain in the M-SERVE respectively. n G is the number of grains and V is the total volume of the M-SERVE. The eigenvalues g α and eigen-vectors v α of the texture tensor I tex , correspond to the texture intensity parameters and material symmetry axes respectively 6 . This tensor accurately represents the overall elastic stiffness for different crystallographic textures in the polycrystalline ensemble, as demonstrated in next section.
(ii) Lattice orientation with respect to material symmetry axes (OMA): To account for the influence of the orientation of slip systems on the homogenized plastic response, the orientation with respect to material symmetry axes (OMA) is introduced as a RAMP. It is defined in terms of the basal and prism Schmid tensors for each grain with respect to the material symmetry axes, and expressed as: Max jv α :s where S ðiÞ 0;basal and S ðiÞ 0;prism are the Schmid tensors for basal and prism 〈a〉slip in i th grain of the polycrystalline ensemble. v α are the unit vectors along the material symmetry axis, derived from the texture tensor I tex . The orientation of the bcc phase in the transformed β colony is determined from the adjacent hcp orientation through the BOR. Consequently, it suffices to characterize the orientation through the OMA functions derived from primary and secondary hcp phase crystallographic orientations. The OMA takes into account the 〈a〉-axes distributions of the grains in the microstructure and is able to represent its influence on the overall plastic slip in the M-SERVE.
(iii) Grain-pair misorientation parameter (A θmis ): The effect of misorientation on the homogenized plastic response is represented using a grain-pair misorientation RAMP A θmis , defined as: This parameter captures the fraction of grain pairs that have smaller than 15 ∘ misorientation angle θ mis between theirĉ axes. It is a measure of the ease of slip transfer across grain boundaries.
(iv) Mean and standard deviation of grain size distribution (D μ g and D σ ln ): Size effect in PHCM is represented using mean D μ g and standard deviation D σ ln of the grain size distribution, which is found to follow a log-normal distribution. The two RAMPs are defined as: Here D i is the equivalent diameter of i th grain in the M-SERVE consisting of n G grains, D μ g is the average grain size and D μ ln and D σ ln are the mean and standard deviation of the log-normal fit to the grain size distribution.
(v) Alpha and Beta lath thickness (l α and l β ): The response of the transformed β phase in the polycrystalline ensemble depends on the distributions of α and β lath thicknesses, in addition to grain size. This is because, the characteristic length in the Hall-Petch relationship of crystal plasticity model in Eq. (7) depends on the orientation of α and β laths with respect to loading direction. Different characteristic lengths are employed for loading along different directions 9 . Therefore the thicknesses l α and l β of the α and β laths are considered as RAMPs in the PHCMs for the transformed β phase. However, as the volume fraction of secondary α in colonies is assumed to be constant (88% in the current work), only l β is considered as an independent RAMP. A value of l α~3 l β has been determined for 88% volume fraction of secondary α in colonies in ref. 9 .
Constitutive equations representing the PHCMs. The general forms of the constitutive equations in PHCMs, representing the homogenized response of polycrystalline microstructures, are chosen to be consistent with microscopic mechanisms of deformation, e.g., anisotropy, tensioncompression asymmetry, cyclic hardening, strain-rate dependency etc. Corresponding to the crystal plasticity relations, the homogenized elasticplastic relations for the primary α and transformed β phases are assumed to be similar, with the difference being in the hardening laws. Thermodynamic consistency of these constitutive equations is demonstrated through relations established in the Supplementary Information. The PHCM constitutive parameters are calibrated from homogenized response of CPFE simulations for different microstructures and loading conditions. The sensitivity of the elastic and plastic response on the microstructural RAMPs is delineated in Table 6. The RAMP-based functional forms of constitutive coefficients are discussed next.
RAMP-dependent functional forms of anisotropic elastic stiffness tensor. The homogenized elastic stiffness tensor for α phase polycrystalline microstructures has been shown to depend only on the underlying crystallographic texture and single crystal (hcp) elastic stiffness in ref. 6 . Similarly, for the transformed β phase, the homogenized elastic stiffness is found to depend on the crystallographic texture, the elastic stiffness of the hcp and bcc phases and their volume fractions. While the hcp and bcc phases are assumed to have transversely isotropic and cubic symmetries respectively, the resulting homogenized elastic stiffness in the transformed β phase may have arbitrary symmetry depending on the crystallographic orientation of the hcp crystals. The bcc laths have a higher stiffness compared to the hcp phase. This leads to an overall increase in elastic stiffness for the transformed β phase over the pure α phase with same crystallographic texture. Orthotropic symmetry has been shown to accurately describe the homogenized elastic stiffness tensor C in ref. 6 . The same symmetry is assumed for both primary α and transformed β phase M-SERVEs in this study. The homogenized stiffness tensor C for the polycrystalline M-SERVE relates the macroscopic elastic Green-Lagrange strain E e to the macroscopic second Piola-Kirchhoff (2nd P-K) stress S through the relation: where the macroscopic elastic Green-Lagrange strain tensor is defined as: The macroscopic elastic deformation gradient F e is obtained from the relation F e ¼ F F pÀ1 , involving the total deformation gradient F and the plastic deformation gradient F p . The texture tensor I tex in Eq. (16) is used to express the crystallographic dependency of the elastic stiffness in material symmetry coordinate system C mat as: For primary α-phase M-SERVEs: For transformed β-phase M-SERVEs: where C Eq ijkl ðv α Þ is the equivalent single crystal elastic stiffness tensor in transformed β colonies. It is obtained from the single crystal stiffness tensor C hcp ijkl ðv α Þ of the hcp phase and C bcc ijkl ðv α Þ of the bcc phase as: The non-zero components of the single crystal stiffness tensors from ref. 13 are given in Table 7, where v 1 is the axis of transverse isotropy of the hcp single crystal.
The accuracy of the proposed parametrization in Eqs. (22) and (23) is assessed by comparing its predictions with elastic stiffness coefficients obtained from CPFE simulations on M-SERVEs having different crystallographic texture distributions. These distributions are represented by their texture intensity parameters as shown in Fig. 8a. The CPFE simulationbased elastic stiffness components are obtained from the derivatives of the homogenized stresses for prescribed perturbations in each component of the strain tensor 6 . The accuracy of the PHCM elastic stiffness can be seen from the plots of the normal components of the stiffness tensor in Fig. 8b and c for primary α and transformed β phases respectively. Excellent agreement is found for all components with less than 2% error. Further, it is seen that the transformed β phases have higher elastic stiffness components compared to the primary α phases for a given crystallographic texture.
Consistent with lattice rotation in CPFE model, the material symmetry coordinate system given by the eigenvectors v α of the texture tensor, is assumed to undergo a deformation-dependent rotation given as: R mat (0) corresponds to the initial material symmetry frame constructed from the initial symmetry axes and expressed in a fixed Cartesian coordinate frame with unit basis vectors e α as: Anisotropic yield function and plastic flow rule in PHCM. The PHCM equations for plasticity must exhibit the following characteristics for consistency with micromechanical observations in CPFEM simulations, discussed in previous sections. They are: • Plastic anisotropy arising from the large difference in slip system resistances for basal/prism and pyramidal slip systems of the α phase.
• Tension-compression asymmetry in the yield surfaces and post-yield behavior arising from tension-compression asymmetry of individual slip systems. Anisotropy and tension-compression asymmetry is more pronounced in transformed β due to the presence of lath structures with soft and hard slip modes.
• Length-scale dependency due to grain or lath size-dependent slip system resistances expressed in Eq. (7). The macroscopic plastic deformation gradient F p is obtained by integrating the macroscopic plastic velocity gradient L p , which is additively decomposed into the rate of deformation D p and plastic spin W p tensors.
Assuming the plastic spin to be negligible as in crystal plasticity model, the plastic velocity gradient is expressed using the flow rule as: where D p 0 and m correspond to the reference strain-rate and the ratesensitivity parameter respectively. N is the normal to the flow surface. The flow rule in Eq. (27) is assumed to be associative, and the unit vector representing the direction of the rate of plastic deformation tensor is expressed in terms of the normalized gradient of the yield function with respect to the second P-K stress as: where Á k k is the norm, defined in Eq. (37). Y is the effective transformed stress for an anisotropic yield function 48 that accounts for tensioncompression asymmetry and is expressed as: where k is the tension-compression asymmetry parameter and a is an exponent that governs the smoothness of the yield surface. λ 1 , λ 2 and λ 3 are the principal values of a macroscopic transformed effective stress P , given as: where L is a microstructure dependent, fourth order transformation tensor containing anisotropy coefficients. These coefficients are expressed in the material symmetry coordinate system defined by the eigen-vectors v α of the texture tensor I tex . The matrix form of the anisotropic tensor is chosen to be 49 : Plastic incompressibility condition leads to the following relation between the components: Evolution of the back stress χ is governed by a modified Armstrong-Frederick type non-linear kinematic hardening law given as: Here C, D, E, G are calibrated constants. The exponential term within brackets manifests a smooth transition from elasticity to plasticity, as observed upon every load reversal in M-SERVEs of the transformed β. The expression within the Macaulay brackets is non-zero only during the transition from elasticity to plasticity in every load reversal, when the back stress and the flow direction have opposing directions.
Deciphering constitutive coefficients in terms of RAMPs through machine learning. The dependency of PHCM coefficients on RAMPs like OMA αβ , D μ g and l β is obtained using machine learning tools. The tools utilize symbolic regression to decipher functional forms of the coefficients in terms of the RAMPs. Unlike traditional regression, where parameters in a given equation are optimized, symbolic regression searches for both the functional form and the coefficients that best fit a given data-set. These functional forms are obtained in this work using a machine learning toolkit Eureqa 30 that is based on genetic programming 50 . The code uses calibrated anisotropy coefficients (α ii , γ ij ) and the RAMPs (OMA αβ , D μ g and l β ) as inputs. Random functional forms for anisotropy coefficients are generated from different combinations of RAMPs and these are represented as tree structures. Fitness values for each of these equations is evaluated using the input data for anisotropy coefficients and the equations are ranked according to their fitness values. Natural selection-based genetic operations, such as mutation, cross over and elitism, are used to generate new functional forms for the next generation. This procedure is repeated for a large number of generations until optimum functional forms that best fit the data are obtained. A Pareto front is generated and the final functional form is chosen based on the complexity and accuracy of the functional form.
The anisotropy tensor L in Eq. (31) is expressed in terms of the RAMP OMA αβ defined in Eq. (17). This dependencies for M-SERVEs of primary α phase have been derived in ref. 6 and summarized in Box 1. For M-SERVEs of transformed β phase, where soft and hard slip systems strongly govern the plastic response, the transformation tensor L depends on the average grain size D μ g and the β lath size l β , in addition to OMA αβ . For developing functional forms of the microstructural dependencies using machine learning, a data-set of M-SERVEs with different crystallographic orientations, grain and lath sizes is generated. Steps in this process are delineated below.  To capture the effect of lath size on plastic anisotropy, the M-SERVE M3 is chosen and two different β lath sizes, viz., l β = 0.8 μm and l β = 3 μm are assigned, while maintaining an approximate ratio lα lβ $ 3. For each of these M-SERVEs, 19 different crystallographic textures used above are assigned, creating a total of 38 different M-SERVEs. The anisotropy coefficients α ii , γ ij in Eq. (31), tension-compression asymmetry parameter k and yield surface exponent a in Eq. (29) are evaluated by calibrating the yield function in Eq. (29) to CPFE simulationbased yield stresses along different directions of the 95 M-SERVEs created above. A total of 30 different uniaxial and biaxial stress-controlled CPFE simulations are performed on each of the above M-SERVEs to extract the initial yield stresses 6 in different directions as shown in Fig. 9. The initial yield stresses are assumed to correspond to an effective plastic strain of 0.3%. Back stresses are assumed to be negligible at this small plastic strain  and are not considered while calibrating the anisotropy coefficients. The calibrated anisotropy coefficients for the M-SERVE with crystallographic texture shown in Fig. 9a, are given in Table 8     of single crystal slip-system hardening parameters. Different aspects of the PHCM hardening models for the primary α and transformed β phases, along with the calibration of their microstructural dependency functions, are established in this section. (i) Hardening laws for the primary α phase: For the primary α phase, the hardening response is characterized by a smooth transition from elasticity to plasticity, followed by a constant hardening slope as shown in Fig. 11a. To account for the transient yield-point phenomenon and subsequent hardening, the flow stress Y 0 in Eq. (27) is represented by a microstructuredependent strain hardening law, which is formulated as a modified Vocetype law 51 as: here e Y 0 is the initial flow stress and ψ is a constant, chosen to be 0.75 in this work.α andβ are microstructure and loading direction-dependent parameters that are used to capture the transient yield-point phenomenon. H represents the microstructure and loading direction-dependent hardening slope. The parameter b controls the rate of cyclic hardening. For Ti alloys, the rate of cyclic hardening is observed to be high for the first~50 cycles, which is followed by a smaller constant-rate. In general, b has been expressed as a function of the size of the plastic strain memory surface 52 . However in the present work, it is taken to be a calibrated constant for simplicity. The effective plastic strain is defined as: where To capture the microstructure and loading direction dependency of the hardening parameters, a dimensionless parameter R Y is defined. This parameter accounts for the amount of 〈a〉 type basal and prism slip occuring in the M-SERVE for a given loading condition. This parameter has been shown to account for the anisotropy in hardening response of primary α M-SERVEs in ref. 6 . It is defined as: Microstructural dependency of R Y in the Eq. (38) comes from the orientation dependency of L. The loading direction dependency is accounted for through the second P-K stress tensor S. For the transformed β phase however, R Y also depends on the mean grain size (D μ g ) and lath thickness (l β ) through L. Since for the primary α phase, only the initial slip system resistance is dependent on tension-compression state, the calibrated tension-compression parameter k is observed to accurately represent the difference in flow stresses in tension and compression.
(ii) Hardening laws for the transformed β phase: The macroscopic hardening response of the transformed β phase depends on the strain-rate and the tension-compression stress-state, in addition to anisotropy of the single crystal hardening parameters. To capture the rate-dependent hardening response, a saturation stress-based hardening law (as in the CPFE model) is assumed for the transformed β phase. The microstructuredependent flow stress is given as: and the evolution of Y 0 is given as: where H is the microstructure-dependent hardening-rate expressed as: H 0 represents the initial hardening-rate and r is an exponent controlling the rate of hardening. Y sat is the saturation stress given as: where Y Ref sat is a reference saturation stress and n is a saturation exponent, which controls the strain-rate dependency of hardening. The PHCM hardening parameters for the transformed β phase are calibrated from CPFE simulations under uniaxial tension and compression loads. The calibration accounts for tension-compression asymmetry of single crystal hardening parameters for the secondary α phase in the CPFE model. Weighted averaging of the calibrated parameters are used for the PHCM hardening parameters subjected to arbitrary loading. The weighting function W t for tensile hardening parameters is determined from the PHCM-based macroscopic stress tensor using Eq. (11), for consistency with tension-compression states in CPFEM.
Machine learning-based determination of hardening parameters in terms of RAMPs. As discussed in previous sections, functional forms of the hardening parameters in Eqs. (36), (39), (41) and (42) in terms of RAMPs are determined using the machine learning toolkit Eureqa 30 that operates on a database created from CPFE simulations of polycrystalline M-SERVEs. The microstructural dependency of different hardening parameters in Eqs. (36), (39), (41) and (42) is obtained by calibrating the PHCM coefficients with homogenized response obtained from CPFE simulations of different M-SERVEs. The calibration process for primary α has been detailed in ref. 6 . Calibration of hardening parameters for the transformed β phase uses the following database of M-SERVEs.  Fig. 11a, b. To calibrate the back stress constants in Eq. (33), strain-controlled cyclic simulations are also performed on these M-SERVEs along three material symmetry axes, as shown in Fig. 11a and b. The sinusoidal load profile is applied for 3/4cycle with a peak strain of 3% and a time period of 120 s.  To determine the strain-rate dependency of transformed β hardening parameters and calibrate the rate sensitivity parameter, a uniform texture M-SERVE is created. Constant strain-rate simulations for three strain-rates, viz., 0.001, 0.01 and 0.1 s −1 are conducted in both tension and compression.
• For the M-SERVE M1 with 19 different crystallographic textures, straincontrolled dwell simulations are conducted along the three material symmetry axes to comprehend stress relaxation behavior with cycles. The imposed dwell loading consists of a hold strain of ϵ 0 = 1.2% with a ramp time of 1 s and a hold time of 120 s, as shown in Fig. 5d. The tensile stress at the end of strain hold and the compressive stress at the end of each cycle are obtained from CPFE simulations. The PHCM calibration process for both the primary α and transformed β phases are plotted in Fig. 11c and d respectively. Since stress relaxation is observed to occur predominantly in the first few cycles of loading, only 10 dwell cycles are simulated for the M-SERVEs. The homogenized stress-strain responses from the constant strain-rate and cyclic load simulations of all the M-SERVEs constitute a representative database. These data-sets are used by the machine learning toolkit 30 for calibrating the hardening parameters in tension and compression. The resulting functional forms of the hardening parameters in terms of RAMPs are given in Box 2.
Uncertainty quantification (UQ) and uncertainty propagation (UP) in PHCM Uncertainty in the multiscale PHCMs may be derived from multiple sources. These may be broadly represented by the following considerations: reduction error inherently exists in the machine learning-generated constitutive parameters, which propagates to the predicted material response variables. e.g. the macroscopic Cauchy stress σ t ð Þ.
• Data sparsity error: A source of uncertainty in the PHCM response is due to the finite size of the calibration data-set. For example, the number of M-SERVEs included in the PHCM calibration data-set is limited by the computational cost of performing the CPFE analyses. This collection of M-SERVEs cover a finite domain in the RAMP-space, and predictions made for microstructures outside this domain contain additional uncertainty as a function of distance from the calibration data points. This uncertain error is referred to as the data sparsity error, and is represented by the Bayesian posterior distribution of model coefficients. Within a Bayesian calibration framework, it may be shown that this type of uncertainty depends on both the number and distribution of the calibration data points in the RAMP-space 14 . For a given microstructure, the total deviation of the PHCM-predicted response from the micromechanics-based response is accounted for by the model reduction and data sparsity errors.

•
Microstructural uncertainty: Microstructural uncertainty may be due to limited experimental characterization data or due to inherent aleatoric uncertainty of the microstructure, which naturally arises from the manufacturing process. This uncertainty necessitates RAMPs to be represented as stochastic variables, rather than deterministic quantities.
During the development of PHCMs, the uncertainties due to model reduction and data sparsity are assessed by computational validation studies following model calibration, by comparing PHCM predictions with CPFEM results. For the dual-phase Ti alloys, these comparisons have been performed in "Results" section. However, the validation studies are performed only once, with a certain set of microstructures and load cases. In general, the prediction error will be unknown to the end-user for the particular microstructure and load case of interest. In addition, quantification of the effect of the microstructural uncertainty, in the model response by conventional sampling-based methods like the Monte Carlo methods requires repeated model evaluations. This is computationally prohibitive for any practical application. It is therefore desirable to have a built-in capability to evaluate the uncertainty in PHCM response, on the fly, as the analysis is taking place within a commercial FE software. To address these needs, stochastic PHCMs with uncertainty quantification and built-in uncertainty propagation capabilities are developed in this paper for dual-phase Titanium alloys. The methods follow a framework that has been introduced in ref. 14 . Assuming the CPFEM simulation-based micromechanical response as the true model, the three main sources of uncertainty discussed above are considered in this framework. A Taylorexpansion-based uncertainty propagation (UP) method calculates the uncertainty in the time and history-dependent material response variables like the macroscopic Cauchy stress, plastic strain and fatigue indicators. A Bayesian framework that is used for the calibration of UQ-enhanced stochastic PHCMs is summarized below, along with the formulation of the built-in UP method.
Formulation of stochastic PHCMs with Bayesian inference. In the stochastic PHCMs, functional forms of the microstructure-dependent constitutive parameters are extended to probabilistic models as: where θ ¼ θ i f g are the microstructure-dependent PHCM constitutive parameters identified in equations given in Boxes 1 and 2. ϕ x R ð Þ ¼ ϕ k x R ð Þ f g , k = 1, . . . , n b are the basis functions identified by machine learning in terms of the microstructure-based RAMPs x R , and C ¼ C ik ½ are random coefficients identified by Bayesian inference from CPFE-based calibration data-set. ε m ¼ ε m;i È É are terms in the model reduction error, which model the disparity between the unknown true values of the constitutive parameters for a given microstructure and the prediction of the functional form. These terms are assumed to be independent Gaussian random variables with zero mean and variances σ 2 m;i , i.e., ε m;i $ N 0; σ 2 m;i : The set of hyper-parameters σ m,i are identified initially by the maximum likelihood estimation, and held as constant during the Bayesian inference of coefficients C, discussed next. The CPFE-generated data-set for PHCM calibration, denoted by D, is used for the Bayesian inference of the stochastic PHCMs. The Bayesian update for the PHCM coefficients C is expressed as: where p C ð Þ and p CjD ð Þ are respectively the prior and posterior distributions of the coefficients. p DjC ð Þis the joint probability of producing the data D, given the PHCM coefficients C, which forms the likelihood function L C; D ð Þ¼p DjC ð Þ on the parameter space of C. Each evaluation of L involves simulations of the input data with stochastic PHCM, i.e., the collection of M-SERVEs under load cases described in previous section. Non-informative, flat priors p C ð Þ / 1 are used to avoid subjective bias of the posterior distribution. The posterior distribution p CjD ð Þ is sampled using the Hamiltonian Monte Carlo technique 53 with No-U-Turn sampler 54 . The uncertainty due to calibration data sparsity is entirely contained in the Bayesian posterior distribution p CjD ð Þ. The statistical moments of the posterior are calculated from the collected samples and are embedded in the stochastic PHCMs for macroscopic response predictions.
The above framework is used for Bayesian inference of the functional forms of the yield surface transformation parameters θ YF ¼ α ii ; γ ij È É of the matrix L in Eq. (31). The yield function calibration data-set D YF is generated in the previous section for (a) α-phase M-SERVE and (b) transformed-β phase M-SERVE. Figure 12 shows the calibrated yield function for two of the calibration microstructures (lines) compared to their CPFEM-based yield points collected under biaxial tests (markers). The black solid line is the mean yield surface and the gray shaded region is the 1σ interval corresponding to the model reduction error. The CPFEM-based yield points from the calibration data-set D YF (markers) are used to identify the hyper-parameters σ 2 m;i of the model reduction error in Eq. (46), as well as to perform Bayesian inference of the form coefficients based on Eq. (48). Figure 13 shows yield surfaces predicted by the Bayesian posterior distribution for two experimentally characterized microstructures, which are not included in the calibration data-set D YF . The black solid line is the posterior mean yield surface and the 1σ interval bounded by the dashed lines corresponds to the total uncertainty in the posterior prediction, including both model reduction and data sparsity errors. The dark-gray shaded region is the 1σ interval corresponding to the model reduction error only. The probabilistic predictions are compared to the CPFEM-based yield points (black markers). Figure 13c shows zoomed-in views of these comparisons.
It is emphasized that in contrast to deterministic PHCMs, the material constitutive parameters θ in stochastic PHCMs, e.g. yield function parameters fα ii ; γ ij g, hardening parameters f e Y 0 ;α;β; Hg, etc. are random variables. This is due to the randomness of the form coefficients C, the model reduction error terms ε m , and the RAMPs x R . Consequently, all material state and response variables like ε p , F p , σ etc. become correlated random variables, whose joint probability distribution evolves with deformation. A method for propagation of the uncertainty through material rate equations is presented next.
Built-in uncertainty propagation in stochastic PHCMs. A Taylor series expansion-based uncertainty propagation (UP) method has been developed in ref. 14 to propagate uncertainty among the PHCM constitutive parameters and evolving material state and output variables. This method is extended to PHCMs for the dual-phase Ti alloys in this paper. First, the operation of the deterministic material time-integration algorithm of PHCM is expressed as: where ξ ¼ ξ i È É is an evolving state vector. It contains all microstructuredependent constitutive parameters θ ¼ θ i f g, as well as the material state variables like the plastic deformation gradient F p , back-stress χ and effective plastic strain ε p . The state vector ξ may be generalized to a stochastic vector with a time-dependent joint probability distribution p ξ; t ð Þ. The mean and covariance of the stochastic state vector ξ are respectively defined as: for the time step t (n) . The non-linear operator I ξ ðnÞ Â Ã is expanded in the neighborhood of ξ (n) = μ (n) , which allows one to calculate the moments (μ (n+1) , ∑ (n+1) ) at time t (n+1) in terms of (μ (n) , ∑ (n) ) at time t (n) . The resulting time-marching rules are expressed using Einstein index notation as: where fI i g ¼ I½μ ðnÞ is the time-integrator evaluated with the mean state μ (n) at time t (n) . The quantities I i k , I i km etc. are partial derivatives of the time-integrator, i.e., evaluated at ξ (n = μ (n) . The fourth order central moment Σ kmrs ðnÞ is approximated using Isserlis' theorem 55 , based on Σ km ðnÞ . In the limit Σ km 0 ð Þ ¼ 0, this scheme is equivalent to deterministic time-integration with μ (n) = ξ (n) , Σ (n) = 0 for all t (n) . Similar expressions, based on series expansions are also derived to propagate the uncertainty from the material state variables to output variables, e.g., Cauchy stress, at every time step. The built-in UP method has been verified by comparison of the calculated stochastic response with Monte Carlo simulation results in ref. 14 .

DATA AVAILABILITY
The datasets generated during and/or analyzed during the current study will be available from the corresponding author on reasonable request.

CODE AVAILABILITY
The codes that are used to generate the results in the current study will be available from the corresponding author in a suitable format upon reasonable request.