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

Aug 5, 2020

### 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 work9,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 co-existing within the β grains31,32,33. These are typically achieved with higher cooling-rates or higher content of β-stabilizing elements34. 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 BOR32. 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 β grains34, and of comparable size to the primary-α grains in the duplex microstructure9, 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,35. The hcp crystal lattice consists of a total of 30 different slip systems divided into five different slip families, viz., basal (leftlangle {bf{a}}rightrangle ({0001}leftlangle 11overline{2}0rightrangle )), prism (leftlangle {bf{a}}rightrangle ({10overline{1}0}leftlangle 11overline{2}0rightrangle )), pyramidal (leftlangle {bf{a}}rightrangle ({10overline{1}1}leftlangle 11overline{2}0rightrangle )), first order pyramidal (leftlangle {bf{c}}+{bf{a}}rightrangle ({10overline{1}1}leftlangle 11overline{2}3rightrangle )) and second order pyramidal (leftlangle {bf{c}}+{bf{a}}rightrangle ({11overline{2}2}leftlangle 11overline{2}3rightrangle )). Similarly, the bcc crystal system consists of a total of 48 different slip systems divided into three different slip families—({{bf{b}}}_{1}({110}leftlangle 111rightrangle )), ({{bf{b}}}_{2}({112}leftlangle 111rightrangle )), ({{bf{b}}}_{3}({123}leftlangle 111rightrangle )). 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)α, ({[1overline{1}overline{1}]}_{beta }| | {[2overline{1}overline{1}0]}_{alpha })36. An example of this relationship is shown in Fig. 1b, which aligns ({{bf{a}}}_{1}([2overline{1}overline{1}0])) slip direction of hcp crystal with the bccb1 (([1overline{1}overline{1}])) slip direction. This constitutes a soft slip mode since plastic slip is easily transmitted at the lath boundary9. On the other hand, significant misalignment exists between α phase ({{bf{a}}}_{2}([overline{1}2overline{1}0])) and β phase b2 slip directions, and also between the ({{bf{a}}}_{3}([overline{1}overline{1}20])) and all ({leftlangle 111rightrangle }_{beta }) 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 a1, a2, a3 basal and prism slips varies significantly13. 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 β grains37 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 β grains13, 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 β grains38,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 direction40. This leads to a slip direction-dependent shear strength for these slip systems, which are considered either soft or hard depending on the slip direction13,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 Fp in the stress-free intermediate configuration and an elastic component Fe as:

$${boldsymbol{F}}={{boldsymbol{F}}}^{{rm{e}}}{{boldsymbol{F}}}^{{rm{p}}}.$$

(1)

The second Piola–Kirchhoff stress S expressed in the intermediate configuration is related to the elastic Green-Lagrange strain Ee as:

$${boldsymbol{S}}={mathbb{C}}:{{boldsymbol{E}}}^{{rm{e}}}quad {rm{where}}quad {{boldsymbol{E}}}^{{rm{e}}}=frac{1}{2}({{boldsymbol{F}}}^{{{rm{e}}}^{{rm{T}}}}{{boldsymbol{F}}}^{{rm{e}}}-{boldsymbol{I}}).$$

(2)

Here ({mathbb{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 Lp is related to the plastic slip-rate ({dot{gamma }}^{alpha }) and the Schmid tensors ({{boldsymbol{s}}}_{0}^{alpha }={{bf{m}}}_{0}^{alpha }otimes {{bf{n}}}_{0}^{alpha }) for a slip system α with slip normal ({{bf{m}}}_{0}^{alpha }) and slip direction ({{bf{n}}}_{0}^{alpha }), given as:

$${{boldsymbol{L}}}^{{rm{p}}}={dot{{boldsymbol{F}}}}^{{rm{p}}}{{boldsymbol{F}}}^{{{rm{p}}}^{-1}}=mathop{sum limits_{alpha = 1}^{{rm{nslip}}}}{dot{gamma }}^{alpha }{{boldsymbol{s}}}_{0}^{alpha }.$$

(3)

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:

$${dot{gamma }}^{alpha }={dot{tilde{gamma }}}^{alpha }{leftlangle frac{| {tau }^{alpha }-{chi }^{alpha }| -{tau }_{{rm{GP}}}^{alpha }}{{g}^{alpha }+{tau }_{{rm{p}}}^{alpha }+{tau }_{{rm{GF}}}^{alpha }}rightrangle }^{frac{1}{m}}mathrm{sign}({tau }^{alpha }-{chi }^{alpha }),$$

(4)

where ({dot{tilde{gamma }}}^{alpha }) 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 law41 as,

$${dot{chi }}^{alpha }=c{dot{gamma }}^{alpha }-d{chi }^{alpha }| {dot{gamma }}^{alpha }|,$$

(5)

with c and d representing direct hardening and dynamic recovery coefficients respectively. ({tau }_{{rm{GP}}}^{alpha }) and ({tau }_{{rm{GF}}}^{alpha }) are the hardening contribution from geometrically necessary dislocations (GND) due to short and long range stresses, given by:

$${tau }_{{rm{GP}}}^{alpha }={c}_{1}^{alpha }{G}^{alpha }{b}^{alpha }sqrt{{rho }_{{rm{GP}}}^{alpha }},quad quad {tau }_{{rm{GF}}}^{alpha }=frac{{Q}^{alpha }}{{c}_{2}^{alpha }{b}^{alpha 2}}sqrt{{rho }_{{rm{GF}}}^{alpha }},$$

(6)

({c}_{1}^{alpha }) and ({c}_{2}^{alpha }) are material constants and Gα, Qα and bα correspond to shear modulus, activation energy and Burgers vector respectively. ({rho }_{{rm{GP}}}^{alpha }) and ({rho }_{{rm{GF}}}^{alpha }) represent GND densities parallel and normal to slip plane and are obtained from the curl of plastic deformation gradient17,28. The slip system hardening resistance gα due to statistically stored dislocations is represented using the following hardening law with stress saturation.

$${dot{g}}^{alpha }(t)={sum limits_{beta = 1}^{{rm{nslip}}}}{q}^{alpha beta }{h}^{beta }| {dot{gamma }}^{beta }| ,text{,},{g}^{alpha }(0)={g}_{0}^{alpha }+{g}_{{rm{HP}}}^{alpha },$$

(7)

where ({g}_{{rm{HP}}}^{alpha }=frac{{K}^{alpha }}{sqrt{{D}^{alpha }}}) is the slip system resistance due to the Hall-Petch effect and Dα is the characteristic length-scale governing the size-effect9. 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,

$${h}^{beta }={h}_{0}^{beta }{left|1-frac{{g}^{beta }}{{g}_{{rm{s}}}^{beta }}right|}^{r}{rm{sign}}left(1-frac{{g}^{beta }}{{g}_{rm{s}}^{beta }}right)quad {rm{and}}quad {g}_{rm{s}}^{beta }={tilde{g}}_{rm{s}}^{beta }{left(frac{{dot gamma}^{beta}}{dot{tilde gamma}}right)}^{n}.$$

(8)

Here ({h}_{0}^{beta }) and ({tilde{g}}_{{rm{s}}}^{beta }) 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,

$${h}^{beta }={h}_{{rm{s}}}^{beta }+mathrm{sec}{h}^{2}left[left(frac{{h}_{0}^{beta }-{h}_{{rm{s}}}^{beta }}{{tau }_{{rm{s}}}^{beta }-{tau }_{0}^{beta }}right){gamma }_{{rm{a}}}right]left({h}_{0}^{beta }-{h}_{{rm{s}}}^{beta }right),$$

(9)

where ({gamma }_{{rm{a}}}=mathop{int}nolimits_{0}^{t}mathop{sum }nolimits_{beta = 1}^{{rm{nslip}}}| {dot{gamma }}^{beta }| mathrm{d}t). The parameters ({h}_{0}^{beta }) and ({h}_{{rm{s}}}^{beta }) are the initial and asymptotic hardening rates, while ({tau }_{0}^{beta }) and ({tau }_{{rm{s}}}^{beta }) 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 ({tau }_{{rm{p}}}^{alpha }) on the basal and prismatic slip systems in hcp crystals, and is represented using the following relationship7.

$${tau }_{{rm{p}}}^{alpha }={tilde{tau }}_{{rm{p}}}^{alpha }exp (-{overline{epsilon }}_{{rm{p}}}/{tilde{epsilon }}_{{rm{p}}}^{alpha }).$$

(10)

Here ({tilde{tau }}_{{rm{p}}}^{alpha }) and ({tilde{epsilon }}_{{rm{p}}}^{alpha }) are the yield-point phenomenon stress and plastic strain magnitudes respectively, and ({overline{epsilon }}_{{rm{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 tension-compression 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 direction9. 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 Wt 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 axis42. The π-plane is divided into 6 regions by pure tensile and compressive stress axes. The weighting parameter Wt 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:

$$begin{array}{l}{W}_{{rm{t}}}=frac{1}{2}(1+cos ,3theta )\ {rm{with}}cos ,theta =frac{{S}_{3}^{{rm{dev}}}}{2sqrt{{J}_{2}}},,sin ,theta =frac{sqrt{3}left({S}_{{rm{2}}}^{{rm{dev}}}-{S}_{{rm{1}}}^{{rm{dev}}}right)}{2sqrt{{J}_{2}}}end{array},$$

(11)

where ({S}_{{rm{1}}}^{{rm{dev}}},{S}_{2}^{{rm{dev}}}) and ({S}_{{rm{3}}}^{{rm{dev}}}) are the principal values of the deviatoric part of the second P-K stress tensor S and J2 is its second invariant. A typical hardening parameter, e.g., ({g}_{0}^{alpha }) is determined as,

$${g}_{0}^{alpha }={W}_{{rm{t}}}{g}_{0}^{alpha }({rm{T}})+(1-{W}_{{rm{t}}}){g}_{0}^{alpha }({rm{C}}),$$

(12)

where ({g}_{0}^{alpha }({rm{T}})) and ({g}_{0}^{alpha }({rm{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:

$${sigma }_{ij}^{{rm{hcp}}}=frac{1}{| {F}_{{rm{e}}}^{{rm{hcp}}}| }{F}_{ik}^{{{rm{e}}}^{{rm{hcp}}}}{S}_{kl}^{{rm{hcp}}}{F}_{jl}^{{{rm{e}}}^{{rm{hcp}}}},{sigma }_{ij}^{{rm{bcc}}}=frac{1}{| {F}_{{rm{e}}}^{{rm{bcc}}}| }{F}_{ik}^{{{rm{e}}}^{{rm{bcc}}}}{S}_{kl}^{{rm{bcc}}}{F}_{jl}^{{{rm{e}}}^{{rm{bcc}}}}.$$

(13)

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,

$${sigma }_{ij}^{{rm{TB}}}={v}_{{rm{f}}}^{{rm{hcp}}}{sigma }_{ij}^{{rm{hcp}}}+{v}_{{rm{f}}}^{{rm{bcc}}}{sigma }_{ij}^{{rm{bcc}}},$$

(14)

where ({v}_{{rm{f}}}^{{rm{hcp}}}) and ({v}_{{rm{f}}}^{{rm{bcc}}}) 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 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., VTB = 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 X-direction. 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 (VTB = 0) and transformed β grains (VTB = 1) in the M-SERVE. The volume-averaged stresses in the M-SERVE for these two limiting cases are denoted as ({overline{sigma }}_{ij}^{{rm{PA}}}) and ({overline{sigma }}_{ij}^{{rm{TB}}}) 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:

$${overline{sigma }}_{ij}=(1-{V}_{{rm{TB}}}){overline{sigma }}_{ij}^{{rm{PA}}}+{V}_{{rm{TB}}}{overline{sigma }}_{ij}^{{rm{TB}}}.$$

(15)

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 condition47 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.

#### 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 (Itex): 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:

$${{boldsymbol{I}}}^{{rm{tex}}}=frac{1}{V}mathop{sum }limits_{i=1}^{{n}_{{rm{G}}}}{hat{{bf{c}}}}^{(i)}otimes {hat{{bf{c}}}}^{(i)}{V}^{(i)},$$

(16)

where ({hat{{bf{c}}}}^{(i)}) and V(i) correspond to the unit vector along the c-axis orientation and the volume of ith grain in the M-SERVE respectively. nG is the number of grains and V is the total volume of the M-SERVE. The eigen-values gα and eigen-vectors vα of the texture tensor Itex, correspond to the texture intensity parameters and material symmetry axes respectively6. 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:

$${overline{mathrm{OMA}}}_{alpha beta }=frac{1}{V} {sum limits_{i = 1}^{{n}_{rm{G}}}},text{Max},left{| {{bf{v}}}_{alpha }.{{boldsymbol{s}}}_{0,{rm{basal}}}^{(i)}.{{bf{v}}}_{beta }| ,| {{bf{v}}}_{alpha }.{{boldsymbol{s}}}_{0,{rm{prism}}}^{(i)}.{{bf{v}}}_{beta }| right}{V}^{(i)},$$

(17)

where ({{boldsymbol{S}}}_{0,{rm{basal}}}^{(i)}) and ({{boldsymbol{S}}}_{0,{rm{prism}}}^{(i)}) are the Schmid tensors for basal and prism 〈a〉-slip in ith grain of the polycrystalline ensemble. vα are the unit vectors along the material symmetry axis, derived from the texture tensor Itex. 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 (({overline{A}}_{{theta }_{{rm{mis}}}})): The effect of misorientation on the homogenized plastic response is represented using a grain-pair misorientation RAMP ({overline{A}}_{{theta }_{{rm{mis}}}}), defined as:

$${overline{A}}_{{theta }_{{rm{mis}}}}=frac{,text{Grain boundary area with},,{theta }_{{rm{mis}}}, < , 1{5}^{circ }}{,text{Total grain boundary area in the SERVE},}.$$

(18)

This parameter captures the fraction of grain pairs that have smaller than 15 misorientation angle θmis between their (hat{{bf{c}}}) axes. It is a measure of the ease of slip transfer across grain boundaries.

(iv) Mean and standard deviation of grain size distribution (({overline{D}}_{{rm{g}}}^{mu }) and ({overline{D}}_{{rm{ln}}}^{sigma })): Size effect in PHCM is represented using mean ({overline{D}}_{{rm{g}}}^{mu }) and standard deviation ({overline{D}}_{{rm{ln}}}^{sigma }) of the grain size distribution, which is found to follow a log-normal distribution. The two RAMPs are defined as:

$${overline{D}}_{rm{g}}^{mu }=frac{1}{{n}_{rm{G}}}{sum limits_{i = 1}^{{n}_{rm{G}}}}{D}_{i}$$

(19a)

$${overline{D}}_{{rm{ln}}}^{sigma }=sqrt{frac{1}{{n}_{rm{G}}}{sum limits_{i = 1}^{{n}_{rm{G}}}}{left[ln({D}_{i})-{overline{D}}_{{rm{ln}}}^{mu }right]}^{2}},quad quad {rm{where}}$$

(19b)

$${overline{D}}_{{rm{ln}}}^{mu }=frac{1}{{n}_{{rm{G}}}}{sum limits_{i = 1}^{{n}_{rm{G}}}}ln({D}_{i})$$

(19c)

Here Di is the equivalent diameter of ith grain in the M-SERVE consisting of nG grains, ({overline{D}}_{{rm{g}}}^{mu }) is the average grain size and ({overline{D}}_{{rm{ln}}}^{mu }) and ({overline{D}}_{{rm{ln}}}^{sigma }) 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 directions9. 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α ~ 3lβ 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, tension-compression asymmetry, cyclic hardening, strain-rate dependency etc. Corresponding to the crystal plasticity relations, the homogenized elastic-plastic 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 (overline{{mathbb{C}}}) in ref. 6. The same symmetry is assumed for both primary α and transformed β phase M-SERVEs in this study.

The homogenized stiffness tensor (overline{{mathbb{C}}}) for the polycrystalline M-SERVE relates the macroscopic elastic Green-Lagrange strain ({overline{{boldsymbol{E}}}}^{{rm{e}}}) to the macroscopic second Piola-Kirchhoff (2nd P-K) stress (overline{{boldsymbol{S}}}) through the relation:

$$overline{{boldsymbol{S}}}=overline{{mathbb{C}}}left({{boldsymbol{I}}}^{{rm{tex}}}right):{overline{{boldsymbol{E}}}}^{{rm{e}}},$$

(20)

where the macroscopic elastic Green-Lagrange strain tensor is defined as:

$${overline{{boldsymbol{E}}}}^{{rm{e}}}=frac{1}{2}({overline{{boldsymbol{F}}}}^{{{rm{e}}}^{{rm{T}}}}{overline{{boldsymbol{F}}}}^{{rm{e}}}-{boldsymbol{I}}).$$

(21)

The macroscopic elastic deformation gradient ({overline{{boldsymbol{F}}}}^{{rm{e}}}) is obtained from the relation ({overline{{boldsymbol{F}}}}^{{rm{e}}}=overline{{boldsymbol{F}}} {{overline{{boldsymbol{F}}}}^{{rm{p}}}}^{-{rm{1}}}), involving the total deformation gradient (overline{{boldsymbol{F}}}) and the plastic deformation gradient ({overline{{boldsymbol{F}}}}^{{rm{p}}}). The texture tensor Itex in Eq. (16) is used to express the crystallographic dependency of the elastic stiffness in material symmetry coordinate system ({overline{{mathbb{C}}}}_{{rm{mat}}}) as:

For primary α-phase M-SERVEs:

$${overline{{mathbb{C}}}}_{{rm{mat}},{ijkl}}({{boldsymbol{I}}}^{{rm{tex}}})=mathop{sum }limits_{alpha =1}^{3}{g}_{alpha }{{mathbb{C}}}_{ijkl}^{{hcp}}({{bf{v}}}_{alpha }).$$

(22)

For transformed β-phase M-SERVEs:

$${overline{{mathbb{C}}}}_{{rm{mat}},{ijkl}}({{boldsymbol{I}}}^{{rm{tex}}})=mathop{sum }limits_{alpha =1}^{3}{g}_{alpha }{{mathbb{C}}}_{ijkl}^{{rm{Eq}}}({{bf{v}}}_{alpha }),$$

(23)

where ({{mathbb{C}}}_{ijkl}^{{rm{Eq}}}({{bf{v}}}_{alpha })) is the equivalent single crystal elastic stiffness tensor in transformed β colonies. It is obtained from the single crystal stiffness tensor ({{mathbb{C}}}_{ijkl}^{{hcp}}({{bf{v}}}_{alpha })) of the hcp phase and ({{mathbb{C}}}_{ijkl}^{{rm{bcc}}}({{bf{v}}}_{alpha })) of the bcc phase as:

$${{mathbb{C}}}_{ijkl}^{{rm{Eq}}}({{bf{v}}}_{alpha })={v}_{{rm{f}}}^{{hcp}}{{mathbb{C}}}_{ijkl}^{{hcp}}({{bf{v}}}_{alpha })+{v}_{{rm{f}}}^{{bcc}}{{mathbb{C}}}_{ijkl}^{{bcc}}({{bf{v}}}_{alpha }).$$

(24)

The non-zero components of the single crystal stiffness tensors from ref. 13 are given in Table 7, where v1 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 simulation-based elastic stiffness components are obtained from the derivatives of the homogenized stresses for prescribed perturbations in each component of the strain tensor6. 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:

$${{boldsymbol{R}}}^{{rm{mat}}}(t)={overline{{boldsymbol{R}}}}^{{rm{e}}}(t){{boldsymbol{R}}}^{{rm{mat}}}(0)quad {rm{where}}quad {overline{{boldsymbol{R}}}}^{{rm{e}}}(t)={overline{{boldsymbol{F}}}}^{{rm{e}}}(t){overline{{boldsymbol{U}}}}^{{{rm{e}}}^{-1}}(t)$$

(25)

Rmat(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:

$${{boldsymbol{R}}}^{{rm{mat}}}(0)={R}_{alpha beta }^{{rm{mat}}}{{bf{e}}}_{alpha }otimes {{bf{e}}}_{beta }={{bf{v}}}_{beta }cdot {{bf{e}}}_{alpha }({{bf{e}}}_{alpha }otimes {{bf{e}}}_{beta }).$$

(26)

#### 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 ({overline{{boldsymbol{F}}}}^{{rm{p}}}) is obtained by integrating the macroscopic plastic velocity gradient ({overline{{boldsymbol{L}}}}^{{rm{p}}}), which is additively decomposed into the rate of deformation ({overline{{boldsymbol{D}}}}^{{rm{p}}}) and plastic spin ({overline{{boldsymbol{W}}}}^{{rm{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:

$${overline{{boldsymbol{L}}}}^{{rm{p}}}approx {overline{{boldsymbol{D}}}}^{{rm{p}}}={dot{overline{{boldsymbol{F}}}}}^{{rm{p}}}{overline{{boldsymbol{F}}}}^{{{rm{p}}}^{-1}}={D}_{0}^{{rm{p}}}{left(frac{Y}{{Y}_{0}}right)}^{frac{1}{m}}overline{{boldsymbol{N}}},$$

(27)

where ({D}_{0}^{{rm{p}}}) and m correspond to the reference strain-rate and the rate-sensitivity parameter respectively. (overline{{boldsymbol{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:

$$overline{{boldsymbol{N}}}=frac{partial Y}{partial overline{{boldsymbol{S}}}}/leftVert frac{partial Y}{partial overline{{boldsymbol{S}}}}rightVert,$$

(28)

where (leftVert cdot rightVert) is the norm, defined in Eq. (37). Y is the effective transformed stress for an anisotropic yield function48 that accounts for tension-compression asymmetry and is expressed as:

$$Y={left[{left(| {lambda }_{1}| -k{lambda }_{1}right)}^{a}+{(| {lambda }_{2}| -k{lambda }_{2})}^{a}+{(| {lambda }_{3}| -k{lambda }_{3})}^{a}right]}^{1/a},$$

(29)

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 (overline{{{sum }}}), given as:

$$overline{{{sum }}}={mathbb{L}}:left(overline{{boldsymbol{S}}}-overline{{boldsymbol{chi }}}right),$$

(30)

where ({mathbb{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 Itex. The matrix form of the anisotropic tensor is chosen to be49:

$${mathbb{L}}=left[begin{array}{cccccc}{alpha }_{11}&{beta }_{12}&{beta }_{13}&0&0&0\ {beta }_{12}&{alpha }_{22}&{beta }_{23}&0&0&0\ {beta }_{13}&{beta }_{23}&{alpha }_{33}&0&0&0\ 0&0&0&{gamma }_{12}&0&0\ 0&0&0&0&{gamma }_{13}&0\ 0&0&0&0&0&{gamma }_{23}end{array}right]$$

(31)

Plastic incompressibility condition leads to the following relation between the components:

$$begin{array}{l}{beta }_{12}=frac{1}{2}({alpha }_{33}-{alpha }_{22}-{alpha }_{11}),,,{beta }_{13}=frac{1}{2}({alpha }_{22}-{alpha }_{33}-{alpha }_{11}),\ {beta }_{23}=frac{1}{2}({alpha }_{11}-{alpha }_{22}-{alpha }_{33})end{array}$$

(32)

Evolution of the back stress (overline{{boldsymbol{chi }}}) is governed by a modified Armstrong-Frederick type non-linear kinematic hardening law given as:

$$dot{overline{{boldsymbol{chi }}}}=C{overline{{boldsymbol{D}}}}^{{rm{p}}}-left(D-Eleft[exp left(leftlangle frac{-overline{{boldsymbol{chi }}}:overline{{boldsymbol{N}}}}{G}rightrangle right)-1right]right)| | {overline{{boldsymbol{D}}}}^{{rm{p}}}| | overline{{boldsymbol{chi }}}.$$

(33)

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 ({overline{OMA}}_{alpha beta }), ({overline{D}}_{{rm{g}}}^{mu }) 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 Eureqa30 that is based on genetic programming50. The code uses calibrated anisotropy coefficients (αii, γij) and the RAMPs (({overline{OMA}}_{alpha beta }), ({overline{D}}_{{rm{g}}}^{mu }) 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 ({mathbb{L}}) in Eq. (31) is expressed in terms of the RAMP ({overline{OMA}}_{alpha beta }) 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 ({mathbb{L}}) depends on the average grain size ({overline{D}}_{{rm{g}}}^{mu }) and the β lath size lβ, in addition to ({overline{OMA}}_{alpha beta }). 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.

• Three different M-SERVEs (M1-M3) with different grain size distributions, having mean grain sizes (({overline{D}}_{{rm{g}}}^{mu })) of 3.5, 25.9 and 70.5 μm, are created. Lath sizes lα and lβ in these M-SERVEs are set to 1 and 0.35 μm respectively. Nineteen different crystallographic textures with different symmetries given in ref. 6 are created for the above three sets of M-SERVEs to capture the effect of orientation on polycrystalline anisotropy. The combination of orientations and mean grain sizes are chosen such that both soft and hard slip systems are activated, thus giving rise to length-scale dependent anisotropy. This set of M-SERVEs consists of a total of 57 different combinations of grain size and crystallographic texture distributions.

• 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 (frac{{l}_{alpha }}{{l}_{beta }} sim 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 simulation-based 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 stresses6 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 for different grain size distributions and lath sizes. The accuracy of the chosen yield function to describe the anisotropy of M-SERVEs of both primary α and transformed β phases is shown in Figs. 9b and Fig. 10a respectively. The values of a and k are calibrated to be −0.1294 and 10 respectively for the transformed β phase. The functional forms of the anisotropy coefficients for M-SERVEs of primary α and transformed β are given in Box 1. The indices used in these equations preserve the objectivity of anisotropy coefficients under arbitrary permutations of material symmetry axes labels6. The machine learning-generated functional forms show interactions between length-scale parameters such as ({overline{D}}_{{rm{g}}}^{mu }), lβ and orientations represented by ({overline{OMA}}_{ij}) in the anisotropy coefficients. Physical insights may thus be gained on the influence of relevant RAMPs on PHCM coefficients representing specific properties.

#### Anisotropic hardening model in PHCM

Macroscopic anisotropic hardening response of polycrystalline microstructures has its roots in the anisotropy 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 Y0 in Eq. (27) is represented by a microstructure-dependent strain hardening law, which is formulated as a modified Voce-type law51 as:

$${Y}_{0}={widetilde{Y}}_{0}+hat{alpha }mathrm{exp}left(-{left(frac{{overline{epsilon }}_{{rm{p}}}}{hat {beta}}right)}^{psi }right)+overline{H}left(1-mathrm{exp}(-b{overline{epsilon }}_{{rm{p}}})right),$$

(36)

here ({widetilde{Y}}_{0}) is the initial flow stress and ψ is a constant, chosen to be 0.75 in this work. (hat{alpha }) and (hat{beta }) are microstructure and loading direction-dependent parameters that are used to capture the transient yield-point phenomenon. (overline{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 surface52. However in the present work, it is taken to be a calibrated constant for simplicity. The effective plastic strain is defined as:

$${overline{epsilon }}_{{rm{p}}}={int nolimits_{0}^{t}}parallel {overline{{boldsymbol{D}}}}^{{rm{p}}}parallel {dt}={int nolimits_{0}^{t}}sqrt{frac{2}{3}{overline{{boldsymbol{D}}}}^{{{rm{p}}}^{prime}}:{overline{{boldsymbol{D}}}}^{{{rm{p}}}^{prime}}}mathrm{d}t,$$

(37)

where ({overline{{boldsymbol{D}}}}^{{{rm{p}}}^{prime}}={overline{{boldsymbol{D}}}}^{{rm{p}}}) is the deviatoric part of ({overline{{boldsymbol{D}}}}^{{rm{p}}}).

To capture the microstructure and loading direction dependency of the hardening parameters, a dimensionless parameter RY 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:

$${R}_{{rm{Y}}}=frac{Yleft({mathbb{L}}:(overline{{boldsymbol{S}}}-overline{{boldsymbol{chi }}})right)}{Yleft({{mathbb{I}}}^{{rm{dev}}}:(overline{{boldsymbol{S}}}-overline{{boldsymbol{chi }}})right)}.$$

(38)

Microstructural dependency of RY in the Eq. (38) comes from the orientation dependency of ({mathbb{L}}). The loading direction dependency is accounted for through the second P-K stress tensor (overline{{boldsymbol{S}}}). For the transformed β phase however, RY also depends on the mean grain size (({overline{D}}_{{rm{g}}}^{mu })) and lath thickness (lβ) through ({mathbb{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 microstructure-dependent flow stress is given as:

$${Y}_{0}={overline{Y}}_{0}left(tright)+hat{alpha }exp left(-{left(frac{{overline{epsilon }}_{{rm{p}}}}{hat{beta }}right)}^{psi }right),$$

(39)

and the evolution of ({overline{Y}}_{0}) is given as:

$${dot{overline{Y}}}_{0}(t)=overline{H}parallel {overline{{boldsymbol{D}}}}^{{rm{p}}}parallel ,{rm{and}},{overline{Y}}_{0}(0)={widetilde{Y}}_{0},$$

(40)

where (overline{H}) is the microstructure-dependent hardening-rate expressed as:

$$overline{H}={overline{H}}_{0}{left|1-frac{{Y}_{0}}{{Y}_{{rm{sat}}}}right|}^{overline{r}}mathrm{sign}left(1-frac{{Y}_{0}}{{Y}_{{rm{sat}}}}right),$$

(41)

({overline{H}}_{0}) represents the initial hardening-rate and (overline{r}) is an exponent controlling the rate of hardening. Ysat is the saturation stress given as:

$${Y}_{{rm{sat}}}={overline{Y}}_{{rm{sat}}}^{{rm{Ref}}}{left(frac{parallel {overline{{boldsymbol{D}}}}^{{rm{p}}}parallel }{{D}_{0}^{{rm{p}}}}right)}^{overline{n}},$$

(42)

where ({overline{Y}}_{{rm{sat}}}^{{rm{Ref}}}) is a reference saturation stress and (overline{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 Wt 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 Eureqa30 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.

• Three M-SERVEs (M1-M3) with different grain size distributions are created. They have mean grain-sizes of ({overline{D}}_{{rm{g}}}^{mu }=3.5 upmu mathrm{m},25.9 upmu mathrm{m}) and 70 μm and standard deviation of 0.11. The α and β lath sizes for these M-SERVEs are taken to be 1 and 0.35 μm respectively. Nineteen different crystallographic textures with different symmetries are assigned to the M-SERVEs. For a few of the M-SERVEs, crystallographic orientations are swapped between grains to obtain three different misorientation distributions. The M-SERVEs are subjected to constant strain-rate, tensile and compressive loading of 0.001 s−1, along three material symmetry axes as shown in 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/4-cycle with a peak strain of 3% and a time period of 120 s.

• For each of the three grain-size distributions, four different M-SERVEs are created with grain-size standard deviations of ({overline{D}}_{{rm{ln}}}^{sigma }=0.12,0.19,0.31,0.39). Constant strain-rate and cyclic simulation tests, discussed above, are performed on these M-SERVEs to extract the dependency of hardening parameters on ({overline{D}}_{{rm{ln}}}^{sigma }).

• For M-SERVEs of transformed β, lath-size dependency of hardening parameters is determined by simulating the M-SERVE M3 with two different lath sizes of lβ = 0.8 and 3 μm, having a fixed lath-size ratio (frac{{l}_{alpha }}{{l}_{beta }}=3). A few of the 19 different crystallographic textures are assigned to these two M-SERVEs. The M-SERVEs are simulated under the constant strain-rate and cyclic loadings.

• 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, strain-controlled 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 toolkit30 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:

• Model reduction error: The functional forms of constitutive parameters in the PHCMs, e.g., those expressed in Box 1 and Box 2, are in essence reduced order models of the true microstructure-dependency that is implicit in the micromechanical model. A certain level of model 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 (overline{{boldsymbol{sigma }}}left(tright)).

• 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-space14. 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 Taylor-expansion-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:

$${boldsymbol{theta }}left({{bf{x}}}_{{rm{R}}}right)={bf{C}}{boldsymbol{phi }}left({{bf{x}}}_{{rm{R}}}right)+{{boldsymbol{varepsilon }}}_{{rm{m}}},,,{rm{or}}$$

(45)

$${theta }_{i}left({{bf{x}}}_{{rm{R}}}right)={sum limits_{k}^{{n}_{rm{b}}}}{C}_{ik}{phi }_{k}left({{bf{x}}}_{{rm{R}}}right)+{varepsilon }_{{rm{m}},{rm{i}}},$$

(46)

where ({boldsymbol{theta }}=left{{theta }_{i}right}) are the microstructure-dependent PHCM constitutive parameters identified in equations given in Boxes 1 and 2. ({boldsymbol{phi }}left({{bf{x}}}_{{rm{R}}}right)=left{{phi }_{k}left({{bf{x}}}_{{rm{R}}}right)right}), k = 1, . . ., nb are the basis functions identified by machine learning in terms of the microstructure-based RAMPs xR, and ({bf{C}}=left[{C}_{ik}right]) are random coefficients identified by Bayesian inference from CPFE-based calibration data-set. ({{boldsymbol{varepsilon }}}_{{rm{m}}}=left{{varepsilon }_{{rm{m}},i}right}) 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 ({sigma }_{{rm{m}},{rm{i}}}^{2}), i.e.,

$${{boldsymbol{varepsilon }}}_{{rm{m}},{rm{i}}} sim {mathcal{N}}left(0,{sigma }_{{rm{m}},{rm{i}}}^{2}right).$$

(47)

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 ({mathcal{D}}), is used for the Bayesian inference of the stochastic PHCMs. The Bayesian update for the PHCM coefficients C is expressed as:

$$pleft({bf{C}}| {mathcal{D}}right) ,text{}propto text{}, pleft({mathcal{D}}| {bf{C}}right)pleft({bf{C}}right),$$

(48)

where (pleft({bf{C}}right)) and (pleft({bf{C}}| {mathcal{D}}right)) are respectively the prior and posterior distributions of the coefficients. (pleft({mathcal{D}}| {bf{C}}right)) is the joint probability of producing the data ({mathcal{D}}), given the PHCM coefficients C, which forms the likelihood function ({mathcal{L}}left({bf{C}};{mathcal{D}}right)=pleft({mathcal{D}}| {bf{C}}right)) on the parameter space of C. Each evaluation of ({mathcal{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 (pleft({bf{C}}right)propto 1) are used to avoid subjective bias of the posterior distribution. The posterior distribution (pleft({bf{C}}| {mathcal{D}}right)) is sampled using the Hamiltonian Monte Carlo technique53 with No-U-Turn sampler54. The uncertainty due to calibration data sparsity is entirely contained in the Bayesian posterior distribution (pleft({bf{C}}| {mathcal{D}}right)). 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 ({{boldsymbol{theta }}}_{{rm{YF}}}=left{{alpha }_{ii},{gamma }_{ij}right}) of the matrix ({mathbb{L}}) in Eq. (31). The yield function calibration data-set ({{mathcal{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 ({{mathcal{D}}}_{{rm{YF}}}) (markers) are used to identify the hyper-parameters ({sigma }_{{rm{m}},{rm{i}}}^{2}) 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 ({{mathcal{D}}}_{{rm{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 ({{alpha }_{ii},{gamma }_{ij}}), hardening parameters ({{widetilde{Y}}_{{rm{0}}},hat{alpha },hat{beta },overline{H}}), etc. are random variables. This is due to the randomness of the form coefficients C, the model reduction error terms εm, and the RAMPs xR. Consequently, all material state and response variables like ({bar{varepsilon }}_{{rm{p}}}), ({overline{{boldsymbol{F}}}}^{{rm{p}}}), (bar{boldsymbol{sigma} }) 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:

$${{boldsymbol{xi }}}_{(n+1)}={mathcal{I}}left[{{boldsymbol{xi }}}_{(n)};{t}_{(n)},{t}_{(n+1)}right],$$

(49)

where ({boldsymbol{xi }}=left{{xi }^{i}right}) is an evolving state vector. It contains all microstructure-dependent constitutive parameters ({boldsymbol{theta }}=left{{theta }_{i}right}), as well as the material state variables like the plastic deformation gradient ({overline{{bf{F}}}}^{p}), back-stress (overline{{boldsymbol{chi }}}) and effective plastic strain ({bar{varepsilon }}^{p}). The state vector ξ may be generalized to a stochastic vector with a time-dependent joint probability distribution (pleft({boldsymbol{xi }};tright)). The mean and covariance of the stochastic state vector ξ are respectively defined as:

$${boldsymbol{mu}}_{(n)}={left{{mu }^{i}right}}_{(n)}quad {text{and}}quad {boldsymbol{Sigma}_{(n)}}={left[{boldsymbol{Sigma} ^{ij}}right]}_{(n)},$$

(50)

for the time step t(n). The non-linear operator ({mathcal{I}}left[{{boldsymbol{xi }}}_{(n)}right]) 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:

$${mu }_{(n+1)}^{i}={overline{{mathcal{I}}}}^{i}+frac{1}{2}{overline{{mathcal{I}}}}_{km}^{i}{Sigma }_{(n)}^{km},$$

(51)

$${Sigma }_{(n+1)}^{ij}={overline{{mathcal{I}}}}_{k}^{i}{overline{{mathcal{I}}}}_{m}^{j}{Sigma }_{(n)}^{km}+frac{1}{4}{overline{{mathcal{I}}}}_{km}^{i}{overline{{mathcal{I}}}}_{rs}^{j}left({Sigma }_{(n)}^{kmrs}-{Sigma }_{(n)}^{km}{Sigma }_{(n)}^{rs}right),$$

(52)

where ({{overline{{mathcal{I}}}}^{i}}={mathcal{I}}[{{boldsymbol{mu }}}_{(n)}]) is the time-integrator evaluated with the mean state μ(n) at time t(n). The quantities ({{mathcal{I}}}_{k}^{i}), ({{mathcal{I}}}_{km}^{i}) etc. are partial derivatives of the time-integrator, i.e.,

$${overline{{mathcal{I}}}}_{k}^{i}={left.frac{partial {xi }_{(n+1)}^{i}}{partial {xi }_{(n)}^{k}}right|}_{{{boldsymbol{mu }}}_{(n)}}quad ,text{and},,{overline{{mathcal{I}}}}_{km}^{i}={left.frac{{partial }^{2}{xi }_{(n+1)}^{i}}{partial {xi }_{(n)}^{k}partial {xi }_{(n)}^{m}}right|}_{{{boldsymbol{mu }}}_{(n)}},$$

(53)

evaluated at ξ(n = μ(n). The fourth order central moment ({Sigma }_{(n)}^{kmrs}) is approximated using Isserlis’ theorem55, based on ({Sigma }_{(n)}^{km}). In the limit ({Sigma }_{left(0right)}^{km}=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.