Stochastic 3D microstructure modeling of twinned polycrystals for investigating the mechanical behavior of γ -TiAl intermetallics

A stochastic 3D microstructure model for polycrystals is introduced which incorporates two types of twin grains, namely neighboring and inclusion twins. They mimic the presence of crystal twins in γ -TiAl polycrystalline microstructures as observed by 3D imaging techniques. The polycrystal grain morphology is modeled by means of Voronoi and –more generally– Laguerre tessellations. The crystallographic orientation of each grain is either sampled uniformly on the space of orientations or chosen to be in a twinning relation with another grain. The model is used to quantitatively study relationships between morphology and mechanical properties of polycrystalline materials. For this purpose, full-field Fourier-based computations are performed to investigate the combined effect of grain morphology and twinning on the overall elastic response. For γ -TiAl polycrystallines, the presence of twins is associated with a softer response compared to polycrystalline aggregates without twins. However, when comparing the influence on the elastic response, the presence of twin grains has a much smaller e ff ect than a statistically di ff erent polycrystalline morphology. Notably, the bulk modulus is almost insensitive to the grain morphology and exhibits much less sensitivity to the presence of twins compared to the shear modulus. The numerical results are consistent with a two-scale homogenization estimate that utilizes laminate materials to model the interactions of twins.


Introduction
At the engineering scale, most polycrystalline metals or alloys consist of a large number of individual crystalline domains, called grains, each typically ranging in size from 1 µm to a fraction of a millimeter.The mechanical behavior of each grain at the local scale, characterized by crystallographic orientation, is anisotropic, while at a larger scale, the mechanical response depends on the spatial arrangement of the grains in the microstructure, i.e., on the morphology of the grain system.The latter arises during fabrication and processing, which typically involves heat treatments and solidification [1].Overall, mechanical properties are influenced by various factors, including inelastic deformation due to the motion of crystalline defects, twinning or phase transformation, as well as brittle or ductile fracture.Even in the purely-elastic regime, the local stress-strain state in each grain under mechanical stress is a complex result of the load distribution within the microstructure, depending on both, the crystallographic texture and the morphology of the grain system.
Numerical approaches to model the behavior of polycrystalline materials traditionally rely on statistical descriptors such as orientation-and misorientation-distribution functions, as seen in the case of nickel-based superalloys [2].Misorientation-distribution functions statistically quantify the correlation between crystallographic orientations of neighboring grains, which is particularly crucial for metals exhibiting crystal twins.Here the symmetries of the crystal lattice determine the relative orientation of adjacent grains, separated by a twin plane, known as habit plane [3].Frequently, grains in polycrystalline materials are in multiple twinning relationships, forming so-called twin-related domains.More advanced statistical descriptors such as "micro-texture functions", that describe clustering with respect to crystal orientations [3], are necessary to model twin-related domains.
A large class of applications involves polycrystals with an overall isotropic mechanical response.This occurs in the important special case where grains exhibit isotropic distribution in shape and crystallographic orientation.Models for realistic polycrystalline materials, therefore, primarily rely on random tessellation of space, with the Laguerre tessellation being a common choice [4,5,6].Note that the Laguerre tessellation is a generalization of the Voronoi tessellation, taking the granulometry, or size-distribution, of the resulting grains into account.To account for curved grain boundaries, strong shape anisotropy or multiscale grain-size distribution, alternative approaches like the Johnson-Mehl model and its anisotropic extensions [7,8] can be used.Another alternative, the generalized balanced power diagram [9,10,11], offers more degrees of freedom compared to the Laguerre tessellation.However, the computational cost associated with numerical optimization of morphological descriptors, especially in 3D [12], limits the use of more general models.A discussion regarding the trade-off between model accuracy and model complexity in terms of morphological descriptors has been provided in [13].Furthermore, in mechanics, numerical results suggest a small effect of the grain morphology on the macroscopic behavior in the elastic regime, whereas crystallography can play an important role, e.g., in the context of thermoelasticity [14].
In linear elasticity, well-established analytical estimates of the self-consistent type [15] are sufficient to estimate the effective properties [16], even in strongly anisotropic cases and in the presence of cracks [17].These estimates may also facilitate accurate reconstruction of histograms for strain and stress tensors [18].Mean-field estimates have been used to account for twin-related domains, modeled as individual grains and their effect on the elasto-plastic strain-stress response [19,20].
Inclusion twins as well as neighboring twins are present in the TiAl intermetallic family, a material of aeronautical interest, characterized by high stiffness, good corrosion resistance, strength at high temperatures and low density.However, its brittleness and low ductility at room temperature hinder its processing and application in modern components.Proper modeling of this material at macroscopic scale requires the understanding of deformation mechanisms at the micromechanical scale, such as dislocation slip and mechanical twinning.These phenomena can be simulated numerically using various techniques, such as finite element analysis and fast Fourier transform (FFT) methods.In this context, the effect on the local mechanical fields of the presence and distribution of twins in the microstructure is not yet fully understood.This motivates the development of stochastic 3D models for the generation of polycrystalline microstructures presenting grains in twinning relations, to be used as input for spectral mechanical simulations.
In this paper, the homogenized elastic properties of stochastically generated representative elementary volumes (RVEs) are numerically investigated and compared to analytical results.Two types of tessellations, Voronoi and Laguerre tessellations, are implemented, and the presence (or absence) of twins in microstructures is explored through three configurations: no twinning, inclusion twinning, and neighboring twinning.The rest of this paper is organized as follows: Section 2 explores the microstructure and material, Section 3 describes the stochastic 3D model for polycrystals containing twins, Section 4 investigates the elastic response of various polycrystalline models and the effect of twins, and Section 5 concludes.

Material microstructure and experimental data
The mechanical properties of binary TiAl alloys can be improved by the addition of heavier elements such as Chromium and Niobium [21] like Ti 48 Al 48 Cr 2 Nb 2 , the alloy considered in the present paper.Depending on heat treatment and precise composition, this material can adopt very diverse types of microstructures [22].However, the focus of the present paper is on the so-called nearly gamma microstructure, consisting of equiaxed grains of the γ-phase, with small α 2 nodules potentially present but considered insignificant for this analysis.The single phase γ-TiAl is obtained through an isothermal hold just above the eutectoid point (∼ 1130 • C) according to the phase diagram established by McCullough et al. [23] and confirmed by several other authors [24,25,26].The obtained γ-phase exhibits an ordered Strukturbericht designation L1 0 structure [21], equivalent to a tetragonal distortion of a face-centred cubic (FCC) structure as shown in Figure 1a.The lattice constants are a = 0.405 nm and c = 0.411 nm as determined by X-ray diffraction measurements [27], i.e., very close to a = 0.400 nm and c = 0.407 nm as found in the same material, for instance, in [28].The ratio c/a = 1.013 indicates its high similarity to the FCC structure.
Electron back-scattered diffraction (EBSD) is employed as a technique of choice to characterize the polycrystalline microstructure of metallic materials in detail.The focused ion beam (FIB) interacts with the specimen surface and the resulting scattered electrons form a specific pattern which is captured by a camera and further indexed to measure the local crystallographic orientation.As a scanning electron microscopy (SEM) technique, it provides a map of the crystallographic orientations of each point on the scanned surface (represented as orientation data such as Euler angles or other orientation representations [29]).Typically, such data is visualized as 2D cross sections with inverse pole figure (IPF) coloring, as illustrated in Figure 1b.In the present study, a detailed reconstruction of the complex γ-TiAl microstructure is performed, see Figure 1c.The data set was obtained by Tribeam tomography performed at the University of California, Santa Barbara [30], using a FEI Versa microscope equipped with both, a FIB column and a femto-second laser.The TiAl sample was prepared with an approximate cross-section of 500 µm × 500 µm, whereas the laser was employed to ablate 1.5 µm of material in the direction perpendicular to the cross-section between each EBSD acquisition.EBSD maps were acquired with a spatial resolution of 1.5 µm to obtain a volume with an isotropic voxel size having an edge length of 1.5 µm.Further details of the tribeam experiment are reported in [27].The volume was processed with Dream3D [31] and several Python routines to obtain the final microstructure with individual segmented grains using an orientation threshold of 5 • .
Stochastic microstructure modeling as described in Section 3 is inspired by the observation and processing of the experimental microstructure.Figure 1b shows the two considered types of twinned grains, neighboring twinned grains, which are side by side and approximately equal in size, as well as inclusion twinned grains, elongated grains, surrounded by their corresponding "parent" grains.Note that both types of twins are annealing twins and identical from the crystallographic point of view.However, they are very different in a morphological sense, which justifies the adopted naming convention.Both types of twins receive specific attention in the model proposed in the present study, considering morphology and the relationships of crystal orientation between neighboring grains, see Section 3.2.
Considering the present γ-TiAl phase with an L1 0 tetragonal crystal structure, closely resembling the FCC structure (often indistinguishable in standard EBSD measurements unless employing high-resolution acquisition with dictionary indexing [32]), it is commonly approximated by an FCC lattice.In this data set, the crystal orientation is indexed accordingly.This prevents to precisely locate the c-axis in each grain but is largely sufficient to determine the type of twin relationship between the grains.
The 3D data set contains 1966 grains with an average-volume equivalent diameter of 31.37 µm.For each grain twinning relations were evaluated, revealing that the vast majority of twins (≈ 90%) are Σ 3 (see Section 3.2 for mathematical details about twinning relationships) as observed previously in [33].

Hierarchical stochastic 3D modeling of polycrystalline microstructures with twinning effects
In this section, we present a stochastic 3D model which has been particularly developed for polycrystalline microstructures with twinning effects as they appear in γ-TiAl intermetallics described in Section 2. The modeling procedure consists of two major steps.First, the initial polycrystalline grain morphology is modeled by random tessellations in the three-dimensional Euclidean space R 3 [34].Second, crystallographic orientations are assigned to the individual grains, including the modeling of neighboring and inclusion twins in the system of grains.
For modeling the initial grain morphology in the first step, we employ Voronoi and -more generally-Laguerre tessellations.These two types of tessellations are defined by a point pattern {s n } and a marked point pattern {(s n , r n )}, respectively, consisting of so-called generators g n = (s n , r n ), which sometimes will be called seed points of the grains, where s n ∈ R 3 and r n ∈ R. Thus, we utilize (marked) random point processes to model random tessellations which subdivide the three-dimensional Euclidean space R 3 into non-overlapping sets, each representing a single grain in the microstructure. 1In the case of Laguerre tessellations, the marks associated with the points of the underlying point pattern allow for controlling the size of individual grains.If all marks are equal, the Laguerre tessellation coincides with the special case of a Voronoi tessellation on the same point pattern.
When modeling crystallographic orientations in the second step, a challenging task is that twin relations between two grains not only depend on their crystallographic orientation, but also on the spatial orientation of their joint grain boundary.More precisely, we have to take into account that the spatial orientation of the grain boundary between a pair of twins coincides with the crystallographic habit plane, see Section 2.

Modeling grain boundaries by random tessellations
We now describe the first step of the modeling procedure, i.e., modeling the polycrystalline grain boundaries.To mimic the structure of polycrystalline materials, random tessellations are used to subdivide the threedimensional Euclidean space R 3 into sets with disjoint interior, which is common approach in literature [6,11,35,36].Note that in the context of stochastic geometry, these sets are usually called cells, but in applications to polycrystalline materials it is more convenient to call them grains.Each grain G n , n ≥ 1, consists of all x ∈ R 3 , which are closer (with respect to a predefined distance function d T ) to its generator g n than to all other generators g m with m n.Formally, the n-th grain G n is defined by Motivated by image data representing the polycrystalline structure of γ-TiAl intermetallics, we assume that all grains, except for parent grains of twin inclusions (see Figure 1b), are convex and the joint boundaries between two neighboring grains are planar.Due to this assumption, for modeling the initial grain morphology we use random Laguerre and Voronoi tessellations, which share the common feature of producing convex grains with planar boundary segments, i.e. convex polyhedra.The randomness of the grain architecture is introduced by a randomization of the (marked) point pattern of generators.For this purpose, Matérn hardcore processes are utilized, as described in Section 3.1.1below.A detailed description of the corresponding random tessellation models is then provided in Sections 3.1.2and 3.1.3.

Modeling of seed points by random point processes
At first we introduce a family of marked point processes which is used in order to model the seed points of grains, i.e., the random generators of a tessellation.In particular, for the spatial locations of generators, we use a Matérn hardcore process in R 3 , denoted by {S n }, with intensity λ h > 0 and hardcore radius r h > 0, see Section 5.4 of [34].In this model, the generators are a thinning of a homogeneous Poisson point process [37], which is performed such that the pairwise distance between generators is larger than the model parameter r h > 0. In the case of Laguerre tessellations, each point of the point process {S n } has to be equipped with a random mark.For this, let R n be a sequence of independent and identically distributed random variables with values in R, which are independent of {S n }.Then, the marked point process {(S n , R n )} induces the random generators of the Laguerre tessellation.It is important to note that a Matérn hardcore process is an isotropic point process and leads to isotropic tessellations.Inducing a Laguerre tessellation by an anisotropic point process instead would lead to an anisotropic tessellation.For more details to (marked) point processes, we refer to [34,38].

Voronoi tessellation
In the case of Voronoi tessellations, the distance function d T defining the grain system, see Eq. ( 1), is given by d Note that the second component r of the generator g = (s, r) is not involved in the definition of the distance d T .To determine the Voronoi tessellation induced by a given point pattern s 1 , s 2 , . . .∈ R 3 , we plug d T into Eq.( 1) and get that Less formally, a point x ∈ R 3 belongs to G n if the Euclidean distance between x and the seed s n is less than or equal to the distance between x and all other seeds s m .The grain boundaries consist of points which are equidistant to at least two seeds, see Figure 2(a) for the visualization of a Voronoi tessellation in 2D.

Laguerre tessellation
In comparison to Voronoi tessellations, Laguerre tessellations allow for more flexibility when modeling polycrystalline materials.In particular, they allow for better controlling the size of the grains through additive weights.The distance function d T leading to the case of Laguerre tessellations is given by d r for all x ∈ R3 and g = (s, r) ∈ R4 .Using this distance function for a given marked point pattern (s 1 , r 1 ), (s 2 , r 2 ), . .., we get that For generating random marked point patterns, we use a random marked point process {(S n , R n )}, where the points S n are modeled by a Matèrn hardcore process and the marks R n are modeled (independently of the points) as independent and identically distributed copies of R, a uniformly distributed random variable on the interval [r − , r + ] ⊂ R. A realization of a Laguerre tessellation in 2D is visualized in Figure 2

Modeling crystallographic orientations and twinning effects
In the present paper, given the γ-TiAl alloy of interest, we restrict ourselves to Σ 3 twin boundaries which constitute the vast majority of twins experimentally observed in the material.Twin boundaries are specific grain boundaries characterized by a certain misorientation (in the axis/angle sense) between adjacent grains. 2 However, a twinning relation depends not only on the misorientation between two neighboring grains.The spatial orientation of the grain boundary, often referred to as the habit plane, is also part of the definition of twins.
The crystallographic nature of the boundary differs depending on the crystal symmetry considered.In the cubic approximation, the rotation is 70.53 • around one of the ⟨110⟩ axes (or equivalently 60 • around one of the ⟨111⟩ axes if the minimal angle representation of misorientations is adopted) [40].In this case the habit plane of the grain boundary is a {111} plane with respect to both grains. 3In the tetragonal case (the actual crystal symmetry of the investigated γ-TiAl), the twin relation is a rotation defined by tan(θ/2) = a/c 4 which, in our case, leads to θ ≈ 69.79 • around either the [100] or [010] direction (which are equivalent with respect to the crystal symmetry).In this case, the habit plane of the grain boundary is either the (011) or the (101) plane of the tetragonal lattice (which are also equivalent) as reported in Table 1.
In Section 2, two different types of twins were described, namely neighboring twins and inclusion twins.Due to their differences in morphology, these two types of twins require different modeling techniques.This leads to the following two-step procedure of twin modeling.First, neighboring twins are generated.To introduce a pair of neighboring twins, two adjacent grains are chosen at random.Subsequently, their joint boundary plane is extracted and the orientations of the two grains are set with respect to the conditions of a twinning relation.After generating all neighboring twins, the remaining grains are assigned with a random crystal orientation.Afterwards, inclusion twins are generated, where the crystallographic and spatial orientation of an inclusion twin depends on the crystallographic orientation of its parent grain, i.e., the grain in which the inclusion is inserted.
The present section is divided into several parts.First, in Section 3.2.1,we restrict the stochastic 3D model to a cubic sampling window V ⊂ R 3 .Subsequently, in Section 3.2.2,we discuss some issues of isotropy.Finally, in Section 3.2.3,we transfer the introduced model onto a discrete grid V ⊂ Z 3 , which mimics the situation of (discrete) experimental EBSD data, described in Section 2. It is important to note that both, the continuous model and its discrete approximation are defined with periodic boundary conditions.This is consistent with the periodic boundary conditions applied for the mechanical simulations in Section 4. Furthermore, for referring to crystallographic planes and directions, we utilize common Miller indexing for the continuous model considered in Section 3.2.1,whereas for its algorithmic realization in Section 3.2.3we use a more practical notation related to coordinates in Euclidean space.
property symbol cubic symmetry tetragonal symmetry twinning angle Table 1: Details of the twining relationships depending on the crystal symmetry considered.

Model description
We now describe the crystallographic aspects of the model, more precisely assigning each grain with a crystallographic orientation on a cubic sampling window V ⊂ R 3 .This procedure also incorporates twin relations between adjacent grains into the virtual grain morphology.
Neighboring twins.At first, pairs of adjacent grains to be neighboring twins are determined.For this purpose, we realize a Poisson distributed random variable N neigh with parameter p neigh •N V /2, conditioned on N neigh ≤ N V /2, where N V denotes the number of grains within the cube V and p neigh ∈ [0, 1] is a model parameter which controls the number of neighboring twins within the virtual microstructure.We then successively choose random pairs of neighboring grains and equip them with a crystallographic orientation, such that a given grain cannot be chosen twice.This implies that just twin pairs, and not so-called twin-related domains, which are "chains" of neighboring twins, are generated in the model.
Thus, the following procedure is repeated until we have either N neigh twin pairs, or no neighboring grains are left, which are not equipped with a crystallographic orientation.At first, we choose a grain G k and one of its neighbors G ℓ at random.Furthermore, let ∂G k,ℓ = ∂G ℓ,k = G k ∩ G ℓ denote the joint grain boundary of G k and G ℓ .Note that the grain boundaries of all grains, derived by the before mentioned tessellations, consist of planar segments, which implies that ∂G k,ℓ is also a planar set (with probability 1).By v ⊥ ∈ R 3 we denote a normal vector of ∂G k,ℓ . 5ecall from Section 2 that for twins, their joint grain boundary is aligned with a habit plane E with respect to the crystal coordinate system of both grains, G k and G ℓ (for instance in the cubic system, the twin boundary is a {111} plane).This property is fulfilled if v ⊥ is aligned with the normal of the habit plane related to G k and G ℓ .To include this into the model, a random orientation O k is chosen for G k , such that the corresponding habit plane normal is aligned with v ⊥ .The orientation of its twin O ℓ is unambiguously determined (except for symmetrically equivalent orientations) by rotating O k by the twinning angle θ around a twinning axis h ∈ Θ, see Table 1.
A 2D visualization of neighboring twins in the case of a cubic symmetry can be found in Figure 3.Note that for a given joint grain boundary ∂G k,ℓ with its normal v ⊥ , the orientations of the two adjacent and twinned grains G k and G ℓ are uniquely defined up to an arbitrary rotation around v ⊥ (and symmetrically equivalent orientations).Inclusion twins.In general, after creating neighboring twins, some grains may not yet be equipped with a crystallographic orientation.They are assigned with orientations, which are sampled from the uniform distribution on the space of rotations SO 3 .Furthermore, independently of possible neighboring twin relations, each grain G can have an inclusion twin I with some probability p incl ∈ [0, 1].If grain G is selected to have an inclusion, the following three-step procedure is performed.Otherwise, if no inclusion is generated, the procedure is skipped and continues with the next grain.
First, the crystallographic orientation O I of I is determined by rotating the orientation O of the parent grain G around a random rotation axis h ∈ Θ at an angle of θ.Furthermore, let H ∈ E denote the corresponding habit plane, which is normal to h for cubic crystal symmetries and aligned with h in the tetragonal case.Recall that H is not only a crystallographic quantity, but (in a morphological sense) parallel to the joint grain boundary between G and I.To realize this morphological behavior, I is modeled by a dilation of the habit plane.
The second step, deriving the location and spatial orientation of I, is visualized in Figure 4(a).To begin with, we choose a random point x on ∂G, the grain boundary of G. Without loss of generality we assume that x = o.Recall that due to the chosen type of tessellations (Voronoi or Laguerre), the boundary of each grain consists of planar boundary segments.Let X ⊂ R 3 denote the boundary segment which contains x6 , and v ⊥ the unit normal vector of X, pointing from x to the interior of G.The direction of inclusion v is derived by an orthogonal projection of v ⊥ onto H, formally expressed by where P H : R 3 → R 3 denotes the orthogonal projection onto H.Note that v is normalized to a length of 1.
In the third step, the shape of I is specified.For this, let diam G (x, v) denote the directional diameter of G at x along v, which is defined by With some probability p split ∈ [0, 1], the grain G is split by I into two separate grains.In the splitting case, the relative depth D of the inclusion equals one.Otherwise, in the case when G is not split, D is sampled from the uniform distribution on the interval [d − , d + ], where 0 < d − < d + ≤ 1. Subsequently, we define the absolute depth as is the absolute thickness of I, depending on the volume equivalent diameter ρ(G) = 3 √ 6/π vol 3 (G), where vol 3 ( • ) denotes the 3-dimensional Lebesgue measure.Finally, I is constructed as the intersection of G with the sets H x and H δ , which are defined as follows.Let P ′ be a plane, normal to v and containing x. Furthermore, let P be a plane which is parallel to P ′ and contains dv. 7hen, H x denotes the half-space resulting from the partition of R 3 by P, containing x. Formally H x is defined by H δ is the dilated habit plane with a thickness of δ, mathematically expressed by The construction of I utilizing the half space H x and the dilated habit plane H δ is visualized in Figure 4(b).Note that after inserting an inclusion I into a grain G, the grain architecture of the initial tessellation is changed, i.e., a new grain is introduced (the inclusion I).In particular, a grain G might (if the inclusion twin spans over the whole grain G) also be split into two new grains (without changing the crystallographic orientation), which corresponds well to the experimentally observed microstructures.

Isotropy
Since an isotropic microstructure is assumed for the mechanical simulations in Section 4.1, it is noteworthy that neither the generation of neighboring twins nor inclusion twins affects the isotropy of microstructure in terms of crystallography and morphology.Recall that the crystallographic orientation of neighboring twins depends on the spatial orientation of their joint grain boundary.Since Voronoi and Laguerre tessellations are isotropic tessellation models 8 , the spatial orientation of the grain boundaries is also isotropic.This implies that the choice of neighboring twins is invariant under rotations and thus isotropic.Grains not assigned with a crystallographic orientation during the neighboring twin procedure, are equipped with a crystallographic orientation sampled from the uniform distribution on the space of rotations SO 3 , which also does not violate the isotropy.Additionally, the crystallographic orientation of inclusion twins is determined on the basis of crystallographic orientation of the parent grain.Consequently, the overall isotropy is not affected either.

Implementation and discretization of the stochastic 3D model
The stochastic 3D model for polycrystalline microstructures introduced in Section 3.2.1 is defined in the threedimensional Euclidean space R 3 .As geometry input for the numerical computation of mechanical properties in Section 4, discrete model realizations are generated on a cube V ⊂ R 3 intersected with the voxel grid, i.e. on V ∩ Z 3 .Thus, in this section, an algorithmic description for the generation of discrete model realizations is provided.Recall that the image data obtained by EBSD measurements is also given as voxelized data, see Section 2.
For each quantity considered in Section 3.2.1 its discrete counterpart is marked by a hat.For instance, a grain G is approximated by its voxel representation G = G ∩ V, where V = V ∩ Z 3 .Furthermore, we use the notation N m (z) for the m-neighborhood of a voxel z ∈ V with m = 18 or m = 125. 9Note that we use periodic boundary conditions and thus each voxel has the same number of neighbors.We also make use of different representations of rotations, namely matrix representations and axis-angle representations.For more details, we refer to [29].Furthermore, we use the convention that a crystallographic orientation O consists of a crystal symmetry (in the present paper either cubic or tetragonal) as well as a rotation R ∈SO 3 .
Neighboring twins.Let G k and G ℓ denote the voxel representations of a pair of arbitrarily chosen neighboring grains G k and G ℓ , which are assumed to share a Σ 3 twin boundary.Recall that the crystallographic orientations of G k and G ℓ depend on the spatial orientation of their joint grain boundary ∂ G k,ℓ as it aligns with its crystallographic habit plane.For Laguerre tessellations in the Euclidean space R 3 , the joint grain boundary ∂G k,ℓ between two neighboring grains G k and G ℓ is normal to the vector connecting their seed points s k and s ℓ .Thus, by defining v ⊥ = s k − s ℓ we obtain the normal v ⊥ of ∂G k,ℓ , which is also used as the normal v⊥ corresponding to the discrete grain boundary ∂ G k,ℓ .
In accordance with the description of neighboring twins G k and G ℓ given in Section 3.2.1, the crystallographic orientation R k of G k is determined by the composition is an orthonormal basis (ONB) of R 3 , chosen such that v 1 and v 2 are spanning ∂G k,ℓ .Because the vectors v 1 and v 2 are deterministically computed 10 based on v ⊥ , the rotation matrix R v has to be randomly rotated around v ⊥ such that R k does not depend on the particular choice of the ONB defining R v .This is accomplished by utilizing the rotation R, the axis-angle representation of which is given by (v ⊥ , ω), where ω is sampled from the uniform distribution on [0, 2π).Subsequently, the composition R • R v is aligned such that v ⊥ coincides with the normal of the habit plane concerning the crystal symmetry by applying R ′ .By composing R ∥ with R k , we derive R ℓ , the underlying rotation of O ℓ .Finally R k and R ℓ are transformed by a random rotation that is equivalent with respect to the crystal symmetry, by applying R equi (R j ), j = k, ℓ.We refer to Table 2 for the specific construction rules for each of these rotations.
face centered cubic symmetry face centered tetragonal symmetry Table 2: Notation used to construct the orientations of neighboring twins.The rotations are presented in the form of matrices, compositions of matrices, and as axis-angle tuples.Note that (• , • , •) indicates a vector in R 3 and not a crystal plane in Miller indices.For the sake of readability, we denote a uniformly distributed random variable X on a set S as X ∼ U(S ).
Inclusion twins.Before generating discretized inclusion twins, recall that each grain G not assigned with a crystallographic orientation in the previous step, is equipped with an orientation sampled from the uniform distribution on SO 3 .For this we use the statistical software package R (version 4.2.2) and the package rotations [41] to sample from the uniform distribution on SO 3 .The precise definitions of quantities related to the crystallographic orientation are not explicitly given in the following paragraph.Instead, these definitions are provided in Table 2 for cubic and tetragonal crystal symmetries.
To model a discretized inclusion I of G, a random twinning axis Θ ′ is simulated relative to the crystal symmetry of G.The orientation O I of I is then determined by the rotation R I = ( Θ, θ), where Θ = R G • Θ ′ with R G , the rotation corresponding to the orientation O G of G. Similar to generating neighboring twins, the orientation of I has to be identified with an equivalent rotation by applying R equi (R I ), see Table 2.Note that the spatial orientation of the discrete habit plane H is directly determined from R I and R G (Table 2).Furthermore, H is represented by its normal ĥ, where R I and ĥ are analytically exact and equal to the corresponding quantities in the continuous model.
After establishing the crystallographic orientation of the inclusion I, the subsequent steps concern its spatial arrangement.Recall that in contrast to neighboring twins, the construction of inclusions requires a modification of the morphology of the grain, which was initially generated by a Voronoi or Laguerre tessellation.Unlike the continuous representation of inclusion twins described in Section 3.2.1,now we are limited to a voxelized representation of the grain architecture, which requires a discrete approximation of the inclusion on the voxel grid.
Let ∂ G = {x ∈ G : there exists a y ∈ N 18 (x) with y G} be the boundary of G, i.e. ∂ G = {x 1 , . . ., x N } ⊂ V for some integer N ≥ 1.Furthermore, let X = x j be a random voxel on ∂ G, where j is sampled from the uniform distribution on {1, . . ., N}.To get the normal v⊥ of ∂ G in X, we perform a principal component analysis on the point cloud x ∈ N 125 X : x ∈ ∂ G .Doing so, we obtain the three components v 1 , v 2 and v⊥ , see e.g.[42].The vectors v 1 , v 2 indicate the two main directions the extension of the point cloud, i.e. the local grain boundary ∂ G in X.The vector v⊥ is normal to v 1 and v 2 .Without loss of generality, we assume that v⊥ is pointing into G.The direction of inclusion v is derived by projecting v⊥ orthogonal onto the habit plane H by v = v⊥ − ĥ, v⊥ ĥ ĥ, where ĥ denotes the normal of H. Without loss of generality, we assume that v is a unit vector of length 1.The approximated diameter of G in X along v is defined by where ⌈ • ⌋ denotes rounding to the closest integer.Furthermore, let ρ( G) = 3 6 # G/π be the volume equivalent diameter, where # G denotes the cardinality of G.The absolute depth d and thickness δ are determined analogously to their continuous counterparts described in Section 3.2.1.Finally, we are able to define I = G ∩ H X ∩ H δ , where face centered cubic symmetry face centered tetragonal symmetry vector representation habit plane normal ĥ = Θ ĥ = R G • ĥ′ Table 3: Notation used to construct the orientations of inclusion twins.The rotations are presented in the form of matrices, compositions of matrices, and as axis-angle tuples.Note, that the for the representation of axes, (• , • , •) indicates a vector in R 3 and not a crystal plane in Miller indices.For the sake of readability, we denote a uniformly distributed random variable X on a set S as X ∼ U(S ).

Specification of model parameters
The study of different microstructures in Section 4 is carried out on the voxel representation V of a cubic sampling window V with an edge length of 512 voxels.Recall that the generated structures adhere to periodic boundary conditions.To generate microstructures according to the stochastic 3D model outlined in Sections 3.1 and 3.2 we have to initially generate a random pattern of seed points by a Matérn hardcore process.This process is defined by a hardcore radius r h = 5 and an intensity λ h = 15360/512 3 (≈ 1.14 • 10 −4 ) which corresponds to an expectation of 15360 seed points within the sampling window V.While a Voronoi tessellation does not need any additional information, the additive weights of the Laguerre tessellation are sampled from the uniform distribution on the interval [−16, 16].
For both tessellation types, we consider three different configurations of twinning parameters.While these configurations are typically not observed in experimental data, they serve as extreme cases of twinning occurring in realistic materials.(i) "Neighboring twins" (NT), where we put p neigh = 1, p incl = p split = 0. Furthermore, the number of neighboring twin pairs N neigh is drawn from a Poisson distribution with parameter where N V denotes the number of grains in V, while ensuring N neigh ≤ N V .Notably, when p neigh ≈ 1 (particularly in this twin configuration) there are less than N neigh twin pairs established in general.The latter occurs due to consistency reasons as it is rather unlikely that such a high predefined number of neighboring grains can be realized in the algorithm.(ii) "Inclusion twins" (IT), where we put p neigh = 0, p incl = p split = 1, i.e., each grain G within the tessellation is assigned with an inclusion twin, effectively splitting G into two distinct parts.(iii) "Without twins" (WT), where we put p neigh = p incl = p split = 0, i.e., the microstructure does not exhibit any twins.
Altogether, a total of 12 scenarios is investigated in Section 4 consisting of two types of tessellations (Voronoi and Laguerre), three twinning configurations (NT, IT and WT) and two crystal symmetries.

Effect of twinning and grain morphology on the elastic response
In this section, mechanical properties of γ-TiAl polycrystals with and without twins are investigated in the elastic regime.To address this problem, numerical computations are carried out utilizing a spectral solver, based on fast Fourier transforms (FFT), proposed by Moulinec and Suquet [43].These methods rely on an implicit integral equation for the strain field and the associated Green's operator.FFT schemes are a class of computationally efficient algorithms designed to compute the physical response (mechanical, thermo-mechanical, conductive) of heterogeneous media under macroscopically periodic boundary conditions.Applying this method to linear elastic theory, strain and stress tensors are defined on a regular grid, i.e. on the voxel grid.Coupled to modern experimental techniques such as 3D microtomography and FIB imaging, this approach is powerful in apprehending not only homogenized properties but also local fields [44] and offering substantial numerical speedup compared to finite element methods [45].
For the sake of simplicity, we adopt infinitesimal strain theory in the present study.Moreover, a discrete Green's operator is applied, which is consistent with centered differences on a rotated grid together with the "direct scheme" presented in [46].Alongside numerical computations, we present self-consistent type estimates for polycrystals containing crystallographic twins.These estimates are explicitly derived and compared to FFT predictions.

Linear elastic problem and FFT-based computations
Consider the following linear elasticity problem: where σ(x) denotes the Cauchy stress tensor, C(x) is the stiffness tensor of the local crystal and ε(u(x)) = grad sym (u(x)), denoted as Green-Lagrange strain tensor, equals the symmetric gradient of u(x).The strain and stress fields are defined on a cubic domain V ⊂ R 3 , and extended, by periodicity, to R 3 .Equivalently, by imposing periodic boundary conditions on the boundary ∂V of V, the strain field ε(u(x)) can be split into a prescribed average E and an unknown periodic term ε(u * (x)) given by where u * denotes the periodic part of the displacement field, and ⟨ • ⟩ spatial averaging over V. Furthermore, we prescribe where # denotes a Q-periodic field and -# an anti-periodic field.The effective properties are defined by the overall stiffness tensor C * defined by Referring to Section 3.2.2, it is important to note that the considered models are statistically isotropic, implying isotropy of the tensor C * [16].For periodic media, this condition is satisfied asymptotically when the domain V surpasses the size of the grains significantly and can be considered as representative volume element (RVE) [47].This was numerically verified for the presented model via FFT computations.
A macroscopic strain field is applied along 6 independent strain directions enabling to compute all 21 components of the stiffness tensor C * .Furthermore, we define the effective bulk modulus κ and shear modulus µ as the closest isotropic elastic tensors in the Euclidean sense [48,49,50].Utilizing single-crystal constants from [51,52], the stiffness tensors of Hooke's law (2) for cubic C C and tetragonal C T crystal symmetries are given by where the corresponding stiffness constants c i j are presented in The numerical computations are carried out on three realizations of each model, namely Laguerre and Voronoi with the three twinning configurations neighboring twins, inclusion twins and without twins (Section 3.3), which results in 18 samples.

Self-consistent estimates
As expected, the Voigt and Reuss bounds, also known as Hill bounds in the context of polycrystals [53], significantly overestimate or underestimate the elastic moduli.Likewise, Hashin-Shtrikman (or Hashin-Shtrikman-Walpole) bounds offer only marginal improvements.For instance, considering the stiffness tensor for tetragonal elastic symmetry class, the Hill bounds provide 109.760GPa ≤ κ ≤ 109.778GPa and 72.3 GPa ≤ µ ≤ 79.1 GPa for the bulk modulus κ and shear modulus µ, whereas Hashin-Shtrikman bounds are slightly more narrow, e.g.72.4 GPa ≤ µ ≤ 79.0 GPa.The determined bulk and shear moduli are shown in Figure 5.
In the following, we instead focus on self-consistent estimates.Closed-form expressions, given as a set of implicit equations, are available for polycrystals with cubic and tetragonal symmetries, see, e.g., [53].To take the effect of twins into account, correlations between neighboring grains have to be taken into account, which is in general not possible for self-consistent estimates based on an isolated inclusion, known as Eshelby's inclusion, embedded in the effective medium.Consider a polycrystal with a laminate substructure, as depicted in Figure 6.This polycrystal comprises grains of two types: monocrystals with uniformly distributed crystallographic orientations (represented in white), and grains exhibiting a laminate substructure, where the layers are in twin relation to each other (blue and red striped grains).
The use of laminates or sequential laminates, for modeling polycrystalline materials is the basis of numerous works, in particular for constructing optimal structures that attain the Voigt and Reuss bounds [54].In another context, Dusthakar et al. [55] make use of a model utilizing laminates for predicting the ferroelectric response of polycrystalline materials with tetragonal crystal symmetry, where laminates with anisotropic layers present notable properties.Goldstein et al. [56] explored notable properties in the context of two-layered plates composed of crystals with cubic crystal symmetry.Their study revealed that for certain misorientations, the Young modulus in the normal direction surpasses that of individual plates.
In the present study, we initially examine the effective elastic tensor of a rank-one laminate consisting of two layers denoted as ℓ 1,2 , where the direction, normal to the layering is denoted as n.The stiffness tensors C ℓ 1 and C ℓ 2 are both given with respect to the reference frame of ℓ 1 .While, this is the canonical choice for C ℓ 1 , the tensor C ℓ 2 has to be transformed by applying a rotation, consistent with the twinning relationship, see Section 3.2.1.Following [57,58] this rotation is obtained by a product of 6 × 6 matrices involving the stiffness tensor, expressed in Voigt notation, and a 6 × 6 "rotation" matrix.The stiffness tensor of the laminate C ℓ is given by [59, Eq. 7.2.7]: 1) , t (1) ) Z(t (2) , t (2) ) Z(t (1) , t (2) ) where (n, t (1) , t (2) ) is an orthogonal basis and f 1 , f 2 = 1 − f 1 are the volume fractions of the layers ℓ 1 and ℓ 2 .The computation of the above formula involves only identifying C ℓ 1 ,ℓ 2 as 6 × 6 matrices, where T ∈ R 6×3 and R ∈ R 3×3 .Specify f 1 = f 2 = 0.5 in Eq. ( 3), the result for the tetragonal case is given by C ℓ,α (1 ≤ α ≤ α max = 2) and for the cubic case, by C C ℓ,α (α = α max = 1).Here α max ∈ N denotes the number of different habit planes with respect to crystal symmetry, which are {111} for cubic, as well as (011) and (101) for tetragonal.Thus, we are left with estimating the effective response of a polycrystal containing three (tetragonal case) or two (cubic case) types of grains, namely the laminates with stiffness tensors C ℓ,α plus the non-twined grains with stiffness tensor C. The orientation of all grains is uniformly distributed on the space of crystallographic orientations and we assume all twining relationships are equiprobable, so that the relative proportions of each of the two types of laminates is fixed to 0.5 for tetragonal symmetry.Subsequently a self-consistent estimate is computed by solving [17] f where f ∈ [0, 1] is the fraction of non-twined grains and the averages are taken over a set of rotations, sampled from U SO 3 , the uniform distribution on the space of rotations SO 3 .The tensor Q ∈ R 3×3 is an isotropic compliant tensor with moduli given by with the shear modulus G and Poisson coefficient ν associated to C * . 11To solve Eq. ( 4), we employ a fixed-point method.The process begins with the mean between the Voigt and Reuss bounds as the initial value, iterating until convergence is achieved.Eq. ( 4) provides a self-consistent estimate for crystals with general triclinic symmetry, where the effective response is statistically isotropic.This feasibility arises because the Eshelby tensor, and consequently Q, remains independent of the elastic tensor within the grains, relying solely on the effective tensor C * .For self-consistent estimates applicable to polycrystals with arbitrary symmetry where the effective response is isotropic, we refer to Kube and Arguelles [60].

Numerical results of FFT-based computations
For cubic crystal symmetry, the bulk modulus is the same for every configuration (κ cubic = 110.333GPa).In contrast, the shear modulus depends on both, the grains morphology and the type of twins (   to both,the self-consistent estimates and FFT predictions, the softest response is provided by polycrystals without twins.Numerical and analytical results for the tetragonal crystal symmetry are presented in Figure 5.The most significant differences between the models (tessellation model and twin configuration) can be observed in shear, with differences up to 0.3 GPa.In contrast, there is approximately no impact on the bulk modulus (≈ 1 MPa).It is noteworthy that the presence of different types of twinning significantly influences the bulk modulus, which can be seen by the alignment of diamonds, circles, and triangles, each lying approximately on a horizontal line.The highest bulk modulus is attained in the absence of twins, while the lowest is achieved for inclusion twins.Note, that the type of tessellation has no significant influence on these values.Table 6 displays the dispersion of the values plotted on Figure 5.For a given twin and tessellation model, the accuracy of the FFT results for the bulk modulus is much lower than the differences observed between the mechanical response of the various models.The standard deviation for the shear modulus, however, is of the same order as the observed differences.
As expected, the predictions of the self-consistent estimate (4) without twins (i.e.f = 1, end of orange line) are very close to the Voronoi and Laguerre model in the absence of twins (red and blue diamonds).The presence of twins (0 < f < 1) lead to a decrease of both the effective bulk and shear moduli, consistently with FFT predictions for either neighboring or inclusion twins.
In order to interpret these results, it is useful to recall the construction of Avellaneda [61], leading to polycrystals with the stiffest properties in terms of bulk modulus.A sufficient and necessary condition is that along all grain boundaries of normal n, the following relationship holds: where C (1) , C (2) are the grains' stiffness tensors, i.e., n is a null eigendirection of I : C (1) − C (2) .This follows from the observation that under hydrostatic loadings, a uniform (across the two grains) strain field results in the Voigt bound.

Conclusion
The stochastic 3D microstructure model for twinned polycrystals developed in this work utilizes Voronoi and Laguerre tessellations of the three-dimensional Euclidean space R 3 for modeling polycrystal grain architectures.The main advantage of these tessellation models is that they provide planar grain boundary segments which are adequate for modeling twin grains.Furthermore, this novel modeling approach allows for incorporating two different types of crystallographic twins, namely inclusion twins, occurring within grains, and neighboring twins, adjacent grains in a twinning relationship.
Additionally, realizations of the aforementioned model served as input for full-field Fourier-based computations to investigate the combined effect of grain morphology and twinning on the overall elastic response.Our analysis showed that the presence of twins in γ-TiAl intermetallics is associated with a slightly softer response, especially in shear.This softening effect of twins has already been observed in the context of plasticity for other polycrystalline materials [19].However, it is important to emphasize that this conclusion depends on the crystal symmetry and the induced twin relationship.On the other hand, it turned out that the type of tessellation (Voronoi or Laguerre) that models the morphology of the grains has less impact on the mechanical response compared to the presence of twins.Finally, we have derived a two-scale homogenization approach that provides predictions as a function of the fraction of neighboring twins, or inclusion twins.The prediction of these analytical estimates are consistent with full-field numerical results obtained by FFT methods.
While the stochastic 3D model introduced in the present paper was utilized to investigate the influence of grain morphology and twinning on elastic response, it holds potential for various applications in virtual materials testing.For instance, in [62] the model was employed to analyze the influence of twins on plastic deformation of polycrystals.Furthermore, beyond enabling the model to incorporate twin-related domains, the development of methods for calibrating the model to grain morphology and crystallographic texture of experimentally measured image data would allow for generating virtual polycrystalline materials with twinning effect that are statistically similar to those observed by 3D imaging.This would come along with reducing the microstructural information to the values of parameters in the stochastic model and opens possibilities for the quantification of process-structure relationships for polycristalline materials with twinning, as it was done in previous studies, e.g., solar cells [63] and battery electrodes [64].

Figure 1 :
Figure 1: Microstructure of the γ-TiAl material: (a) L1 0 tetragonal crystal structure very close to FCC with a measured ratio c/a of 1.013; (b) EBSD measurement (IPF coloring) of a polished surface of the material showing the numerous twins present (the material has been indexed as cubic here); (c) 3D microstructure reconstruction (IPF coloring) obtained by stacking several hundred EBSD scans after a serial sectioning experiment [27]. (b).

Figure 2 :
Figure 2: Comparison of different tessellation models in 2D.The black dots indicate the seed points s n , n ∈ {1, . . ., 10} of the tessellations.(a) Each seed point of a Voronoi tessellation generates a convex grain.Boundaries between two neighboring grains are exactly in the middle between two seeds.(b) A Laguerre tessellation has also convex grains with planar boundary segments.However the joint grain boundary associated with two neighboring seed points is not necessarily between the two seeds, see the dark green and lower yellow grains.Moreover, not each generator g = (s, r) induces a grain (no blue grain).

Figure 3 :
Figure 3: 2D sketch of a neighboring twin relation for cubic crystal symmetries.The crystallographic orientations of G k and G ℓ are rotated by π/3(= 60 • ) around an axis v ⊥ (red) to each other.In a morphological sense, v ⊥ is normal to the joint grain boundary (gray) of the twinned grains G k and G ℓ .

Figure 4 :
Figure 4: 2D sketch of the construction of an inclusion twin relation.(a) Let R 3 be decomposed into two half spaces by the plane P, which is normal to v and contains x + dv.H x is defined as the half space containing x.(b) H δ is derived by dilating the habit plane H to a thickness δ.Intersecting G, H δ and H x defines the inclusion I. (c) The parent grain G and its inclusion I are in a twin relation at their two main joint grain boundary segments.

Figure 5 :
Figure 5: Effective elastic response of polycrystals with tetragonal elastic symmetry class.FFT results: Laguerre (L, blue symbols) and Voronoi models (red symbols), with neighboring twins (circles), inclusion twins (triangles), and without twins (diamonds).The orange line represents self-consistent estimates for varying twinned volume fraction.

Figure 6 :
Figure 6: Schematic example of a laminate polycrystalline microstructure.Red and blue striped grains represent the stiffness obtained for different twinning configurations, where the habit plane of the boundary either (011) or (101).White grains represent grains without twinning.

Table 4 :
Stiffness constants in GPa for cubic and tetragonal crystal symmetries.

Table 5 )
. According are 6 GPa and 12 GPa, which are far from the lowest eigenvalue.

Table 6 :
Standard deviation values (in GPa) for bulk modulus κ and shear modulus µ obtained for the three realizations of both tessellation models of polycrystals with and without twinning for tetragonal symmetry class.