Mean-field elasto-plastic self-consistent (EPSC) modeling of polycrystalline metallic materials has a long history since the original introduction. The EPSC formulation was first theorized by Hill [1] and then implemented in a computer code by Turner and Tomé [2] as a small strain formulation. Much later, a large strain formulation of the model was developed in [3] and the implicit numerical aspects of the model were improved in [4]. The model had been used to interpret and predict a wide range of mechanical behaviors of polycrystalline metals including strength [5], texture [6], twinning [7], phase transformations [8], residual stress [9], elastic lattice strains [3] and damage [10] under monotonic and strain-path-change deformation conditions. To facilitate modeling of all these phenomena, a number of sub-models were implemented in EPSC including dislocation density-based hardening, twinning and detwinning [11], backstress [12], latent hardening [13], and stress and strain induced phase transformation kinetics [8,14,15].
The EPSC formulation was embedded in the implicit finite element method (FEM) to model geometrical shape changes under complex deformation boundary conditions such as Taylor impact [16], deep drawing [17,18], and extrusion [19]. Compared with the more computationally demanding but more accurate crystal plasticity finite element (CPFE) [20][21][22][23] and elasto-viscoplastic fast-Fourier transform (EVPFFT) [24][25][26][27] full-field models, the mean-field EPSC model offers a favorable balance between accuracy and computational efficiency [28,29].
At very large plastic strains, polycrystalline microstructures develop high intragranular stress -strain fields and misorientation spreads. These are phenomena not typically considered in mean-field models. Instead, the models are based on average mechanical fields per grain/inclusion. As a result, deformation textures at large plastic strains are predicted sharper than measured [30]. To relax these issues, a higher order formulation of the mean-field viscoplastic self-consistent (VPSC) framework considering spatial distributions of properties over inclusions was developed first [31], and then adapted in EPSC [32]. Although both EPSC and VPSC employ self-consistent homogenization and solve the Eshelby inclusion problem, the codes are fundamentally different. EPSC considers elasticity and does not use the commonly employed viscoplastic power-law for slip activation. The differences required re-derivation of the higher order terms from the perspective of EPSC. The present work builds upon the recent higher order formulation to create an advanced strain gradient (SG)-EPSC model capable of calculating geometrically necessary dislocations (GND) to influence hardening along with backstress and resulting size sensitive response of metallic materials. GNDs represent an extra storage of dislocations required to accommodate the lattice curvature and misorientations that arises whenever there is a non-uniform plastic deformation.
Higher order formulations of mean-field models were first established for composites [33][34][35], then adapted to models for polycrystalline metals [35][36][37]. An algorithm for the implementation of the second moments of stress in the VPSC model was presented first in [38]. With the second moments of stress available in VPSC, algorithms for obtaining intragranular fluctuations in lattice rotation rates and misorientation trends were presented in [31]. The trends were validated by comparisons with corresponding full-field calculations. With misorientation spreads available in VPSC, the subsequent works developed a grain fragmentation model, which improved predictions of texture evolution during deformation [39] and a recrystallization model, which enabled predictions of texture evolution during recrystallization [40]. These developments were unified in a comprehensive field fluctuations (FF)-VPSC model [41]. Adapting the field fluctuation formulation from FF-VPSC to FF-EPSC laid the foundation for the present development of the SG-EPSC model.
SG plasticity (SGP) theories originate from the works reported by Fleck and Hutchinson [42,43] over the work of Arsenlis and Parks [44] that introduced a new formulation of Nye's dislocation tensor to the fundamental papers of [45] and Gurtin [46] for phenomenological plasticity and [47] for strain gradient crystal plasticity. Moreover, the two-part companion papers by Gurtin and Anand [48,49] detailed the SGP theory in both small strain and finite deformation regimes. Some issues in the basic constitutive equations pertaining to the original formulation of SGP and advances for non-proportional loading scenarios have been pointed out by Fleck et al. [50]. A comprehensive review of the computational methods and comparisons for Nye's dislocation tensor and the underlying GND density is given by Das at al [51]. Arora and Acharya [52,53] developed a mesoscopic field dislocation mechanics theory considering SGP. An efficient implementation of the theory involving both energetic and dissipative higher-order effects in finite element frameworks was presented in Panteghini and Bardella [54]. While full field models like finite elements, CPFE, or EVPFFT were typically used as base models to implement SG crystal plasticity formulations [55], a novel modeling strategy is put forward in the present work to develop a mean-field SG formulation that accounts for the evolution of GNDs in EPSC for the first time. The SG-EPSC formulation merits some focus because of its computational efficiency providing incentives for scaling up to the component level via embedding into the FEM simulation tools. Moreover, the existence of several sub-models in EPSC facilitates multi-level complex simulations of multiple phenomena using a single code.
In the core of the novel SG-EPSC formulation is the calculation of intragranular orientation spreads in the grains of a polycrystalline aggregate. These spreads are based on the second moments of the stress fields governing fluctuations in the lattice rotation rates. The intragranular orientation spreads are converted into an intragranular field of elastic deformation gradients for spatial derivatives to obtain the Nye dislocation tensor and associated GNDs at the grain scale resolution. The calculated GNDs are incorporated into an existing hardening law based on statistically stored dislocations (SSDs) and an advanced composite grain (CG) model for handling twinning in EPSC. Moreover, a backstress law to influence the activation of slip systems based on GNDs is implemented into the overall SG-EPSC formulation. Importantly, the GND calculations enable the model to be verified experimentally in addition to facilitating the understanding of the underlying mechanisms governing size sensitive behavior, kinematic hardening, and strain localizations leading to failure of metallic materials. To this end, the predicted and measured grain size sensitive response of α-titanium (α-Ti) and evolution of GNDs during three strain-path-change deformations of aluminum alloy (AA) 6016-T4 are compared. To examine the GND-based backstress formulation, a cyclic deformation of AA6016-T4 is modeled and compared with corresponding experimental measurements. The successful validation of the model demonstrated the advantages and versatility of the first implementation of SG-EPSC, which is an important step in advancing the mean-field self-consistent formulations. The model formulation and these results are presented and discussed in the paper.
The layout of the paper is as follows. An overview of the EPSC formulation followed by the field fluctuation equations is presented in Section 2. The details pertaining to the novel SG formulation as well as the enhanced dislocation density-based hardening law and twinning model along with a new slip system level backstress law are described in Section 3. Sensitivity of GNDs on selected grids including spatial distributions, sizes, and grid point densities is explored in Section 4. Simulation results are presented in Section 5. Concluding remarks and a summary of the main findings are given in Section 6.
In this paper, the symbols ⋅ and ⊗ denote a dot product and a tensor product, respectively. Tensors are bold letters. Tensor components and scalars are italic and not bold. Moreover, the slip families / modes are enumerated by the superscript α; the slip systems are denoted by the superscript, s, while the twin families / modes are enumerated by the superscript β; the twin systems are denoted by the superscript, t.
The EPSC model originated from the initially small strain formulation proposed by Hill [1], which was implemented by Turner and Tomé [2], and later extended to large strains by Neil et al. [3]. The particular version advanced in the present work is from [4]. In EPSC, a polycrystalline aggregate is represented as a collection of single crystal inclusions with ellipsoidal shapes and having a crystal orientation and weight. The inclusions are embedded in a homogeneous effective medium (HEM) during self-consistent homogenization. It is important to note that all EPSC formulations are incremental in nature. The stress rate of a single crystal grain can be given by:
where σc is the Cauchy stress rate for a grain c, C c is the fourth rank elastic stiffness tensor, D c is the total rate of deformation tensor, representing the strain rate tensor in this work, and D pl,c is the plastic strain rate stretching tensor. The plastic stretching tensor is:
where γs,c is the shear strain rate for a slip system, s, in a grain, c, m s,c is the symmetric Schmid tensor defined by the slip plane normal unit vector, n s,c , and a dimensionless unit vector indicating the slip direction aligned with the Burgers vector direction, b s,c . The slip systems also include pseudo-slip modes such as twinning and phase transformations. The last term in Eq. ( 1) is the rate of elastic dilatation that is typically small, especially at small strain levels.
Eq. ( 1) can be conveniently expressed as a linear relationship between single crystal stress rate and strain rate as:
and the same relationship at the aggregate level is:
where σ is the polycrystal stress rate, D is the polycrystal strain rate, and L is the polycrystal instantaneous elastoplastic stiffness to be calculated in an iterative process described shortly. The single crystal instantaneous elastoplastic stiffness tensor, L c , is derived to be:
where I is a second-rank identity tensor and X ss ʹ is a square matrix of size corresponding to the number of active systems as:
with h ss ʹ ,c = ∂τ s c ∂γ s ʹ ,c derived as the hardening rate of active system. Note that the active systems can include pseudo-slip systems such as twinning. h ss ʹ bs is the contribution to the hardening matrix from the similar derivatives of backstress. The derivatives will be given in Section 3.2.
To account for the rotation of single crystals during plastic deformation with respect to the polycrystal specimen frame, the lattice rotation rate tensor or spin of a single crystal is:
where the applied macroscopic spin, W c,app , is a sum of an overall macroscopic spin rate applied over the polycrystal W and an additional spin, Π c , that arises from solving the Eshelby inclusion problem for embedded inhomogeneity inside the HEM, subject to applied boundary conditions [56]:
where P c is the anti-symmetric Eshelby tensor. Lastly, W pl,c is the plastic spin calculated from shear strain rates and the anti-symmetric Schmid tensors, q s,c :
The large strain deformation using the Jaumann rate of Cauchy stress is now invoked to account for rotations of the crystal lattice with respect to the global frame. Eqs. ( 1) and (3) are then reformulated as:
Similarly, Eq. ( 4) for the polycrystal is also corrected:
The Jaumann and Cauchy stress rates on the single crystal and polycrystal levels are related in a common frame, considering the spins via:
In this paper, we use the bar above a quantity to indicate a volume average obtained by summing the products of the per grain quantities and the relative weights of the grains, w c .
The equilibriums between the single crystal and polycrystal levels of stress and strain rates are solved via the interaction equation:
where L c * is the interaction tensor:
that defines a localization tensor:
which is used to calculate the polycrystal instantaneous elastoplastic stiffness tensor from the single crystal quantities as:
Since the right-hand side of Eq. ( 16) is also a function of L via Eq. ( 15), these equations must be solved iteratively until convergence with a chosen tolerance of 0.001.
In EPSC, a slip system in a grain becomes active when the resolved shear stress reaches the slip resistance i.e. critical resolved shear stress (CRSS), satisfying the following conditions:
where τ s,c c is the CRSS. τ s bs influences the driving force for activation. Backstress was previously modeled using phenomenological models in EPSC [13,57], but in the current work, it is derived from the GNDs, as will be described shortly. The first condition indicates that a slip system is active when its resolved stress state is on the yield surface. The second condition then ensures that the stress state remains on the yield surface during plastic deformation. The CRSS evolves with plastic deformation via:
Lastly, the shear strain rate is calculated as:
The model summarized in Section 2.1 simulates large strain mechanical response of polycrystalline aggregates based on average properties of constituent inclusion. Recent works have shown that the effects of intra-inclusion field fluctuations can be approximated via calculating the second moment quantities such that the spatial variation of the properties are considered. The present work builds upon these extensions, specifically upon the field fluctuations (FF) formulation initially implemented in VPSC [40] and then adopted to EPSC [32]. For completeness, we summarize the FF formulation of EPSC in this section to provide the basis for the strain gradient formulation. The detailed expressions of the terms used in the FF formulation can be found in Appendix A.
The second moment formulation summarized below is the culmination of several previous works [31,[38][39][40][41][58][59][60]. Two sources of stress fluctuations: intragranular variation of mean grain properties within the polycrystal [38] and intragranular variation of orientation [39] are considered. The second moment of stress rate fluctuations is:
The angle brackets 〈〉 indicate a second moment quantity. M and M c are the polycrystal and single crystal compliances calculated by taking the inverse of the respective elastoplastic stiffness tensors from Eqs. ( 5) and ( 16) [61], and 〈 σ ⊗ σ〉 c is the total fluctuation of stress rate:
In the above equation, σt,c is the mean stress rate at time t for a grain,
the second moment of stress rate induced by the average grain stress rate fields, and 〈 σt,δr ⊗ σt,δr 〉 c is the second moment of stress rate fluctuations induced by the fields of misorientations, δr. 〈 σt,q ⊗ σt,δr 〉 c and 〈 σt,δr ⊗ σt,q 〉 c are transpose of each other and describe the cross-covariance between the two fluctuation sources.
Having the second moments of stress rate, the second moments of lattice rotation rates can be calculated at the same time step, t:
ij , is the dual vector of the spin tensor from Eq. (7). 〈w t, σ ⊗ w t, σ 〉 c and 〈w t,δr ⊗ w t,δr 〉 c are the second moments of lattice rotation rate terms caused by the intragranular stress rate and misorientation fluctuations, respectively. 〈w t, σ ⊗ w t,δr 〉 c and 〈w t,δr ⊗ w t, σ 〉 c are the cross-covariance terms that are transpose of each other. Finally, the second moments of misorientation from the time step, t, to the current time step, τ, are:
where δr t inc ≈ Δt 2 w t, σ ,c is the vector part of the incremental rotation, while δr t,rot = R t,c inc δr t is the vector part of the rotated misorientation of the misorientation vector, δr t , using the increment in rotation, R t,c inc , written in the rotation matrix representation. Note that in these calculations, the second rank orientation and misorientation tensors are expressed as quaternions. Then the quaternions are reduced to their vector parts. The derivations of each terms in Eqs. ( 21)- (23) are summarized in Appendix A. The second moments of misorientation also enable a formulation of a grain fragmentation model, which is provided in Appendix B.
Relying on the second moments of misorientation, a strain gradient (SG)-EPSC is formulated to calculate the evolution of geometrically necessary dislocations (GNDs). The Nye tensor based on the wellestablished theory is [62]:
where F c,t e is the elastic deformation gradient of a single crystal, c, at a simulation time, t. F c,t e is made up of the rotation, R c,t e , and elastic stretching, U c,t e parts. α c,t is the grain's Nye tensor that has been treated as an equivalent to GND density via [42,46,47]:
where b c is the magnitude of the Burgers vector. It is apparent that the Nye tensor is not possible to calculate in the present state of EPSC since the misorientation spreads are not associated with spatial points, so the curl operation cannot be applied. Therefore, a procedure to generate spatial distributions is developed.
The second moments of misorientation defined in Eq. ( 23), 〈δr ⊗ δr〉 c , describe the intragranular orientation spreads in grains from their mean grain orientation. A set of discrete δr rotation vectors can be sampled from a given misorientation spread. The steps involved in the sampling were described in Zecevic et al. [40]. The same procedure is adopted here. The Cholesky decomposition of the second moment is:
where L is the decomposed lower triangular matrix that transforms random vectors, v, with a mean of zero and unit variance to random vectors with a mean of zero and a variance following the second moment to give the vectors δr:
The random vectors, v, are generated using the gasdev subroutine from Press et al. [63]. Given the discrete δr vectors belonging to the given spread, the corresponding increment in rotation tensors, ΔR c e , can be readily calculated using:
where the tensor components are the components of δr. The number of sampled δr vectors can be any desired number.
To obtain the elastic deformation gradient at the sampled vector point needed for the calculation of the Nye tensor, the elastic stretch tensor, ΔU c e , needs to be calculated. The elastic stretch incremental tensor is obtained from the elastic strain increment tensor, ΔE c e , following the logarithmic strain definition:
The elastic strain tensor increments, and the total elastic strain tensor are both available in the code.
Although the transformed vectors follow the variance of the distribution, they are not associated with any spatial information to allow for the calculation of the spatial derivatives. Therefore, in the next step, we distribute the sampled rotations and associated elastic stretches into spatial points on a 3D grid to obtain a function of increments in deformation gradients, ΔF c e (x), where x is the position vector of the spatial point. Examples of 3D grids of 5 × 5 × 5 spatial points are shown in Fig. 1. In these examples, 125 discrete ΔF c e (x) are needed to populate the grid points.
As shown in Fig. 1, the physical size and shape of the grid per ellipsoid follows the length scale of the ellipsoidal grain in which the grid is embedded. The grid arrangement updates during plastic deformation to reflect the evolving shape of the ellipsoidal grain with plastic strain. As a result, the GND density calculations in EPSC are sensitive to the shape and size effects. The embedded grid is either a cube in equiaxed spherical grains or a rectangular parallelopiped in ellipsoidal grains. The effects of distributing ΔF c e (x) over the intragranular grids on GNDs is discussed in Section 4.
With a spatial discrete distribution of the elastic deformation gradients per grain, a differentiable functional can be defined. To this end, the tricubic interpolation procedure is adopted. Each component in the tensor ΔF c e (x) is interpolated separately to obtain a continuous tensor field for calculating the curl necessary to get the Nye tensor. The tricubic interpolation stems from trilinear interpolation, which is a multivariate procedure utilizing a 3D grid. It generates a cubic functional form of the values on the grid and approximates the values at intermediate positions between the grid points. The cubic function of ΔF c e,11 is:
Expanding all the terms, the following function is obtained:
involving a total of 64 terms with unique coefficients. A system of equations to represent Eq. ( 31) can be constructed:
where B is a column vector containing values at the grid points. For the grids shown in Fig. 1, there are 125 grid points, so the vector is of 125 × 1 size. X is a 64 × 1 column vector containing the 64 coefficients:
Finally, the matrix A has a dimension 125 × 64 and contains the position information:
The coefficients vector X in Eq. ( 32) is solved using QR decomposition to arrive at the functional form of Eq. ( 31), and the spatial derivatives can be readily evaluated:
In this work, the curl of a tensor field is defined as Das, Hofmann and Tarleton [51]:
Combining Eqs. ( 24), ( 35) and ( 36), the Nye tensor increment is:
It should be noted that the Nye tensor increment is not a physical quantity by itself in terms of density of GNDs because the concept of density of GNDs, in continuum theories, is strictly related to the concept of incompatibility requiring a total quantity. The increments are summed to obtain the physical quantity of GNDs. The GND densities at each grid point can then be calculated using Eq. ( 25). Then, the average GND density increment of the grain is the average over the grid points:
The quantity calculated using Eq. ( 38) is the evolved GND density in the current simulation increment. The hardening law in our model requires dislocation densities to be defined per slip system, as will be defined in the next section. To this end, we assume that the GNDs are developed on active slip systems due to shearing. Therefore, we distribute the GND density increments in each grain over the active slip systems using:
The relative slip activity term on the right is based on the shear strain rates of the active systems. The calculated GND density per active slip system can then be added to the statistically stored dislocation (SSD) density per system for the evolution of hardening, as will be described shortly. The GND density on a slip system in a grain begins to evolve from a specified initial value per slip system, ρ s,c,t=0 GND = ρ s,c GND,0 , as:
where the current GND density at τ = t + Δt, ρ s,c,τ GND , is the sum of the GND density accumulated up to t, ρ s,c,t GND , and the current GND density increment, Δρ s,c GND . With the GND density determined per slip system, a GNDbased backstress formulation from Kapoor, Yoo, Book, Kacher and Sangid [64] is adopted to calculate the backstress on each slip system:
where k bs is a scaling parameter and μ s is the shear modulus. The contribution of backstress to the hardening matrix in Eq. ( 6) is given by:
where the derivative is:
Eq. (43a) defines the diagonal terms, while Eq. (43b) defines the latent hardening off-diagonal terms contributing to the hardening matrix.
It is important to note the backstress is passively accumulated in the opposite directions of the active slip systems during forward loading. Upon load change or load reversal (reverse loading), the accumulated backstress lowers the yield requirements (the Bauschinger effect).
This section describes a dislocation density-based hardening law and a twinning model advanced for the first time to include the GND density effects. The resistance to slip in Eq. (17a) consists of four terms [65][66][67]:
where τ α 0 is an initial slip resistance per deformation mode as a fitting constant, which can be sensitive to temperature and strain rate [16]. τ s HP is a Hall-Petch (HP)-like barrier term due to initial grain size or due to the formation of twin lamellae [68]:
In the above Eq. ( 45) for twin free grains, H a is the Hall-Petch fitting parameter and d g is the grain size. When twinning is active and twin lamellae nucleate, Eq. ( 45) is replaced with Eq. (46). Then the grain size is replaced by the mean-free path, defined as d s mfp =
(1-f PTS )dc sin (λ) , with the angle, λ, between the slip plane and that of a predominant twin system (PTS) in the grain [69]. The constant, d c = dg n lamellas , is the distance between twin lamellae calculated using the grain size, and a chosen number of lamellae developing per grain, n lamellas = 4. The barrier factor, B = f PTS -f PTS,0 f PTS,max -f PTS,0 includes two constants, f PTS,max = 1.0 and f PTS,0 = 0.01, corresponding to the maximum and minimum fractions that twins could occupy in a grain, while f PTS is the current fraction of the PTS in that grain.
The dislocation density-based hardening law available in EPSC was based solely on the statistically stored dislocation densities (ρ s for -forest dislocation density and ρ deb -debris dislocation substructure dislocation density), which evolve with the slip rates. Such law corresponds to the dissipative constitutive law [47]. On the other hand, the driving force for slip influenced by backstress represents the energetic constitutive law due to GNDs. As in a variety of previous analyses based on Gurtin (e. g. [70]), the effect of GNDs was only incorporated in the energetic term.
The presumption was that the main contribution to the dissipative hardening comes from the statistical dislocations. In principle, the effect of GNDs can be incorporated into the dissipative hardening, as shown by [71]. Inclusion of such terms involves no additional computational complexity. It is emphasized that there is no fitting parameters involved in including the effects of GNDs. Such advanced formulation should be appropriate to capture phenomena associated with strain-path changes such as non-linear unloading, change in yield strength (e.g. Bauschinger effect if the loading is reversed), and changes in hardening or softening.
The dissipative hardening terms have several physically based parameters, as described next.
τ s DD is the strain hardening term accounting for the effects of forest dislocation density and GND density on hardening. Previously, the term only considered statistically stored forest dislocations [72] but now the evolved GNDs are also included:
where χ is an interaction parameter set at 0.9 for all simulations and L ss ʹ is a latent hardening interaction matrix. For FCC, it is defined with parameters a 0 to a 5 [73,74], and for HCP, L ss ʹ is an identity. ρ s ʹ for is the forest dislocation density, and ρ s ʹ GND,tot is the total GND density.
τ debris accounts for the strain hardening with evolving debris dislocation density, ρ deb :
The total dislocation density is:
where ρ s forw is the forward forest dislocation density on s, while ρ s + rev and ρ s - rev are the reversible dislocations on the s + and s -directions. The consideration of reversible dislocation populations makes the hardening law strain path sensitive. The evolution of the forest densities with shearing strains is [7,[74][75][76][77]:
The initial conditions are:
while
is specified as an input to the model. In the above equations, k a 1 is one of the fitting parameters controlling the rate of generation of dislocations, k a 2 is a temperature (T) and strain rate sensitive parameter driving dynamic recovery, p is a reversibility parameter set to 0.2, g ss ʹ is the dislocation interaction matrix set to g ss = 1 and g ss ʹ = 1 [74,78,79], m is a constant governing the rate of dislocation recombination set to 0.5 [80], and ρ s 0 is the dislocation density at the point of local load reversal [76].
The coefficient k a 2 is:
where, k B is the Boltzmann constant, ε0 = 10 7 is a reference strainrate, g α is a fitting parameter representing the effective activation enthalpy, and D α is the drag stress. Lastly, the debris / substructure dislocation density evolves with:
where q α = 4.0 is a fitting parameter determining the number of dislocations that become debris / substructure. The debris dislocation density evolves from a very small initial value of 0.1 m -2
The twinning volume fraction increment is [81]:
where Δγ t is the shear strain increment on the twin system, t, and S β is the intrinsic twinning shear. The twinning volume fraction of a twin slip reaching the critical value of 0.01 satisfies the nucleation criteria. Nucleated twins inherit state variables of their parents. Multiple systems can nucleate twins, but only the first variant reaching f pts = 0.03 is selected as the PTS. The twinning resistance evolves as:
Fig. 2. Magnitudes of the 125 discrete δr vectors sampled from the 〈δr ⊗ δr〉 c second moment of a randomly selected grain/inclusion, c, of AA6022-T4 organized into the 5 × 5 × 5 regular grids, depicted in Fig. 1a. The grids of corresponding misorientation angles from the mean orientation, δθ, are also provided in degrees. The same datapoints are spatially organized in three different distributions designated as: linear, concave, and convex. The second moment for the grain/inclusion of AA6016-T4 is taken after biaxial tension in RD and TD to an equivalent strain of 0.2, followed by unloading, and then final simple tension along RD to a total equivalent strain of 0.3.
In the above equation, τ β 0 is an initial friction term, τ t 0,HP is a barrier term, and τ b slip is a term that couples the twin resistance with slip. τ b 0 is calculated as:
where τ β crit and τ β prop are the nucleation and propagation stresses, set to be equal in the current work so Eq. ( 59) reduces to τ β 0 = τ β prop . The initial barrier term is:
After nucleation when twin volume on the PTS reaches the critical value, the barrier effect evolves to:
The twin-slip interaction term is:
where C αβ is the twinning hardening matrix and ρ s tot is the total dislocation density from Eq. (49).
The derivative of Eq. ( 44) gives the hardening matrix used in Eqs. ( 6) and ( 18):
Where the components are:
The Hall-Petch derivative reflecting the grain size is zero, ∂τ s HP ∂γ s ʹ = 0. When twins are present in the grain, the derivatives of Eqs. ( 46) and ( 61) with respect to shear strain contribute to the hardening matrix as: ∂τ t HP ∂γ t,pts = ∂τ t HP ∂f pts ∂f pts ∂γ t,pts
Each of the above derivative can represent the slip-twin and twin-twin interactions, respectively, with the partial derivative ∂f pts ∂γ t,pts calculated as:
Unlike full-field models in which spatial relationships between material points exist, mean-field models do not associate spatial coordinates to material points. A procedure for embedding spatial grids of material points into single crystals in EPSC is described in Section 3.1. This section evaluates the sensitivity of calculated GNDs on grids representing the elastic deformation gradient fields over inclusions. Spatial arrangement, sizes, and grid point density are all explored. The results are for the initial texture of AA6016-T4 subjected to biaxial tension in the rolling direction (RD) and transverse direction (TD) to an equivalent strain of 0.2, followed by unloading, then reloading in uniaxial tension in RD to 0.3 strain. The SG-EPSC simulated intragranular fields are presented first and then the effects of these factors on the final predicted GND evolution are compared at the end of this section.
To assign each of the 125 sampled δr vectors to the 5 × 5 × 5 grid points, the vector magnitudes are quantified. Additionally, corresponding misorientation angle, δθ, between each orientation from the spread and the mean grain orientation of the given grain is calculated. Note that δθ is not a misorientation increment. The misorientation angles are obtained from the 125 rotation tensors belonging to the given spread, ΔR c e (x), relative to the mean orientation of the given grain. Three spatial distribution/arrangement methods, named "linear", "concave", and "convex", to facilitate the spatial assignment of the fields are implemented in the SG-EPSC code. Fig. 2 shows the plots of the three assignments for a randomly selected grain. Every grain of the polycrystalline aggregate being modeled follows the same spatial distribution method throughout the entire simulation. The 5 × 5 × 5 grids per grain get populated with ΔF c e (x) every strain increment to calculate GNDs using one of the selected spatial distribution methods.
Orientation spreads vary from grain to grain in a polycrystalline aggregate. Moreover, the orientations spreads evolve differently/heterogeneously in grains with plastic strain. A given orientation spread per grain belonging to a polycrystalline aggregate is represented using 125 orientations, which are sorted by magnitude. Out of many ways to associate the 125 orientations to 125 grid points, we follow one of the three assignments. The linear distribution exhibits directionality in the intragranular spread. In the example shown in Fig. 2, the rotational magnitude is lowest at the corner point with the coordinates of (-12.5, 12.5, 12.5), and highest at the corner point with the coordinates of (12.5, -12.5, -12.5). The diagonal connecting these two corners is the axis corresponding to the axis of largest plastic deformation in this grain. The fields are plotted on the original undeformed grid for visualization, while the actual simulated grid will be distorted subject to the size and shape changes. The remaining 123 points of the 3D grid receive assignment to the grid points layer-by-layer based on the distances from the corner point (12.5, -12.5, -12.5) and in the descending order of the sampled vectors' magnitudes.
To generate the concave and convex distributions on the 5 × 5 × 5 grid, the sorted grid points are separated into three layers. The first layer contains the 8 corner points of (-12.5, -12.5, -12.5), (-12.5, -12.5, 12.5), (-12.5, 12.5, -12.5), (-12.5, 12.5, 12.5), (12.5, -12.5, -12.5), (12.5, -12.5, 12.5), (12.5, 12.5, -12.5), (12.5, 12.5, 12.5). The middle layer contains 56 points that are adjacent to the corner points. The final layer contains the remaining 61 points placed onto the three XY, XZ, and YZ planes that go through the center point (0,0,0). For the concave distribution, the first layers' 8 corners points receive the 8 largest rotation angles. The next 56 largest rotations are assigned to the second layer. Finally, the remaining 61 rotations go to the third layer. The points in the layers are randomly shuffled after assignment to avoid sharp local gradients due to numerical loops that iterate through the Fig. 5. Contour plots of ΔF c e,11 (x) after tricubic interpolation over the grids shown in Fig. 2. The data to form the 5 × 5 × 5 grids with 25 µm edge for the interpolations was taken for a randomly selected grain of the AA6016-T4 polycrystal after biaxial tension to a strain of 0.025. points during assignment, but some local gradients in the magnitude are still present. For the convex distribution, the steps involved in creating the concave distribution are reversed. The largest 61 rotation vector magnitudes are assigned to the third layer, while the smallest 8 rotation vector magnitudes are assigned to the corner points in the first layer.
Next, the effect of grid size on calculated GNDs is evaluated. In general, GNDs are defined based upon the length scale of interest meaning that the GND density is expected to vary with spacing between the grid points [82]. The initial texture provided to EPSC as an input from a file contains a set of Euler angles and a relative weight for each grain. The code calculates the actual sizes of every grain based on the relative weights and a specified average grain size of the material being simulated, resulting in every grain having unique sizes and weights. Three SG-EPSC simulations are performed in which only the average grain size input parameter is changed from 15 μm to 25 μm to 50 μm, while the rest of the setup remains the same in every simulation. In Fig. 3, the visualized grains from every simulation case are those with their size closest to the average grain size. The range of field values for the three grain sizes are similar as they are from the same deformed state.
Finally, the role of grid point density in predicting GNDs is evaluated. We maintain an odd number of grid points per dimension, so a physical center point exists as needed by the implemented distribution methods. The number of grid points evaluated was limited by the memory usage of a Fortran application compiled on Windows and the available memory of the computer. A maximum of 7 × 7 × 7 was possible to evaluate so only two grid densities are simulated since a grid smaller than 5 × 5 × 5, i.e. 3 × 3 × 3, would distribute only 27 sampled statistics, which cannot sufficiently describe the intragranular spreads. Fig. 4 compares the distributed fields of the same 25 μm grain as in Fig. 3b with the same fields represented using a 7 × 7 × 7 grid. The 343 vectors are samples from the same distribution as the 125 vectors. The 7 × 7 × 7 grid relies on more sampled points from the same distribution and shows slightly smoother field, but the overall magnitude and range of the values are about the same.
After the grids are established based on distributing rotation vectors to the grid points, the corresponding elastic deformation gradients, ΔF c e (x), are calculated and assigned. The different distributions are expected to influence GNDs. Before showing the sensitivity of calculated GNDs on the selected grids, the tricubic interpolation method is illustrated by showing three slices of ΔF c e,11 (x) fields. Fig. 5 shows the XY, XZ, and YZ slices of the 3D spatially interpolated fields from the three types of distributions from Fig. 2. The coordinate position of the third dimension for each slice is labeled above each figure. The figure shows that the tricubic interpolation procedure successfully creates the continuous fields from the discrete ΔF c e (x) values at the grid points for the calculation of spatial derivatives. As a side note, trilinear, triquadratic, and tricubic interpolation methods were explored to represent the fields. While the trilinear interpolation method produced high degree of errors, especially at higher strain levels, the triquadratic interpolation method achieved similar interpolation accuracy to the tricubic interpolation method at a minimal computational advantage. Nevertheless, the tricubic interpolation method is used for the case studies presented in this paper owing to the highest accuracy.
Fig. 6 shows the sought sensitivities on the predicted GND evolutions. The three different spatial distributions and grid point densities predicted similar GND trends. However, the effect of grain sizes is significantly more pronounced. Evidently, the smaller grain size results in higher GND density. The grid point density effects are minor meaning that the 125 sample points of the 5 × 5 × 5 grid are sufficient for calculating the derivatives. Interestingly the jumps in GND density at an equivalent strain of ~0.2 are predicted owing to the strain path change from biaxial tension to simple tension. The predicted jump occurs because of the sudden change in active slip systems accommodating plastic strains and reorientation trends at the strain path change. Therefore, the SG-EPSC model is size and strain path sensitive. More details on the changing strain path cases of AA6016 are given in Section 5.2.
While taking an average of multiple spatial arrangements may lead Fig. 6. Sensitivity of the predicted evolution of GNDs with equivalent strain on selected grids for a polycrystal of AA6016-T4: (a) grids of different spatial distribution, Fig. 2, (b) grids of different size, Fig. 3, and (c) grids of different point density, Fig. 4. The straining was biaxial tension in RD and TD to an equivalent strain of 0.2, followed by unloading, and simple tension along RD to a total equivalent strain of 0.3.
Model parameters for α -Ti adopted from [83]. .0) H α 0.26 (0.15) 0.28 (0.1) 0.36 (0.2) H αβ=T 0.24 (0.14) 0.26 (0.16) 0.28 (0.16) H αβ=C 0.23 (0.11) 0.33 (0.21) 0.26 (0.2) k α 1 [10 7 m -1 ] 1.5 (4.5) 7.2 (8.7) 5.7 (7.5) g α 0.041 (0.03) 0.088 (0.04) 0.046 (0.03) D α [MPa] 600 (300) 800 (400) 800 (400) ρ α forw,0 [10 12 MPa] 1.0 (0.5) 1.0 (0.5) 1.0 (0.5) b a [10 -10 m] 2.95 5.54 2.95 b β [10 -11 m] S β τ β 0 [MPa] H β [MPa ̅̅̅̅ m √ ] C αβ [10 3 ] α = P α = π α = B Tensile Twin β = T 3.0174 0.175 49.2 0.4 7.11 6.05 4.08 Compression Twin β = C 2.725 0.218 77.7 0.4 6.13 6.03 5.25
to more accurate predictions of GND density, the upper bound concave spatial arrangement is selected for the results presented in the next section of the paper for the case studies. The default setting in the code is therefore a concave 5 × 5 × 5 grid.
The high purity α-Ti dataset used in the present work is taken from Savage et al. [83]. The dataset consists of flow curves, texture evolution, and twin volume fractions measured during compression along TD for two grain sizes: 25.5 μm and 150.4 μm. Specimens to measure the data were taken from a straight-rolled plate, which was annealed to produce the two grain sizes. The initial and deformed textures were measured by electron backscatter diffraction (EBSD) and the twin volume fractions were determined via the analyses of the scans [84].
The EPSC modeling parameters for the material are shown in Table 1. The parameters are also taken from Savage et al. [83]. Like in Savage et al. [83], three families of slip modes of prismatic, first order pyramidal 〈c + a〉, and basal are considered in the present simulations.
Their notation is as follows: α = P, {1100}〈1120〉, α = π, {1011}〈2113〉, and α = B, {0001}〈2110〉, respectively. The two twinning modes considered are the tensile and compression twins: β = T, {1012}〈1011〉 and β = C, {1122}〈1123〉, respectively. Note that since the α-Ti data did not involve changing deformation paths, the backstress law was not active in the simulations. The predicted and measured α-Ti stress-strain curves using the parameters are shown in Fig. 7a. As expected, the fits using the same parameters are reasonably good as the dislocation density-based hardening law is nearly identical to that in Savage et al. [83]. In the first simulation setup (Fig. 7a), the novel SGP formulation in the model was not needed and the grain size effects are captured solely by the Hall-Petch effect. Thus, the first simulation setup is labeled as "HP only." This acts as the reference setup to demonstrate the benefits of the novel SG-EPSC model.
Two additional simulation setups relying on the novel SG-EPSC are introduced to predict the same grain size sensitive mechanical responses of α-Ti. The second setup involved suppressing the Hall-Petch (HP) effect and relying on GNDs to govern the size sensitivity. The Hall-Petch coefficients were set to zero. Therefore, the second simulation setup is labeled as "GND only". The predicted stress-strain curves are shown in Fig. 7b Importantly the grain size effect is predicted reasonably well. The SG-EPSC with GND model alone predicts the stronger mechanical response for the smaller grain size, but not enough to capture the measured response. However, the hardening parameters can be adjusted to improve the fit.
The third simulation setup involved both HP and GND effects and some adjustment of the hardening parameters was needed to achieve good fits of the data. Otherwise, the model would over predict the stressstrain responses. The adjusted parameters are listed in Table 1 as the values in parentheses. Half of the initial SSD density is assumed to be the initial GND density, as shown in Table 1. The size effects arise from a combination of GNDs and HP with the former governing more of the effects than the latter. As a result, more of the size sensitive mechanical responses is driven from the physics-based formulation. Only twins contribute to the HP size effect. Fig. 7c shows the predictions demonstrating the best fits achieved using the third setup, especially at larger strains. The decrease of strain hardening rates at very larger compressive strains is owing to the loss of uniaxial stress state in the specimens because of the occurrence of barreling due to friction and onset of localizations unavoidable during compression testing. Since these damage aspects are not modeled, the model is increasingly deviating from the measured data at very large compressive strains.
Comparison of measured and predicted texture evolution for the two α-Ti materials of different grain size is presented in Fig. 8. Pole figures of the initial and deformed textures are shown. The initial materials have typical rolled texture for α-Ti exhibiting the split of the c-axes for most of the grains in the TD. The 150 μm grain size material had a sharper initial texture than the 25 μm grain size material. The initial textures were represented using 325 weighted orientations for the 25 μm grain size material and 272 weighted orientations for the 150 μm grain size material. The measured data had about 500,000 indexable points by EBSD.
The measured textures were compacted for the simulations using the procedure described in [85,86]. The deformed textures were measured at 0.2 strain for both materials [83]. The predicted textures are compared for the three simulation setups. First, the predicted texture evolution agrees reasonably well with the measured data for all three setups for both grain sizes. Second, poles showing the predicted texture evolution using SG-EPSC show a bit weaker intensity texture with slightly larger spreads of orientations owing to the second moments and any underlying grain fragmentation taking place. The strain of 0.2 is not sufficient to cause a great deal of grain fragmentation taking place. Similar prediction were reported in Zecevic et al. [39] where weaker textures were attributed to the second moments and underlying fragmentation of grains.
Given the good predictions of the flow behavior and texture evolution indicate that the predicted relative activities of deformation modes accommodating the plastic strains under compression should be inferred correctly because the mechanical responses, texture evolution, and twinning evolution are directly the consequences of the activities. The relative activities are compared in Fig. 9. Minor differences in the activities are predicted using different setups. Evidently, prismatic slip predominates during deformation. As the prismatic slip hardens, there is an increase in the activity of basal slip. Pyramidal slip is the least active. The coarser grain size material shows more compressive twin activity than the smaller grain size material. Larger grain size is known to promote more twinning activity. Interestingly, a little less tensile twinning operates in the coarser grain size material than in the smaller grain size material.
Lastly, the predicted SSD and GND evolutions are shown for the three setups in Fig. 10. The smaller grained material is predicted to develop more SSDs (forest dislocations) and GNDs to rationalize the grain size sensitive mechanical response shown in Fig. 7.
A dataset for modeling of GNDs and flow behavior of AA6016-T4 is taken from [89]. Specimens of the alloy were subjecting to three strain paths: biaxial tension to uniaxial tension (BTUT), plane-strain tension to uniaxial tension (PSTUT), and uniaxial tension to uniaxial tension (UTUT). The equivalent strain levels upon the pre-strains were approximately 0.2. The specimens were unloaded and then reloaded in uniaxial tension along the RD. Flow response during reloading in uniaxial tension was measured. GND density was measured via high resolution EBSD (HREBSD) during pre-straining and during reloading. The measured flow behavior and the evolution of GND density during these strain path changes are modeled using the novel SG-EPSC.
The stress-strain experimental data of the alloy was previously modeled by [90,91] using an older version of EPSC. The novel SG-EPSC model has the same strain-path sensitive dislocation density-based hardening law as the older version of the model. Therefore, the model parameters established the previous works are used as a good initial guess [92]. After including the contribution of GNDs into the hardening law, some tunings of the parameters were necessary. The SG-EPSC model parameters for AA6016-T4 are provided in Table 2. The previous work also reported the measured initial GND density that is used in the present work to properly initialize the SG-EPSC model. It should be noted that the Hall-Petch barrier effects were not considered in the AA6016-T4 simulation cases presented in Sections 5.2 and 5.3 because the alloy was not characterized for multiple grain sizes.
The initial texture of the alloy is shown in Fig. 11a. The texture was Table 2 Model parameters established for the evolution of slip resistance, latent hardening, and GND-based backstress for the octahedral slip, {111}〈110〉, of AA6016-T4. The invariant single crystal elastic constants taken for the alloy are: C 11 = 108.2 GPa, C 12 = 61.3 GPa, and C 44 = 28.5 GPa [95]. τ 0 [MPa] k1 [10 8 m - 1 ] g D [MPa] b [10 -10 m] ρ forw,0 [10 12 m -2 ] ρ GND,0 [10
12 m -2 ] 30 1.0 0.07 400 2.86 4.9 5.2083 a0 a1 a2 a3 a4 a5 k bs 0.068 0.068 0.0454 0.625 0.137 0.122 0.02
measured by neutron diffraction and compacted to 100 weighted orientation using the procedure developed in [85,93,94] for simulations. The predicted mechanical responses of the alloy during the three strain paths are presented in Fig. 11. The predicted and measured flow curves can only be compared during the reloading in UT. It should be noted that backstress played a secondary role during these strain path changes because of no load reversals present during straining upon unloading. The backstress parameter listed in Table 2 was adjusted using the cyclic curves presented shortly in Section 5.3. Fig. 12 compares the measured and predicted evolution of GND density for the three strain paths. The predicted GNDs are of comparable magnitude to the measured data. Both BTUT and PSTUT show jumps in the evolution upon the strain path change. The observed phenomenon is predicted reasonably well using the SG-EPSC model. The jump arises within the relatively short elastic range during reloading upon unloading. Interestingly, the predicted SSD evolution also shows a more continuous jump at the strain path change. The jumps are owing to the activation of different slip systems and different texture reorientation trends after the strain path changes relative to those at the end of the pre-strain paths.
Relying on the hardening parameters presented in the previous section, the physically based backstress formulation driven by GNDs along with the strain-path sensitive hardening law is used to simulate several load reversal deformation conditions measured for AA6016-T4. The experimental data is taken from [90]. Fig. 13 shows the model predictions with and without backstress. Evidently, the predictions are slightly more accurate with backstress. It should be noted that the hardening law is strain path sensitive because it includes the reverse glide of dislocations. Therefore, the predictions are relatively good even without backstress. The original hardening law formulation as well as the contribution of the reversible dislocations mechanism to predicting load reversals was discussed in [57]. Fig. 13 only isolates the contribution of backstress to the predictions as the novel aspect of the present work. When the backstress contribution to the slip system activation is included to obstruct dislocation motion, the level of flow stress to cause the activation increases, as evident from the comparison in Fig. 13. Upon a load reversal, the backstress action reduces the needed driving force to yield governing the so called Bauschinger effect. Upon yielding in the reverse direction, the backstress builds up in the opposite direction causing the increase in the flow stress again. Therefore, the cyclic curves reach lower stress levels as well as displays a slightly sharper transition at the re-yielding points upon load reversals when the backstress formulation is disabled. These observations are consistent with those observed in [57] employing a phenomenological backstress model in EPSC.
Additionally, Fig. 14 shows the evolution of backstress versus accumulated shear strain on the mist active slip system in a randomly selected grain. The backstress values between the directions are modeled to exhibit symmetry, consistent with the available phenomenological model in EPSC [90]. The phenomenological backstress law has four fitting parameters, while the more physical one implemented in the present work has only one. The purpose is to provide insights into the profile and magnitude of backstress given the GNDs and fitting parameter. The magnitude of backstress in the given grain approaches 5 MPa. It should be noted that some grains have it slightly higher, while some lower. Also, multiple slip systems active in every grain increase the overall intragranular magnitude of backstress fields. The magnitude is adjustable with the backstress scaling parameter, k bs . Since the parameter is adjusted to provide good fit of the data, the magnitude should be reasonable. ) and predicted (Sim.) GNDs during pre-straining in biaxial tension (BT), plain strain tension (PST) and uniaxial tension (UT) followed by a comparison between measured and predicted GNDs in UT for AA6016-T4 (left). Corresponding predictions of SSDs (right).
This work developed a new strain-gradient elastoplastic selfconsistent (SG-EPSC) model integrating a dislocation density-based hardening law, an advanced composite grain twinning sub-model, and a backstress sub-model. The new development constitutes a major improvement of the EPSC model since its original formulation because it addresses the fundamental shortcoming of the mean-field self-consistent schemes where any spatial relationships of material properties and intragranular spreads are represented as averages. The development was made possible with the recent uplift of the EPSC model into the fieldfluctuations (FF)-EPSC model estimating the second moment of stress, lattice rotation rates, and resulting misorientation spreads. The SG-EPSC formulation is therefore a natural extension of the FF-EPSC formulation to estimate geometrically necessary dislocations (GNDs) from the misorientation spreads. To this end, a procedure for creating spatial arrangements of the orientation spreads in grains was developed to obtain a functional form of the rotation tensors for spatial derivatives. Furthermore, the right stretch tensors for orientations belonging to the granular spreads were attained to form the polar decomposition of the elastic deformation gradients. The spatial derivatives involving the "curl" operation were then taken to obtain the Nye dislocation tensor from the elastic deformation gradients for calculating the GNDs in individual grains. The spatial derivatives were strongly influenced by grid size. In contrast, the derivatives were slightly dependent on the spatial arrangement and density of grid points.
The model was successful in predicting and interpreting the grain size sensitive strength of α-Ti. Predicted texture evolution and twin volume fractions were in good agreement with corresponding experimental measurements. The good predictions were a consequence of the hardening law and the twinning sub-model being appropriately modified to incorporate the effects of GNDs. The model also predicted the evolution of strength and GNDs of AA6016-T4 during simple tension of pre-strained AA6016-T4 sheets in simple tension, biaxial tension, and plain strain tension. Finally, the load reversal flow curves of AA6016-T4 were predicted to agree well with the measurements verifying the backstress sub-model formulated based on GNDs.
The developed SG-EPSC model opened the door for a plethora of future model advancements and applications. Distinguishing between edge and screw dislocation types and evaluating merits is left for future works [44]. A possible approach can be taking gradients of slip rates like described in [96]. Future works will utilize other sub-models available in EPSC such as the martensitic phase transformations and recrystallization to evaluate the benefits of SG-EPSC. Additionally, the SG-EPSC will be applied to more benchmark studies such as shear/torsion including a wider range of grain sizes for further validation. Moreover, the SG-EPSC formulation is expected to improve the predictions of elastic lattice strains, residual stress, and diffraction peak broadening. The improved predictions of these phenomena will improve the modeling of ductility and geometrical shape changes under complex deformation boundary conditions in the FEM-SG-EPSC component-level modeling framework of processes like rolling, recrystallization, and deep-drawing while predicting microstructural evolution. U.S. National Science Foundation and the Center for Extreme Events in Structurally Evolving Material under U.S. Army CC-APG-RTP Division contract number W011NF-23-2-0073.
For completeness of the field fluctuations formulation, we provide below a few necessary partial derivatives specific to EPSC. Expressing the second rank 3 × 3 tensors as 6 × 1 vectors via Voigt notation, Eqs. (17b) and ( 18) are combined to take a derivative with respect to the misorientation vector, δr m , as: The partial derivatives of shear strain rate with respect to Schmid tensor and stress rate can both be derived from Eq. (19). Taking the partial derivative with respect to Schmid tensor gives:
Finally, the partial derivative with respect to stress rate is:
] .
(A4)
The detailed derivations for each term in Eqs. ( 21)-( 23) can be found in earlier works of [31,[39][40][41]58]. The partial derivatives re-derived for EPSC-specific formulations are provided in Appendix A.1. This section gives a summary of the expressions used to calculate the second moments.
The fluctuations due to mean orientation is calculated at each time step:
The second moment from Eq. (A5) can be decomposed into: 〈 σt-Δt,q ⊗ σt-Δt,q 〉 c = T t-Δt (T t-Δt ) T , (A6a)
where T t-Δt and T t are the lower triangular matrices from the decomposition of the second moment of stress rate fluctuation at a previous time step, t -Δt, and the current time step, t. Then, the linear map, Z t-Δt , of the stress rate fluctuations can be constructed as:
The fluctuations between the time steps are then:
σt,q = Z c σt-Δt,q . (A8)
The linear map from each step is then used to update the cross-covariance term as [39]:
where Y δr,c and Y σ ,c are:
Where R inc is the mean increment in rotation.
The partial derivatives of the spin with respect to stress rate and misorientation vector are calculated from the above-described partial derivatives via:
∂w ∂ σ = -∂w pl ∂ σ = -∑ s t s ⊗ δγ s ∂ σ , (A11a) ∂w ∂δr = -∂w pl ∂δr = -∑ s [ t s ⊗ ∂γ s ∂m s ∂m s ∂δR ∂δR ∂δr + γs ∂t s ∂δr ] . (A11b)
where t s k = -1 2 ε ijk q s ij is the dual vector of the anti-symmetric Schmid tensor, q s . The second moment of lattice rotation rate terms from Eq. ( 22) are:
σ ⊗ w t, σ 〉 c = ∂w t,c ∂ σt,c 〈 σt,q ⊗ σt,q 〉 c ( ∂w t,c ∂ σt,c ) T , (A12a) 〈w t,δr ⊗ w t,δr 〉 c = ∂w t,c ∂δr t,c 〈δr t ⊗ δr t 〉 c ( ∂w t,c ∂δr t,c ) T , (A12b) 〈w t, σ ⊗ w t,δr 〉 c = ∂w t,c ∂ σt,c 〈 σt ⊗ δr t 〉 c ( ∂w t,c ∂δr t,c ) T , (A12c) 〈 σt ⊗ δr t 〉 c = 〈 σt,q ⊗ δr t 〉 c + ∂ σt,c ∂δr t,c 〈δr t ⊗ δr t 〉 c . (A12d) Finally, the second moment of misorientation terms from Eq. (23) are: 〈 δr t inc ⊗ δr t inc 〉 c = Δt 2 4 〈w t, σ ⊗ w t, σ 〉 c , (A13a) 〈δr t,rot ⊗ δr t,rot 〉 c = R t,c inc 〈δr t ⊗ δr t 〉 c R t,c inc T , (A13b) 〈 δr t,rot ⊗ δr t inc 〉 c = R t,c inc 〈δr t ⊗ δ σt 〉 c ( ∂w t,c ∂ σt,c ) T Δt 2 + R t,c inc 〈δr t ⊗ δr t 〉 c ( ∂w t,c ∂δr t,c ) T Δt 2 . (A13c)
The grain fragmentation model in EPSC stems from the second moments of misorientation spreads describing the intragranular distribution of local orientations per grain deviating from the grain's average orientation. From the second moments of misorientation, eigenvalues and eigenvectors, 〈δr t ⊗ δr t 〉 c , are calculated, sorted in the descending order, and denoted as
. The first components, λ 1 and v 1 , are associated with the direction with the largest variation. The misorientation angle in each grain is [40,41]:
When the misorientation angle reaches the critical value of 15 • , the grain fragments into two equal halves, labeled here as f 1 and f 2 . The mean misorientations of the fragments from the mean orientation of the fragmenting grain are obtained using [97,98]:
Then, the fragments' initial misorientation quaternions are:
δq c f1 = ⎧ ⎨ ⎩ ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅ 1δr c f1 ⋅δr c f1 √ δr c f1 ⎫ ⎬ ⎭ (B3a) δq c f2 = ⎧ ⎨ ⎩ ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅ 1δr c f2 ⋅δr c f2 √ δr c f2 ⎫ ⎬ ⎭ (B3b)
Finally, the mean orientations of the fragments are calculated using the misorientations of the fragments and the mean orientation of the grain:
q c f1 = δq c f1 q c (B4a)
The fragments are assigned their initial volumes, representative grain size, and shape following the deformed state of the parent grain. Moreover, the fragments inherit the second moments and other state variables from the parent. The fragments are then treated as independent grains in the same phase that continue to evolve their own properties and can fragment further if the criteria for fragmentation are fulfilled.