Accessing Colony Boundary Strengthening of Fully Lamellar TiAl Alloys via Micromechanical Modeling

In this article, we present a strategy to decouple the relative influences of colony, domain and lamella boundary strengthening in fully lamellar titanium aluminide alloys, using a physics-based crystal plasticity modeling strategy. While lamella and domain boundary strengthening can be isolated in experiments using polysynthetically twinned crystals or mircomechanical testing, colony boundary strengthening can only be investigated in specimens in which all three strengthening mechanisms act simultaneously. Thus, isolating the colony boundary strengthening Hall–Petch coefficient KC experimentally requires a sufficient number of specimens with different colony sizes λC but constant lamella thickness λL and domain size λD, difficult to produce even with sophisticated alloying techniques. The here presented crystal plasticity model enables identification of the colony boundary strengthening coefficient KC as a function of lamella thickness λL. The constitutive description is based on the model of a polysynthetically twinned crystal which is adopted to a representative volume element of a fully lamellar microstructure. In order to capture the micro yield and subsequent micro hardening in weakly oriented colonies prior to macroscopic yield, the hardening relations of the adopted model are revised and calibrated against experiments with polysynthetically twinned crystals for plastic strains up to 15%.


Introduction
After decades of academic and industrial research, γ-based fully lamellar titanium aluminides (TiAl) outperform most competing high-temperature lightweight materials up to temperatures around 800 • C [1][2][3]. Their exceptional properties originate from the dense arrangement of three types of microstructural boundaries, namely lamella, domain and colony boundaries [4,5]. Despite advances in understanding the three corresponding Hall-Petch effects, their relative contributions to the strength of fully lamellar TiAl is not yet consistently quantified. Micromechanical modeling helps to separate and thus quantify these effects as shown in the following.  The tetragonal L1 0 lattice of the γ phase exhibits four twinning systems 1 6 112], four ordinary slip systems 1 2 110] and eight super slip systems, namely all systems involving a c component. Plastic deformation of the hexagonal D0 19 lattice of the α 2 phase is carried by prismatic, pyramidal and basal slip. Table 1 lists respective deformation mechanisms for both phases.

Gamma Gamma
Al Ti [0001] 11 20 1126 [001] [010]  Table 1. Slip and twinning systems in γ and α 2 phase with morphological classification according to [6]. The lattices of the α 2 and γ lamellae in each colony are aligned following the orientation relation (111) γ (0001) α 2 allowing coexistence of six orientation variants of the γ lattice that manifest in the domain structure, shown in Figure 1 (see, e.g., [7] for details). Separation distances of lamella, domain and colony boundaries are denoted by λ L , λ D and λ C throughout this paper.

Influence of Microstructural Boundaries on Strength
There is general agreement that colony, lamella and domain boundaries all give rise to Hall-Petch strengthening [4,5,[8][9][10][11][12][13][14][15][16][17]. Thus, the yield strength of fully lamellar TiAl is a function of three Hall-Petch slopes K and the inverse square roots of corresponding microstructural lengths: The Hall-Petch slopes K L and K D for lamella and domain boundary strengthening can be separately determined from experiments with polysynthetically twinned crystals (i.e., specimens that only contain parallel lamellae of one specific orientation) [8,11,18,19] or by micromechanical testing [20,21]. Colony boundary strengthening, however, can exclusively be investigated in specimens in which naturally all three strengthening mechanisms act simultaneously. Isolating the corresponding Hall-Petch slope K C thus requires a sufficient number of specimens with different colony sizes but (at least nearly) the same lamella thickness and domain size, hard to produce even with sophisticated alloying techniques [13,22]. Thus, only few of the reported K C values [10,13,14,16,17,23,24] were determined meeting this demand. In fact, the domain size λ D was given in neither of the cited references, but, since it is suspected to be correlated to lamella thickness λ L [12,25], specimens with the same λ L should exhibit similar values of λ D . However, the K C values determined for (nearly) constant λ L [13,14,16]-and thus constant λ D -still differ significantly for different values of λ L and λ D . Understanding colony boundary strengthening as the dislocation pile-up stress, required to activate slip/twinning in an adjacent colony (cf. Figure 3), supports the idea of a functional relation between colony boundary strength and strength of respective slip/twinning systems as suggested in [12]. Figure 3. Illustration of dislocation pile-up stress at a colony boundary possibly activating slip/twinning in the adjacent colony. For λ I I C = λ I C but λ I I L < λ I L and λ I I D < λ I D , the strength of respective systems in the right image will be higher than in the left one, requiring a higher pile-up stress to be activated and thus making colony boundary strengthening a function of λ L and λ D .
Consequently, this implies that • K C is a function of λ L and λ D , since strengths of slip/twinning systems in adjacent colonies are determined by lamella and domain boundary strengthening and • experimentally determined K C values are only valid for the given combination of λ L and λ D , rendering identification of the functional relation K C = f (λ L , λ D ) unreasonably labor-intensive.
With these implications, reported experimental values for the Hall-Petch slope K C may be regarded as a few points in a space spanned by K C , λ L and λ D , insufficient to reveal the functional relation K C = f (λ L , λ D ).

Scope of Present Paper
Other than in experiments, the choice of any specific combination of colony size, domain size and lamella thickness is not restricted in numerical simulations. Therefore, a well designed micromechanical model can be used to reveal K C = f (λ L , λ D ), finally enabling microstructure sensitive prediction of macroscopic yield stress of fully lamellar TiAl alloys.
While capturing selected micromechanical effects like, e.g., the plastic anisotropy of polysynthetically twinned crystals, the micromechanical models set up in the past (e.g., [13,[26][27][28][29][30][31][32][33][34][35][36][37][38][39]) are limited to isothermal conditions and/or temperatures near room temperature. Furthermore, the model parameters in these formulations are generally found for one specific combination of temperature and microstructural lengths, necessitating recalibration of the model whenever the mechanical behavior shall be analyzed at different temperatures or for a different set of microstructural parameters. Therefore, we set up a thermomechanically coupled crystal plasticity model of a polysynthetically twinned crystal that incorporates lamella and domain boundary strenghtening as well as the typical yield stress temperature anomaly to enable yield stress prediction between room temperature and 900 • C [18].
In the present paper, we aim to adapt this model to a representative volume element (RVE) of a polycolony fully lamellar microstructure in order to identify the functional relation K C = f (λ L , λ D ). Since a pronounced microyield has been observed in weakly oriented colonies prior to macroscopic yield (cf. digital image correlation (DIC) analyses in [40]), the hardening relations of the constitutive model [18] are revised to capture hardening of single colonies/polysynthetically twinned crystals up to several percent of plastic strain. The revised hardening model is subsequently calibrated against the uniaxial compression experiments with polysynthetically twinned crystals reported in [9].

Modeling
In this section, the thermomechanically coupled crystal plasticity model from [18] is revisited and extended in order to capture above mentioned microhardening as well as colony boundary strengthening.

Kinematics and Stress Measures
With the multiplicative split of the deformation gradient F into an elastic part F E and a plastic part the elastic representation of the right Cauchy-Green tensor C E can be defined as strain measure via and the second Piola-Kirchhoff stress S reads Here, J E = det F E is the Jacobian and σ the Cauchy stress. Correspondingly, the Mandel stress reads

Thermomechanics and Temperature Evolution
The procedure of thermomechanical coupling follows the line of arguments in, e.g., [41][42][43] and is thus just briefly recalled here. The Helmholtz free energy ψ is assumed to be a function of C E , absolute temperature θ and plastic internal variables q n . It is additively split into an elastic and a plastic part, reading As derived in, e.g., [41], introducing the relation between Helmholtz free energy ψ, internal energy ε and entropy η in the form ψ = ε − ηθ and inserting the corresponding time derivativė ψ =ε −ηθ − ηθ and subsequently the balance of internal energy into the Clausius-Duhem inequality, yields the following relations for stress and entropy: with ρ 0 being the density in reference configuration. Furthermore, the dissipation D reads In this, L P denotes the plastic velocity gradient and Q is the heat flux vector described via Fourier's law Q = −κ∇ 0 θ. The first two terms in Equation (9) represent the mechanical dissipation D mech during plastic deformation, while the last term denotes the thermal dissipation D therm due to heat conduction.
To obtain the temperature evolution equation, the time derivative of the Helmholtz free energy is inserted into the energy balance together with Equations (7) and (8). With the specific heat capacity with r being the external heat supply per unit mass. With evolving plasticity (q n > 0), D mech is lowered, i.e., energy is stored in the plastic internal variables q n instead of being dissipated as heat. If the plastic internal variables decrease, e.g., due to thermally activated recovery (q n < 0), the before stored energy is released in the form of heat.
Following the ideas from [44], mechanical twinning in the γ phase is treated as a unidirectional shear mechanism acting in the twinning plane and obeying Schmid's law, i.e., Equation (11). This leads to the following definition of the plastic velocity gradient L P (cf. [44]) Herein, shear rate of slip system α and the twinning rate of twinning system β are denoted by ν α and g β , the total twinned volume fraction is given by f and γ T denotes the material specific twinning shear. Other than in the original formulation of L P from [44], lattice reorientation and subsequent slip in twinned regions are neglected since twins in γ lamellae are generally very narrow [5].

Slip
The flow rule relates the resolved shear stress τ α to the resulting shear rate ν α on slip system α via the viscoplastic powerlaw (cf. [45]) In this, τ Y α is the current critical resolved shear stress and ν 0 and n are the reference shear rate and the strain rate sensitivity exponent.

Twinning
To ensure that twinning is unidirectional and that the twinned volume fraction does not exceed the theoretical limit of f = 1.0, the relation of resolved shear stress τ β and twinning rate g β is modeled via (cf. [44]) Herein τ T β denotes the current twinning resistance of twinning system β and the reference twinning rate is determined by dividing the reference shear rate ν 0 by twinning shear γ T [44,46].

Defect Density Evolution
In crystal plasticity, plastic deformation is represented by the accumulated shear of underlying deformation mechanisms instead of discretely resolving defects like dislocations or twins. The defect density evolution is, however, directly correlated to the stored energy of cold work, enables a consistent definition of dissipation and provides a physics-based way to describe work hardening and thermal recovery processes as, e.g., investigated in [15,47]. Therefore, we introduce the twinned volume fractions f β for twinning systems β and the dislocation densities ρ dis α for slip systems α as well as the total dislocation density ρ dis given by ρ dis = ∑ N sl α ρ dis α and the total twinned volume fraction f (see Equation (12)), which is determined via f = ∑ N tw β f β ≤ 1.0. Although not specifically intended here, a micromechanical model capable of tracking the evolution of these defect densities in TiAl alloys may, e.g., grant valuable insight into forming processes in which the defect density dictates necessary forming forces and the number of annealing steps [47].

Dislocation Density Evolution
The dislocation densities ρ dis α are assumed to evolve with the generation/recovery formulation (cf. [48][49][50] The first term on the right side represents the dislocation generation due to shear on slip system α (|ν α | > 0). Herein, A α is described by the saturation relation [48,49] where A α,0 is the reference accumulation coefficient, ρ dis α,sat is the saturation value for dislocation density and p α > 0 is a constant.
The second term in Equation (15) describes static thermal recovery and follows an Arrhenius type law (cf. [47,49]) where R α,0 denotes the reference recovery rate, Q R is the activation energy for recovery, k B denotes the Boltzmann constant, ρ dis α,min is the minimum dislocation density for recovery to take place and ρ dis ref is a reference dislocation density. The exponent q α > 0 is a constant.

Twin Evolution
In the context of thermomechanical modeling, it has to be noted, that-other than in conventional metallic materials-twinning is not only a room temperature mechanism in TiAl alloys. In fact, twinning is rather pronounced at elevated temperatures and even plays a significant role in creep deformation [5,51]. Therefore, twinning is not restricted to room temperature in this model. The twinned volume fractions f β are assumed to evolve directly with the corresponding twinning rates, i.e.,ḟ

Critical Resolved Shear Stresses
The slip and twinning system strengths τ Y α and τ T β from Equations (13) and (14) can be written as follows: Here, τ Y α,0 and τ T β,0 are the temperature and micro structure dependent initial slip and twinning system strengths and ∆τ Y α and ∆τ T β denote their evolution due to plastic deformation, i.e., represent work hardening.

Initial Slip/Twinning System Strength
As mentioned earlier, the yield strength of fully lamellar TiAl alloys is mainly determined by three coexisting Hall-Petch effects. Depending on its orientation with respect to the lamellae, a slip/twinning systems shear plane does either cross the lamella or the domain boundaries. Respectively, either λ L or λ D is the determining microstructural length for Hall-Petch strengthening [6,8,19]. The colony size λ C , however, has the same influence on all slip/twinning systems independent of their orientation. The initial slip/twinning system strengths can thus be written as with τ R being the lattice resistance to slip/twinning. The here introduced Hall-Petch coefficients k are defined on the slip/twinning system level and are not directly comparable with the measured Hall-Petch coefficients K. In polysynthetically twinned crystals, however, k D and k L may be determined from measured K D and K L values using the Schmid factors of selectively activated slip/twinning systems in certain orientations of the lamella plane with respect to load (cf. e.g., [18]). The last terms in Equations (21) and (22) are obviously only meaningful when modeling polycolony microstructures and thus have to be set to zero for simulations of polysynthetically twinned crystals. The temperature dependence of τ Y α,0 and τ T β,0 of the γ phase-especially the typical yield stress temperature anomaly-is incorporated in the parameters in Equations (21) and (22) as shown in [18]. The temperature dependent initial critical resolved shear stresses of the α 2 phase are modeled according to [18], thus not involving the Hall-Petch strengthening shown in Equation (21).

Evolution of Slip System Strength
Once the resolved shear stress τ α exceeds τ Y α,0 , the slip systems strength increases due to dislocation interaction and nucleation of twins. In the reported crystal plasticity models of fully lamellar TiAl, work hardening was usually treated in a simplified way using, e.g., linear [26,29,[31][32][33]35,37] or hyper secans [36,39,46] hardening laws, mostly neglecting the interaction with and between evolving twins. This proved sufficient for modeling the yield point of polysynthetically twinned crystals since, in this case, the plastic strains remain small. The micro yield in weakly oriented colonies, however, involves considerable local plastic strains. With the dislocation densities and twinned volume fractions, introduced earlier, we may thus set up a more physical, defect density-based work hardening description.
In general, hardening ∆τ Y α of slip system α can be expressed by where ∆τ Y α,s|s denotes strengthening due to dislocation interactions and ∆τ Y α,s|t denotes strengthening of slip system α by twin activity.
The slip|slip interaction ∆τ Y α,s|s is best described in terms of the dislocation density via the well-known relation ∆τ Y α,s|s = aGb α ρ dis .
Herein, G = E 2[1+ν] denotes the shear modulus determined in terms of Young's modulus E and Poisson's ratio ν and b α is the length of the Burgers vector of system α. The constant coefficient a ≈ 0.5 in TiAl alloys [5].
The boundaries of evolving twins β act as strong barriers to dislocation motion (cf. e.g., [52,53]), thus reducing the free path length of non-coplanar (ncp) slip systems α (i.e., n α ∦ n β ). Since the representation of twins by their volume fraction neither provides information about their number nor their thickness, the resultant strengthening of non-coplanar slip systems cannot be modeled as function of the inverse square root of the actual free path length as done in e.g. [54,55]. Thus, several authors introduced formulations to account for this source of Hall-Petch type strengthening as function of twinned volume fraction [27,52,56]. Based on the formulation from [27], we introduce the following relation for the corresponding term ∆τ Y with h αβ being a coefficient for hardening of slip system α due to twinning on non-coplanar system β.

Evolution of Twinning System Strength
The strength τ T β of twinning system β increases with nucleation of non-coplanar twins (i.e., n β ∦ n β ) and with interaction of twinning dislocations with the slip dislocation network [55]. In a general form, this can be written as introducing ∆τ T β,t|t as the strengthening due to nucleation of non-coplanar twins and ∆τ T β,t|s to account for interaction of twinning dislocations with slip dislocations.
Hall-Petch strengthening of twinning system β by non-coplanar twins β is modeled as for slip systems, reading Herein, h ββ is again a hardening coefficient. As described in [55], twinning dislocations may interact with slip dislocations, although this effect is not as severe as the other strengthening mechanisms. The corresponding formulation from [55] reads where C βα accounts for interaction between twinning system β and slip system α. Figure 4 qualitatively shows the evolution of slip/twinning system strength.  The hardening relations are not dependent on temperature, which is a reasonable assumption for the modeled defect structures, since they cannot be overcome by thermal activation [5].

Helmholtz Free Energy
As stated earlier, the Helmholtz free energy is assumed to be a function of the elastic Cauchy-Green strain tensor C E , temperature θ and plastic internal variables q n . The dislocation densities ρ dis α are well suited plastic internal variables, since they are directly related to the stored energy of cold work, enabling a physically meaningful representation of mechanical dissipation and temperature evolution.
While the stored energy of cold work increases linearly with dislocation density, the energy stored in the form of twin/matrix interfaces is given by the fixed interface energy per twin boundary (i.e., the stacking fault energy) and, thus, does not change with twin growth [5]. Therefore, the discrete number of twins has to be known to correctly assess their contribution to stored energy of cold work. Since resolving the discrete number of twins is not possible with the here presented continuum theory, we neglect the corresponding contribution to stored energy of cold work by assuming that the Helmholtz energy is no function of the twinned volume fraction.
With these assumptions, the Helmholtz free energy from Equation (6) is rewritten as The elastic part of the Helmholtz free energy ρ 0 ψ E is assumed to follow a Neo-Hookean behavior (cf. [18]) where µ = E 2(1+ν) and λ = νE (1+ν)(1−2ν) are the Lamé constants, α t denotes the thermal expansion coefficient, K = E 3−6ν is the bulk modulus and S 0 is the absolute entropy density. The plastic part of the Helmholtz free energy ρ 0 ψ P reads [48,49] With this definition of the Helmholtz free energy and relations (7) and (10), the stress and the temperature evolution read

RVE Generation and Discretization
In the following, two representative volume elements (RVEs) are set up to serve as geometries for the intended finite element analysis-one for a polysynthetically twinned crystal and one for a polycolony fully lamellar microstructure.

RVE of a Polysynthetically Twinned Crystal
As shown in [18,33], it is sufficient to discretize only seven lamellae to represent the geometry of a polysynthetically twinned crystal, i.e., one α 2 lamella and one for each γ orientation variant. Figure 5 schematically depicts the chosen representative volume element. This RVE is subjected to periodic boundary conditions and its rotation with respect to the uniaxial load is realized as described in [57]. The geometry is meshed using linear hexahedral elements.

RVE of a Polycolony Microstructure
Due to the high aspect ratio of the lamellae and their large number per colony, a one-to-one discretization of a polycolony fully lamellar microstructure is computationally highly inefficient. Therefore, numerical homogenization schemes and/or geometrical simplifications are inevitable.
We set up an RVE with a reduced number of lamellae per colony but ensure that the volume fractions of α 2 and γ lamellae are correctly reflected in each colony. The RVE is shown in Figure 6. The shapes of the 36 colonies in this RVE are based on a randomized 2D Voronoi diagram. Periodicity is achieved by repeating the diagram in a matrix pattern. The 2D Voronoi diagram is extruded in the z-direction, resulting in a columnar 3D representation of colonies. Subsequently, the lamella boundaries are introduced by Boolean intersection of appropriately oriented planes with the single columnar Voronoi cells. The orientation of the lamella plane in each colony is uniquely defined by only one angle ϕ i rotating around the z-axis. To minimize the influence of texture, the orientations of the colonies are evenly distributed between ϕ i = 0 • and ϕ i = 360 • (i.e., ϕ i = i 360 • n col for i = 1, 2, ..., n col ). Furthermore, the orientation of the γ phase along the lamellae is not altered, i.e., each lamella contains only one γ orientation.
This RVE is subjected to periodic boundary conditions and meshed using linear wedge elements.

Parameter Identification
The introduced defect density-based hardening model is best calibrated against experiments with polysynthetically twinned crystals or micropillar compression tests since both allow for investigating the micromechanical behavior of a single colony without the influence of neighboring colonies and colony boundaries. Although micropillar compression has the clear advantage of enabling the analysis of single colonies within the actual microstructure [20,21], currently more data is available on the plastic deformation of polysynthetically twinned crystals. In [9], the plastic deformation of polysynthetically twinned crystals is analyzed for eight different load angles between 0 • and 90 • and up to plastic strains of about 15% yielding a good data base for the here intended model calibration.

Morphological Classification
As stated in [6], all slip and twinning systems in fully lamellar TiAl can be uniquely classified according to their morphology to be either • longitudinal (s lamellar plane; n ⊥ lamellar plane), • mixed (s lamellar plane; n ⊥ lamellar plane) or • transversal (s ∦ lamellar plane; n ⊥ lamellar plane). systems (cf. Table 1). As shown in [6] and later confirmed by other authors, e.g., [6,18,26,35,57], assigning the same model parameters to all slip and twinning systems of a morphological class effectively reproduces the plastic deformation behavior of the lamellar compound while leaving only three parameter sets that have to be identified instead of one per individual slip/twinning system.

Hall-Petch Strengthening by Evolving Twins
Due to their different orientations with respect to the lamella plane, evolving longitudinal and transversal twins in the γ lamellae presumably contribute to a different extent to the strength of non-coplanar slip and twinning systems. Transversal twins evolve as many thin needles that subdivide the γ lamellae [58] and, thus, strongly reduce the free path length of longitudinal and mixed deformation systems. Other transversal deformation systems do, however, not necessarily have to cross these twins (cf. Figure 7), consequently being strengthened by them to a lesser extent. Little is known about the evolution of longitudinal twins since they are harder to capture in experiments because they cannot be distinguished from original lamellae after plastic deformation. However, it is very unlikely that several longitudinal twins nucleate within a single lamella since this is energetically unfavorable. Longitudinal twins will more likely grow from an already existing lamella boundary slightly reducing the lamella thickness but most probably not causing excessive strengthening of non-coplanar deformation systems (cf. Figure 7). Therefore, we distinguish strengthening by longitudinal and transversal twinning in the here derived parameter set.

Initial Critical Resolved Shear Stress
Since it has been frequently observed experimentally that super slip systems are less active than ordinary slip systems, they are suspected to exhibit a (slightly) higher critical resolved shear stress [4]. Unfortunately, the isolated critical resolved shear stress of super slip systems could not yet be accessed experimentally. Therefore, the contribution of super slip to the deformation of fully lamellar TiAl and appropriate modeling approaches have been discussed to some extent in the past. Some authors did, e.g., not allow any super slip in their models [59,60], introducing a pronounced tension/compression asymmetry [33] (due to the unidirectionality of twinning) that was not observed in experiment [61], while others scaled the strength of super slip systems by a factor Q so ≥ 1 with respect to the strength of ordinary slip systems [33,37]. Since there is no agreement on the actual strength of super slip systems, we choose a different approach here. To keep the benefit of Lebensohns morphological classification [6], the critical resolved shear stresses of all deformation mechanisms of a morphological class are assumed to be the same, neglecting the potentially higher stresses necessary to activate super slip.

Taylor Hardening
Lattice restoring super slip requires full 011] translations and may, e.g., be achieved by two identical 1 2 011] super partial dislocations separated by an antiphase boundary [5]. Thus, Taylor hardening of super slip systems can be described using the Burgers vector of either the two individual super partial dislocations 1 2 011] or the perfect super dislocation 011] in corresponding Equation (24). In order to account for the collective movement of the two super partial dislocations, bonded by the antiphase boundary, we use the Burgers vector of the perfect super dislocation. This results in a stronger hardening of super slip systems and, thus, reflects the fact that both super partial dislocations have to pass the forest of dislocations together. In consequence, ordinary slip is preferably activated in most crystallographic orientations as observed in experiments.

Recovery
Since, for the moment, we are only interested in (temperature independent) work hardening behavior of polysynthetically twinned crystals, the recovery related terms in Equations (15) and (17) are neglected in the following.

Model Parameters
Some parameters of the presented crystal plasticity model can be determined directly from experiments or have well-established standard values from literature. The Hall-Petch coefficients k D and k L as well as the parameters for the viscoplatic flow rules (Equations (13) and (14)) are taken from our previous work [18]. Respective parameters are gathered in Table 2. The α 2 content is not reported in [9], but, since the composition of the tested specimens is Al rich (Ti-49.3at %Al), we assume the α 2 content to be as low as 2 Vol. % for the calibration process.

Onset of Yield
With the Hall-Petch coefficients from Table 2, solely lattice resistance τ R , domain size λ D and lamella thickness λ L remain to be identified for determining the initial critical resolved shear stresses τ Y α,0 and τ T β,0 defined in Equations (21) and (22). Unfortunately, neither λ D nor λ L are reported in [9]. Furthermore, τ R depends on multiple factors like e.g., composition and impurity of a given alloy and can thus not be identified universally. To calibrate the hardening model against the given data, we thus set τ R = 25 MPa (based on experimental findings from [11]), λ D = 50 µm and λ L = 1 µm, which reproduces the initial yield of the experimental results from [9] reasonably well (see Figure 8).

Dislocation Accumulation and Hardening Interaction
With above introduced constitutive assumptions and the fixed initial critical resolved shear stresses, the remaining parameters for dislocation accumulation (Equation (16)) and hardening (Equations (25), (27) and (28)) are found by successively adjusting them in trial simulations until the post yield behavior of the polysynthetically twinned crystals tested in [9] is reproduced well. Respective model parameters are gathered in Table 3.

Results
Comparison of the simulation results with experimentally determined stress strain curves from [9] is shown in Figure 8.
Given the high number of deformation systems involved in the plasticity of polysynthetically twinned crystals-namely six γ orientation variants with 12 slip systems and four twinning systems each as well as 12 slip systems for the α 2 phase-the agreement between simulation and experimental results is very good. The qualitative features of the experimental curves are well reproduced, which is not possible with simple linear hardening relations as in [35,57]. We refrain from a quantitative improvement of the agreement between simulation and experiment by, e.g., numerical optimization methods, given the fact that in [9] only one specimen was tested per orientation.

Determining the Hall-Petch Coefficient for Colony Boundary Strengthening
With the calibrated defect density based model of a polysynthetically twinned crystal, the actual aim of the present paper is tackled: isolating the Hall-Petch coefficient for colony boundary strengthening K C and proving its dependence on lamella thickness λ L and domain size λ D .

Calculation Scheme
The Hall-Petch relation in its original form σ Y = σ 0 + KD −0.5 requires knowledge of σ 0 and Hall-Petch slope K to predict the yield stress σ Y of a polycrystalline material from its grain size D. In this, σ 0 represents the (theoretical) yield stress of a polycrystalline alloy with infinitely large grains, i.e., grain size D → ∞ and thus D −0.5 → 0. Transferring this idea to colony boundary strengthening, σ 0 represents the (theoretical) yield stress of a fully lamellar alloy with infinitely large colonies, i.e., without the influence of colony boundaries. Although such a microstructure can obviously not exist in reality, a model representation of it can still be created by directly applying the presented constitutive model of a polysynthetically twinned crystal-which exactly represents colonies without the influence of colony boundaries-to the polycolony RVE shown in Figure 6. Simulation results from such a model yield σ sim 0 as function of λ L and λ D and quasi represent a homogenization of a polysynthetically twinned crystals anisotropic yield stress. Being able to calculate σ 0 (λ L , λ D ), we may rearrange the Hall-Petch relation as follows: in order to determine the Hall-Petch slope K C (λ i L , λ i D ) individually for any combination of λ i D and λ i L from simulated σ sim 0 (λ i L , λ i D ) and corresponding experimentally determined yield stress σ exp Y (λ i L , λ i D , λ i C ). This is illustrated in Figure 9.
Repeating this for different combinations of λ L and λ D reveals K C = f (λ L , λ D ).
This interpolation scheme allows for determining K C (λ i L , λ i D ) for a certain combination of λ i D and λ i L by only one simulation and one experiment (although more experiments would obviously help to reveal the scatter of mechanical behavior) and if repeated for different combinations of λ D and λ L reveals the relation K C = f (λ L , λ D ).

Simulation Results
For determining σ 0 (λ L , λ D ) as described above, we assume τ R to be 30 MPa in the γ phase. Since the domain size is reported in neither of the references gathered in Table 4 but is suspected to be correlated to the lamellar thickness [12,25], we further assume λ D = 50λ L . This assumption seems reasonable evaluating lamella/domain size ratios reported for polysynthetically twinned crystals [8,19]. Figure 10a shows K C values determined via the above introduced scheme using the experimental results from [10,13,14] (cf. Table 4) plotted over λ L .
; open symbols: experimentally determined in [12][13][14]16]. K C value from [13] was determined from experiments with the indicated range of lamella thicknesses λ L . (b,c): comparison of experimentally determined σ exp Y (λ L , λ D , λ C ) [10,13,14] and simulated yield stresses σ sim Y (λ L , λ D , λ C ). In simulations, colony boundary strengthening was incorporated by introducing the interpolation from Figure 10a to Equations (21) and (22). Despite the scattering of K C values, the simulated yield stresses reproduce the experimental results very well.
The improvable correlation of the determined K C values in Figure 10a most likely results from the fact that neither domain sizes λ D nor α 2 volume contents were completely reported for the underlying experiments [10,13,14]. Still, the simulations successfully demonstrate that the colony boundary strengthening coefficient K C must be a function of lamella thickness λ L and domain size λ D as it was suspected in [12]. In order to finally incorporate colony boundary strengthening in the model via Equations (21) and (22), a functional relation for K C has to be established. Since usually λ D λ L and, thus, 1 , we expect that K C is considerably more sensitive to changes in λ L than in λ D . Thus, we assume in the following that K C is a function of λ L only. With this simplification, we choose the following interpolation of the calculated K C values in Figure 10a: In this, K C,0 corresponds to colony boundary strengthening in globular γ TiAl alloys (i.e., λ L = ∞) and is, thus, assumed to be K C,0 = 1 MPa √ m [12]. The linear coefficient K C,λ L is determined to K C,λ L = 4.5 × 10 −4 MPa[ √ m] 2 . The chosen functional relation K C (λ L ) has to be resolved to the slip and twinning systems in order to be used in Equations (21) and (22). Since no common Schmid factor can be determined for slip/twinning systems of arbitrarily oriented colonies, we use a factor of 0.3 to map K C (λ L ) to k C (λ L ). This yields: Including definition (36) into Equations (21) and (22) and rerunning the simulations with the reported colony sizes λ C from Table 4 yields a good qualitative and reasonable quantitative agreement with experimentally determined yield stresses as depicted in Figure 10b,c. These results indicate that the predicted yield stress σ sim Y (λ L , λ D , λ C ) is not too sensitive to the deviations between calculated K C values and used interpolation function (35) (cf. Figure 10a), which is evident since the inverse square root λ −0.5 C of colony size takes comparatively small values of 25 m −0.5 to 200 m −0.5 .

Conclusions
Most reported micromechanical models of polysynthetically twinned crystals were designed to capture their yield point and, thus, do not reproduce their hardening behavior for higher plastic strains. As opposed to previously reported models, the here presented hardening relations and interactions are able to reproduce all distinct features experimentally observed in [9]. This is associated with the incorporated strengthening effect through non-coplanar twins, which triggers transition between dominant deformation systems, causing changes in slope of the apparently linear hardening behavior. The presented hardening model is therefore helpful for further investigation of defect density evolution in polysynthetically twinned crystals and polycolony microstructures.
With the presented modeling approach, one can determine the Hall-Petch-coefficient K C to a good degree from very few experiments. Of course, the more data available, the better the computational determined value for K C . Although the accuracy of the here determined K C values suffers from missing domain sizes and unknown α 2 volume contents in used literature experimental data (cf. Table 4), we showed that K C depends on the lamellar microstructure and needs to be modeled at least as a function of the lamella thickness λ L . Depending on the desired degree of detail, the domain size λ D may be added as well, but can also be disregarded as its influence is much smaller.
In summary, for the introduced approach, only a very small number of lamellar microstructures with arbitrary combinations of λ L , λ D and λ C needs to be tested experimentally. This significantly simplifies the determination of the Hall-Petch coefficient K C = f (λ L , λ D ).