Crystal plasticity based constitutive model for deformation in metastable β titanium alloys

Due to attractive mechanical properties, metastable β titanium alloys have become very popular in many industries including aerospace, marine, biomedical, and many more. It is often the complex interplay among the different deformation mechanisms that produces many of the sought-after properties, such as enhanced ductility, super-elasticity, and shape memory effects. Stress induced martensitic transformation is an important deformation mechanism for these alloys. Understanding of it and the influence it has on the microstructural evolution of materials is of great importance. To this end we have developed a crystal plasticity based constitutive model which accounts for both martensitic phase transformation and slip based plasticity simultaneously in metastable β titanium alloys. We present a new formulation for the evolution of martensite transformation, based on physical principles and crystal plasticity theory. To understand and demonstrate this feature of the model, a parametric assessment of the newly developed constitutive model is conducted. This is followed by first of its kind analyses of stress induced martensitic transformation in metastable β titanium alloys. We firstly present validations against uniaxial loading experiments for different metastable β titanium alloys exhibiting stress induced martensite transformation. As part of this, single crystal simulations in metastable β titanium alloys are used for the first time to investigate the interaction of individual transformation systems during unconstrained transformation. This study shows good agreement between the experimental and simulated responses during all stages of deformation in which elastic, transformation and finally the slip stage are exhibited. Relatively ‘strong’ and ‘weak’ orientations for transformation are observed, consistent with experimental studies. The work done here demonstrates the ability of this crystal plasticity finite element method to capture physical mechanisms while bringing new insight about the interaction of different deformation mechanisms in metastable β titanium alloys.


Introduction
Titanium alloys find wide-ranging application in many industries including aerospace, marine, biomedical, and many more.One class of titanium alloys, metastable β titanium alloys, have become very popular due to the superior mechanical properties and workability [1].These metastable β titanium alloys can display several reversible or irreversible deformation mechanisms such as dislocation slip, twinning, the formation of stress-induced martensite (SIM), and a combination of these.The deformation mechanisms are largely a function of the β phase stability [2].It has been shown that SIM formation in this class of alloy can produce an improved compromise between strength and ductility [3][4][5], while the transformation from body centred cubic (BCC) β to orthorhombic α ′ ′ martensite can produce shape-memory and super-elasticity effects in titanium alloys [6][7][8][9].Understanding of martensitic transformations in this context is therefore of great importance.It is well understood that the α ′ ′ martensite can form in six equivalent crystallographic variants [10].However, the exact deformation conditions for SIM, including the influence of different alloying elements, microstructural evolution, the interaction between SIM and elastic/plastic deformation, and failure of the material are not fully understood.Numerous experimental studies have been conducted on metastable β titanium alloys for this reason [3][4][5][6].Earlier investigations mainly concentrated on the conditions to induced SIM [4,5,[11][12][13][14][15][16][17] and the macroscale implication of SIM, such as improved strength and ductility [5,17,18], fracture toughness [3,19,20], work hardening [13,18], super-elasticity [21], and shape memory effects [6,22].The conditions to induced SIM are usually quantified with the trigger stress.This is the point at which SIM is initiated in metastable β titanium alloys [11,12,17,21,23].Trigger stress was found to rise linearly with logarithmic strain rate [16,25,26].A dependence on the loading direction has also been reported [27,28].While recent studies have focused more on the microstructural evolution of the material and the detailed interaction with SIM [10,[28][29][30][31][32][33][34][35][36][37][38].Studies have shown, initially SIM is found to nucleate at β grain boundaries or at hard inclusions in the β grain matrix such as α phase [10,28,29,31,34].A strong martensite variant selection is often noted with initially formed martensite [28,29,37].Furthermore, when SIM is a dominant deformation mechanism many authors have had success predicting the observed, initially nucleated, martensite variants using available work and Schmid factor type analyses [10, 28, 33-35, 37, 38].However, such methods provide poor justifications for the variants formed in the later stages of deformation or when SIM is not dominant [10,31,33,34].This is attributed to interactions between transformation products and complex local stress fields caused by SIM [33][34][35].Being able to computationally model the entire deformation process, including martensitic transformation, is critical in the analysis of these materials.Allowing for improved understanding of physical mechanisms observed during microstructural evolution.
There have been many constitutive models proposed for martensitic transformations in different metals and alloys, which are mostly limited to transformation induced plasticity (TRIP) steels and shape memory alloys (SMAs) (for details please see references [39][40][41][42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57]).It is evident from these references that the development of constitutive models for martensitic transformation has largely been for TRIP steels and SMA.That being said, Grujicic and Sankaran [58] developed a constitutive model for metastable β phase in Ti-1023 which takes into account the basic thermodynamic and kinetic aspects of deformation-induced martensitic transformation and slip plasticity.However, the transformation in this formulation is modelled in a strongly phenomenological nature, as a lattice-less continuum with no account for the crystallography of martensitic transformation.Grujicic and Zhang [59] attempted to address this, developing a single crystal constitutive model combining the basic theories for crystallography, thermodynamics and kinetics of martensitic transformation in metastable β titanium alloy Ti-1023.The model was limited to martensitic transformation only and does not include any plastic deformation mechanisms.Hence, there is a need for a combined thermodynamically consistent model which accounts for the crystallography of martensitic transformation and slip based plasticity for metastable β titanium alloys.Furthermore, modern metal forming processes often utilise finite element methods (FEMs) as an optimisation tool [60,61].Optimisation of process parameters in metal forming is a vital task to reduce manufacturing cost; and is dependent on understanding the materials deformation behaviour during metal forming.The correctness of the simulated material behaviour is governed by the accuracy of the constitutive equations [61,62].With the continuing popularity of metastable β titanium alloys the need for constitutive models to represent the various deformation processes is of ever more importance.This is precisely the aim of this work.In this study, a new crystal plasticity FEM (CPFEM) formulation including SIM transformation along with conventional slip based plastic deformation is presented.The formulation takes account of elastic anisotropy of both the β and α ′ ′ phases, and crystal orientational effects.A novel implementation of the SIM evolution is derived here in a comprehensive manner, allowing for calibration of the model against accessible experimental data.With the intention of allowing the investigation of complex interaction between SIM, slip-based plasticity, damage nucleation and growth in metastable β titanium alloys.Furthermore, providing new insights into the influence SIM has on the bulk material response and microstructural evolution of metastable β titanium alloys during deformation.
The outline for this paper is as follows.In section 2, the kinematics and thermodynamics of local deformation mechanisms inside a representative volume element are presented.A parametric assessment of the constitutive equations for SIM is then conducted in section 3.Sections 4 gives model calibration and validation against available experimental data.Finally, a summary of the work in section 5.

Model formulation
We start the model description with the kinematics, decomposing the deformation gradient into elastic, transformation, and plastic components.Following this, the thermodynamic framework is outlined and dissipated energy and driving forces are defined.We introduce a new phenomenological representation for the transformation evolution based on classical crystal plasticity theory.
The main feature of this model is that the inelastic deformation of SIM takes place by a transformation of the parent material with the body centred cubic crystal structure into the orthorhombic structure of martensite.The model accounts for the basic crystallographic lattice structure of the martensite transformation which can form in several crystallographicequivalent variants [63].In addition, the material may also undergo elastic stretch and rotation, as well as conventional slip-based plastic deformation.
The notation used throughout this section is as follows.Scalar quantities are represented by lowercase letters/symbols.Boldface lowercase letters represent vectors, the dyadic product of two vectors is a ⊗ b and the contraction of a second-order tensor with a vector as Aa.Secondorder tensors are written as boldface capital letters.Tensor operations between two secondorder tensors are represented as AB for the inner product (a second-order tensor), A ⊗ B for the dyadic product (a fourth-order tensor), and A:B for the scalar product.The transpose, A T , of a tensor A is denoted by a superscript 'T', and the inverse, A −1 , by a superscript '−1'.Finally, fourth-order tensors are scripted capital letters (i.e.A, B).The contraction operation over two indices between a fourth-order tensor and a second order tensor is denoted as A : B. Any deviations or additional notation/operation used will be indicated in the text or made clear from the context.

Kinematics
The kinematics of crystalline materials can be expressed using the total deformation gradient, F. This can be decomposed into the elastic, plastic, and SIM transformation parts where F p and F tr are the conventional (slip based) plastic and transformation deformation gradients respectively.Please note that both slip based plasticity and transformation are nonrecoverable as reported previously in the literature [52,62].The elastic deformation gradient (F e ) is decomposed here into the left symmetric elastic stretch tensor V e and rigid body rotation R e .It must be noted that breaking down the material's behaviour (decomposition) does not necessarily correspond to the exact sequence in which different subdomains of the material actually deform under stress (loading path).Also, we have assumed that equation (2.1.1)captures the average effect of all the smaller-scale elastic, plastic, and transformation deformations happening within the mate`arial (see [52,62] and references there in).
Figure 1 shows a visual kinematic representation of the deformation process present in this formulation.It shows the material element moving from the reference state B 0 to its final state B, via three intermediate configurations.These intermediate configurations represent the contribution to the total deformation gradient from the individual component of equation (2.1.1).The relationship between the gradients of deformation and velocity is given as which can be expressed as, Alternatively, equation (2.1.3)can be expressed in the intermediate configuration 3 as where L * describes both the plastic flow due to crystallographic slip and martensitic transformation, and the rotation of the lattice.Therefore, L * can lead to This equation is related to the plastic and transformation parts in intermediate configurations 2 and 1 (figure 1) In equation (2.1.6),the vectors s(i) and n(i) are unit vectors describing the slip direction and the normal to the slip plane of the slip system (i ), respectively.γ(i) is the effective shear strain rate on the same slip system and N s is the total number of active slip systems.In equation (2.1.7),Ltr is the velocity gradient due to transformation at intermediate configuration 1 and ξ(α) is the rate of transformation.In metastable β titanium alloys, the parent β phase (BCC) can be transformed to one of six possible orthorhombic α ′′ 'variants' during martensitic transformation.The lattice correspondence for each of the six martensite variants are given in table 1.Furthermore, it is common in these alloys to observe an internally twinned structure [10,28,34].This pairwise arrangements of two different twin-related α" variants, is termed a transformation system.The theory of martensitic transformations from [64][65][66][67][68][69][70] is used to determine the number and properties of the twinned transformation systems possible.A total of 96 transformation systems are possible for this cubic to orthorhombic transformation.In this current formulation, we adopt the assumption that any or all 96 transformation systems are potentially viable.However, it is essential to acknowledge that future refinements, to improve the computational efficiency, might be made to better align with the observed preferential transformation systems found in these alloys.In a manner akin to what is typically done for NiTi SMA's [49,50,71,72].In this formulation, b(α) 0 and m(α) 0 are the martensite shape strain and habit plane normal vector for each transformation system respectively, and are determined from [65,66].N t is the total number of transformation systems, a total of 96 of which are possible from BCC to orthorhombic transformation.The procedure in the work of Hane and Shield [66] is used to determine these transformation systems.The following conditions must be obeyed for transformation, Here ξ α is the volume fraction for transformation systems (α = 1, 2, . . .N t ) and ξ A denotes the volume fraction of the parent phase measured in the reference configuration.
The velocity gradient can now be broken down into symmetric and skew parts (L = D + W), where D = sym (L) and W = skew (L).D is interpreted as the rate of deformation in the final configuration (2.1.9) while W is interpreted as the rate of spin in final configuration To obtain the rate of the deformation tensor in intermediate configuration 3, D, the definition of L (2.1.4)is combined with D = sym( C e L) to obtain D = V eT D V e , where D is from equation (2.1.9).This is expanded to get Following the work of Marin [73] one can also show where

Thermodynamics
This formulation uses the state variable model of crystal plasticity in the context of thermodynamics.The mechanical version of the Clausius-Duhem inequality per unit volume in the unloaded configuration B is given by Free energy is defined to be a function of applied elastic strain Ẽe and set of internal state variables ϵ s and ξ (α)   ψv = ψv Ẽe , ϵ (i )  s , ξ (α) , (2.2.6) where ϵ s and ξ (α) represent the internal structure of the material during evolution from plastic slip and the SIM transformation.These variables are incorporated to represent the state of the evolving internal structure of the material generated during deformation.The internal variable ϵ (i ) s is a strain-like variable assumed to be related to the elastic lattice strain fields generated on the slip systems by dislocations.This will effectively model the isotropic hardening effects induced by the interaction among dislocations and the variable is connected to the density of dislocations ρ (i ) s by the relation ϵ s where b is the Burger's vector.While ξ (α)  is volume fraction of martensite from the αth transformation system contained in the parent lattice structure.Thus, we can give the material time derivative of ψv as Substituting this back into (2.2.5), rearranging and collecting like terms gives, In view of the decomposition of the total dissipation, the dissipation inequality can be written as and it is further assumed that the dissipation inequality holds for the plastic deformation (D p ) and phase transformation (D tr ) independently, i.e.
D p ⩾ 0 and D tr ⩾ 0. (2.2.10) Dissipation due to slip plasticity and phase transformation are given as The shear rate on slip system (i) is determined using the well-known power law relationship [74,75] where γ0 and m s are the reference shear strain rate and rate sensitivity parameter, respectively.We further assume that dislocations are only generated in the parent (BCC) phase.The effects of transformation on the plastic deformation is considered by weighting the effect of slip by the remaining volume fraction of the parent (untransformed) phase, ξ A [52] The simplified hardening law is taken from Kocks and Mecking [76] is used in this formulation where all systems harden at the same rate.Hence the following can be given:

.2.16)
where h 0 is the initial hardening, κ s,0 is the initial critical resolved shear stress (CRSS) and κ s is the current CRSS.Lastly, κ s,S is the saturation strength given by the relationship of Follansbee and Kocks [77], with κ s,S0 , γ(i) S0 , and m ′ all being material parameters.The internal variable ϵ s taken account of by adopting the following relationship, and the total elastic lattice strain field generated by dislocations (ε s ) is related to the average hardening on each slip system, as postulated by Clayton [78]: We now turn our attention to the elastic contribution to the Helmholtz free energy: and as elastic deformation is reversible energy dissipation is absent.Therefore where Ẽe is the elastic Green-Lagrange strain, conjugate to the second Piola-Kirchhoff stress in the second intermediate configuration.Furthermore, the effective elasticity tensor C is determined by an averaging scheme taking the contributions from the elastic stiffness of both the parent and martensitic phases (C A and C (α) ) [51,52,79,80] where δ tr = b (α) • m (α) is the volumetric growth associated with each transformation system and J tr is the Jacobian of F tr J tr = detF tr4 .Integrating this with respect to Ẽe yields the quadratic form of the elastic strain energy per unit volume in the final intermediate configuration, in which C is an integration constant.It must be noted that the set of internal state variables are incorporated to represent the state of the evolving internal structure of the material generated during plastic slip and transformation, respectively.From above this constant can be assumed to be a function of ϵ (i ) s and/or ξ (α) .Assuming it is a function of both, equation (2.2.23) can be written as and C is replaced with ψ1 (ϵ s , ξ (α) ).Let us turn our attention to the elastic lattice strain fields generated by dislocations.In this model, similar to the work of Tjahjanto et al [52], the kinetics of dislocations is not resolved explicitly and the elastic lattice strain energy density is also assumed to be described quadratically from ϵ ψsp ϵ (i )  s , ξ where µ (ξ ) is the effective shear modulus and c k is a scaling material parameter for the elastic lattice strain energy.Like equation (2.2.22) the effective equivalent shear modulus from equation (2.2.25) is expressed as the homogenised average equivalent shear moduli of the parent and the martensitic phases [52].The equivalent shear moduli is expressed as where µ A and µ (α) are the isotropic shear moduli of the parent material and martensitic transformation system α, respectively.The isotropic shear moduli are calculated from the anisotropic elasticity tensors for both the phases, using the Voigt approximation [81] where C ij is the ith and jth component of the 6 × 6 Voigt notation elasticity tensor.We can now write the partial derivative of ψv [ Ẽe , ϵ s , ξ (α) ] with respect to ϵ Following the work of [51,52] we have two equivalent expressions for the Helmholtz energy Differentiating the above with respect to Ẽe and ϵ ψm Finally, integrating the above results where ψ3 is dependent on ξ (α) .Once again following the example of [51] the Helmholtz free energy is now written as In this model, ψ3 ξ (α) is interpreted as the internal energy barrier resisting transformation accounting for all energy contributions associated with SIM.Returning to the dissipation inequality, the transformation dissipation in equation (2.2.12) can be written as (2.2.38) separating them out gives (2.2.40) Defining the driving force for transformation as and combining this with equations (2.2.23) and (2.2.25) gives It is assumed that the driving force for transformation is related to the total elastic lattice strain fields generated by dislocations from all slip system contributions (ϵ s ).Therefore, the driving force for each transformation system is, Finally, the transformation dissipation can be expressed as As stated, the second term here is interpreted as an internal energy barrier resisting transformation.With an analogy from slip-based crystal plasticity, this term gives the stored work due to the accumulation of the martensite phase within the BCC lattice structure.In ∂ξ (α)   represents the critical mechanical driving force (f (α) cr ) required for transformation.This is an internal energy barrier for transformation and is made up from various contributions such as the nucleation of parent/martensite interfaces, internal frictional resistance due to parent/martensite boundary movement, and the interaction among individual martensite plates/systems [14,39,59,82].Other authors have attempted to include some or all contributions explicitly [49,52,59,80].We have not attempted to explicitly accounted for these mechanisms.Instead, we collated such contributions into f (α) cr and its evolution.The following simple phenomenological formulation for transformation resistance is deduced and used throughout in this work in which, v 0 is a normalising resistance parameter (MPa) and R 0 is the initial resistance (MPa); R s is saturation resistance (MPa).In the second term of equation (2.2.46), R h and η are all material parameters used to increase the resistance to transformation once it has begun due to the interaction between transformation systems.R h is the transformation system interaction parameter and represents the magnitude of the transformation system interaction, while η controls the critical amount of transformation for the interaction to become significant.The initial critical value for driving force is assumed same for all systems.
Figure 2 shows a graphical representation of the transformation resistance model.The figure presents 3 stages of resistance evolution, each controlled by different parts of equation (2.2.46).This allows for better control of the transformation process at different stages of the deformation, depending on the material type and SIM evolution.This relationship can be further simplified and fine-tuned through experiments.The evolution of martensite ( ξ(α) ) is modelled using the following power based, rateindependent, kinetic law similar to the work of Anand and Gurtin [49] where ξ0 is the reference transformation rate, m is a rate sensitivity parameter, f (α) is the driving force for transformation, f cr is the critical value and ξ (α) is the total volume fraction of each system.
Finally, this formulation is implemented using a small elastic strain assumption as a user material subroutine (UMAT) in ABAQUS finite element software.Detailed numerical implementation can be obtained upon request to the corresponding author and are not provided for brevity.

Parametric assessment of the model
In this section, we present the effect of key transformation model parameters on the stressstrain response and martensite volume fraction progression.The parameters of Ti-10-2-3 alloy given in tables 2 and 3 are used as the baseline for this assessment; the transformation parameter values used in this section are as given in table 2 unless stated otherwise.It should be noted that the complete deformation process comprises the elastic, plastic, and transformation parts due to the deviatoric and volumetric contributions.The current constitutive model reflects this, and material parameters are divided into three categories: elastic parameters, slip based plasticity parameters and transformation parameters.Since the purpose of this parametric study is to assess the transformation parameters, the elastic and plastic material parameters are kept constant throughout and are given in table 3. The value of the martensite elastic constant scaling parameter, χ in table 3 was also kept constant.This parameter is a scaling parameter for the martensite elastic constants and is used to elevate the lack of experimental data, see section 4.This study has been conducted for the case of uniaxial loading conditions.All the tests have been performed by varying the following material constants: f cr , ξ0 , R 0 , R s , v 0 , R h , η.To link these individual parameters with physical mechanisms, the influence of each parameter on stress-strain response and martensite volume fraction progression are presented below.

Transformation initiation
The point at which transformation begins is controlled by the critical driving force (f cr ).This point at which SIM is initiated in metastable β titanium alloys is termed the trigger stress [11,12,17,21,23].The main factors which affect the trigger stress are: • β phase stability, • β phase grain domain size, shape, and orientation, • loading conditions.
The stability of the β phase is intrinsically determined by the alloying elements and their volume, i.e. the alloy composition.However, it is also controlled by the thermo-mechanical processing and heat treatments used.It is well known that externally applied stress may alter the driving force for transformation, favouring the formation of variants aligned with the applied stress [10,34,35,63,83].This is accounted for in the derivation of equation 2.2.43.The contribution of this assists the chemical driving force for transformation of the β phase structure; determined by the inherit stability of the alloy and the testing temperature.However, the model presented in section 2 is isothermal and hence the testing temperature does not change during the deformation.The effect of the testing temperature on the trigger stress is therefore not considered.Therefore, except for the loading conditions, f cr comprehensively incorporates all these factors into one value.A value that may be calibrated from experimental data.Figure 3 shows the effect of raising f cr from 3-30 MPa.It is found that increasing f cr delays the transformation (increasing the trigger stress) up to the point where the transformation is incomplete at the end of the simulation (f cr = 30 MPa).As we are considering a full β phase alloy under isothermal conditions, f cr is effectively modelling the β phase stability and the additional mechanical driving force to trigger transformation.
Figure 4 shows that the reference transformation rate is very sensitive in the range ξ0 = 1 × 10 −5 to 5 × 10 −5 .Above this range the change in response is far more subtle.This parameter is effectively the maximum rate of transformation for each system.Figure 5 shows the effect of rate sensitivity (m) on the dependence of transformation rate against the applied driving force for transformation.The transition up to the maximum transformation rate occurs over a range of the driving force for each system, f (α) , up to f cr .The transition is more rapid for smaller m values.In a totally rate independent model, this transition does not exist, and transformation will only occur when f (α) ⩾ f cr .From figure 5 we can see this behaviour is approximated at m ≈ 0.005.

Transformation resistance
The initial resistance to nucleation of martensite is given by R 0 .Figure 6 shows the influence of this parameter of R 0 on the transformation.It can be seen in figure 6 that it controls the initial stage of transformation.Increasing this parameter raises the initial resistance during Stage I, slowing down the initial nucleation of new transformation systems.R s from equation (2.2.46) is termed the saturation resistance (giving the point of maximum rate of transformation for  any system).The two resistance values are related by the volume fraction of the individual transformation system and the normalising parameter v 0 .This results in the interfacial resistance decreasing as the volume fraction of the one transformation system increases, ensuring that the interfacial resistance is at a minimum as a transformation system begins to occupy an entire area.The point at which the saturation resistance (R s ) is reached is also significantly influenced by v 0 (figure 7).The results presented in figure 6 are given using v 0 of 5.0, this to depict the influence more clearly.Effectively, v 0 determines how fast the transformation system moves from the initial resistance R 0 to the saturation resistance R s (R 0 → R s ).This can be inferred from figure 7 where the stress-strain response reaches the stress plateau region faster.For the curves in figure 7, R 0 has been made equal to 500.This is to exaggerate the response of v 0 and better show its influence.The plateau region has been 'shifted' upwards, or Stage I is extended, because of v 0 increasing.As stated, R s gives the minimum resistance to transformation for all individual transformation systems.Figure 8 shows that increasing R s increases the rate of transformation in Stage II.This is seen to equates to a flatter or more prominent stress plateau region.
Finally, the interaction of α ′ ′ plates with each other [14,34,35] and additional deformation mechanism become dominant features of the deformation process.Here, R h and η are material parameters which control the resistance to transformation based on the interaction among systems of martensite.R h represents the magnitude of the interaction resistance (the strength of the interaction among systems) while η controls how quickly this interaction begins, or the point at which the interaction in the crystal becomes significant.This resistance due to  interaction among individual martensite systems is governed by the total volume fraction of martensite, second term in equation (2.2.46).
Figure 9 shows the influence of R h .The point at which the interaction begins remains the same, however the strength of the interaction increases with R h .This is evident in figure 9, with the total volume fraction of martensite increasing more slowly and the stress plateau region ending sooner. Figure 10 shows how increasing η delays the effects of transformation system interaction.Reducing η to 0.5 induces a significant interaction early in the transformation process.This is seen to reduce the rate of transformation (increase the resistance).
The approach presented facilitates the stress plateau and double yield characteristic of SIM deformation in metastable β titanium alloys.

Numerical examples and experimental validations
As explained in the previous sections, the novel transformation model presented is based on the physical observations seen in experiments and hence parameters must be calibrated against experimental data.This has been done by comparing the resulting stress-strain curve from the proposed model with outcomes from room temperature uniaxial tensile experiments [6,21] (figure 12).The transformation parameters so identified are given in table 2. In this work, manual fitting of the experimental stress-strain data determined the material parameters, with the methodology described in the following section below.A single ABAQUS C3D8R element is used in conjunction with a Taylor model homogenisation scheme to represent 200 randomly orientated grains (figure 11).MTEX toolbox [84] is used along with the BCC crystal symmetry to produce the random orientations used.The particular aggregate here consists of 200 initially randomly oriented crystals chosen from a uniform orientation distribution with the range for the randomness was 0 • -360 • for Euler angles.For this case the Taylor model homogenisation produces the volume-averaged Cauchy stress σ T at a continuum material point as, Table 3.
Elastic and plastic parameters.
Elastic parameters Martensite Scaling Parameter (a)   , χ : 0.5 Plastic deformation parameters (BCC Slip) System ms γ0  where, N g is the number of grain orientations considered, and σ (i) is the Cauchy stress of each orientation.The applicability of this method for polycrystalline materials has been previously shown by [50,73,75,[85][86][87].
Elastic constants for both the BCC and orthorhombic α ′ ′ phases are presented in table 3. The BCC values are obtained from data for similar materials in the literature [88][89][90][91][92][93]; which are then adjusted iteratively to match experimental results (keeping within the range identified in the literature).For all the materials, fitting the anisotropic elastic constants of the α ′ ′ phase is a challenge.With the metastable nature of the α ′ ′ phase it is very difficult to experimentally determine these constants for most alloys.In this study we have used the values given by the first principle calculations of Sun et al [91].The properties are further scaled in the formulation to help fit the experimental data.Similar approaches have been used for the SMA CPFEM where the scaled down β elastic constants have been used to represent martensite [47,49,71,94].
As discussed earlier, the {110}<111>, {112}<111> and {123}<111> slip families are all possible in the β phase [95][96][97].The plastic deformation parameters for these were iteratively fitted until a close match to experimental results was achieved.However, the order of activation is held consistent with the literature with {110}<111> and {112}<111> being more dominant and {123}<111> activating later (for details on this please see refs [95][96][97]).This trend is kept consistent and the CRSS (κ s,i ) values are fitted to the data for the individual alloy.The parameters identified for both elastic and plastic (slip) deformation are presented in table 3. The parameter κ 0 is the CRSS for each slip system when shear strains due to plastic slip are zero, while h 0 , κ s,0 , κ s,S0 , m ′ , γS0 are material parameters from equations (2.2.16) and (2.2.17).For c k , in equation (2.2.28), a similar approach to that of [52,54] is adopted and the value is kept constant at 0.5.
The result of the inverse modelling in figure 12 show good agreement between the experimental and simulated responses during all stages of deformation, i.e. elastic, transformation, and slip, successfully capturing the double yield effect so often referred to in the literature for SIM in metastable β titanium alloys [10,15,26].The model has been further validated through single crystal simulations to evaluate the response of a free-standing unconstrained single crystal.Three different directionalities, representing the extreme orientations in a standard stereographic projection plot [8], are analysed.The loading direction and corresponding single crystal Euler angles (ψ , θ, ϕ ) are given in table 4; with a visual representation shown in figure 13.The parameters identified for alloy Ti-10-2-3 and β-Cez are used.All figures in this section are related to Ti-10-2-3 unless stated otherwise.Uniaxial tension tests, followed by volumetric expansion and contraction of the model is presented.From here on, the orientations will be referred to as either [001, 011], or [111].

Effect of uniaxial loading on transformation and slip based plasticity, and their interaction
The true stress-strain curve for the three orientations can be seen in figure 14.For each of the curves presented, to varying degrees one can observe a plateau of stress.This plateau is due to martensitic transformation and would be the first yield point for the material and later plastic deformation would constitute the second yield [10,15,26].
The [011] direction shows the most distinct plateau, both in gradient and extent, suggesting that this direction has a larger transformation strain contribution.This is in keeping with the literature for metastable β titanium alloys [8,98].Kim et al [8] carried out theoretical and experimental calculations of the transformation strain in Ti-Nb alloys.The [011] transformation strain was calculated as almost twice that of the other directions, 4.2% compared to 2.1% in [001] and 1.5% in [−111], while [6,14] reported the maximum transformation strain of c.+7.2% in the [110] direction for T-10-2-3.Additionally, Chen et al [10] found that martensitic variants with < 110> β // < 001> α ′ ′ orientation relationship produce the maximum transformation strain along the tensile direction and are activated first.Pinilla Ducreux et al [98] reported lattice strains for β and α ′ ′ phases during uniaxial deformation.It was found that the (011) β phase plane loaded in the axial direction showed the least amount of lattice strain upon commencement of transformation.The transformation in this direction accounted for nearly all the lattice strain during that stage of deformation.Figure 14 further shows that there can be relatively 'strong' and 'weak' directions for transformation.The strongest resistance to transformation is exhibited by [111], requiring a far larger stress to initiate transformation.On the other hand, the [011] direction is the most favourable and transforms first.The result for the [001] direction in figure 14 lies between the two others, although slightly closer to the response of [011].This is an observation also made experimentally [10,99,100].Cai et al [100] reported experimentally that the [011] orientation was the first to transform in Ti-Nb-Fe alloys; while [99] found that the (111) plane experienced virtually no transformation at room temperature for tensile tests of Ti-27.96Nb-11.97Ta-5.02Zr%wt alloy.
Figure 14 also shows the total volume fraction of martensite for each direction.The first to transform is [011] and has the slowest rate of martensite volume fraction increase.This is a product of a high transformation strain in this direction; requiring less transformation to accommodate the deformation.Finally, figure 14 also shows the interaction of SIM and slip on the three orientations.It is clear for this case, where SIM and slip are largely separated, the interaction is minimal.The [001] and [011] directions show a clear and distinct region of elasticity between the two deformation modes.On the other hand, the [111] direction shows a far less segmented response.Figure 15 shows the true stress-strain and total volume fraction curves in the [011] and [111] orientations for the β-Cez alloy.In this alloy SIM and slip mechanisms are not so separate.This material shows the same broad trends as the ones described by figure 14.There is a strong plastic mismatch between the [011] and [111] orientations, with [011] activating at c. 225 MPa before [111] responds.The interaction between SIM and slip is far greater in this alloy than for Ti-10-2-3, as presented in figure 14.There is no distinct double yield in the [111], indicating that SIM and slip are happening simultaneously.This highlights the 'strength' of this [111] direction for transformation.If we consider polycrystalline response, for an alloy like β-Cez the plastic mismatch and strength of [111] may result in no transformation for this direction.This is possibly the case reported by Maghsoudlou et al [99] where they report that 'virtually' no transformation is seen on the (111) plane during tensile deformation.
The individual transformation systems for each direction are depicted in figures 16 and 17. Figure 16 shows the systems activating in the [001] and [011] directions.As explained in section 2.1, all transformation systems in this context are composed of two variations.Within this pair, each system exhibits one variant as its primary component (dominant variant).For instance, in the case of Ti-10-2-3, the twin combination is established with a distribution of approximately 75% for one variant and 25% for the other 5 .Consequently, we can categorize all the systems into six groupings based on their dominant variant.Table 5 provides the primary (dominant) variant for each of the 96 transformation systems.
Direction [001] shows four distinct groups of transformation systems, the transformation is fairly ordered with the first to activate, having the largest impact on the result.This direction has the greatest number of active systems, 32, compared to 16 in [011], and 22 in [111].The maximum contribution of each transformation system is less than that of the other directions.This suggests that in this direction there are more equally viable transformation systems available.
Direction [011] shows only two groups of systems, the transformation is again fairly ordered with the first to activate, having the largest impact on the result.This direction appears  to have the clearest set of viable systems to accommodate the transformation, i.e. having the lowest number of systems active in only two groups.Inspection of figure 14 and the lattice correspondence reveal that martensite 'wants' to form in the [011] direction.This direction shows less competition for transformation system activation and a clearer path for a single dominant variant (table 5).Comparing figure 16 and table 5 we can see that only transformation systems dominated by variant 6 are activated.The other directions need a combination of expansions in differing directions to achieve the required deformations.
On the other hand, [111] shows a far more complex and chaotic order of transformation.It is clear by comparing the y-axes of figures 16 and 17 that this direction has the largest contribution from any one system.However, there is no clear grouping of systems.There is also a large amount of cross-over between systems, showing that initially favoured systems do not necessarily dominate the final response.This behaviour also suggests that the act of transformation changes the state of the crystal, negatively impacting continued activation of some active systems and enhancing the response of others.Nissen et al [34] commented on not being able to predict the transformation systems that developed later in the deformation process, suggesting that SIM transformation changed the stress state of the grain to an extent that new, previously less responsive transformation systems become activated.However, this was for a polycrystal and changes in stress state could have been due, in part, to intergranular interactions.Observations of Gao et al [28] support our initial statement; they describe the lattice strains formed during transformation having the effect of controlling variant selection for subsequent transformation.This is shown more clearly in the magnified image of figure 17, where transformation system 61 is initially one of the most responsive systems.However, upon its initial transformation its volume fraction begins to saturate at c.0.005.At this point transformation system 26 is activated and becomes the most prominent system in the rightmost magnified plots of figure 17-[111] Direction Magnified.This orientational effect, to the authors' best knowledge, has not been previously investigated in metastable β titanium alloys.

Volumetric response of transformation
Volumetric response of the transformation is investigated in this section.In figures 18 and 19 the volume change is defined as a ratio between the current volume (V) and the initial volume (V ini ).Figures 18 and 19 show that transformation is only present with volumetric expansion.All systems are activated equally.There is no transformation for volumetric contraction.The martensitic transformation of BCC to α ′ ′ in Ti-1023 is accompanied by an expansion and it is found that all systems have the same volumetric expansion (see [65] for details).Hence, there is no transformation in pure contraction and all the systems are activated equally.This provides validation of the model and its implementation.With this new CPFEM formulation for metastable β titanium alloy we investigated the single crystal response of Ti-10-2-3 and β-Cez alloys.This study shows good agreement between the experimental and simulated responses during all stages of deformation in which elastic, transformation and slip are exhibited, successfully capturing the 'double yield effect' to differing degrees for different materials.Furthermore, the analysis shows there are relatively 'strong' and 'weak' directions for transformation.The strongest resistance to SIM is exhibited by [111], requiring the largest stress to initiate.The progression of SIM is found to be strongly dependant on orientation in this metastable β alloy, presenting a strong plastic mismatch between orientations.The [011] direction is the weakest and transforms first with the largest transformation strain.The lattice correspondence (table 1) between β and α ′ ′ phases  shows this to be a 'natural' direction for martensite to form, producing the large transformation strain seen in this direction, an experimental observation that has been remarked upon numerous times [6,8,10,14,98].The [111] orientation shows a far more complex and chaotic order of transformation.Significant cross over among transformation systems suggests that transformation changes the elastic and inelastic (slip and/or SIM) properties of the crystal.Altering the subsequent transformation systems required to accommodate the continued deformation.
This highlights the significance of orientation effects on the progression of SIM in single crystals.This study brings new insights regarding transformation, crystal orientation and slip interactions, while results are also in line with experimental observations.

Conclusions
A new CPFEM based model for SIM transformation and slip based plasticity, in metastable β titanium alloys, has been developed.The model incorporates a fully comprehensive framework for the evolution of martensitic transformation in the material.Microstructural information taken from literature is included using a series of averaging/homogenising techniques.A detailed parametric assessment has been performed to understand and clarify the effect of individual transformation parameters on material response during deformation.Numerical validation of the model has been presented by the manual fitting of experimental stress-strain curves and single crystal analysis.Notable outcomes are as follows: • Good agreement between the experimental and modelled responses during all the stages of the deformation where elastic, transformation and the slip stage was shown.Successfully capturing the 'double yield effect' to differing degrees for different alloys.SIM was found to be strongly dependant on orientation in metastable β alloys, presenting a strong plastic mismatch between orientations.• Less favourable orientations experienced significant cross over between the transformation systems.Indicating the act of transformation changes, the elastic and plastic properties of the crystal; altering the required transformation systems to accommodate the subsequent deformation.
The model presented here has shown capable of capturing physical mechanism associated with SIM, slip based plasticity and their interactions in metastable β titanium alloys at microscale.The need for accurate constitutive models to represent the various deformation processes in these alloys in increasingly important with their continuing popularity.However, the formulation has potential for expansion to incorporate thermal effects during deformation, as these can significantly impact SIM transformation and β phase stability, particularly influenced by temperature.This expansion could pave the way for research into high strain rate and cyclic thermal loading scenarios commonly found in hot metal forming operations.Furthermore, refining the transformation model is necessary; for instance, re-evaluating the assumption that all transformation systems share the same f cr value would enhance the formulation.Additionally, the formulation currently use simplistic approach to the influence of transformation on slip plasticity.Experimental observations suggest that SIM transformation induces a work hardening effect by reducing the mean free path for dislocations [26,97,101,102].However, this dynamic hardening or dynamic path-hall effect is not accounted for in the current model formulation.Engaging in dedicated experimental studies could enrich the transformation evolution law and attribute explicit physical meaning to the transformation parameters.Furthermore, there is a need to establish a more robust, nuanced interaction between SIM transformation and the flow strength of slip-based plasticity.
The next steps in our research involve conducting comprehensive polycrystalline investigations in 3D.One of our main objectives is to perform a direct comparison with experimental data using a complete 3D analysis, which will offer quantitative validation at the polycrystalline level.This forward-looking approach also aims to further enhance the accuracy and reliability of the formulation presented here.In a numerical setting, this rotation matrix is used to transform to the slip system vectors (s 0 ) and the fourth order tensor of elastic moduli C e 0 to the sample coordinate system in configuration B. Hence, The initial orientation C 0 is prescribed as part of the crystal initial state and is numerically updated during deformation process using C = R e C 0 (R e is the rotation tensor from the kinematics).
Martensitic transformation is also regarded as a mode of deformation, one which causes a change in the crystalline structure [12].At the lattice scale, we can describe the transformation from parent lattice to martensitic lattice as a deformation as there is no diffusion.Suppose we have parent lattice vector e p and a martensite lattice vector e m , we can define a matrix to describe the homogeneous deformation such that, This matrix U i is called the 'transformation matrix'.In metastable β titanium alloys, the parent β phase (BCC) can be transformed to one of six possible orthorhombic α ′′ 'variants' during martensitic transformation.Hence, there are six different transformation matrices (U i , i = 1, 2, . . ., 6).Martensitic transformation here consists of a stretch and a rotation of the parent lattice.The size of the stretch is dependent on the lattice parameter (a, b, c) of the material.The lattice correspondence for each of the six martensite variants are given in table 1.
The equation (B.6) are used to calculate the lattice stretches [91,93], where the m subscripted lattice parameters belong to martensite, and a 0 is the parent BCC phase.Equation (B.7) give the transformation matrices for all six possible variants of martensite which can arise in cubic to orthorhombic transformation It is common in metastable β titanium alloys to observe an internally twinned structure [13][14][15].This pairwise arrangements of two different twin-related α ′ ′ variants, is termed a transformation system.The theory of martensitic transformations from [16,[17][18][19][20][21][22], is used to determine the number and properties of the twinned transformation systems possible.A total of 96 transformation systems are possible for this cubic to orthorhombic transformation.In this current formulation, we adopt the assumption that all 96 transformation systems are potentially viable.However, it is essential to acknowledge that future refinements, to improve the computational efficiency, might be made to better align with the observed preferential transformation systems found in these alloys.In a manner akin to what is typically done for NiTi SMA's [23][24][25][26].At the microscale the transformation systems can be described by the habit plane normal vector m (α) 0 and the shape strain vector b (α) 0 .The link between these transformation vectors and the kinematics of transformation can be realised when looking at the habit plane equation from [18], where F α tr and F β tr are the transformation deformation gradient of the α ′ ′ martensite and β parent phase, respectively.The theory of martensitic transformation given in [16,[17][18][19][20][21][22], is developed assuming that the parent phase is stress free.Hence, F β tr = I in (B.8) is a part of the theory of [16,[17][18][19][20][21][22].We assume that the main characteristics of the transformation system remain under applied mechanical loading [27].The rotations R(α) ij and R (α) ij are the habit plane rotation and variant rotation (describing the rotation between variants), respectively.Finally, λ i is the volume fraction of the ith variant in the system.All this information is given by the procedure outlined by [18].The reader is directed to [27] for an excellent explanation of the kinematics of martensitic transformation in this context.
Moving to elasticity, the number of constants needed to specify the anisotropic elasticity tensor C e 0 depends on the crystal structure, where Ce is the effective elasticity matrix in the sample coordinate system.With respect to the crystal axis C e 0 , for both phases, can be written as For the parent BCC phase (3 material constants required), and, For the orthorhombic martensitic phase (9 material constants required).The elastic constants for the parent BCC phase can be taken from literature (table 1).
The elasticity tensor for the martensitic phase is determined using the same method as [27], where, each transformation system will have its own unique elasticity tensor.As stated, the procedure from [18] is used to determine the transformation systems.Firstly, the martensite elasticity matrix in the crystal axis, C e 0,m , must be transformed to match the lattice correspondence for each variant given in table.We denote this rotation as R β * i .Each of these rotations can be combined with the necessary rotations from the transformation systems.This determines the rotations required for each variant in each of the 96 systems.The rotations are as follows, The additional rotation in Q (α) i is to move variant 'i' to the same axis as 'j' before rotating both to the habit plane axis.As stated, the procedure used also provides the volume fraction of each variant in the twinned transformation systems (λ i ).This information is combined with elastic constants for the orthorhombic martensitic phase to determine an effective elasticity tensor for the transformation system, which is the volume averaged combination of each variant in the twin system, The effective elasticity matrix Ce is then determined by equation (2.2.15) of the main text.The constants for C e 0,m can also be taken from literature.However, due to the metastable nature of the α ′ ′ phase there is considerably less data available.Table 2 presents a collection of elastic constants for this phase from the available literature.
We consider Voigt (vector) notation for the stress and strain tensors.Using this notation, one can write the elasticity law in matrix form as, With this definition of elasticity, we also use the vector representation of deviatoric stress and strain, devτ and devϵ e respectively [37,38].

B.2 Constitutive integration scheme
The crystal constitutive model can be considered as a set of coupled first order ordinary differential equations for the variables τ , R e , κ . The time integration of these evolutionary equations is carried out in the sample axis and proceeds by discretisation of the deformation history in time and numerically integrating the equations over each time step.For this purpose, we consider the configurations of the body at time t n and t n+1 with t n+1 = t n + ∆t.Accordingly, variables evaluated at t n and t n+1 will be denoted with the subscripts n and n + 1 respectively.This integration scheme is developed assuming that: i. the crystal deformation represented by L n+1 is given, ii. the variables τ n , R e n , κ i s,n , f α cr,n are known, and iii. the time-independent slip system vectors, transformation system vectors, the elasticity tensors, the initial crystal orientation, and the plasticity/transformation material parameters are input.
The integration of the model will then give the updated values τ n+1 , R e n+1 , κ i s,n+1 , f α cr,n+1 .The numerical integration proceeds as follows,  Turning our attention to the evolution of the rotation tensor, slip hardness and transformation resistance.These quantities are determined, at time t n+1 , using an exponential map [37,38] for R e n+1 and backward Euler schemes for the slip hardness and transformation resistance, as shown below where ).

B.3 Specifics for transformation
Driving force for transformation where using the result of equation (2.2.42) and altering it to account for initial slip hardness, the elastic lattice strain fields generated by dislocations ϵ s,n+1 can be given as: The volume fraction of martensite is: Note that in this formulation we are solving for τ n+1 , R e n+1 , κ i s,n+1 , f α cr,n+1 .By solving these quantities, the volume fraction is also updated consistently.cr,n+1 is also dependant on the volumetric response.We focus our attention to solving equations (B.37)-(B.41).For this purpose, using these equations we can formally write the residuals, A staggered scheme is used in the formulation instead of full Newton Raphson (N-R) scheme.First R e n+1 , κ cr,n+1 are estimated and then first instance of devτ n+1 and p τ,n+1 are solved.Of course, R 1 and R 2 could be incorporated together to solve for τ n+1 , but have been presented this way for transparency and consistency with (B.14).A N-R solution is used to solve τ n+1 at the first instance.The solution for this requires calculation of a local Jacobian matrix.This is derived using a first-order forward finite-difference numerical differentiation scheme of the general form: where, in this general form k, i and j represent the primary variable in X (or residual number), and the components of R k and X k , respectively.Furthermore, h is a perturbation parameter used to disturb the solution on the jth component of unit vector e j .The unit vector e j being of equal size to X k ; the perturbation parameter was taken as 10 −8 [39][40][41].Analytical derivation of the Jacobian was also considered.However, was abandoned due to the difficulties and uncertainties associated with the extensive derivatives required in the calculations.
Once the solution for τ n+1 is obtained, an N-R solution for κ cr,n+1 is determined while keeping (τ n+1 , R e n+1 , κ (i ) s,n+1 ) constant and finally R e n+1 can be found.The initial guess for τ n+1 in this scheme is based on the viscoplastic solution i.e., the elastic strain is taken as zero at the first instance.The guess for the other quantities of κ i s,n+1 , f α cr,n+1 , R e n+1 estimated explicitly, evaluated at time t n , using a forward Euler scheme for κ i s,n+1 and f α cr,n+1 and an exponential map scheme for R e n+1 .Convergence for this two-level iterative scheme is based on the changes in the norm of residuals 1-5.Iterations are carried out until changes in the norm are less than a user defined limit.In the present work the limit is taken as 10 −8 .
The Cauchy stress is then calculated as where The main appeal of an implicit scheme is the larger steps towards a solution which can be taken.As shown in the previous sections martensitic transformation is a complex problem.A possible 96 different systems can potentially be activated, each of which influencing the plastic and elastic behaviour of the crystal.Maintaining large implicit steps on such a process can be challenging for convergence.To preserve the larger global steps for the global model a step stepping algorithm has been used.This allows complex transformation to be discretised while keeping the large step size for other areas of the model (such as areas not transforming).This algorithm takes the global step time ∆t, deformation rate d n+1 , and spin rate w n+1 and divides into several sub-steps.The deformation and spin rate are discretised assuming a linear variation in ∆t.The number of sub-steps taken will depend on the convergence of the global step.Initially the step is broken into 2 sub-steps.If these sub-steps still prove too large for convergence the total number of sub-steps is increased to 4. The total sub-steps is doubled until convergence of the global step or until some user defined max number of sub-steps is reached.

B.6 Consistent elastoplastic tangent moduli
In implicit finite element procedures, the global tangent modulus is used to preserve the rate of convergence during the equilibrium iterations.The tangent moduli here can be defined as As has been stated by other authors, the deriving an analytical expression for the global tangent modulus can be challenging [39].Using the same numerical method as used for the residuals, equation (B.51) can be expressed as Here, i and j represent the components of σ n+1 and d n+1 , respectively.h is again the perturbation parameter used to disturb the solution on the jth component of unit vector e j .This method gives a convenient way to determine the global tangent modulus; the main benefit being the simplicity of the numerical implementation.
However, the main drawback is the additional computational time required to solve each timestep.Using this method, the constitutive equations must be solved an additional six times.When combined with the sub stepping procedure solve time increases drastically.
An effort has been made, therefore, to develop an approximate analytical solution for the global tangent modulus.This approximate closed form solution is developed by selective linearisation of the constitutive equations.The linearisation does not consider the rotation tensor R e n+1 , the Jacobian J n+1 , the hardness, and the transformation resistance or the change in elastic constants due to the change in stress.Hence, it is an approximate solution not the actual moduli.Starting with the same definition as equation (B.54) and expanding to make it about the Kirchhoff stress, Recognizing that the martensitic transformation is mainly connected to the deviatoric response of the crystal (i.e. the shear is far larger than the accompanying dilation).Using this we can linearise devτ n+1 and p τ,n+1 only considering the influence deviatoric stress on ξ(α) n+1 .Furthermore, the linearisation of devτ n+1 and p τ,n+1 assumes that R e n+1 is constant.As such P  The methods presented here represent approximations for the tangent moduli.The numerical method (NM) gives an all-inclusive approximate, taking into consideration all factors which influence the deformation in the model.On the other hand, the approximate analytical method (AAM) gives a selective or exclusive approximation.A comparison between the NM and AAM is given in figure B.1.This figure clearly shows both solutions appear indistinguishable upon inspection.The maximum difference in the solutions is c.0.2%.It can be concluded that either method is suitable.
Inspection of the iteration output, from ABAQUS, for both methods show the NM can reach the solution in less iterations.This is due to the more accurately calculated tangent modulus.The AAM must use smaller time increments to solve the same solution.However, speed at which the AAM can solve each iteration is far superior to the NM.This is reflected in table 3 showing the AAM is the clear choice, with no difference in result and a significant benefit in computational time.

B.7 Procedure overview
The following procedural overview (given below) highlights the numerical procedure that is used to solve the formulations constitutive equations.

== 13 )
sym C e Ω e + R e Dp R eT , W * = skew C e Ω e + R e Wp R eT (Ṙe R eT and Dp = sym Ce L P + sym Ce L tr Wp = skew Ce L P + skew Ce L tr , Ce is the elastic right Cauchy-Green tensor.The rate of the deformation and spin tensors due to plastic deformation are Dp and Wp respectively.They are both dependent on the slip activation and the amount of transformation occurring in the crystal structure. rate type of E e based on the elastic spin of lattice, Ω e .The same procedure can be followed through to get the spin tensor in the intermediate configuration 3. The final relations for D * can be found by combining equations (2.1.6),(2.1.7),(2.1.11),and (2.1.12):D * = sym C e Ω e + Ns i=1 γ(i) sym Ce Z (i) + Nt α=1 ξ (α) sym Ce K (α) , (2.1.17)W * = skew C e Ω e + Ns i=1 γ(i) skew Ce Z (i) + Nt α=1 ξ (α) skew Ce K (α) , (2.1.18)

Figure 2 .
Figure 2. Graphical representation of transformation resistance model.

Figure 3 .
Figure 3. Parametric study of critical driving force fcr.Showing stress-strain response (left) and martensite volume fraction (right).

Figure 5 .
Figure 5.The effect of rate sensitivity on the dependence of transformation rate against the applied driving force for transformation.

Figure 9 .
Figure 9. Parametric study of transformation parameter R h , showing stress-strain response (left) and martensite volume fraction (right).

a
Martensite elastic values are kept the same for both materials.

Figure 13 .
Figure 13.Visual representation of crystal orientation in relation to loading (Note: This is not the finite element mesh; grid is for improved visualisation only).

Figure 14 .
Figure 14.The true stress-strain (left) and total martensite volume fraction (right) curves for the three orientations.This figure shows significant plastic mismatch between the orientations.Also all orientations exhibit the double yield effect.

Figure 15 .
Figure 15.The true stress-strain (left) and total martensite volume fraction (right) curves in the [011] and [111] orientations for β-Cez alloy, showing a large plastic mismatch between orientations and a distinct double yield for the [011] direction only.

Figure 16 .
Figure 16.Individual transformation systems for the [001] (left) and [011] (right) directions.Clear groupings of the transformation systems can be seen in both directions.

Figure 17 .
Figure 17.Individual transformation systems for the [111] direction (left) and magnified view of shaded section (right).This direction shows no clear grouping and a large amount of cross-over among systems.

Figure 18 .
Figure 18.Volumetric expansion results: pressure vs. volume change (left) and total volume fraction vs. volume change (right).This figure shows a plateau for pressure as volume increases (left) and all transformation system are equally active (right).

Figure 19 .
Figure 19.Volumetric contraction results: pressure vs. volume change (left) and total volume fraction vs. volume change (right).This figure shows a linear increase in pressure against volume change (left) and all transformation system have zero martensite volume fraction (right).

τ = Ce ϵ e , (B. 13 )
For computational convenience, the elastic law given by equation (B.13) is written as,τ = devτ + p τ 1 (B.14)where devτ and p τ are the deviatoric and volumetric parts of τ , respectively.Given as: devτ = Ce d : devϵ e + H Ce d is the fourth order deviatoric elasticity tensor, He is the second order deviatoricvolumetric elastic coupling tensor, and M e is the elastic volumetric coefficient.Pd = Ĩ − 1 3 1 ⊗ 1 is a fourth order deviatoroc projection tensor, with Ĩ and 1 being fourth and second order unit tensors.Note that, Ce d = C ⊗ − C : C e d,0 : C ⊗ − C 0 C T , Me = M e .(B.18)

. (B. 20 )
R e ∂ ∂t R eT ϵ e R e R eT = d − D pThe previous is integrated using a backward Euler scheme, giving: ∆R e is the incremental elastic rotation tensor.The plastic strain rate Dp n+1 can he expressed as

1 : 1 ::M e n+1 1 :Figure B. 1 .
Figure B.1.Plot showing the difference in solution between the NM and the AAM, for 50 grain simulation.

Table 1 .
Lattice correspondence for possible β

Table 4 .
Euler angles and loading direction for single crystal simulations.

•
The single crystal analysis provided further experimental comparison and validation to the model.Single crystal simulations showed that there are relatively 'strong' and 'weak' directions for transformation.The strongest resistance to SIM is given by [111], requiring the largest stress to initiate.•The [011] direction is the weakest and transforms first with the largest transformation strain.

Table B1 .
Elastic constants for BCC β phase taken from literature.

Table B3 .
Comparison of CPU time between NM and AAM for various number of grains.