A Review on Capturing Twin Nucleation in Crystal Plasticity for Hexagonal Metals

Abstract

Owing to its ability to incorporate Schmid’s law at each integration point, crystal plasticity has proven a powerful tool to simulate and predict the slip behavior at the grain level and the ensuing heterogeneous stress/strain localization and texture evolution at the macroscopic level. Unfortunately, notwithstanding substantial efforts during the last three decades, this remarkable capability has not been replicated for materials where twinning becomes a noticeable deformation mechanism, namely in the case of low-stacking fault energy cubic, orthorhombic, and hexagonal close-packed structures. The culprit lies in the widely adopted unphysical pseudo-slip approach for capturing twin formation. While the slip is diffuse, twinning is a localized event that occurs as a drastic burst of a confined number of partial twinning dislocations establishing an interface that pursues growth through a thread of perfect twinning dislocations in the sense of bicrystallography. Moreover, at earlier stages, twin nucleation may require atomic diffusion (Shuffling) and faceting, generally demanding higher stress levels not necessarily on the twin shear plane, while triaxiality at adequate sites might be needed or preferred such as lower grain boundary misorientations or other twin boundaries. Identifying a mathematical framework in the constitutive equations for capturing these twin formation sensitivities has been a daunting challenge for crystal plasticity modelers, which has stalled ameliorating the design of key hexagonal materials for futuristic climate change-related industries. This paper reviews existing approaches to incorporating twinning in crystal plasticity models, discusses their capabilities, addresses their limitations, and suggests prospective views to fill gaps. The incorporation of a new physics-based twin nucleation criterion in crystal plasticity models holds groundbreaking potential for substantial progress in the field of computational material science.

1. Introduction

Twinning plays a prominent role in the plastic deformation and recrystallization of low symmetry crystals such as those showing hexagonal close-packed (hcp), orthorhombic, body-centered cubic and even face-centered cubic crystal structures holding a very low stacking fault energy, e.g., twinning induced plasticity steels. In hcp metals such as magnesium (Mg), twinning can be profuse depending on the alloy chemical composition, texture, strain rate, and temperature, while concomitant slip leads to complex interactions with existing twins in the deforming lattice, thereby exacerbating hardening and ultimately damage. Mitigating these interactions in materials design paradigms to ameliorate their properties for many futuristic applications has been a daunting challenge. For instance, wrought Mg alloys have not yet penetrated the automotive market to replace their aluminum and steel counterparts, notwithstanding intense research studies. This incites the need to develop a rigorous understanding of all the factors surrounding the origin of twinning, how it thrives in the lattice, and interacts with slip dislocations and other twins and defects. Crystal plasticity models, whether mean-field or full–field, have proven a good test for the depth of our understanding of the mechanisms driving plastic deformation in polycrystals [1,2,3,4,5,6,7,8,9,10,11,12,13,14]. A comparison of crystal plasticity predictions of the macroscopic and microscopic (case of full–field) behaviors under various conditions of external loading against experimental results, whether testing or digital image correlations, proved to be very helpful in identifying gaps in our understanding and calling for more small–scale modeling and characterization to reveal the associated missing mechanisms. Over the last decades, molecular dynamic (MD) simulations played a crucial role in identifying these mechanisms, which should be captured up in the crystal plasticity models [15,16,17,18,19,20,21,22]. This paper provides a review and a futuristic view of various scholastic works made to enhance crystal plasticity models in their prediction of the mechanical behavior of materials with twins. We emphasize efforts to incorporate twin nucleation in these models, which remains the most challenging event in this entire endeavor.

1.1. Twinning Mechanisms

Twinning in hcp metals is primarily due to the inability of the material to provide “easy” 〈c〉–axis slip deformation with both a low temperature and high strain rate. For instance, in hcp metals, the critical resolved shear stress of glissile pyramidal 〈c + a〉 slip mode is so high at room temperature that it exceeds that of $\left\{10\overline{1}2\right\}$ twinning, which then readily triggers when the external load requires tension/extension or compression/contraction of the 〈c〉–axis for $c/a$ ratio less or larger than $\sqrt{3}$, respectively.

A crystallographic description of twinning was originally given by Christian and Mahajan [23]. A twin is formed mainly by four symmetry operations:

• reflection in $K_1$;
• rotation of ${180}^{\circ }$ about ${\eta }_1$; and
• reflection in the plane normal to ${\eta }_1$; and
• rotation of ${180}^{\circ }$ about the direction normal to $K_1$.

Twinning, or more precisely compound twinning, in general, is the re-orientation of a fraction of the original lattice to the mirror image position due to a simple shear of the lattice. The plane about which this mirroring occurs is the $K_1$ plane. Figure 1 depicts the four crystallographic twinning elements. $K_1$ and ${\eta }_1$ are the invariant twin plane and the shear direction, respectively. $K_2$ is the conjugate twin plane, while P is the plane of shear. The intersection of $K_2$ and P results in the conjugate direction ${\eta }_2$. The shear plane P contains ${\eta }_1$, ${\eta }_2$, and the normals to $K_1$ and $K_2$. Based on these crystallographic elements, hcp twins can be classified as Type I ($K_1$ and ${\eta }_1$ are rational) or Type II ($K_2$ and ${\eta }_2$ are rational) even if all are rational sometimes. In hcp metals, a simple shear may not be enough for generating appropriate twinned crystal structures [23], and hence, an additional small atomic displacement, called shuffle [24], in a direction different from the lattice shear may be required. El Kadiri et al. [25] developed the first comprehensive crystallographic theory that provides with a unique sets of equation for the amount of shear and shuffle displacement for all compound twins in hcp lattices.

Based on the formation mechanism, twins can also be classified into 3 groups: (i) growth twins, (ii) transformation twins, and (iii) deformation twins [23]. Growth twins are formed during nucleation and growth processes, such as crystal growth from the vapor or liquid phases, and belong to special kinds of grain boundaries (GBs) [23,26]. Transformation twins are formed during solid state transformation, while deformation twins are formed due to the application of external stress and/or internal stresses, such as residual stresses. This paper focuses on the plasticity due to deformation twinning or, as termed by Christian, glide twinning.

Current twin nucleation models tuned for hcp metals can be divided into homogeneous and heterogeneous nucleation models. The homogeneous model relies on the nucleation of a twin domain through high stress concentration in a pristine parent domain [27], whereas the heterogeneous model accounts for the pre-existing twin defects as a cause for twin nucleation [27,28,29,30,31,32,33]. The experiments by Beyerlein et al. [33], Wang et al. [34] suggest that the latter type of twin nucleation models are superior. This review discusses both types of nucleation models and their applications.

1.1.1. Homogeneous Twinning

The first homogeneous twin nucleation model was introduced by Orowan [35], where he derived an expression for the minimum stable size of a twin nucleus in Zn. Price [36] studied differences between twinning in bulk crystals and in dislocation free whiskers and platelets and construed that homogeneous twinning requires higher stresses than heterogeneous twinning. Christian and Mahajan [23] calculated energy for homogeneous nucleation of twins in Zn and found a value of 75 eV, which was very large. This led to the conclusion that homogeneous twinning is hardly possible unless there is a very high stress and very low surface and strain energies. Lebensohn and Tomé [1] also found in their calculations that very high energies are required for homogeneous twinning. Such energy can not be provided by thermal fluctuations, thereby eliminating the chance of homogeneous twin nucleation. Using MD simulations, Barrett et al. [37] showed that $\left\{10\overline{1}2\right\}$ twinning is impossible to nucleate in pure Mg unless a void is introduced to the simulation box, which lowers the energy barrier for nucleation. In return, in a pristine lattice, $\{11\overline{2}1\}$ twinning was the only active twin mode, which is obviously not experimentally observed in Mg. The absence of atomic shuffling in $\left\{11\overline{2}1\right\}$ twinning suggests that the energy barrier for homogeneous nucleation is exacerbated by atomic shuffling. In fact, since shuffles correspond to the transport of atoms through diffusion mechanisms, a sufficient level of the gradient of hydrostatic pressure would be required to form a $\left\{10\overline{1}2\right\}$ embryo [38]. A void promotes triaxiality and, as such, the required gradient of hydrostatic pressure.

1.1.2. Heterogeneous Twinning

Bell and Cahn [28] found in their experiments that stresses required for twinning decreased as the specimen was slightly bent or contained an existing twin lamella. In another study, Bell and Cahn [39] mentioned that crystals indented under a stress twin immediately. Price [36] conducted experiments on Zn whiskers and platelets and found that in large crystals, twins are initiated in bursts and grow at fluctuating stresses and high speeds, while in nearly perfect crystals, they grow at almost constant stress and a controllable rate. Lebensohn and Tomé [1] stated that any microstructural feature acting as a stress concentrator or having an associated extra energy improves the chances of inducing twin nucleation. Experiments performed by Lay and Nouet [40] showed that for small twins, few defects are contained in the interfaces, and for large twins, many dislocations are generally found in the twin walls.

1.2. Twin Nucleation

Based on the experiments that showed nucleation of twin is easier with pre-existing defects in Zn crystals, Bell and Cahn [28] hypothesized that there are three stages in twinning: creation, spreading, and thickening of twins and these processes require stresses in decreasing order. Bilby and Entwisle [41] examined the stress fields around two types of inhomogeneities arising during plastic deformation (bounded slip bands and kink bands) and their influence on the formation of mechanical twins. The presence of plastic deformation is required for twin formation. Orowan et al. [35] derived a formula for the minimum stable size of a twin nucleus in zinc crystal by considering energies of dislocation loops and work effected by the applied stress. Bell and Cahn [39] proposed an argument that a twin nucleus must be created as a whole by a locally homogeneous shear of the lattice and not by progressive motion of dislocations. Chyung and Wei [42] studied the effects of stress concentrations on twin nucleation in Zn bicrystals and concluded that strongly locked dislocations, such as the leading dislocation of a pile-up against a barrier, may become a potential twin nucleus. A twin embryo can grow to a stable sized twin nucleus when sufficient stress is applied to it. Figure 2a shows the ∼ 30 nm thick $\left\{10\overline{1}2\right\}$ twin nucleus as observed in Mg AZ31 alloy using transmission electron microscope (TEM) [43]. The new twin nucleus has a lenticular shape with pointed twin tips that can be further observed in the high–resolution TEM images in Figure 2b,e. Efforts were made by Meyers et al. [44] and Fischer et al. [45,46,47,48] to quantify twin morphologies and twin stresses using thermodynamics and other energy-based approaches so that this information can be bridged to higher scale crystal plasticity models.

Hooshmand et al. [49] studied the interactions between 〈a〉, 〈c〉, and 〈c + a〉 dislocations with $\left\{\overline{1}011\right\}$ and $\left\{\overline{1}013\right\}$ twin boundaries in Ti using MD simulations and found that these interactions are twin growth mechanisms. They also found that $\left\{10\overline{1}2\right\}$ twin embryo could nucleate from $\left\{10\overline{1}3\right\}$ twin boundaries. Wang et al. [50] studied the nucleation of twins in hcp structures via MD simulations and identified twin nucleation through a mechanism they termed “pure shuffle” and growth via a glide-shuffle mechanism. Ghazisaeidi and Curtin [51] studied $\left\{10\overline{1}2\right\}$ twin nucleation in Mg using atomistic simulations and inferred that dislocation–assisted mechanisms for twinning in Mg are feasible, apart from twin nucleation at grain boundaries. Wang et al. [34] used atomistic simulations to study twin nucleation in Mg and found that twins are most likely to nucleate at grain boundaries—preferably more at low misorientation angle GBs than high angle GBs. Kumar et al. [52] studied $\left\{10\overline{1}1\right\}$, $\left\{10\overline{1}2\right\}$, and $\left\{10\overline{1}3\right\}$ coherent twin boundaries (CTB) and coherent basal prismatic (CBP)/coherent prismatic basal (CPB) boundaries in six hexagonal metals using the VASP code [53,54,55]. Several important conclusions were drawn, such as the creation of excess volume during twin formation, the solubility of solute atoms being better at twin boundaries than in the bulk, and the formation energy of $\left\{10\overline{1}2\right\}$ CTB being higher in metals with $c/a<\sqrt{8/3}$ and lower in metals with $c/a>\sqrt{8/3}$ than the formation energies of $\left\{10\overline{1}1\right\}$ and $\left\{10\overline{1}3\right\}$ CTBs. Another conclusion was that $\left\{10\overline{1}2\right\}$ CTB has higher energy than the CBP/CPB boundary for all metals. Wang et al. [17] conducted DFT calculations in Mg to study $\left[10\overline{1}1\right]$$\left\{10\overline{1}2\right\}$ twin nucleation in hcp metals. A $\left\{10\overline{1}2\right\}$ twin in an hcp lattice was found to have a minimum thickness of 17 crystallographic $\left\{\overline{1}012\right\}$ planes in their simulation. A two-layer twin is unstable, and the formation energy of twin boundaries decreases with increasing thickness of $\left\{10\overline{1}2\right\}$ twins. In another study, Wang et al. [18] studied $\left\{10\overline{1}2\right\}$ twinning mechanisms in hcp crystals using topological analysis and atomistic simulations. Normal and zonal twinning mechanisms were studied using DFT and empirical potentials for $\left\{10\overline{1}2\right\}$ twin nucleation. For the normal twinning mechanism, the stable twin nucleus consists of three twinning dislocations corresponding to a thickness of six $\left\{10\overline{1}2\right\}$ planes. For the zonal twinning mechanism, the stable twin nucleus consists of n (n is determined from DFT calculations) twinning dislocations corresponding to $(2n+1)$ crystallographic $\left\{\overline{1}121\right\}$ planes. Xu et al. [56] performed MD simulations of Mg to study twin nucleation and twin growth and argued that $\left\{\overline{1}\overline{1}21\right\}$–type twinning is the most preferential mode as their growth is led by the pure shear process. It was observed that increasing the stress can lead to secondary $\left\{\overline{1}\overline{1}22\right\}$ twinning, which migrates by pyramidal slip. All of these mechanisms and properties associated with twinning demonstrate the complexity of understanding and predicting the mechanical behavior of materials implicating twinning.

1.3. Twin Propagation

The formation and propagation of twins are greatly affected by the grain size and the texture in a polycrystalline material [57,58,59,60,61,62]. Antonopoulos et al. [63] conducted TEM observations in Zn specimens to observe the deformation twin modes. It was observed that dislocations with a weaker contrast are visible only along a segment of a twin boundary and responsible for twin growth. The effect of secondary twinning on twin growth was studied by Barnett et al. [64] using EBSD analysis on an Mg AZ31 alloy, who found no evidence supporting an effect. Out of the four $\left\{10\overline{1}1\right\}$ –$\left\{10\overline{1}2\right\}$ double twinning variants, the one with the misorientation relation of 37.5 degrees $<1\overline{2}10>$ was found to display high rate of lateral growth, and it is due to the tendency of this twinning variant to minimize the compatibility strain.

Beyerlein and Tomé [10] mentioned that twinning dislocations of the predominant twin system (PTS) are responsible for twin growth as they are often the most active twin dislocations. On the other hand, slip dislocations provide resistance to twin growth and eventually suppress it. Capolungo and Beyerlein [27] developed a 3D model of twinning in hcp materials. The chances of twinning dislocations breaking away and contributing to the formation of twins or twin growth are higher if twins are at a larger equilibrium distance from the reaction center. In other words, a twin can grow by the coalescence of twin nuclei or by the expansion of twin fault loops, or it can also grow by the creation of twinning dislocation dipoles. Capolungo et al. [65] investigated the interactions between slip dislocations and $\left\{10\overline{1}2\right\}$ and $\left\{11\overline{2}2\right\}$ twins. It was observed that a slip dislocation interacting with a $\left\{10\overline{1}2\right\}$ twin boundary helps twin growth if it dissociates into glissile twin dislocations. Twin nucleation and growth can also occur by dislocation pile-ups. TEM observations performed by Wang and Agnew [66] show interactions between dislocation and $\left\{10\overline{1}2\right\}$ twin. The interfacial interactions with an abundant amount of transmuted dislocations along the twin boundary, as shown in Figure 3, results in twin advancement and rapid hardening. Figure 2e shows various features, such as dislocation pileups at twin boundary, transmuted dislocations, and stacking faults observed in Figure 3b–d, respectively. Capolungo et al. [65] also performed experiments in strongly textured pure Zr samples under different loading conditions of temperature and compression. A conclusion was made that pre-existing $\left\{10\overline{1}2\right\}$ twin boundaries and prismatic dislocations do not have a significant effect on $\left\{10\overline{1}2\right\}$ twin growth, but these twin boundaries and dislocations can hinder the propagation of $\left\{11\overline{2}2\right\}$ twins. Meanwhile, using the climbing image nudged elastic band method, Tang et al. [67] showed the diffusive nature of the extension twin boundary and its athermal characteristics above critical stress. The faceted boundary plays a significant role in twin embryo nucleation and twin boundary mobility [16]. The twin may grow on boundaries other than the twin plane boundary and overcome obstacles or lower interface energy by means of faceting. Transmission electron microscope (TEM) experiments performed by Zhang et al. [43], as shown in Figure 2, show the microstructural features of a new twin nucleus with four representative areas shown by red boxes and the high–resolution TEM images in Figure 2b–e, respectively. The image in Figure 2d confirms the $\left\{10\overline{1}2\right\}$ twin orientation relationship with the matrix. The lenticular-shaped twin has a twin boundary that consists of straight terraces and step-like structures from CPB/CBP boundaries.

Upon the application of stress, the atoms move in an organized way to create defects such as dislocations, disclination dipoles, twinning, faceting, etc. At the continuum scale, these defects make it complex to formulate models that can capture or extrapolate those small–scale phenomena. Twinning mechanisms, such as nucleation, propagation, and growth, play a significant role in determining the final microstructure of the material. Site-specific embryonic twin nucleation in a homogeneous stress state is not yet deterministic, although Barnett et al. [60], Paudel et al. [68], Siska et al. [69] have done some work to estimate the nucleation of the next twin based on the stress state of the first one. Giri et al. [70,71] tried to predict twin formation susceptibility based on energy criteria using the nudged elastic band (NEB) method. A discrete dislocation-based modeling of a twin is performed by Lloyd [72] to study the stress state of a twin and predict the twin transfer at grain boundaries. The continuum-based crystal plasticity method does not directly comprehend the discrete information obtained from MD, and hence, it is necessary to bridge these two scales for an appropriate flow of information to better understand the material behavior. The next section discusses the CP methods and the evolution of CP models to incorporate the above-mentioned features of twinning in CP models.

2. Crystal Plasticity Modeling

2.1. Early Rate-Independent Crystal Plasticity Models

Numerical methods based on crystal plasticity theories have taken a pragmatic approach to understanding the macroscopic and microscopic behavior of different material systems. The early crystal plasticity theory can be traced back to Taylor and Elam [73,74] and Taylor [75,76], where a quantitative description of plastic flow due to crystallographic slip is developed. From the standpoint of modern continuum mechanics, Mandel [77] and Hill [78] formulated constitutive equations for elastoplastic behavior in single ductile crystals. These crystal plasticity theories were extended to finite deformations by Rice [79] and Hill and Rice [80]. Van Houtte [5] performed, for the first time, a rate-independent simulation to predict deformation textures considering the role of mechanical twinning along with crystallographic slips. All of these models use a rate-independent method for crystal plastic deformation, which has three main issues [81]:

1.

determining the slip systems that are active;

2.

determining the increments of shear on the active slip systems; and

3.

selecting slip systems required to produce an arbitrary deformation increment, which is not necessarily unique.

Peirce et al. [82] implemented the first two-dimensional numerical scheme for rate-independent elastic–plastic single crystals with non-homogenous deformation. This strategy lacked the determination of active slip systems and the amount of slip on these systems, leading towards the path of rate-dependent schemes.

2.2. Rate-Dependent Crystal Plasticity Models

Peirce et al. [83] developed a rate–sensitive CP model able to predict unique constitutive responses for arbitrary deformation histories. Asaro and Needleman [84] extended the rate-dependent crystal plasticity model to predict deformation textures and large-strain hardening behavior with varying stress–strain histories. For plastic deformation that occurs by crystallographic slip, a comprehensive constitutive relation was established by Asaro and Needleman [84] to model the evolution of crystallographic textures and the anisotropic stress–strain responses. These constitutive models used a semi–implicit time integration scheme developed by Peirce et al. [82,83]. Kalidindi et al. [85] developed a fully implicit time-integration procedure that rigorously predicted the stress–strain behavior and the texture evolution for both homogenous and non-homogenous deformation conditions.

2.3. Incorporation of Twinning

All of these crystal plasticity models are confined to single–phase, high stacking fault energy (SFE), cubic metals, in which plastic deformation occurs predominantly by slip. As mentioned in Section 1, for a number of other cubic metals and non-cubic metals, the deformation is dominated by deformation twinning that changes the stress–strain response [23,86,87,88] as well as the texture evolution [89,90,91]. Van Houtte [5], Chin et al. [92], and Tomé et al. [6] faced challenges of successfully incorporating deformation twinning in a crystal plasticity model. The major challenge was the proliferation of orientations created by the twinned regions, which are different from the untwinned regions of the grain. For instance, if each grain produces one twin at each time step, it will produce $2^n$ grain orientations after n time-steps. Van Houtte [5] resolved the problem by tracking the volume fraction of activated twins for each grain and selecting single orientations at random, where the probability of a given parent grain and twin variant is scaled by the current volume fraction. At the end of each time step, the number of grains remained the same. Tomé et al. [6] showed that this scheme has two main disadvantages: (i) the statistical criterion used in the scheme requires a large number of grain orientations and (ii) the reorientation decision at every time step does not consider the deformation history. To resolve this problem, Tomé et al. [6] proposed two different approaches to track the twin orientation at different time steps. The first approach is “the volume fraction transfer (VFT)” scheme. The VFT approach specifies a discretized space to represent an orientation space and tracks the weighted volume to represent the texture of the polycrystal. The volume fraction transfer between these discrete regions, which is proportional to shear accommodation, determines the crystal reorientation during deformation. The second approach is the predominant twin reorientation (PTR) scheme. Like Van Houtte’s scheme for tracking different orientations, the PTR scheme chooses the crystallographic orientation of either twinned region or untwinned region. However, the PTR scheme compares the accumulated volume fraction represented by the twin-reoriented grains with accumulated twinned fraction associated with the twinning shear as a criterion for making the reorientation choice.

2.4. Kalidindi’s Lagrangian Method of Incorporating Twin

Kalidindi [7] introduced a new constitutive framework incorporating deformation twinning in a polycrystal plasticity model together with an efficient time-integration scheme [85]. Unlike approaches by Van Houtte [5] and Tomé et al. [6], the calculation occurs in the relaxed configuration (see Figure 2) in which the lattice orientation of the twinned and the untwinned regions are pre-defined based on the initial lattice orientation of the crystal. This constitutive framework is the most comprehensive and amendable technique in the crystal plasticity finite element method (CPFEM). Therefore, this technique is widely used as a framework to construct new CPFEM models.

Staroselsky and Anand [93] proposed constitutive models to compute both texture and twin volume fraction, where a model based on the total Lagrangian approach to predict twin volume fraction and texture evolution was advanced by Kalidindi et al. [85], Kalidindi [94].

Wu et al. [95] presented a Taylor–type crystal plasticity model to predict the evolution of crystallographic texture and an anisotropic stress-strain curve during large plastic strains in Titanium. Clausen et al. [96] simulated the Van Houtte [5] scheme with grain proliferation since the technological advancement allowed to track all grains. In the study, Clausen et al. [96] maintained the relationship between parent and twin material although the original, precise twin crystallographic relationship between the two may no longer be strictly enforced. The stress field in a matrix–twin composite after the formation of a deformation twin has been analyzed by Knezevic et al. [97] using Kalidindi’s CPFE framework to study the conditions affecting the mobility of the twin-matrix interface. Lévesque et al. [9] developed a new framework incorporated with in-house codes by Inal et al. [98] to simulate large deformation behavior in hcp metals that considered the crystallographic slip and deformation twinning as principal deformation mechanisms. In the study, the crystallographic framework considered the dislocations inside the twinned region, which will be discussed in Section 3. Izadbakhsh et al. [99] proposed a new crystal plasticity constitutive model that incorporated plastic deformation due to primary extension twins, primary contraction twins, secondary (double) twins, as well as basal and non-basal slip system in all parent grain, primary twins, secondary twins, and double twins.

Kalidindi’s CPFE approach [7,85] was used by Abdolvand and Daymond [100] to investigate the internal strain and texture development based on various assumptions of twin–parent interactions. Abdolvand and Daymond [101], Abdolvand et al. [102] used the same model to simulate the effect of grain boundary geometry and texture on twin inception and propagation, as well as the 3D stress development in parent and twin pairs in Zircaloy-2. In addition, FEM with polycrystal constitutive description in each integration point is employed by [103,104,105,106] to investigate the deformation behavior of hcp metals. Wu et al. [12] developed a new twin nucleation, propagation, and growth (TNPG) model motivated by the studies of [107,108] on deformation mechanisms in amorphous glassy polymers. Similarly, Cheng and Ghosh [109] developed a physics-based crystal plasticity FE model of dislocation-mediated heterogeneous deformation and tensile twin nucleation in single-crystal Mg and the polycrystalline alloy Mg AZ31. They expanded their FE model with a new propagation approach based on the glide velocity of twin partial dislocation through shear and shuffling mechanisms [110].

Ardeljan et al. [111,112,113] used a 3D CPFE framework incorporated with a deformation twin through a novel approach, where twin lamella is automatically inserted across its parent grain when the twin volume fraction reaches the threshold value. The propagation process was implemented by Qiao et al. [114], who used the CPFE framework with Wu et al. [12]’s TNPG model and Tomé et al. [6]’s PTR scheme. Qiao et al. [114] noticed that stress relaxation from twin nucleation in the TNPG model played an important role for twin propagation across a grain boundary. Lévesque et al. [115] extended the framework developed in [9,98] to consider separate rigid body rotation of the twinned region and parent matrix. A stress relaxation associated with twin propagation accounted in the TNPG method was incorporated with the TDT model of Wang et al. [11,116] in a CPFE framework by Qiao et al. [114]. Figure 4 shows the evolution of the twinned area fraction during fatigue deformation behavior of a 2% precompressed Mg AZ31 sample observed by Hong et al. [117]. The microstructural changes during the 5/4 fatigue cycle show that the twinning in precompressed Mg AZ31 specimen accommodates deformation by the detwinning process. Reversed loading to compression causes twinning again, making an addition to the preexisting twins. The detwinning starts to operate from the first compressive unloading and runs through the second tensile loading with a final twinned area fraction of 9.7% reduced from initial 22% at the precompressed condition, as seen in Figure 4.

Several researchers Proust et al. [118], Knezevic et al. [119] have used the composite grain (CG) approach to determine the deformation twin in their crystal plasticity models. All these CPFE models for crystal plasticity of twinning are based on the framework developed by Kalidindi (1998), where twinning is considered as a pseudo-slip mechanism.

2.5. Fast–Fourier Transfer (CPFFT) Based Method

Finite element methods (FEM) possess limitations related to meshing (requiring fine meshes) and difficulties related to a large number of degrees of freedom, limiting the complexity and size of the microstructures that can be investigated. Similarly, the prediction of intracrystalline states using homogenization techniques requires n–site self–consistent approaches [120,121] instead of classical one-site formulation [122,123] in which ideal crystals are considered to deform embedded in a homogeneous medium with average properties [3]. These limitations of both the small–scale FEM and n-site self–consistent models can be overcome with a novel approach. As an alternative to the CPFE approach, a full–field approach using the crystal plasticity fast-Fourier transform (CPFFT) method was developed by Lebensohn [3], Lebensohn et al. [124], Prakash et al. [125] motivated by the work of Moulinec and Suquet [126], Michel et al. [127] on composite materials. Initially, a self–consistent elasto-viscoplastic model was developed by Hutchinson [128], Molinari et al. [129] and evolved to account for large anisotropic viscoplastic deformation by Lebensohn and Tomé [130], Lebensohn et al. [131]. Based on the FFT algorithm, Moulinec and Suquet [126,132] developed an iterative method to compute the effective properties and local response of elastic and elastoplastic composites. This method works with better numerical performance than CPFE for periodic heterogeneous microstructures to provide an exact solution of the governing equations. However, the convergence of this method is limited to few initial iterations for low rate sensitivity and strongly anisotropic materials [133]. This limitation was overcome by an application of an augmented Lagrangian method [127] on the FFT formulation for isotropic composites with a high contrast of properties. Since the CPFFT method is limited to periodic boundary conditions, CPFE models are suitable for more complex boundary conditions.

In 2008, Lebensohn et al. [124] used visco-plastic fast Fourier transform (VPFFT) to model the subgrain texture evolution in polycrystalline copper. Stress gradients and orientation gradients were captured using VPFFT, showing stress localizations near grain boundaries [134]. Extensions of the viscoplastic self–consistent (VPSC) model have been proposed by Proust et al. [4] and Knezevic et al. [119]. Lebensohn and Needleman [135] implemented a non-local polycrystal plasticity theory using FFT on three-dimensional FCC polycrystals to analyze the mechanical response of polycrystalline aggregates accounting for size dependency that arose from plastic strain gradient. In recent years, Barnett et al. [60], Siska et al. [69], and Paudel et al. [68] showed that the stress relaxation due to twin formation is significant during modeling of twin formation and growth in a crystal plasticity model. The effects of different microstructural aspects such as misorientation, low angle boundaries, autocatalysis, and twin aspect ratios during the growth and propagation of deformation twins are studied by Kumar et al. [69,136,137,138,139,140] using the crystal plasticity fast Fourier transform (CPFFT) method.

2.6. Coupled Crystal Plasticity and Phase-Field Model (CP-PFM)

In recent years, crystal plasticity models have been coupled with phase-field (PF) methods to study the twinning mechanisms [141,142,143,144,145], where twinning is modeled as a phase transformation, and plastic deformation is computed using crystal plasticity methods. Phase-field models, first started by Khachaturian [146], have been successful to predict austenite-martensite transformation in steels [147] and tetragonal–monoclinic transformation in zirconia [148]. In hcp Mg, Pi et al. [149] used the phase-field model to study the tensile twins using twin interfacial energy as a driving force. Later, the coupled CP-PF models are successfully used to study deformation behavior associated with slip induced plasticity, and twinning and detwinning [141], twin nucleation and thickening [142,143], twin morphologies [144], and twin interactions with other slip systems, twins, and grain boundaries [142,145].

All of these models have evolved with time to incorporate various deformation mechanisms observed at mesoscale and microscale levels. The main challenge for recent material scientists is to model physics-based twin nucleation criteria to incorporate into CPFEM or CPFFT models. In the following sections, different CP approaches as well as the current techniques developed to incorporate twins in these CP models are discussed.

3. Crystal Plasticity Finite Element Method (CPFEM)

The mathematical formulation for the crystal plasticity finite element method is based on the theory of continuum mechanics and the decomposition property of a deformation gradient. In this section, crystal plasticity formulation by Lévesque et al. [115] is summarized that uses an updated Lagrangian approach to incorporate the twinning mechanism. Throughout the formulation, bold letters and symbols are used for tensors and vectors unless defined otherwise in the text. Similarly, $\left(\right)·\left(\right)$ represents dot product, and $\left(\right)\otimes \left(\right)$ represents tensor product.

3.1. Total Deformation

The deformation gradient, F, can be decomposed into the elastic and the plastic part (see Figure 5).

$$\mathbf{F}={\mathbf{F}}^{*}·{\mathbf{F}}^P,$$

(1)

where F* is the elastic deformation. Rigid body rotation and crystallographic slip and twinning are included in ${\mathbf{F}}^P$. The velocity gradient can be written as:

$$\mathbf{L}=\dot{\mathbf{F}}·{\mathbf{F}}^{-1}={\mathbf{L}}^{*}+{\mathbf{L}}^P$$

(2)

where ${\mathbf{L}}^{*}={\dot{\mathbf{F}}}^{*}·{{\mathbf{F}}^{*}}^{-1}$, ${\mathbf{L}}^P={\mathbf{F}}^{*}({\dot{\mathbf{F}}}^P·{{\mathbf{F}}^P}^{-1}){{\mathbf{F}}^{*}}^{-1}$. The velocity gradient has the symmetric stretching part, D, and asymmetric part spin tensor, $\mathsf{\Omega }$, which can be decomposed into elastic and plastic parts.

$$\mathbf{L}=\mathbf{D}+\mathsf{\Omega },\mathbf{D}={\mathbf{D}}^{*}+{\mathbf{D}}^P,\mathrm{and}\mathsf{\Omega }={\mathsf{\Omega }}^{*}+{\mathsf{\Omega }}^P$$

(3)

where ${\mathbf{D}}^{*}$ and ${\mathsf{\Omega }}^{*}$ are elastic strain rate and rigid lattice rotation rate, respectively, and ${\mathbf{D}}^P$ and ${\mathsf{\Omega }}^P$ are plastic strain rate and plastic spin rate, respectively. A constitutive equation by Lévesque et al. [9], which is based on Asaro and Needleman [84] and Kalidindi [7], is presented in this subsection. The plastic strain rate and plastic spin for the crystal are written as

$${\mathbf{D}}^P=\left(1-\sum _{\beta =1}^{N^T}{f}_{\beta }\right)\left. [\sum _{\alpha =1}^{N^S}{\mathbf{P}}_{\alpha }{\dot{\gamma }}_{\alpha }+\sum _{\beta =1}^{N^T}{\mathbf{P}}_{\beta }{\dot{f}}_{\beta }{\gamma }^{\mathrm{tw}}\right]+\sum _{\beta =1}^{N^T}\left. [f_{\beta }\sum _{\alpha =1}^{N^S}{\mathbf{P}}_{\alpha }^{tw\beta }{\dot{\gamma }}_{\alpha }^{tw\beta }\right]$$

(4)

and

$${\mathsf{\Omega }}^P=\left(1-\sum _{\beta =1}^{N^T}{f}_{\beta }\right)\left. [\sum _{\alpha =1}^{N^S}{\mathbf{W}}_{\alpha }{\dot{\gamma }}_{\alpha }+\sum _{\beta =1}^{N^T}{\mathbf{W}}_{\beta }{\dot{f}}_{\beta }{\gamma }^{\mathrm{tw}}\right]+\sum _{\beta =1}^{N^T}\left. [f_{\beta }\sum _{\alpha =1}^{N^S}{\mathbf{W}}_{\alpha }^{tw\beta }{\dot{\gamma }}_{\alpha }^{tw\beta }\right]$$

(5)

where ${\mathbf{P}}_{\alpha }$ and ${\mathbf{W}}_{\alpha }$ are symmetric and skew-symmetric tensors for each slip systems $\alpha$, $\dot{{\gamma }_{\alpha }}$ is the shear rate on the slip system $\alpha$, $f_{\beta }$ is the volume fraction of the grains that have been twinned to the system $\beta$, ${\gamma }^{\mathrm{tw}}$ is the amount of shear associated with twinning, and $\dot{f_{\beta }}$ is the rate of twinning of the system. The parameters associated with the twinned region of the twin system $\beta$ are indicated by superscript $tw\beta$. $N^T$ and $N^S$ represent the number of slip systems and twinning systems, respectively. ${\mathbf{P}}_{\alpha }$ and ${\mathbf{W}}_{\alpha }$ for each slip system $\alpha$ can be obtained using the lattice vectors ${\mathbf{s}}_{\alpha }$ and ${\mathit{m}}_{\alpha }$ stretched and rotated to ${\mathbf{s}}_{\alpha }^{*}={\mathbf{F}}^{*}{\mathit{s}}_{\alpha }$ and ${\mathit{m}}_{\alpha }^{*}={\mathit{m}}_{\alpha }{{\mathbf{F}}^{*}}^{-1}$:

$${\mathbf{P}}_{\alpha }=\frac{1}{2}\left. [{\mathit{s}}_{\alpha }^{*}\otimes {\mathit{m}}_{\alpha }^{*}+{\mathit{m}}_{\alpha }^{*}\otimes {\mathit{s}}_{\alpha }^{*}\right]$$

(6)

$${\mathbf{W}}_{\alpha }=\frac{1}{2}\left. [{\mathit{s}}_{\alpha }^{*}\otimes {\mathit{m}}_{\alpha }^{*}-{\mathit{m}}_{\alpha }^{*}\otimes {\mathit{s}}_{\alpha }^{*}\right]$$

(7)

3.2. Crystallographic Slip and Twinning

Power law equations are used to describe the slip rate (${\dot{\gamma }}_{\alpha }$) and the rate of twinning ($\dot{f_{\alpha }}$) as follows:

$$\dot{{\gamma }_{\alpha }}=\dot{{\gamma }_{\circ }}sign\left({\tau }_{\alpha }\right){\left|\frac{{\tau }_{\alpha }}{g_{\alpha }}\right|}^{\frac{1}{m}}$$

(8)

$$\dot{f_{\alpha }}=\frac{\dot{f_{\circ }}}{{\gamma }^{\mathrm{tw}}}sign\left({\tau }_{\alpha }\right){\left|\frac{{\tau }_{\alpha }}{g_{\alpha }}\right|}^{\frac{1}{m}}$$

(9)

where $\dot{{\gamma }_{\circ }}$ and $\dot{f_{\circ }}$ are the reference shear rate of the slip system, ${\gamma }^{\mathrm{tw}}$ is a twinning shear, m is a rate sensitivity parameter, ${\tau }_{\alpha }$ is the resolved shear stress on the slip system $\alpha$, and $g_{\alpha }$ is the hardness of the slip system, whose rate is defined by hardening law:

$${\dot{g}}_{\alpha }=\sum _{\beta }{h}_{\alpha \beta }\left|{\dot{\gamma }}_{\beta }\right|$$

(10)

where $h_{\alpha \beta }$ is the hardening modulus. This modulus has the form

$$h_{\alpha \beta }=\sum _{\alpha =1}^{N^S}{q}_{\alpha \beta }{h}_{\beta }$$

(11)

where $h_{\beta }$ is a hardening rate of the slip system, and $h_{\alpha \beta }$ is the latent hardening matrix. The hardening rate for the slip and the twinning systems is represented by a power-law equation as follows:

$$h_{\alpha }=h_{\circ }{\left. [\frac{h_{\circ }{\gamma }_{\alpha }}{{\tau }_{\alpha }n}+1\right]}^{n-1}+h_1$$

(12)

where $h_{\circ }$ is the initial hardness of the slip systems, ${\gamma }_{\alpha }$ is the sum of the accumulated slip all time, and n is the hardening exponent.

$${\gamma }_{\alpha }={\int }_0^t|{\dot{\gamma }}_{\alpha }|dt$$

(13)

Hardness caused by the formation of twin boundaries is incorporated into the formulation with the parameter $h_1$. The hardening parameter $h_1$ depends on the density of the twin boundaries within a grain.

$$h_1=\left\{\begin{array}{cc}{h}_{TB}{\left(1-\left|\frac{\zeta -{\sum }_{\beta }{f}_{\beta }}{\zeta }\right|\right)}^e,\hfill & \mathrm{if}\sum _{\beta }{f}_{\beta }<=\zeta \hfill \\ \\ h_{TB}{\left(1-\left|\frac{\zeta -{\sum }_{\beta }{f}_{\beta }}{1-\zeta }\right|\right)}^e\hfill & \mathrm{if}\sum _{\beta }{f}_{\beta }>\zeta \hfill \end{array}\right.$$

(14)

where $\zeta$ represents the volume fraction of twins with maximum twin boundary density, and $h_{TB}$ is a coefficient accounting for the effect of twin boundaries on the hardening of slip systems. The model accounts for softening from lattice reorientations of twins by computing the resolved shear stress (RSS) on slip planes inside the twins using a new orientation. The reorientation of the twinned region is reflected on the elastic modulus using the transformation law:

$${\mathbf{E}}_{ijkl}^{tw\beta }={\mathbf{Q}}_{im}{\mathbf{Q}}_{jn}{\mathbf{E}}_{mnop}^{mt}{\mathbf{Q}}_{ok}{\mathbf{Q}}_{pl}$$

(15)

where $\mathbf{Q}$ is transformation matrix and is defined in the work of Van Houte [5] as follows:

$${\mathbf{Q}}_{ij}=2{\mathit{m}}_i{\mathit{m}}_j-{\delta }_{ij}$$

(16)

where $\mathit{m}$ is the twin plane normal.

Finally, the Cauchy stress, $\mathbf{\sigma }$, over a grain is obtained by averaging the Cauchy stresses in the twinned region (${\mathbf{\sigma }}^{mt}$) and un-twinned regions (${\mathbf{\sigma }}^{tw\beta }$) of the grain:

$$\mathbf{\sigma }=\left(1-\sum _{\beta =1}^{N^T}{f}_{\beta }\right){\mathbf{\sigma }}^{mt}+\sum _{\beta =1}^{N^T}{f}_{\beta }{\mathbf{\sigma }}^{tw\beta }$$

(17)

The Cauchy stresses, ${\mathbf{\sigma }}^{mt}$ and ${\mathbf{\sigma }}^{tw\beta }$, are computed using the elastic constitutive equation specified by

$$\stackrel{\mathsf{\Delta }}{\mathbf{\sigma }}=\mathcal{L}\mathbf{D}-{\dot{\mathbf{\sigma }}}^0-\mathbf{\sigma }tr\mathbf{D},$$

(18)

where $\stackrel{\mathsf{\Delta }}{\mathbf{\sigma }}$ is the Jaumann rate of Cauchy stress, $\mathcal{L}$ is the elastic modulus, and ${\dot{\mathbf{\sigma }}}^0$ is a visco-plastic type stress rate [150].

4. Fast Fourier Transform Method

In this section, we present the fast Fourier transform based elasto-viscoplastic formulation (EVPFFT) developed by Lebensohn et al. [151]. Based on the crystal plasticity approximation by Asaro and Needleman [84], the elastic strains are comparatively negligible to the plastic strains and the constitutive relation between the viscoplastic strain rate ${\dot{\mathbf{ϵ}}}^P\left(\mathit{x}\right)$ and stress $\mathbf{\sigma }\left(\mathit{x}\right)$ at a single-crystal material point is given by:

$$\begin{array}{cc}\hfill {\dot{\mathbf{ϵ}}}^P\left(\mathit{x}\right)& =\sum _{s=1}^{N^S}{\mathbf{m}}^s\left(\mathit{x}\right){\dot{\gamma }}^s\left(\mathit{x}\right)\hfill \\ \hfill & ={\dot{\gamma }}_{\circ }\sum _{s=1}^{N^S}{\mathbf{m}}^s\left(\mathit{x}\right){\left(\frac{|{\mathbf{m}}^s\left(\mathit{x}\right):\mathbf{\sigma }\left(\mathit{x}\right)|}{{\tau }_{\circ }^s\left(\mathit{x}\right)}\right)}^{n}sgn\left({\mathbf{m}}^s\left(\mathit{x}\right):\mathbf{\sigma }\left(\mathit{x}\right)\right)\hfill \end{array}$$

(19)

where ${\dot{\gamma }}^s\left(\mathit{x}\right)$, ${\tau }_{\circ }^s\left(\mathit{x}\right)$, and ${\mathbf{m}}^s\left(\mathit{x}\right)$ are, respectively, the shear rate, the critical resolved shear stress (CRSS), and the Schmid tensor, associated with slip system (s) at point $\mathit{x}$ and n is the stress exponent (inverse of the rate-sensitivity exponent). ${\tau }_{\circ }^s(\mathit{x},{\mathbf{ϵ}}^p\left(\mathbf{\sigma }\left(\mathit{x}\right)\right))$ is a function of the accumulated plastic strain of a crystal, which, in turn, is a function of the stress due to the strain-hardening in slip systems. To solve the EVP problem for small strain conditions, an Euler implicit time discretization scheme is used along with Hooke’s law in order to obtain the stress $\mathbf{\sigma }\left(\mathit{x}\right)$ at time $(t+\mathsf{\Delta }t)$:

$${\mathbf{\sigma }}^{t+\mathsf{\Delta }t}\left(\mathit{x}\right)=\mathbf{C}\left(\mathit{x}\right):{\mathbf{ϵ}}^{e,t+\mathsf{\Delta }t}\left(\mathit{x}\right)=\mathbf{C}\left(\mathit{x}\right):\left({\mathbf{ϵ}}^{t+\mathsf{\Delta }t}\left(\mathit{x}\right)-{\mathbf{ϵ}}^{p,t}\left(\mathit{x}\right)-{\dot{\mathbf{ϵ}}}^{p,t+\mathsf{\Delta }t}(\mathit{x},{\mathbf{\sigma }}^{t+\mathsf{\Delta }t})\mathsf{\Delta }t\right)$$

(20)

where $\mathbf{\sigma }\left(\mathit{x}\right)$ is the Caucy stress tensor, $\mathbf{C}\left(\mathit{x}\right)$ is the elastic stiffness tensor, $\mathbf{ϵ}\left(\mathit{x}\right)$, ${\mathbf{ϵ}}^e\left(\mathit{x}\right)$, and ${\mathbf{ϵ}}^p\left(\mathit{x}\right)$ are the total, elastic, and plastic strain tensors, and ${\dot{\mathbf{ϵ}}}^p$ is the plastic strain-rate tensor. They are related as follows:

$$\begin{array}{c}\hfill \mathbf{\sigma }\left(\mathit{x}\right)=\mathbf{C}\left(\mathit{x}\right):\left(\mathbf{ϵ}\left(\mathit{x}\right)-{\mathbf{ϵ}}^{p,t}\left(\mathit{x}\right)-{\dot{\mathbf{ϵ}}}^p(\mathit{x},\mathbf{\sigma })\mathsf{\Delta }t\right)\\ \hfill \mathbf{ϵ}\left(\mathit{x}\right)={\mathbf{C}}^{-1}\left(\mathit{x}\right):\mathbf{\sigma }\left(\mathit{x}\right)+{\mathbf{ϵ}}^{p,t}\left(\mathit{x}\right)-{\dot{\mathbf{ϵ}}}^p(\mathit{x},\mathbf{\sigma })\mathsf{\Delta }t\end{array}$$

(21)

where the supra-indices $t+\mathsf{\Delta }t$ are omitted, and only fields corresponding to the previous time step $\left(t\right)$ are explicitly indicated. [1,3,124,125,126,131,135,138,152,153,154,155,156] have evolved this method to simulate various deformation behavior in hcp metals.

4.1. Green’s Function Method

A stress–strain relation can be expressed in terms of linear reference medium with stiffness ${\mathbf{C}}_{ijkl}^0$ as:

$${\mathbf{\sigma }}_{ij}\left(\mathit{x}\right)={\mathbf{\sigma }}_{ij}\left(\mathit{x}\right)+{\mathbf{C}}_{ijkl}^0{\mathit{u}}_{k,l}\left(\mathit{x}\right)-{\mathbf{C}}_{ijkl}^0{\mathit{u}}_{k,l}\left(\mathit{x}\right)$$

(22)

where ${\mathit{u}}_{k,l}\left(\mathit{x}\right)$ is the displacement-gradient tensor with ${\mathbf{ϵ}}_{kl}=\frac{1}{2}\left({\mathit{u}}_{k,l}\left(\mathit{x}\right)+{\mathit{u}}_{l,k}\left(\mathit{x}\right)\right)$. The terms from the above equation can be regrouped to form a polarization field such that:

$$\begin{array}{c}\hfill {\mathbf{\phi }}_{ij}\left(\mathit{x}\right)={\mathbf{\sigma }}_{ij}\left(\mathit{x}\right)-{\mathbf{C}}_{ijkl}^0{\mathit{u}}_{k,l}\left(\mathit{x}\right)={\mathbf{\sigma }}_{ij}\left(\mathit{x}\right)-{\mathbf{C}}_{ijkl}^0{\mathbf{ϵ}}_{i,l}\left(\mathit{x}\right)\\ \hfill {\mathbf{\sigma }}_{ij}\left(\mathit{x}\right)={\mathbf{C}}_{ijkl}^0{\mathit{u}}_{k,l}\left(\mathit{x}\right)+{\mathbf{\phi }}_{ij}\left(\mathit{x}\right)\end{array}$$

(23)

Implementing the equilibrium condition ${\mathbf{\sigma }}_{ij,j}\left(\mathit{x}\right)=0$,

$${\mathbf{C}}_{ijkl}^0{\mathit{u}}_{k,lj}\left(\mathit{x}\right)+{\mathbf{\phi }}_{ij,j}\left(\mathit{x}\right)=0$$

(24)

The auxiliary problem solving the differential Equation (23) for a periodic unit cell under an applied strain $\mathit{E}=<\mathbf{ϵ}\left(\mathit{x}\right)>$ using Green’s function method is given by replacing the polarization term by a unitary body force acting at a given point and along a given direction. Thus, we obtain:

$${\mathbf{C}}_{ijkl}^0{\mathbf{G}}_{km,lj}(\mathit{x}-{\mathit{x}}^{\prime })+{\delta }_{im}(\mathit{x}-{\mathit{x}}^{\prime })=0$$

(25)

where ${\mathbf{G}}_{km}\left(\mathit{x}\right)$ is the Green’s function associated with the displacement field ${\mathit{u}}_k\left(\mathit{x}\right)$. Assuming that the Green’s function is known, the displacement and displacement gradient fields can be written as a convolution integral:

$$\begin{array}{c}\hfill {\mathit{u}}_k\left(\mathit{x}\right)={\int }_{{\mathbf{R}}^3}{\mathbf{G}}_{ki,j}(\mathit{x}-{\mathit{x}}^{\prime }){\mathbf{\phi }}_{ij}\left({\mathit{x}}^{\prime }\right)d{\mathit{x}}^{\prime }\\ \hfill {\mathit{u}}_{k,l}\left(\mathit{x}\right)={\int }_{{\mathbf{R}}^3}{\mathbf{G}}_{ki,jl}(\mathit{x}-{\mathit{x}}^{\prime }){\mathbf{\phi }}_{ij}\left({\mathit{x}}^{\prime }\right)d{\mathit{x}}^{\prime }\end{array}$$

(26)

The convolution becomes simply a multiplication in Fourier space; thus, the strain field can be calculated in Fourier space and then transformed back to real space as follows:

$${\mathbf{ϵ}}_{ij}\left(\mathit{x}\right)={\mathbf{E}}_{ij}+F{T}^{-1}\{sym\left({\widehat{\mathsf{\Gamma }}}_{ijkl}^0\left(\mathbf{k}\right)\right){\widehat{\mathbf{\phi }}}_{kl}\left(\mathbf{k}\right)\}\left(\mathit{x}\right)$$

(27)

where the hat symbol represents the Fourier transform (FT), $\mathbf{k}$ is a frequency in Fourier space, and ${\mathbf{E}}_{ij}$ is the strain prescribed on the unit cell. The Green’s operator in Fourier space, which is only a function of the reference stiffness tensor and the frequency, is given by

$$\begin{array}{cc}\hfill {\widehat{\mathsf{\Gamma }}}_{ijkl}^0\left(\mathbf{k}\right)& =-k_j{k}_l{\widehat{\mathbf{G}}}_{ik}\left(\mathbf{k}\right)\hfill \\ \hfill {\widehat{\mathbf{G}}}_{ij}\left(\mathbf{k}\right)& ={\left. [{\mathbf{C}}_{ijkl}^0{\mathbf{k}}_l{\mathbf{k}}_j\right]}^{-1}\hfill \end{array}$$

(28)

The grid size ${\mathbf{k}}^{\mathbf{d}}$ in the Fourier space is same as the grid size $(N_1\times N_2\times N_3)$ in Cartesian space. The above equations will directly give the displacement and strain fields if the polarization field ${\mathbf{\phi }}_{ij}\left(\mathit{x}\right)$ is known. The model guesses the initial polarization field ${\mathbf{\phi }}_{ij}\left(\mathit{x}\right)$ to calculate the strain field ${\mathbf{ϵ}}_{ij}\left(\mathit{x}\right)$. Using the strain field, the stress fields, and the slip rate fields are updated to re-evaluate the polarization field, and the process is repeated until the polarization field is converged within the given tolerance value.

4.2. Algorithm

Lebensohn et al. [151] used a robust and faster convergence technique called the augmented Lagrangian scheme adapted from [127,154] along with the Newton–Raphson method to solve nonlinear equations. The iterative procedure used by Lebensohn et al. [151] is shown below:

1.

Assume that ${\mathbf{\lambda }}_{ij}^{\left(i\right)}$ and ${\mathit{e}}_{ij}^{\left(i\right)}$ are auxiliary guess stress and strain fields at iteration $\left(i\right)$.

2.

Compute the polarization field at iteration $\left(i\right)$:

$${\mathbf{\varphi }}_{ij}^{\left(i\right)}\left(\mathit{x}\right)={\mathbf{\lambda }}_{ij}^{\left(i\right)}\left(\mathit{x}\right)-{\mathbf{C}}_{ijkl}^0{\mathit{e}}_{kl}^{\left(i\right)}\left(\mathit{x}\right)$$

(29)

3.

A new guess for strain field at iteration (i+1) is then given by:

$${\mathit{e}}_{ij}^{(i+1)}\left(\mathit{x}\right)={\mathbf{E}}_{ij}+F{T}^{-1}\left(sym\left({\widehat{\mathsf{\Gamma }}}_{ijkl}^0\left(\mathbf{k}\right)\right){\widehat{\varphi }}_{kl}^{\left(i\right)}\left(\mathbf{k}\right)\right)$$

(30)

The above equations can be combined to avoid the calculation of the polarization field, as given by Michel et al. [154].

$${\mathit{e}}_{ij}^{(i+1)}\left(\mathit{x}\right)={\mathbf{E}}_{ij}+F{T}^{-1}\left({\widehat{\tilde{\mathit{e}}}}_{ij}^{\left(i\right)}+sym\left({\widehat{\mathsf{\Gamma }}}_{ijkl}^0\left(\mathbf{k}\right)\right){\widehat{\lambda }}_{kl}^{\left(i\right)}\left(\mathbf{k}\right)\right)$$

(31)

4.

The strain field ${\mathit{e}}_{ij}^{(i+1)}$ in above equation is used to calculate the stress ${\mathbf{\sigma }}_{ij}^{(i+1)}$ using augmented Lagrangian scheme. At every material point $\mathit{x}$, a residual $\mathit{R}$, which is a function of ${\mathbf{\sigma }}^{(i+1)}$, is defined by:

$${\mathit{R}}_{ij}\left({\mathbf{\sigma }}^{(i+1)}\right)={\mathbf{\sigma }}_{ij}^{(i+1)}+{\mathbf{C}}_{ijkl}^0{\mathbf{ϵ}}_{kl}^{(i+1)}\left({\mathbf{\sigma }}^{(i+1)}\right)-{\lambda }_{ij}^{\left(i\right)}-{\mathbf{C}}_{ijkl}^0{\mathit{e}}_{kl}^{(i+1)}=0$$

(32)

Here, ${\mathbf{ϵ}}^{(i+1)}$ is a function of ${\mathbf{\sigma }}^{(i+1)}$ as they are related by constitutive relations.

5.

The Newton–Raphson (NR) method is implemented to solve the nonlinear Equation (32), where ${(j+1)}^{th}$ guess for stress field ${\mathbf{\sigma }}_{ij}^{(i+1)}$ is calculated using:

$${\mathbf{\sigma }}_{ij}^{(i+1,j+1)}={\mathbf{\sigma }}_{ij}^{(i+1,j)}-{\left({\left.\frac{\partial {\mathit{R}}_{ij}}{\partial {\mathbf{\sigma }}_{kl}}\right|}_{{\mathbf{\sigma }}^{(i+1,j)}}\right)}^{-1}{\mathit{R}}_{kl}\left({\mathbf{\sigma }}^{(i+1,j)}\right)$$

(33)

The derivative in Equation (33) can be obtained using the expression:

$${\left.\frac{\partial {\mathit{R}}_{ij}}{\partial {\mathbf{\sigma }}_{kl}}\right|}_{{\mathbf{\sigma }}^{(i+1,j)}}\approxeq {\delta }_{ik}{\delta }_{jl}+{\mathbf{C}}_{ijpq}^0{\mathbf{C}}_{pqkl}^{-1}+\left(\mathsf{\Delta }tn{\dot{\gamma }}_0\right){\mathbf{C}}_{ijpq}^0\sum _{s=1}^N\frac{{\mathbf{m}}_{pq}^s{\mathbf{m}}_{kl}^s}{{\tau }_0^s\left({\mathbf{\sigma }}^{(i+1,j)}\right)}{\left(\frac{{\mathbf{m}}_{kl}^s:{\mathbf{\sigma }}_{kl}}{{\tau }_0^s\left({\mathbf{\sigma }}_{kl}^{(i+1,j)}\right)}\right)}^{n-1}$$

(34)

The convergence on ${\mathbf{\sigma }}^{(i+1)}$ (and thus on ${\mathbf{ϵ}}^{(i+1)}$) is achieved with a provided tolerance.

6.

The new guess for the auxiliary stress field ${\lambda }^{(i+1)}$ is obtained by

$${\lambda }^{(i+1)}\left(\mathit{x}\right)={\lambda }^{\left(i\right)}\left(\mathit{x}\right)+{\mathbf{C}}^0:\left({\mathit{e}}^{(i+1)}\left(\mathit{x}\right)-{\mathbf{ϵ}}^{(i+1)}\left(\mathit{x}\right)\right)$$

(35)

The iteration continues until the normalized average differences between the stress fields, $\mathbf{\sigma }\left(\mathit{x}\right)$ and $\lambda \left(\mathit{x}\right)$, and the strain fields, $\mathbf{ϵ}\left(\mathit{x}\right)$ and $\mathit{e}\left(\mathit{x}\right)$, are below the provided tolerance.

5. Twinning Criterion in CP Models

In this section, we will discuss the schemes found in the literature to incorporate deformation twinning into CP models. Most of the crystal plasticity models treat deformation twinning as a pseudo-slip deformation mechanism, which defines the volume fraction of the twin, f, as the ratio of accumulated shear ${\gamma }^{\mathrm{acc}}$ to the characteristic shear of the respective twinning ${\gamma }^{\mathrm{tw}}$ mode, as shown in equation below [99]:

$$f=\frac{{\gamma }^{\mathrm{acc}}}{{\gamma }^{\mathrm{tw}}}$$

(36)

Figure 6 shows the chronological order in which different twinning criteria have been developed in recent years and the twinning behaviors these criteria intend to implement in the crystal plasticity models.Early models, such as PTR and VFT schemes, tend to predict the texture evolution associated with twinning and recent models incorporate various mechanisms of twinning such as detwinning, twin propagation, and twin nucleation. Here, we will discuss each of these models in brief.

5.1. Predominant Twin Reorientation (PTR) Method

The PTR method, proposed by Tomé et al. [6] using Van Houtte [5]’s scheme, is one of the most widely used twinning criterion in CP models. This method tracks the twin volume fraction in each grain, determines the twin system with dominant twin volume fraction using statistical criterion (based on the volume fraction of twinned regions in the grain and entire polycrystal aggregate), and reorients the entire grain into a dominant twin. Shear strain $\mathsf{\Delta }{\gamma }^{n,t_i}$ contributions by each twinning system $t_i$ in the grain n are used to compute the volume fraction of the grain $\mathsf{\Delta }{g}^{n,t_i}$ using Equation (36). $f^n$ represents the weighted volume fraction of the grain in a polycrystal. Thus, the overall twin volume fraction in each grain is calculated by:

$$\mathsf{\Delta }{f}^{n,t_1}=f^n\frac{\mathsf{\Delta }{\gamma }^{n,t_i}}{{\gamma }^{\mathrm{tw}}}=f^n\mathsf{\Delta }{g}^{n,t_i}$$

(37)

After each deformation step, the twin volume fraction in each grain accumulates as:

$$g^{n,t_i}=\sum _{\mathrm{steps}}\mathsf{\Delta }{g}^{n,t_i}$$

(38)

The total twinned volume fraction in a polycrystal is a weighted aggregate of the twin volume fraction in each grain, $g^{n,t_i}$.

$$F_R=\sum _n{f}^n\sum _{t_i}{g}^{n,t_i}$$

(39)

After reorientation of grain into a twinned region, the effective polycrystal twinned fraction is calculated by taking the sum of the grain volume fraction of reoriented grains.

$$F_E=\sum _n{f}^n$$

(40)

The criterion for the reorientation of a grain is measured by the twin fraction accumulated in each twinning system for each grain at every incremental step, $g^{n,t_i}$, compared to the “threshold" fraction $F_T$ defined as:

$$F_T=0.25+0.25\frac{F_E}{F_R}$$

(41)

The main disadvantage of the PTR scheme is that it only accounts for the most active twinning system in each grain. However, in the case of deformation by twinning, more than one twinning system can be active inside the grains. This problem can be mitigated by using a volume fraction transfer scheme, which considers all reorientation from all twin systems [6].

5.2. Volume Fraction Transfer (VFT) Method

This method, proposed by Tomé et al. [6], discretizes the Euler space (space containing all Euler angles representing grain orientations) into identical equiaxed cells of volume $(\mathsf{\Delta }{\varphi }_1\times \mathsf{\Delta }\mathsf{\Phi }\times \mathsf{\Delta }{\varphi }_2)$. Each of these cells represents an orientation defined by its center: ${\varphi }_1^n,{\mathsf{\Phi }}^n,{\varphi }_2^n$ (the Bunge convention) and contains the initial volume fraction $f^n$. After each deformation step, a reorientation of $(\delta {\varphi }_1,\delta \mathsf{\Phi },\delta {\varphi }_2)$ is experienced by each grain. Assuming the discrete cells are small, the cell shifts by the reorientation vector, as shown in Figure 7a, and overlaps with the neighboring cells. The volume fraction, $f^{nm}$, associated with the overlapped region is transferred from the cell m to the neighboring cell n.

The volume fraction of twinning active within the grain is obtained from Equation (36). As a result, the orientation of the twinned region, which is different from the parent grain, is implemented by transferring a finite and non-incremental displacement into the Euler space. As shown in Figure 7b, the fraction $\mathsf{\Delta }{f}^{n,t_1}$ is transferred from the cell m to a non-contiguous cell n in the Euler space.

5.3. Total Lagrangian Approach

The major drawbacks of the PTR and VFT schemes are that the twinned regions are treated as new grains, which can further deform by slip and twinning, just like untwinned regions. However, Asgari et al. [88] showed that twinned regions are harder than the parent matrix. These drawbacks were addressed by Kalidindi et al. [85], where he considers primary and secondary slip systems in addition to primary twin systems in their crystal plasticity models. Kalidindi et al. [85] used the decomposition of the total deformation [157] and incorporated deformation twinning into the plastic deformation with an assumption of the same deformation gradient in the twinned and the untwinned regions, as seen in Figure 5. The incorporation of twinning in the evolution of plastic deformation was done using the constitutive equations shown below:

$${\dot{\mathbf{F}}}^p={\mathbf{L}}^p·{\mathbf{F}}^p$$

(42)

with

$${\mathbf{L}}^p=\left(1-\sum _{\beta }^{N^T}{f}_{\beta }\right)\sum _{\alpha }^{N^S}{\dot{\gamma }}_{\alpha }{\mathbf{S}}_{\alpha }^{o-\mathrm{sl}}+\sum _{\alpha }^{N^T}{\dot{f}}_{\alpha }{\gamma }^{\mathrm{tw}}{\mathbf{S}}_{\alpha }^{o-\mathrm{tw}}+\sum _{\beta }^{N^T}{f}_{\beta }\left(\sum _{\alpha }^{N^{S-T}}{\dot{\gamma }}_{\alpha }{\mathbf{S}}_{\alpha }^{\mathrm{tw}-\mathrm{S}}\right)$$

(43)

The three terms in the right hand side of Equation (42) represents the velocity gradient due to slip in the untwinned region, twins in the untwinned region, and slip in the twinned region, respectively. The evolution of twinning, determined by twin volume fraction, was assumed to depend on the resolved shear stress and a twin system resistance.

$${\dot{f}}_{\alpha }={\dot{f}}_{\alpha }({\tau }_{\alpha },{\mathrm{s}}_{\alpha }^{\mathrm{tw}},\dots ),f_{\alpha }\ge 0,\sum _{\alpha }^{N^T}{f}_{\alpha }\le 1,{\dot{f}}_{\alpha }\ge 0$$

(44)

where ${\mathrm{s}}_{\alpha }^{\mathrm{tw}}$ is associated with the resistance of the twin system. The same approach was exploited by Salem et al. [8] with different resistance evolution for primary slip and twins. The evolution of slip resistance was described by:

$${\dot{\mathrm{s}}}_{\alpha }^{\mathrm{slip}}=h^s\left(1+C\sum _{\beta }{f}_{\beta }\right){\left(1-\frac{{\mathrm{s}}_{\alpha }}{{\mathrm{s}}_{\alpha }^s}\right)}^a\sum _k^{N^s}{\dot{\gamma }}_k$$

(45)

$${\mathrm{s}}_{\alpha }^s={\mathrm{s}}_0^s+{\mathrm{s}}^{pr}{\left(\sum _{\beta }{f}_{\beta }\right)}^{0.5}$$

(46)

$${\dot{\mathrm{s}}}_{\beta }=h^{\mathrm{tw}}{\left(\sum _m{f}^m\right)}^b\sum _k^{}{\gamma }^{\mathrm{tw}}{\dot{f}}_k+h^{\mathrm{tw}-\mathrm{sl}}{\left(\sum _{\alpha }{\gamma }_{\alpha }\right)}^d\sum _k{\dot{\gamma }}_k$$

(47)

The above equations for slip and twin resistances consider some of the important behaviors seen due to deformation twinning, such as hardening of all slip systems, the Hall–Petch mechanism (controlled by the parameter ${\mathrm{s}}^{pr}$), the Basinski–hardening mechanism (controlled by the parameter C), and the saturation of twin volume fraction [8]. Wu et al. [95] further improved this model by incorporating secondary slip and introducing a grain fragmentation technique for accurate simulation of $\alpha$-Titanium. In this model, the criterion for the grain fraction was determined by comparing the twin volume fraction in a grain $\left({\sum }_{\beta }^{}{f}_{\beta }\right)$ with a predetermined saturation value $\left(f_{\mathrm{sat}}\right)$.

5.4. Updated Lagrangian Approach

Lévesque et al. [9], Lévesque et al. [115], Izadbakhsh et al. [158], Izadbakhsh et al. [99] used rate-dependent CP models that incorporated primary slip and twins, secondary slip and twins, as well as tertiary slip for plastic deformation of hcp metals. A detailed discussion on the incorporation of twinning in Lévesque et al. [9] is presented in Section 3.

5.5. Composite Grain (CG) Model

The CG model was proposed by Proust et al. [4] to address the limitation of the PTR scheme in VPSC polycrystal code. The PTR scheme considers either the initial orientation or the twin reoriented grain, which hardly accounts for the directional barrier effect that twins have upon the propagation of dislocation [4]. The CG model takes care of parent and twin reoriented fractions, their interaction, and their evolution with deformation. It adopts the PTR scheme to define the volume fraction of predominant twin system (PTS), the twin system with maximum grain volume fraction.

$$f^{\mathrm{PTS}}=\frac{\mathsf{\Delta }{\gamma }^{\mathrm{PTS}}}{{\gamma }^{\mathrm{tw}}}$$

(48)

Here, PTS is the twin system that forms into a parallel twin lamellae in the grain, as shown in Figure 8. The CG model introduces two parameters: (i) $d^c$ the separation distance of the center planes and (ii) the maximum volume fraction of twin inside the grain $f_{\max }^{\mathrm{PTS}}$. The CG model assumes that the twins with thickness of $d^{\mathrm{tw}}$ form inside a grain of size $d^g$ in an equally spaced manner. The twins can grow up to a maximum value of $d_{\max }^g=f_{\max }^{\mathrm{PTS}}{d}^c$. The evolution of twin and parent thickness is determined by the volume fraction.

$$\begin{array}{cc}\hfill d^{\mathrm{tw}}& =f^{\mathrm{PTS}}{d}^c\hfill \\ \hfill d^{\mathrm{mat}}& =\left(1-f^{\mathrm{PTS}}\right)d^c\hfill \end{array}$$

(49)

The volume fraction of each twin system is tracked until it exceeds the threshold value (typically $5\%$ of grain volume) and defines PTS as the twin system that reach the threshold value first. Within the CG model, the evolution of slip and twin resistances are determined by three mechanisms: the evolution of statistical dislocations with strain, the evolution of geometrically necessary dislocation (GND), and a directional Hall–Petch (HP) effect associated with twin interfaces.

$${\tau }_{\alpha }={\tau }_{\alpha }^{\mathrm{STAT}}+{\tau }_{\alpha }^{\mathrm{GND}}+{\tau }_{\alpha }^{\mathrm{HP}}$$

(50)

The first term in Equation (50) is defined using a classical saturation Voce law for statistical dislocation as well as a latent hardening effect coupled with shear increments $\mathsf{\Delta }{\gamma }_{\beta }$ in the system $\beta$ with the increase in strength in system $\alpha$:

$$\mathsf{\Delta }{\tau }_{\alpha }^{\mathrm{STAT}}=\frac{\partial {\widehat{\tau }}_{\alpha }}{\partial \mathsf{\Gamma }}\sum _{\beta }{h}_{\alpha \beta }\mathsf{\Delta }{\gamma }_{\beta }$$

(51)

where

$${\widehat{\tau }}_{\alpha }\left(\mathsf{\Gamma }\right)={\tau }_0+{\tau }_1\left(1-exp\left(-\frac{\mathsf{\Gamma }{\theta }_0}{{\tau }_1}\right)\right)$$

(52)

Here, $\mathsf{\Gamma }$ is the accumulated shear in the grain. The effect of twin interfaces, which is central to the CG model, comes from a mean free path for each slip and twin system. The second term for Equation (49) is defined by Karaman et al. [159] and Kok et al. [160]:

$$\mathsf{\Delta }{\tau }_{\alpha }^{\mathrm{GND}}=\frac{h_{\alpha }^m}{d_{\alpha }^{\mathrm{mfp}}({\tau }_{\alpha }^{\mathrm{STAT}}+{\tau }_{\alpha }^{\mathrm{GND}})}\mathsf{\Delta }{\gamma }_{\alpha }$$

(53)

$d_{\alpha }^{\mathrm{mfp}}$ represents the mean free path introduced due to specific orientation of twin plane (A) with respect to the slip and twin planes of the other active systems and it is defined as:

$$d_{\alpha }^{\mathrm{mfp}}=\left\{\begin{array}{cc}{d}^{\mathrm{tw}}/sin\left(A\right)\hfill & \mathrm{in}\mathrm{the}\mathrm{twin}\mathrm{and}\hfill \\ d^{\mathrm{mat}}/sin\left(A\right)\hfill & \mathrm{in}\mathrm{the}\mathrm{matrix}\hfill \end{array}\right.$$

(54)

The Hall–Petch–type mechanism is determined using:

$${\tau }_{\alpha }^{\mathrm{HP}}=\frac{h_{\alpha }^{\mathrm{HP}}}{\sqrt{d_{\alpha }^{\mathrm{mfp}}}}$$

(55)

Here, the twin resistance due to the Hall–Petch effect is ${\tau }_{\alpha }^{\mathrm{HP}}$, the Hall–Petch hardening parameter is $h_{\alpha }^{\mathrm{HP}}$, and the mean free path is denoted by $d_{\alpha }^{\mathrm{mfp}}$.

5.6. Twinning Detwinning (TDT) Model

The twinning criteria discussed above mainly describe the reorientation due to twinning and plastic hardening caused by twinning. The CG model employed the empirical treatment of twin nucleation and twin propagation. However, none of these models were able to describe the detwinning process observed by Roberts [161] during a subsequently reversed loading. The TDT model proposed by [11,116] addresses the detwinning behavior of twinning in the crystal plasticity models. A schematic of the TDT model is shown in Figure 9. For a pristine grain shown in Figure 9a, a TDT process has mainly four mechanisms:

(a)

the nucleation of a twin (child) grain with mirrored symmetry with the parent grain (shown in Figure 9b);

(b)

propagation of twin into the parent grain (Figure 9c);

(c)

propagation of the parent grain into the twin/child grain (Figure 9d); and

(d)

activation of twinning dislocation inside twinned region (Figure 9e).

In the TDT model, the resolved shear stress, ${\tau }^{\alpha }=\mathbf{\sigma }:{\mathbf{P}}^{\alpha }$, is assumed to be the driving force for both twinning and detwinning. However, the utmost importance should be given to choosing an appropriate resolved stress. For instance, the twin nucleation process should be driven by the stress state of the parent grain, as no twin grain has formed yet. Likewise, detwinning in Figure 9c is driven by the stress state of the twinned region. However, twinning and detwinning mechanisms in Figure 9c,d are driven by the twin–parent interface. Constitutive equations addressing the TDT model are presented in the following paragraphs.

The shear rates during each TDT mechanism are defined separately.

$$\begin{array}{cc}\hfill {\dot{\gamma }}^{\alpha A}& =\left\{\begin{array}{cc}{\dot{\gamma }}_0{\left|{\tau }^{\alpha }/{\tau }_{cr}^{\alpha }\right|}^{1/m}\hfill & {\tau }^{\alpha }={\tau }^{\alpha P}>0\hfill \\ 0\hfill & {\tau }^{\alpha }={\tau }^{\alpha P}\le 0\hfill \end{array}\right.,\hfill \\ \hfill {\dot{\gamma }}^{\alpha C}& =\left\{\begin{array}{cc}-{\dot{\gamma }}_0{\left|{\tau }^{\alpha }/{\tau }_{cr}^{\alpha }\right|}^{1/m}\hfill & {\tau }^{\alpha }={\tau }^{\alpha P}<0\hfill \\ 0\hfill & {\tau }^{\alpha }={\tau }^{\alpha P}\ge 0\hfill \end{array}\right.\hfill \end{array}$$

(56)

Operation A and C are twinning and detwinning process and cannot be activated simultaneously as operation C is only possible when a twin exists. The evolution of the twin volume fraction for each operation differs as:

$${\dot{f}}^{\alpha A}=\frac{\left|{\dot{\gamma }}^{\alpha A}\right|}{{\gamma }^{\mathrm{tw}}},{\dot{f}}^{\alpha C}=\frac{-\left|{\dot{\gamma }}^{\alpha C}\right|}{{\gamma }^{\mathrm{tw}}}$$

(57)

The shear rates for twinning system $\alpha$ inside a twin are from operations B and D:

$$\begin{array}{cc}\hfill {\dot{\gamma }}^{\alpha B}& =\left\{\begin{array}{cc}-{\dot{\gamma }}_0{\left|{\tau }^{\alpha }/{\tau }_{cr}^{\alpha }\right|}^{1/m}\hfill & {\tau }^{\alpha }={\tau }^{\alpha C}<0\hfill \\ 0\hfill & {\tau }^{\alpha }={\tau }^{\alpha C}\ge 0\hfill \end{array}\right.,\hfill \\ \hfill {\dot{\gamma }}^{\alpha D}& =\left\{\begin{array}{cc}{\dot{\gamma }}_0{\left|{\tau }^{\alpha }/{\tau }_{cr}^{\alpha }\right|}^{1/m}\hfill & {\tau }^{\alpha }={\tau }^{\alpha C}>0\hfill \\ 0\hfill & {\tau }^{\alpha }={\tau }^{\alpha C}\le 0\hfill \end{array}\right.\hfill \end{array}$$

(58)

Again, operations B and D represent twinning and detwinning and cannot act simultaneously. Their activation is possible only in the presence of the twinned region. The evolution of the twin volume fraction is given by equations similar to Equation (56).

$${\dot{f}}^{\alpha B}=\frac{\left|{\dot{\gamma }}^{\alpha B}\right|}{{\gamma }^{\mathrm{tw}}},{\dot{f}}^{\alpha D}=\frac{-\left|{\dot{\gamma }}^{\alpha D}\right|}{{\gamma }^{\mathrm{tw}}}$$

(59)

The overall evolution of twin volume fraction is determined by the sum of weighted volume fraction evolution in the parent $(f^0=1-f^{\mathrm{tw}})$ and twins $(f^{\mathrm{tw}}={\sum }_{\alpha }{f}^{\alpha })$.

$${\dot{f}}^{\alpha }=f^0({\dot{f}}^{\alpha A}+{\dot{f}}^{\alpha C})+f^{\alpha }({\dot{f}}^{\alpha B}+{\dot{f}}^{\alpha D})$$

(60)

A threshold twin volume fraction for the termination of twinning is determined using the accumulated twin fraction $V^{\mathrm{acc}}$ and effective twin fraction $V^{\mathrm{eff}}$.

$$V^{\mathrm{th}}=min\left(1.0,A_1+A_2\frac{V^{\mathrm{eff}}}{V^{\mathrm{acc}}}\right)$$

(61)

where $A_1$ and $A_2$ are two materials constants.

Wu et al. [162] have used this twinning criterion technique in order to study rapid hardening and exhaustion behavior of twinning in Mg alloys.

5.7. Twin Nucleation, Propagation, and Growth (TNPG) Method

Wu et al. [12] proposed a new physics-based constitutive model to describe twin nucleation, propagation, and growth, mainly to overcome the twin propagation mechanisms missing in the methods mentioned above. A schematic of a twinning process and evolution of twin resistances are shown in Figure 10a. The TNPG method assumes that the twin nucleates at a grain boundary under certain conditions and propagates on the twinning plane along the twin shearing direction, forming a twin band. Figure 10b represents the evolution of twin resistances ${\tau }^{\alpha }$ as a function of twin volume fraction $f^{\alpha }$.

In the constitutive model for Elastic Viscoplastic Self–Consistent (EVPSC) polycrystal model by Wang et al. [163], the evolution of the CRSS ${\tau }_{cr}^{\alpha }$ is given by

$${\tau }_{cr}^{\alpha }=\frac{d{\tau }^{\alpha }}{d\mathsf{\Gamma }}\sum _{\beta }^{}{h}^{\alpha \beta }\left|{\dot{\gamma }}^{\beta }\right|$$

(62)

where the shear rate due to twinning ${\dot{\gamma }}^{\beta }$ is defined by a piecewise function due to its polar nature.

$${\dot{\gamma }}^{\alpha }=\left\{\begin{array}{cc}{\dot{\gamma }}_0{|{\tau }^{\alpha }/{\tau }_{cr}^{\alpha }|}^{\frac{1}{m}}\hfill & {\tau }^{\alpha }>0\hfill \\ 0\hfill & {\tau }^{\alpha }\le 0\hfill \end{array}\right.$$

(63)

Finally, the threshold stress that characterizes twinning with the incorporation of stress relaxation associated with twin propagation is defined using the following equation.

$${\tau }^{\alpha }=\left\{\begin{array}{cc}{\tau }_0^{\alpha }-\frac{{\tau }_0^{\alpha }-{\tau }_g^{\alpha }}{f_g^{\alpha }}{f}^{\alpha }\hfill & f^{\alpha }\le f_g^{\alpha }\hfill \\ {\tau }_g^{\alpha }+\left({\tau }_1^{\alpha }+h_1^{\alpha }{\mathsf{\Gamma }}^{\alpha }\right)\left(1-exp\left(-\frac{h_0^{\alpha }}{{\tau }_1^{\alpha }}{\mathsf{\Gamma }}^{\alpha }\right)\right)\hfill & f^{\alpha }>f_g^{\alpha }\hfill \end{array}\right.$$

(64)

Here, the accumulated shear of the twinning ${\mathsf{\Gamma }}^{\alpha }$ is calculated after $f^{\alpha }>f_g^{\alpha }$. Like PTR and TDT models, a threshold twin volume fraction to terminate twinning is also defined in the TNPG model based on two statistical variables: the accumulated twin fraction $V^{\mathrm{acc}}$ and the effective twinned fraction $V^{\mathrm{eff}}$, as shown in Equation (60).

5.8. Dislocation Density Based Model

In this section, we will discuss the constitutive model based on dislocation densities for slip and twinning by Beyerlein and Tomé [10] in order to address the effect of deformation histories involving changes in temperature and strain rate. Each slip mode has its own dislocation evolution law that includes a thermally-activated recovery process that leads to either annihilation or debris formation. The hardening law for twin propagation accounts for temperature effects through its interaction with slip dislocations. The dislocation density-based model uses the CG approach for twin grain formation with the PTS scheme with dislocation density-based slip and twin resistances.

The slip resistance is defined as a summation of friction stress ${\tau }_0^{\alpha }$ (depends on Peierls stress and initial dislocation content), a barrier term ${\tau }_{0,\mathrm{HP}}^{\alpha }$ (depends on initial grain size), and forest and debris interaction stresses that evolve with strain ${\tau }_{\mathrm{for}}^s,{\tau }_{\mathrm{for}}^{\alpha },and{\tau }_{\mathrm{deb}}^{\alpha }$ (depends on a spatially random and ordered distribution of stored dislocations, respectively) [10,164]:

$${\tau }^s={\tau }_0^{\alpha }+{\tau }_{0,\mathrm{HP}}^{\alpha }+{\tau }_{\mathrm{for}}^s+{\tau }_{\mathrm{deb}}^{\alpha }$$

(65)

The friction stresses, ${\tau }_0^{\alpha }$, are the function of temperature and strain rate and defined as:

$${\tau }_0^{\alpha }(\dot{\mathbf{ϵ}},T)=A^{\alpha }(1+B^{\alpha }log\left(\dot{\mathbf{ϵ}}\right))exp\left(-\frac{T}{C^{\alpha }}\right)\left(1+D^{\alpha }Lognormal({\mathbf{\sigma }}^{\alpha },{\nu }^{\alpha })\right)$$

(66)

where $A^{\alpha },B^{\alpha },\mathrm{and}{C}^{\alpha }$ are constants and $D^{\alpha }$ and lognormal terms, ${\mathbf{\sigma }}^{\alpha }$ and ${\nu }^{\alpha }$ are used to capture the saturation of the resistances at high temperatures. T and $\dot{\mathbf{ϵ}}$ are the current temperature and strain rate, respectively. The barrier term, ${\tau }_{0,\mathrm{HP}}^{\alpha }$, is given by:

$${\tau }_{0,\mathrm{HP}}^{\alpha }=\left\{\begin{array}{cc}{\mu }^{\alpha }\left(T\right)H{P}^s\sqrt{\frac{b^{\alpha }}{d_g}},\hfill & \mathrm{without}\mathrm{twins}\hfill \\ \\ \frac{f^{\mathrm{PTS}}-f_0^{\mathrm{PTS}}}{f_{\max }^{\mathrm{PTS}}-f_0^{\mathrm{PTS}}}{\mu }^{\alpha }\left(T\right)H{P}^s\sqrt{\frac{b^{\alpha }}{d_{\mathrm{mfp}}^s}},\hfill & \mathrm{with}\mathrm{twins}\mathrm{present}\hfill \end{array}\right.$$

(67)

where $b_g$ is the magnitude of burgers vector, $d_g$ is the initial grain size, ${\mu }^{\alpha }\left(T\right)$ is the effective shear modulus (same for all directions), and $H{P}^{\alpha }$ is the Hall–Petch parameter. $f_0^{\mathrm{PTS}}$ is the minimum volume fraction identifying the PTS.

The later two terms in Equation (64) evolve with the strain governed by the forest dislocation density ${\rho }_{\mathrm{for}}^s(\dot{\mathbf{ϵ}},T)$ and debris dislocation density ${\rho }_{\mathrm{deb}}(\dot{\mathbf{ϵ}},T)$, as shown by:

$$\begin{array}{cc}\hfill {\tau }_{\mathrm{for}}^s(\dot{\mathbf{ϵ}},T)& =b^{\alpha }{\mu }^{\alpha }\left(T\right)\sqrt{{\chi }^{s{s}^{\prime }}{\rho }_{\mathrm{for}}^{s^{\prime }}(\dot{\mathbf{ϵ}},T)},\hfill \\ \\ \hfill {\tau }_{\mathrm{deb}}(\dot{\mathbf{ϵ}},T)& =k_{\mathrm{deb}}{\mu }^{\alpha }\left(T\right)b^{\alpha }\sqrt{{\rho }_{\mathrm{deb}}(\dot{\mathbf{ϵ}},T)}\times log\left(\frac{1}{b^{\alpha }\sqrt{{\rho }_{\mathrm{deb}}(\dot{\mathbf{ϵ}},T)}}\right)\hfill \end{array}$$

(68)

where ${\chi }^{s{s}^{\prime }}$ is a dislocation interaction matrix accurately characterized through specialized experiment and simulation [165]. The evolution of forest density, ${\rho }_{\mathrm{for}}^s(\dot{\mathbf{ϵ}},T)$, is given by the following relations [10,166,167]:

$$\begin{array}{cc}\hfill & \frac{\partial {\rho }_{\mathrm{for}}^s}{\partial {\gamma }^{s^{\prime }}}=\frac{\partial {\rho }_{\mathrm{gen},\mathrm{for}}^s}{\partial {\gamma }^{s^{\prime }}}-\frac{\partial {\rho }_{\mathrm{rem},\mathrm{for}}^s}{\partial {\gamma }^{s^{\prime }}}=k_1^s\sqrt{{\rho }_{\mathrm{for}}^s}-k_2^s(\dot{\mathbf{ϵ}},T){\rho }_{\mathrm{for}}^s,\hfill \\ \\ \hfill & \frac{k_1^s}{k_2^s}=\frac{1}{\sqrt{{\left({\chi }^{-1}\right)}^{s^{\prime }s}{\left(\frac{{\tau }_{\mathrm{sat}}^s}{b^{\alpha }{\mu }^{\alpha }}\right)}^2}},{\tau }_{\mathrm{sat}}^s=\frac{D^{\alpha }{\left(b^{\alpha }\right)}^3{g}^{\alpha }{\mu }^{\alpha }}{D^{\alpha }{\left(b^{\alpha }\right)}^3-kTln\left(\frac{\dot{\mathbf{ϵ}}}{{\dot{\mathbf{ϵ}}}_0}\right)}.\hfill \end{array}$$

(69)

where $k_1^s$ is the rate-insensitive coefficient that accounts for dislocation storage by statistical trapping of mobile dislocations, and $k_2^s$ is the rate-sensitive coefficient for dynamic recovery through thermally activated mechanisms, respectively. $k,{\dot{\mathbf{ϵ}}}_0,g^{\alpha },\mathrm{and}{D}^{\alpha }$ are the Boltzmann constant, a reference strain rate, an effective activation enthalpy, and a drag stress, respectively. The evolution of debris dislocation density is obtained by:

$${d\rho }_{\mathrm{deb}}=\sum _s{q}^{\alpha }{b}^{\alpha }\frac{\partial {\rho }_{\mathrm{rem},\mathrm{for}}^s}{\partial {\gamma }^{s^{\prime }}}\left|d{\gamma }^{s^{\prime }}\right|$$

(70)

where $q^{\alpha }$ defines a fraction of dislocations leading to debris formation and the remaining fraction that annihilates.

Similarly, resistance for twinning is defined by three different terms: (i) a temperature-independent friction term ${\tau }_0^{\beta }$, (ii) a Hall–Petch–like term ${\tau }_{0,\mathrm{HP}}^{\mathrm{tw}}$, and (iii) a latent hardening term coupling slip and twin systems.

$${\tau }_c^{\mathrm{tw}}={\tau }_0^{\beta }+{\tau }_{0,\mathrm{HP}}^{\mathrm{tw}}+{\tau }_{\mathrm{slip}}^{\beta }$$

(71)

where tw and $\beta$ in superscripts refer to the twin variant and twinning mode, respectively. These terms are defined by:

$$\begin{array}{cc}\hfill & {\tau }_0^{\beta }={\tau }_{\mathrm{prop}}^{\beta }+({\tau }_{\mathrm{crit}}^{\beta }-{\tau }_{\mathrm{prop}}^{\beta })exp\left(-\frac{\sum {\rho }_{\mathrm{for}}^s}{n^{\alpha }{\rho }_{\mathrm{sat}}^s}\right)\hfill \\ \hfill & {\tau }_{0,\mathrm{HP}}^{\mathrm{tw}}=\frac{H{P}^{\beta }}{\sqrt{d_g}}\hfill \\ \hfill & {\tau }_{\mathrm{slip}}^{\beta }={\mu }^{\beta }\left(T\right)\sum _s{C}^{\beta \alpha }\left(\dot{\mathbf{ϵ}}\right)b^{\beta }{b}^{\alpha }{\rho }_{\mathrm{for}}^s\hfill \end{array}$$

(72)

where ${\tau }_{crit}^{\beta }$ and ${\tau }_{\mathrm{prop}}^{\beta }$ are the stresses required for nucleation and propagation of twin $\beta$, whose contribution to twin activation is leveraged by the probability term $\frac{\sum {\rho }_{\mathrm{for}}^s}{n^{\alpha }{\rho }_{\mathrm{sat}}^s}$ [10]. ${\mu }^{\beta }$ and $C^{\alpha \beta }$ are the elastic shear modulus on the system and latent hardening matrix as a function of strain rate.

This model assumes twin resistance as a shear stress $\tau$ resolved in the twin direction and twin plane. This assumption applies best for the twin propagation and not for nucleation [10]. The CP model for dislocation-driven deformation of hcp metals producing equivalent results from the dislocation density simulation was presented by Messner et al. [168].

5.9. Probabilistic Nucleation Method

Based on the numerical observation from atomistic simulations by Wang et al. [34] and experimental statistical analyses by Beyerlein et al. [33] and Capolungo et al. [169], a probabilistic nucleation framework was presented by Beyerlein and Tomé [32] that is able to predict the dispersion in the onset of twinning and variant selection with respect to crystallographic orientation. The resolved shear stress (RSS) $\tau$ in the direction of twinning dislocation is an important quantity to describe twin nucleation events. In hcp crystals, the RSS ${\tau }^{\alpha }$ on twin variant $\alpha$ due to average stress in the grain ${\mathbf{\sigma }}_{ij}^g$ is given by

$${\tau }^{\alpha }={\mathbf{P}}_{ij}^{\alpha }{\mathbf{\sigma }}_{ij}^g,$$

(73)

where ${\mathbf{P}}_{ij}^{\alpha }$ is defined by Equation (6), and the RSS ${\tau }^{\alpha }$ for the grain is calculated using the average properties of the grains. In this nucleation model, the twins are assumed to nucleate from the grain boundaries. Thus, the inhomogeneities in RSS near the grain boundaries and/or interior of the grain are defined by $\mathsf{\Delta }{\tau }_{ij}\left(x\right)$. The RSS at any point inside the grain, ${\tau }^{\alpha }\left(x\right)$, can be defined using:

$${\tau }^{\alpha }\left(x\right)={\tau }^{\alpha }\pm \mathsf{\Delta }{\tau }^{\alpha }\left(x\right)={\mathbf{P}}_{ij}^{\alpha }\left({\mathbf{\sigma }}_{ij}+\mathsf{\Delta }{\mathbf{\sigma }}_{ij}\left(x\right)\right),$$

(74)

The effects of temperature, strain rate, deformation mode, deformation histories, boundary misorientation, alloying contents, and other sources of heterogeneity are accounted by ${\tau }^{\alpha }\left(x\right)$.

The next step is defining the grain boundary that is most favorable to nucleate a twin variant. The selection of a grain boundary is quantified by its areal fraction, $a^{*}$, which is assigned to every grain (see Figure 11a). One of the possible choices for $a^{*}$, which is compatible with mean-field polycrystal approaches, is the area of a spherical cap $\left(a_{\mathrm{cap}}\right)$ defined by the area of the intersection between a cone and the spherical grain, as shown in Figure 11b. This area can be computed using the equation

$$a^{*}=a_{\mathrm{cap}}=\frac{\pi }{2}{d}_g^2\left(1-cos\left(\theta \right)\right),$$

(75)

where $\theta$ is the angle of the cone. The value of $\theta$ varies between 30 and 80${}^{\circ }$ [31,33].

A twin nucleation criterion using the probabilistic nucleation method depends on critical nucleation stresses, $S_{\min }\mathrm{and}{S}_{\max }$, as defined below:

$$n_{\mathrm{twins}}=\left\{\begin{array}{cc}0,\mathrm{no}\mathrm{nucleation}\hfill & \tau <S_{\min }+\frac{\mathsf{\Delta }S}{3}\hfill \\ \\ n^{*}\mathrm{twins}\hfill & \tau >S_{\min }+\frac{\mathsf{\Delta }S}{3}\hfill \end{array}\right.$$

(76)

where $\mathsf{\Delta }S=S_{\max }-S_{\min }$. The critical nucleation stresses are related to defect structures within the grain boundary $a^{*}$ and calculated by using equations below:

$$\begin{array}{c}\hfill S_{\min }={\tau }_c{\left(-\frac{ln(1-Y)}{n^{*}}\right)}^{\frac{1}{\mathbf{\phi }}}\\ \hfill S_{\max }={\tau }_c{\left(-ln(1-Y^{\frac{1}{n^{*}}})\right)}^{\frac{1}{\mathbf{\phi }}}\end{array}$$

(77)

where ${\tau }_c={\tau }_0{\left(\frac{a_0}{a_c}\right)}^{\frac{1}{\mathbf{\phi }}}$, Y is a randomly picked number from uniform distribution $U(0,1)$, and ${\tau }_0$ is a characteristic stress (scalar) that corresponds to a grain boundary surface area $a_0$. $\mathbf{\phi }$ is a material property. Here, the number of twins with characteristic area of $a_c\approx 0.5d^g$ is given by:

$$n^{*}=\frac{a^{*}}{a_c}$$

(78)

Once nucleated, the propagation of $n_{\mathrm{twins}}$ twins is proportional to the shear rate ${\dot{\gamma }}^{\alpha }$ and described by the power law:

$${\dot{\gamma }}^{\alpha }={\dot{\gamma }}_0{\left(\frac{{\tau }^{\alpha }}{{\tau }_{\mathrm{prop}}}\right)}^n,$$

(79)

where ${\tau }_{\mathrm{prop}}$ is a characteristic propagation stress.

5.10. Explicit Incorporation of Twin

An explicit incorporation of twinning in the CPFE model was formulated by Ardeljan et al. [111], where a dislocation-based hardening model was used. In order to incorporate a twin explicitly in a microstructure, a certain procedure was followed that is represented in Figure 12.

Based on the distribution of grain size and grain orientations, a voxel-based microstructure and surface meshes for each individual grain is generated using DREAM.3D software. A voxel density/resolution in the model is defined to achieve the desired number of surface finite elements. Starting from the surface mesh, which is the bridge between a voxel-based model and a volumetric mesh, a 3D solid meshing of individual grains is performed to ensure mesh conformance between grain boundaries, as seen in Figure 12b. This “conformal” grain boundaries is critical to modeling grain-to-grain interactions and is not found in other spatially resolved full-field techniques [111].

After each deformation step, the procedure starts by obtaining the surface mesh of deformed grains (Figure 12c) regardless of whether a twin already exists, which twin variants are first formed, or a pre-existing twin variant thickens. The extraction of surface mesh, finding the intersecting points between two cutting planes (that form twin), and the 3D mesh of the parent grain are performed by a Python script, as seen in Figure 12d,e. The script also contains the commands to generate and export surface meshes for each intersection plane. The surface meshes of the neighboring grains are also adjusted to preserve the mesh conformity of the grain boundaries of neighboring grains touching the twins (see Figure 12f). At this point, the model that is made of each grain is assembled, and a surface mesh cleanup is performed. Then, an Abaqus input file containing the new mesh, state variables, and boundary conditions for the subsequent step is generated. After meshing, the state variables ${\mathbf{F}}^P$, dislocation densities, and crystal orientation are mapped onto the new mesh.

The criterion used to insert a twin lamella is based on the threshold value of $1\%$ for the twin volume fraction $f^{\mathrm{PTS}}$. In each deformation step, a twin volume fraction for each twin variant is calculated. When the twin volume fraction of a PTS reaches the threshold value, the deformation is interrupted to insert the twin lamella. The twins are nucleated at the grain boundaries (places with high von Mises stresses) and placed across the parent grain with twin thickness such that the total volume fraction is matched. The plastic portion of the deformation gradient in a freshly formed twin is given by:

$${\mathbf{F}}^{P^{-1}}\left(\tau \right)=\left\{\mathbf{I}-{\mathbf{S}}^{tw,PTS}\mathit{m}{ß}_0^{tw,PTS}\right\}$$

(80)

The plastic deformation in the parent grain is updated to account for the strain accommodation due to twin formation and growth as:

$${\mathbf{F}}^{P^{-1}}\left(\tau \right)={\mathbf{F}}^{P^{-1}}\left(t\right)\left\{\mathbf{I}+f^{tw,PTS}{\mathbf{S}}^{tw,PTS}{\mathit{m}}_0^{tw,PTS}\right\}$$

(81)

Simultaneously, Kumar et al. [139,140] used FFT-based explicit twin incorporation of a twin into a single crystal and polycrystal materials to study the localized twinning behavior.

5.11. Energy-Based Micro-Twin Nucleation Model

Based on energy-partitioning following the dislocation dissociation process, micro-twin nucleation criteria were developed by Cheng and Ghosh [109,110]. For a heterogeneous type of nucleation, a dislocation-assisted mechanism for $\left\{10\overline{1}2\right\}$ twins involves non-planar dissociation of a sessile 〈c + a〉 lattice dislocation into a residual sessile stair-rod partial dislocation and n layers of glissile twinning dislocations. The occurrence of such a dissociation process, shown in the equation below, at a sufficiently large applied shear stress on a 2nd order pyramidal 〈c + a〉 dislocation system $\left\{\overline{1}2\overline{1}3\right\}$ $\langle 1\overline{2}12\rangle$ is presented by Ghazisaeidi and Curtin [51].

$${\mathit{b}}_{\mathrm{ini}}\to {\mathit{b}}_{\mathrm{r}}+{\mathit{b}}_{\mathrm{tw}}$$

(82)

where ${\mathit{b}}_{\mathrm{ini}}$ is the Burgers vector of the initial 〈c + a〉 dislocation, ${\mathit{b}}_{\mathrm{r}}$ is a Burgers vector for the residual stair-rod dislocation, and ${\mathit{b}}_{\mathrm{tw}}$ is a Burgers vector for twin partial dislocations. For 〈c + a〉 2nd order pyramidal dislocation with ${\mathit{b}}_{\mathrm{ini}}=\frac{1}{3}$ $\langle \overline{1}2\overline{1}3\rangle$, the twin partial dislocation is $\langle \overline{1}0\overline{1}1\rangle$, where n is the number of twin layers, and s is the magnitude of shear on each layer [109]. Figure 13 shows $\left\{10\overline{1}2\right\}$ twin nucleation where dislocation on $(11\overline{2}\overline{2})$ plane with the line direction along $\left[4\overline{2}\overline{2}3\right]$ dissociates onto the $\left(1\overline{1}02\right)$ plane and a schematic of 〈c + a〉 dislocation dissociation on the $\left\{10\overline{1}2\right\}$ twin plane.

The initial energy before the dissociation of dislocation is given by:

$$E_{\mathrm{ini}}=\frac{L}{4\pi }\left. [K_{\mathrm{ini}}^e{\left(b_{\mathrm{ini}}^e\right)}^2+K_{\mathrm{ini}}^s{\left(b_{\mathrm{ini}}^s\right)}^2\right]ln\frac{R}{r_0}$$

(83)

where edge and screw dislocations are represented by superscript e and s, and $K_{\mathrm{ini}}$ is the elastic energy coefficient for dislocation, for which the values are calculated using integral methods [170] and tabulated in [109]. R and $r_0$ are outer and inner radii of the dislocation core, respectively. The energy in the system after dissociation is given as:

$$E_F=E_{\mathrm{tw}}+E_r+E_{\mathrm{inter}}+E_{\mathrm{fault}}-W_{\mathrm{ex}}$$

(84)

where energies on the right side are associated with self–energy of the twinning dislocation loop $E_{\mathrm{tw}}$, self–energy of the stair-rod dislocation $E_r$, interaction energy between the twinning dislocation and stair-rod dislocation $E_{\mathrm{inter}}$, stacking fault energy from twinning $E_{\mathrm{fault}}$, and $W_{\mathrm{ex}}$ represents the work done by external stress. The self–energy of twinning dislocation $E_{\mathrm{tw}}$ is evaluated as the sum of the line energy of the front segment, the transverse segment, and the interaction energy between the transverse segments as:

$$\begin{array}{cc}\hfill E_{\mathrm{tw}}=& \frac{1}{4\pi }\left(L\left. [K_{\mathrm{tw}}^{s,ft}{\left(b_{\mathrm{tw}}^{s,ft}\right)}^2+K_{\mathrm{tw}}^{e,ft}{\left(b_{\mathrm{tw}}^{e,ft}\right)}^2\right]+2d\left. [K_{\mathrm{tw}}^{s,tr}{\left(b_{\mathrm{tw}}^{s,tr}\right)}^2+K_{\mathrm{tw}}^{e,tr}{\left(b_{\mathrm{tw}}^{e,tr}\right)}^2\right]\right)ln\frac{R}{r_0}\hfill \\ \hfill & -d\left(K_{\mathrm{tw}}^{e,tr}\left. [\frac{({\mathit{b}}_{\mathrm{tw}}^{e,tr}\otimes {\mathsf{\xi }}^{tr}).({\mathit{b}}_{\mathrm{tw}}^{e,tr}\otimes {\mathsf{\xi }}^{tr})}{2\pi }ln\frac{L}{r_0}+\frac{[({\mathit{b}}_{\mathrm{tw}}^{e,tr}\otimes {\mathsf{\xi }}^{tr}).\mathbf{L}][({\mathit{b}}_{\mathrm{tw}}^{e,tr}\otimes {\mathsf{\xi }}^{tr}).\mathbf{L}]}{2\pi L^2}\right]\right)\hfill \\ \hfill & -d\left(K_{\mathrm{tw}}^{s,tr}\frac{{\left(b_{\mathrm{tw}}^{s,tr}\right)}^2}{2\pi }ln\frac{L}{r_0}\right)\hfill \end{array}$$

(85)

The superscripts $fr$ and $tr$ denote the front and transverse parts of the twin partial dislocation loop. L is the length of the front segment, which is equal to the length of the initial 〈c + a〉 dislocation $\mathbf{L}$. The transverse segment of twinning dislocation has a length of the dissociation distance d, and ${\mathsf{\xi }}^{tr}$ represents the unit vector of the transverse portion of twin partial dislocation.

The self–energy for stair-rod dislocation is given as:

$$E_r=\frac{L}{4\pi }\left. [K_r^e{\left(b_r^e\right)}^2+K_r^s{\left(b_r^s\right)}^2\right]ln\frac{R}{r_0}$$

(86)

The interaction energy $E_{\mathrm{inter}}$ is given as:

$$\begin{array}{cc}\hfill E_{\mathrm{inter}}=& -L{K}_{\mathrm{tw}}^{s,ft}\frac{({\mathit{b}}_{\mathrm{tw}}^{s,ft}.{\xi }^{ft})({\mathit{b}}_r^s.{\xi }^r)}{2\pi }ln\frac{d}{r_0}\hfill \\ \hfill & -L{K}_{\mathrm{tw}}^{e,ft}\left(\frac{({\mathit{b}}_{\mathrm{tw}}^{e,ft}\otimes {\xi }^{ft})({\mathit{b}}_r^e\otimes {\xi }^r)}{2\pi }ln\frac{d}{r_0}+\frac{[({\mathit{b}}_{\mathrm{tw}}^{e,ft}\otimes {\xi }^{ft}).\mathit{d}][({\mathit{b}}_{\mathbf{r}}^{\mathit{e}}\otimes {\xi }^{\mathbf{r}}).\mathit{d}]}{2\pi d^2}\right)\hfill \end{array}$$

(87)

Similarly, the stacking fault energy and work done due to external shear stress ${\tau }_{\mathrm{tw}}$ are given by:

$$\begin{array}{cc}\hfill E_{\mathrm{fault}}& ={\nu }_{\mathrm{tw}}dL\hfill \\ \hfill W_{\mathrm{ex}}& ={\tau }_{\mathrm{tw}}{b}_{\mathrm{tw}}dL\hfill \end{array}$$

(88)

where ${\nu }_{\mathrm{tw}}$ is the twin boundary energy.

The stable separation distance $d_s$ after the dissociation of dislocation is calculated by minimizing the energy, $E_F$, such that:

$$\frac{\partial E_F}{{\partial }_d}=0,\frac{{\partial }^2{E}_F}{\partial d^2}\ge 0,d=d_s$$

(89)

The minimum stable separation distance $d_s$ is assumed to be $d_s>2r_0$. The criterion for twin nucleation is determined by comparing initial self–energy of 〈c + a〉 dislocation to the sum of self–energies for twinning dislocations and stair-rod dislocation for each twin variants:

$$\begin{array}{c}\hfill E_{\mathrm{ini}}\ge E_{\mathrm{tw}}(d=0)+E_r\mathrm{and}\\ \hfill E_{\mathrm{ini}}>E_F(d_s,{\tau }_{\mathrm{tw}})\forall d_s>2r_0\end{array}$$

(90)

This work is extended to include the twin propagation model by Cheng and Ghosh [110], which incorporates both shear and shuffle processes.

5.12. Thermal Activation Based Propagation Model

Cheng and Ghosh [110] expanded their energy-based nucleation model to incorporate the thermal activation-based propagation model. The velocity of a twin partial dislocation gliding on a plane, ${\nu }_{\mathrm{glide}}$, is provided by Keshavarz and Ghosh [171] as a function of the thermal activation energy barrier, which also accounts for the shear-shuffle process.

$${\nu }_{\mathrm{glide}}=f_{\mathrm{shuffle}}{\lambda }_{\mathrm{shear}}\left. [exp\left(-\frac{\mathsf{\Delta }F-\tau A_p{b}_{\mathrm{tw}}}{K_{B}T}\right)\right]$$

(91)

where the frequency of shuffle is $f_{\mathrm{shuffle}}$, the shear distance is ${\lambda }_{\mathrm{shear}}$, and the probability of gliding in the presence of the energy barrier $\mathsf{\Delta }F$ is $exp\left(-\frac{\mathsf{\Delta }F-\tau A_p{b}_{\mathrm{tw}}}{K_{B}T}\right)$. The effective resolved shear stress on the twin plane is $\tau$, and the shearing area during the plastic deformation is given by $A_p$.

Further, the twin growth along the twin plane perpendicular is determined as a function of the growth velocity, which is obtained from the stimulated slip model by Yu et al. [172]. The twin growth, the transfer of twin-partial dislocations from current twin plane to the adjacent twin plane, is promoted by the dislocation dipoles. The rate of encountering a promoter is a function of the total dislocation density $\left({\rho }_{\mathrm{tot}}\right)$, the fraction of dislocations that penetrates the twin plane $\left(P_{\mathrm{promoter}}\right)$, the glide velocity $\left({\nu }_{\mathrm{glide}}\right)$, and the length of twin $\left(l_{\mathrm{tw}}\right)$.

$$R=P_{\mathrm{promoter}}{\rho }_{\mathrm{tot}}{\nu }_{\mathrm{glide}}{l}_{\mathrm{tw}}$$

(92)

The twin growth is determined by the velocity of twin partial dislocations crossing twin planes, ${\nu }_{\mathrm{grow}}$, which is given as the distance between two planes divided by the time interval:

$${\nu }_{\mathrm{grow}}=\frac{d_{\mathrm{tw}}}{\mathsf{\Delta }{t}_{\mathrm{tw}}}=d_{\mathrm{tw}}\times R=d_{\mathrm{tw}}{P}_{\mathrm{promoter}}{\rho }_{\mathrm{tot}}{l}_{\mathrm{tw}}{\nu }_{\mathrm{glide}}$$

(93)

Here, the time interval $\mathsf{\Delta }{t}_{\mathrm{tw}}$ is inversely proportional to the rate R. The time-averaged plastic shear rate, ${\dot{\gamma }}_{\mathrm{tw}}$, due to micro-twinning is derived from the Orowan equation as a function of twin partial dislocation ${\rho }_{\mathrm{tw}}$ and its Burgers vector $b_{\mathrm{tw}}$:

$${\dot{\gamma }}_{\mathrm{tw}}={\rho }_{\mathrm{tw}}{b}_{\mathrm{tw}}{\nu }_{\mathrm{glide}}={\dot{\gamma }}_{0,\mathrm{tw}}exp\left(-\frac{\mathsf{\Delta }F-\tau A_P{b}_{\mathrm{tw}}}{K_{B}T}\right)$$

(94)

where ${\dot{\gamma }}_{0,\mathrm{tw}}={\rho }_{\mathrm{tw}}{b}_{\mathrm{tw}}{f}_{\mathrm{shuffle}}{\gamma }_{\mathrm{shear}}$. Further manipulation of Equation (93) and using the fact that there is no thermal activation at 0 $\mathrm{K}$ temperature, the shear rate for twin can be written as:

$${\dot{\gamma }}_{\mathrm{tw}}={\dot{\gamma }}_{0,\mathrm{tw}}{\left|\frac{\tau \left(T\right)}{\tau (T=0)}\right|}^{\frac{\mathsf{\Delta }F}{K_{B}T}}sign\left(\tau \left(T\right)\right)={\dot{\gamma }}_{0,\mathrm{tw}}\left|\frac{\tau \left(T\right)}{{\tau }_{\alpha }^{\mathrm{tw}}}\right|sign\left(\tau \right)$$

(95)

The rate of shear resistance for the twin system $\alpha$ due to interactions with dislocation slips is given by the phenomenological hardening law:

$${\dot{\tau }}_{\alpha }^{\mathrm{tw}}=\sum _{\beta =1}^{N^s}{h}_{\alpha \beta }\left|{\dot{\gamma }}_{\beta }\right|$$

(96)

The hardening of slip systems due to twinning is expressed as:

$${\dot{\tau }}_{\alpha }^{\mathrm{sl}-\mathrm{tw}}=\sum _{\beta =1}^{N^{\mathrm{tw}}}{h}_{\alpha \beta }\left|{\dot{\gamma }}_{\beta }^{\mathrm{tw}}\right|$$

(97)

The numerical implementation of the twin propagation model is based on the criterion determined by glide and growth velocity.

$$\begin{array}{cc}\hfill {\nu }_{\mathrm{glide}}& \ge \frac{l_{\mathrm{glide}}}{\mathsf{\Delta }{t}_{\mathrm{twin}}}\hfill \\ \hfill {\tau }^{crit}& \ge \frac{ln\left(\frac{l_{growth}}{\mathsf{\Delta }{t}_{\mathrm{twin}}{f}_{\mathrm{shuffle}}{\lambda }_{\mathrm{shear}}}\right)K_{B}T+\mathsf{\Delta }F}{A_p{b}^{\mathrm{tw}}}\hfill \end{array}$$

(98)

Likewise, the twin growth criterion is a function of the growth velocity for which the critical stress requirement is given by:

$$\begin{array}{cc}\hfill {\nu }_{\mathrm{grow}}& \ge \frac{l_{growth}}{\mathsf{\Delta }{t}_{\mathrm{twin}}}\hfill \\ \hfill {\tau }^{crit}& \ge \frac{ln\left(\frac{l_{\mathrm{glide}}}{\mathsf{\Delta }{t}_{\mathrm{twin}}{f}_{\mathrm{shuffle}}{\lambda }_{\mathrm{shear}}{d}_{\mathrm{tw}}{P}_{\mathrm{promoter}}{\rho }_{\mathrm{tot}}{l}_{\mathrm{tw}}}\right)K_{B}T+\mathsf{\Delta }F}{A_p{b}_{\mathrm{tw}}}\hfill \end{array}$$

(99)

Equations (97) and (98) above provide the stress criterion to propagate a twin to a neighboring material point in longitudinal and transverse directions of a twin plane, respectively. Ghosh and Cheng [173] used subcycling–augmented CPFEM, which uses an adaptive multi–time subcycling algorithm, to accelerate the simulation time required for the prediction of a discrete twin formation and the resulting heterogeneous deformation.

5.13. Phase-Field Twinning Model

In a phase field model, twin and parent regions are regarded as two different phases and described by an order parameter $\eta$: $\eta =0$ for parent region and $\eta =1$ for twinned region, and $0<\eta <1$ represents the interface. While the evolution of the twin order parameter ($\eta$) is driven by a discrete twin interface energy, the twin nucleation mechanism is often incorporated using the stochastic method. The Ginzburg–Landau equation describing the twin evolution ($\eta$) Liu et al. [142] is shown below

$$\dot{\eta }=M\left(\mathrm{Div}\mathbf{\kappa }\nabla \eta +\mathsf{\Delta }f(2\eta -1)-{\delta }_{\eta }{\varphi }_{CP}-\mu \right)$$

(100)

where M is a mobility parameter, $\kappa$ is a symmetric second–order tensor describing anisotropic twin interface energy, $\mathsf{\Delta }f$ is the height of the GIBBS energy barrier between the twin and its parent, and $\mu$ is the sub–differential of the indicator function. The mechanical contribution to the free energy density is represented by ${\varphi }_{CP}$ and derived using the crystal plasticity formulation [142].

5.14. Others

A CP model has also been incorporated within the Marciniak–Kuczynski (M–K) approach for forming limit curve predictions [174,175], which is not discussed in this paper. A thermodynamics based twin nucleation model was proposed by Mareau and Daymond [156]. Similarly, Zhang et al. [176,177] have used FE methods to perform a parametric study of stress state development during twinning. Recently, Grilli et al. [178] incorporated discrete twin nucleation and growth using criterion driven by stress-induced nucleation as well as an atomic driving force based on energy barrier. They used CPFE with the discrete twin model to study crack nucleation near the twin tip. A summary on each twinning criterion discussed above is tabulated in

Table 1. A summary of twinning criteria used in crystal plasticity models.

SN	Method	Ref.	Capabilities
1	PTR scheme	[5,6,125,138,163,174,179,180,181,182,183,184]	- stress response for the materials system that involves deformation from twinning - predicts texture evolution using predominant twin system
2	VFT scheme	[6,130,185,186,187]	- discrete orientation space and texture evolution driven by volume fraction transfer - capable of predicting more accurate texture evolution
3	Lagrangian method	[7,8,188,189,190]	- accurate prediction of texture evolution - considers volume fraction evolution for all twin variants
4	Updated Lagrangian method	[9,99,115,158,183]	- includes the effects primary, secondary, and/or tertiary twins in constitutive relations - study effect of extension twinning on dynamic recrystallization in Mg AZ31 alloy
5	Composite Grain model	[4,118,191]	- incorporates physical twin in a CP model - both twinning and detwinning can be implemented
6	Twinning detwinning moddel	[11,116,162,163,192,193,194]	- incorporates detwinning mechanism in a CP model - simulate stress–strain hysteresis due to twinning and detwinning - describe the tensile and compressive yield asymmetry, anisotropy, and strain hardening behavior along arbitrary directions [194]
7	Dislocation density based model	[10,119,168,191,195,196,197,198,199]	- includes dislocation transmutation and twin accommodation effects to predict plastic anisotropy of texture Mg alloys - Reproduce stress–strain response, anisotropic hardening, and evolution of twin volume fraction and texture - incorporate directionality of dislocation glide and de-twinning and their effects on the evolution of stored dislocation density populations and texture evolution
8	Probabilistic nucleation method	[32,33,109,137,142,188,200,201]	- bridge between grain boundary defect interactions at the atomic scale, localized stress concentrations, intergranular stresses at the grain sclae and effective material response at polycrstalline aggregate level - effect of grain size and stress on twin volume fraction, fractional twin length and the fraction of twin contact - effect of grain size on the yield stress - effect of grain sie on the general shape of the stress–strain curve at low strains
9	Twin nucleation, propagation, and growth model	[12,114]	- stress relaxation effect for twinning - explicitly show a twin propagation in the parent grain
10	Explicit incorporation method	[101,111,112,136,139,140,190,202,203,204]	- simulates local stress fields produced by twins - relate spatially resolved fields of stress, strains with microstructural changes during twin formation and thickening
11	Energy-based micro-twin nucleation model	[109,110,173,178,205,206]	- study on large grain aggregates revealing the critical role of crystallographic orientation and grain boundaries on micro-twin formation - heterogeneous twin formation with strain localization - twin variant selection
12	Thermal activation based twin propagation method	[110,173,205]	- predict explicit twin formation and heterogeneous deformation in hcp materials such as Mg - simulate heterogenous deformation and strain localization due to twinning effects
13	Phase field twinning model	[141,142,143,144,145]	- couples twin evolution as a phase transformation - predict twin nucleation, propagation, and growth - study twin morphologies, twin thickness, as well as interaction of twins with other grain defects

. It should be noted that the PTR scheme or the VFT scheme is incorporated along with other twinning criterion so as to include a physical twin in a simulation domain.

Figure 14 shows how these twinning criteria fit with incorporating different mechanisms of twinning, such as twin nucleation, twin propagation, or twin growth, or how the twins are represented such as physical twin (with actual reorientation), or tracking of twin volume fraction. The schematic shows that more efforts have been made to incorporate the twinning mechanisms associated with twin growth. Only in recent years, researchers have turned their focus to twin propagation and twin nucleation. However, it is essential to incorporate all twin nucleation, propagation, and growth mechanisms with the crystal plasticity framework to capture the damage in hcp metals.

6. Gaps and Future Directions

6.1. Preface

Despite appreciable efforts over the last three decades by various prominent authors, a physically motivated model that reflects the early stages of embryonic twin nucleation and lengthwise burst is lagging. This can be clearly understood when considering the bases of the phenomenological criteria employed to capture twinning, exemplified by the pseudo-slip approach, the absence of site-specific nucleation, and the lack of subsequent evolution of the observed twin spacing, including the case of multivariant twinning. Likewise, the evolution of deformation twins with strain and ultimately calculations of constitutive response along with a proper description of twin interfaces and twin boundaries are yet to be incorporated into the CP models [31]. The dynamics of twin expansion have not been treated realistically either. The Schmid factor based on far-field stress turned out to be inadequate to describe variant selection [31].

While atomistic simulations and experimental results provide great insight into the mechanisms underlying twin nucleation and propagation, there is no widely accepted mechanism for twin nucleation and twin expansion. Continued progress in atomic–scale studies can provide the needed insight to develop physics-based criteria for nucleating and propagating twins in the CPFE and FFT frameworks.

The following subsections pinpoint some key gaps and propound ideas that might have the potential to fill these gaps for what they are worth.

6.2. Embryonic Twin Nucleation Mechanisms and Their Crystal Plasticity Rendering

One of the major challenges in crystal plasticity modeling efforts for twinning is to render twin nucleation events and their localized effects. While nucleation of the twin itself is a great focus of study in the materials community, major advances on the understanding of twin nucleation events have been made through experimental observations and atomistic simulations. In the following subsections, we will discuss major features/mechanisms that assist in twin nucleation and how these mechanisms can be incorporated into crystal plasticity models.

6.2.1. Faceting vs. Pure Shuffle for Nucleation

Twinning is mediated through shear and shuffle at the atomic level. While some researchers [207,208] may argue that $\left\{10\overline{1}2\right\}$ twin mode in magnesium is predominantly mediated through pure shuffle without the action of twinning dislocations, [209,210,211,212] showed that the purely “shuffle-dominated” mechanism for twin nucleation is inaccurate. Serra et al. [213] presented how the core width of active twinning dislocations is affected by the complexity of atomic shuffles, whereas El Kadiri et al. [25] derived the values for shear and the net shuffles required for each atom at the intermediate planes for the $\left\{11\overline{2}2\right\}$ and $\left\{10\overline{1}2\right\}$ twins. In particular, shear displacements and net shuffles for ${\mathit{b}}_2$ twinning dislocation for the $\left\{10\overline{1}2\right\}$ twinning mode at the intermediate planes are shown in

Table 2. Shear displacements and net shuffles for $b_2$ twinning dislocation for the $\left\{10\overline{1}2\right\}$ twinning mode.

Plane n	Shear Displacement	Net Shuffle
$0_a$	$\frac{1}{12}{\mathit{b}}_2$	$-\frac{1}{12}\left[10\overline{1}1\right]$
$0_b$	$-\frac{1}{12}{\mathit{b}}_2$	$\frac{1}{12}\left[10\overline{1}1\right]$
$1_a$	$\frac{7}{12}{\mathit{b}}_2$	$\frac{1}{12}\left[\overline{1}011\right]$
$1_b$	$\frac{5}{12}{\mathit{b}}_2$	$-\frac{1}{12}\left[\overline{1}011\right]$

. The net shuffles for $\left\{10\overline{1}2\right\}$ twinning in Table 2 depict that the shuffle directions switch from one K1 plane to another, and the signs of the net shuffles are opposite in the same plane.

Li and Ma [207] argued for shuffle-mediated twinning by explaining the large deviations of twin boundaries from the twinning planes using the shuffle-dominated mechanism. Similarly, Tang et al. [208] showed that the shuffle-induced twin nucleation energy requirement is on the same level or even less than the energy required for a partial pyramidal dislocation when the c–axis stress/strain is sufficiently large, which resulted in an initial 90° misorientation. This idea is debunked in the works by [16,210,211,213] as faceting is responsible for these large deviations in twin boundaries. Barrett and El Kadiri [16] showed that faceting mechanics are vital to understand the nucleation of twin embryos and mobility of twin boundaries. They also showed that faceting promotes the nucleation of an octahedral faceted twin from a defect inside of a single crystal magnesium, and the mobility of twinning dislocation along these facets, disclination dipoles, and cross-faceting, among others, aid in the nucleation and growth of twin embryos into a stable twin.

Crystallographic evidence has accumulated that shuffling plays a major role in the nucleation events for hcp materials. Atomic shuffling is, in fact, required for all twin modes with the exception of $\left\{11\overline{2}1\right\}$ shear twinning. Atomic shuffle pertains to atomic motion by diffusion. Thus, the associated driving force must be connected to the gradient of the hydrostatic pressure and not only the deviatoric stress, which is traditionally utilized for plasticity. This is why triaxiality should also play an important role in the embryonic stages of twin nucleation, consistent with the observed needs of defects to trigger $\left\{10\overline{1}2\right\}$ twins and higher stress levels tallied for nucleation than propagation. Consequently, any twin nucleation constitutive model should engage both the resolved shear stress and hydrostatic pressure, and this is yet to be implemented in the flow rule.

Corroborating evidence for the need to including hydrostatic pressure in the flow rule is supported by the results recently reported by Russell et al. [214] on in situ EBSD characterization of tension and compression of an AM30 extruded sheet showing a double in–plane basal fiber. It was found that $\left\{10\overline{1}1\right\}$ twinning occurs easier under contraction than the compression of the 〈c〉–axis with nearly an order of magnitude difference in the CRSS levels. In a contraction $\left\{10\overline{1}1\right\}$ twin, <c> axis compression happens as a result of $\left\{10\overline{1}0\right\}$ tension, thus resulting in different stress states than that of actual <c> axis compression [214]. The difference of CRSS in compression and contraction $\left\{10\overline{1}1\right\}$ twins was explained by the higher triaxiality naturally induced along the 〈c〉–axis under contraction than compression. It confirms the importance of the hydrostatic pressure in twin nucleation, especially for twinning disconnections requiring substantial shuffles at their core of $\left\{10\overline{1}1\right\}$ $b_4$ disconnection. The effect of hydrostatic pressure during $\left\{10\overline{1}2\right\}$ twinning is further elaborated in the work of Serra and Bacon [21]. In the paper, the local fluctuation in hydrostatic pressure at the twin interface is demonstrated. Wang et al. [50] suggested a process where a twin embryo nucleates from lattice dislocation pile-up at a symmetric tilt wall. This now acts as stress-concentrator motivating nucleation through “pure-shuffle” rather than the dissociation of lattice dislocations into twinning dislocations. All of these twin nucleation mechanisms exhibit the necessity for the current crystal plasticity models to incorporate the combined effects of shear and shuffle during twin nucleation events.

6.2.2. Grain Boundary Energy and Defects

Another effect of considerable importance to polycrystalline materials is the sensitivity of twin nucleation to the grain boundary atomic structure. Generally, grain boundaries with low misorientation angles and high energy have been observed to contain more twins. However, this could be an artifact of the autocatalysis phenomena [215] or the fact that a highly misoriented grain is inherently associated with a neighboring grain, which has a low Schmid factor, so few twins would nucleate at this grain boundary. Generally, twins nucleate at the free surface of the tested sample, and if the texture is sharp enough, they find their way to punch through boundaries and continue propagating by autocatalysis. In the case of weak textures (e.g., rare-earth containing Mg alloys), however, grain boundary nucleation is needed in grains that are not completely surrounded by the grains with a high Schmid factor for the same twinning mode. Therefore, the sensitivity of twinning to grain boundary atomic structure must be considered.

Twins often nucleate from grain boundaries. Misorientation [33], interface angle, defects, segregated elements, and strain rate are some of the major nucleation criteria. Studies show that symmetric tilt grain boundaries with the $\langle \overline{2}110\rangle$ axis are able to spontaneously twin under a relatively small loading and a pile-up of basal dislocations [50]. Boundaries with this axis and low angle misorientation have high energy and hence are likelier to encourage twin nucleation by grain boundary energy relaxation.

Twin nucleation, a heterogeneous event, depends upon the embryo’s interaction with existing defects. The effect of dislocations, vacancies, tilt axes, and tilt angles have been studied by Giri et al. [71,216] using atomistic simulations. Their studies show that while the twin nucleation energy barrier increases with a higher number of vacancies, twin nucleation eventually becomes easier as there is more energy relaxation. Through atomistic simulations on $\left\{10\overline{1}2\right\}$ twin nucleation on a number of $\left[\overline{2}110\right]$ axis-symmetric tilt grain boundaries, the study showed that there is a simple relation between grain boundary energy and applied stress that has the potential to be implemented in larger scale plasticity models. The observation of the most stable $\left\{10\overline{1}2\right\}$ twin embryo on the $\left\{10\overline{1}1\right\}$ grain boundary depicts the significance of the misorientation angle, apart from the grain boundary energy, on the nucleation of a twin.

Figure 15 shows that 〈a〉 tilt grain boundaries demonstrate a strong correlation of grain boundary energy with the stress required for exothermic twin nucleation [216]. This is consistent with the experimental observation that shows the correlation of the non–Schmid twin nucleation behavior with the neighboring grains and the grain boundary effects [214,217,218]. Basu et al. [218] showed that the twin nucleation was favored for the twin variant that has a higher local Schmid factor. Furthermore, highly ordered low-energy boundaries nucleate twins more easily than disordered boundaries, and the addition of vacancies into the boundaries encouraged twin nucleation [216]. This correlation between nucleation stresses, grain boundary energy, and grain boundary defects shows a major role in the twin nucleation mechanism that should be generalized and inform a deterministic site-specific twin nucleation model for the robust prediction of twin nucleation sites, double twinning, or secondary twinning.

6.2.3. Lengthwise Shooting Mechanism

Atomistic simulations performed by Barrett and El Kadiri [16] showed that a stable twin embryo nucleated through faceting acquires an octagonal diamond–shaped morphology, as shown in Figure 16a. Experimental observations present a propagating twin as a sharp pointed lenticular twin. In metals, the deformation twin propagates at the speed of sound [219]. The crystallographic explanation for rapid lengthwise propagation of twins has been a challenge, though Oberson and Ankem [220] attempted to explain why some deformation twins do not propagate at the speed of sound.

Using the interfacial defect theory, Barrett and El Kadiri [16] proposed an explanation on why embryonic twin propagates at a rapid speed. They used faceted boundary mechanics to understand the nucleation of twin embryos and the mobility of twin boundaries. A twin may readily use boundaries other than the twin plane boundary and overcome obstacles or lower interface energies through faceting. In the particular case of $\left\{\overline{1}012\right\}$ twins, twinning disconnection pile–up is converted into new facets on the BP, PB, or $\left\{10\overline{1}2\right\}$ twin boundaries. The BP and PB boundaries that are fully coherent contain long–range strains proportional to their length, making them unstable except over short distances. This explains the octagonal shape for a stable twin nucleus, as shown in Figure 16a. After the critical size of the embryonic twin, BP and PB boundaries are highly unstable; thus, they transform into a kinematically equivalent pile–up of highly mobile twinning partials, which, finding themselves outside the boundary, immediately shoot over a grain or several grains at the speed of sound (see Figure 16b). This possible explanation for rapid twin propagation adheres with the fact that high stress concentrations develop ahead of a shooting twin [68].

6.2.4. Computation of Adequate Twin Spacing

Imagine a single crystal or a grain in a polycrystal undergoing a very high resolved shear stress on a twin plane. Though every single plane of the grain parallel to this plane is equally highly stressed, only a few localized planes will yield to twinning. The result is a single or a few twin lamellae burst with a characteristic spacing. It would be hard to imagine that there might be a way to exactly capture the sites on which these twin lamellae nucleate unless the free surface or grain boundary experiences large stress fluctuations due to the presence of defects (e.g., case of intergranular precipitates), which is only a particular case and not a general rule. However, twin spacing must be properly captured to adequately simulate the microstructure evolution. So, perhaps instead of site–specific nucleation, a sound crystal plasticity model should accurately reproduce how many twins should first nucleate per unit grain boundary area as well as their spacing. There should also be a critical strain where no more nucleation events at the grain boundary are needed, paving the way to homogeneous twin nucleation to begin in spaces that are confined by existing twin lamellae from various variants [221].

Figure 17a shows one such scenario where the multiple twin variants nucleate, grow, and engulf a parent matrix showing a pattern of multiple twins nucleation with characteristic spacing [221]. Later, Paudel et al. [68] provided the insight into the mechanisms that drive twin spacing. In fact, they developed a micromechanics–based analytical model that was able to suggest the minimal distance at which a second twin could simultaneously nucleate in an idealistic single variant twinning scenario within a single crystal. The analytical solution to the stress field was obtained for a spherical oblate twin domain with uniform eigenstrain in homogenous, isotropic, and half–infinite space (see Figure 17b). The model was validated through in situ three–point bending characterization of a highly textured AZ31 sheet, which behaves similarly to a single crystal. Under these conditions, twinning occurs by localized bands separated by a critical distance during nucleation, which conformed well with model predictions. Furthermore, the model justifiably suggested that this critical distance depends upon the height of the twin domain, while the maximum stress relaxation depends upon the thickness of the twin domain.

6.3. Twin Thickening

6.3.1. Thickening of Different Variants and Effect of Grain Boundaries

The selection of variant in $\left\{10\overline{1}1\right\}$–$\left\{10\overline{1}2\right\}$ double twins are found to be strongly dependent upon the orientation and loading condition (tension or compression) [64]. For primary twins in Mg, the variant selection for twins does not fully obey Schmid’s law [222]. Jonas et al. [222] observed that 5% of the twins had very low Schmid’s factors, ranging from 0.03 to 0.15, and more interestingly, many potential “high Schmid factor” twins did not form. They also found that twin variant selection was associated with the “easy" accommodation of strain in highly plastic anisotropic Mg grains. Figure 18 shows the $\left\{10\overline{1}1\right\}$ contraction twin boundaries and variants of double twin boundaries after tensile strains of 0.08 and 0.15 strains. The statistical analysis of twin variants, as observed in Figure 18, shows that two of the secondary twin variants ( 38 ${}^{\circ }$ and $69.9$ ${}^{\circ }$) are more favorable than the others. The former twin variant is favored due to the least amount of strain accommodation requirements, whereas the latter twin variant is favored by CRSS (obeys Schmid’s law) [223].

Furthermore, there could be an absence of potential variants for primary, secondary, and tertiary twins in a deformed Mg alloy [224]. The selection process is based on the amount of work required for strain accommodation in most of the neighboring grains [222,224,225]. The statistical analysis on primary $\left\{10\overline{1}1\right\}$ contraction twins performed by Jonas et al. [222] showed that Schmid factors ranging from 0 to 0.5 are associated with the observed twins (see Figure 19). Their results suggested that the value of the Schmid factor itself does not determine the criterion for twin variant selection. Figure 19d reveals that ∼45% of the twins had Schmid factor ranks of 3 or 4 ($0.15<SF<0.3$). Another study showed that the selection of the twin variant for extension twinning in Mg is determined by a geometric compatibility factor $\left(m^{\prime }\right)$ between slip and twin [226]. This type of non–Schmid’s behavior for twin variant selection have also been observed for titanium alloys [227,228].

It is not clear yet, however, whether this behavior is due to grain–grain interactions so it could be systematically captured by a sound full–field crystal plasticity model or due to the grain boundary atomic structure, which would then truly pertain to a non–Schmid’s behavior that warrants considerations of stresses other than the shear on the twin plane.

6.3.2. Effect of Twin–Twin Interactions

Twin–twin interactions leading to secondary twinning and nanovoids [37,229,230] are abundantly evidenced through experiments and atomistic simulations. These twin–twin interactions are responsible for stress localization and damage [37,231]. It has been observed that twin nucleation is quickly saturated with the activation of a single twin variant, where the strain accommodation is provided by twin thickening; however, the twin propagation rate is reduced with the activation of two twin variants in Mg AM30 [221]. Although twin nucleation stress is higher than twin propagation stress, twin–twin interaction in grains with two twin variants results in profuse twin nucleation, which also challenges the classical assumptions of twinning stress sensitivity to grain size [221]. The twin–twin interactions have been correlated with an effect on hardening [221,229,232].

The thickening of a given twin in a crystal plasticity model, though can be handled with PTR, shows discrepancy if the twin has an excessive source of twinning disconnections compared to other engulfed variant as observed in the experiment [214]. This is the case of an early twin nucleating at a grain boundary segment. Twin variants that can thicken along the grain boundary segments showed growth within a grain compared to other variants that are more likely to get engulfed away through twin–twin interactions [214]. A way for crystal plasticity to capture such an effect is to consider slightly higher CRSS for twin nucleation other than grain boundaries. Furthermore, twins that start meeting with other twins will have their thickening become sluggish, allowing for twin nucleation in core regions. At the macroscale level, this effect is observed on the twin band formation where nucleation of a second twin variant obstructs the thickening of initial twin bands observed during bending of highly basal textured Mg AZ31, as seen in Figure 20.

6.3.3. Detwinning and Pseudoelasticity

A pseudoelastic behavior was observed upon stress removal during a three–point bending experiment [233], which is associated with untwinning. Experiments have shown that detwinning is not strictly limited to reverse loading but can also occur during unloading [11,192,234,235,236,237,238,239]. Twinning and detwinning processes have been associated with the pseudoelastic behavior in hcp materials, such as zirconia [148], nitinol [240], and magnesium [238]. The final microstructure from the three–point bending experiment shown in Figure 20 shows the role of the second twin variant in accommodating higher strains during the maximum loading. The residual twins indicate the untwinning of these variants upon unloading that resulted in pseudoelastic behavior. This extensive untwinning, not seen in simple compression or tension, is associated with the significant fraction of stopped elastic twins in Garber’s notation. It is because autocatalysis phenomena are stopped at the neutral axis, whereupon accommodation effects are not complete. Localized twins form during a three–point bending load. However, since they do not reach the free surface, unresolved twins vanish with unloading.

6.4. Crystal Plasticity Rendering

Current crystal plasticity models are limited in terms of physics-based site–specific twin nucleation or crystal plasticity rendering of twin nucleation events. Any CP model attempting to capture the phenomenon discussed in previous subsections must use a size threshold factor compatible with the maximum separation distance of the disclination dipole associated with the coherent BP boundary before its disintegration into twinning dislocations. Thereafter, the lattice regions laying ahead of the embryo should be reoriented by simple geometrical considerations to make it grow till the next grain boundary. Recently, Jin et al. [201] used a probabilistic approach to explicitly incorporate deformation twinning in a crystal plasticity model, as seen in Figure 21a. Likewise, Cheng et al. [205] used explicit twin nucleation based on the non-planar dissociation process of sessile pyramidal dislocations, which are assumed to follow a Weibull–type probability distribution function. Multiple twins within a grain, see Figure 21b, nucleate with a characteristic spacing as a result of stress field from the initial twin. In contrast, Paudel et al. [241] used a characteristic twin spacing parameter, which is a function of twin length, as established by the micromechanics-based analytical solution [68], to explicitly introduce multiple twins inside of grain (see Figure 21c). Within a CPFFT framework, Paudel et al. [241] demonstrated a way to include the effect of both shear and shuffle on site–specific nucleation. Figure 21c shows the multiple twin nucleation in a soft grain of a tri–crystal system and autocatalysis of different twin variants in another twin. Combined with the dislocation density-based hardening mechanisms, hydrostatic stress gradient driven criterion is incorporated for a site–specific twin nucleation criterion, rendering the twin microstructure observed in the experiment, as seen in Figure 17a. The inclusion of hydrostatic pressure as a criterion for twin nucleation renders the effect of initial grain boundary defects and grain boundary energy into the crystal plasticity model. The site–specific twin nucleation criterion used by the authors can nucleate a non-Schmid twin, where hydrostatic stress is favorable in the nucleation of the twin. Figure 22 shows the site–specific twin nucleation based on both shear and shuffle contribution, where the overall volume fraction of the twin is determined using a shear calculation, but the site selection is determined with hydrostatic stress gradient.

Figure 22a,b show a Schmid factor of the most favorable twin variant for each grain and hydrostatic stress gradients in a seven crystal system, whereas Figure 22c shows the twin nucleation based on the shear and hydrostatic stress gradient. The twin system favorable inside of the grain that crosses a threshold volume fraction of twin ($f^{thres}$) at the twin nucleation site is shown by the red arrow in Figure 22b. For the site–specific twin nucleation, the code chooses the region at the grain boundary with a high hydrostatic stress gradient. The results show that the combination of shear and shuffle can produce an accurate twin microstructure evolution for HCP metals. Although more work is required to fine–tune the algorithm, authors believe that the crystal plasticity models should refocus on the inclusion of hydrostatic stress and shuffle mechanisms in the rendering of site–specific twin nucleation. In addition to the hydrostatic stress gradient, one can quantify the grain boundary energy based on the defects, dislocation density, and grain misorientation [216] to make a robust twin nucleation criterion.

7. Conclusions

This paper reviews existing accomplishments with the incorporation of deformation twinning into crystal plasticity models and the challenges they faced to break away from a pseudo-slip approach, which has been consistently discounted by many phenomena related to twin nucleation. In light of recent discoveries achieved by novel in situ experimental techniques and powerful simulations at the atomistic level, we emphasized a few key mechanisms for twin nucleation and suggested methods to reflect those mechanisms in a crystal plasticity model:

1.

Just like shear banding, twinning is a localized event, but one which challenges the ability to deterministically pinpoint a favorable nucleation site. Due to the shuffling requirement in hexagonal closed packed metals (with the exception of $\left\{11\overline{2}1\right\}$ twinning), the choice of a nucleation point must account for the local hydrostatic pressure gradients because twins always nucleate at a defect structure, most notably, the sample–free surface or otherwise grain or phase boundaries. This incites a radical modification of the flow rule to include the trace of the stress tensor in addition to its deviatoric part. Furthermore, the minimum twin thickness, which does not seem to be in the order of a few plane layers, number of twin variants inside a grain, and an adequate twin spacing are factors that could affect each other but would need to be quantified in order to predict the microstructure evolution during twinning and the ensuing local deformation.

2.

Variant selection for primary twins does not always obey Schmid’s law due to the atomistic structure of the grain boundary, which seems to affect the adopted plan for twin transformation. This is quite a challenge for crystal plasticity, and it is still unclear how the choice of a variant is made during the nucleation stages at the grain boundary. Experimental evidence has accumulated, however, that grain boundaries offering a good combination between low misorientation angles and high energies are the best candidates for twin nucleation.

3.

Profuse twin nucleation and rapid lengthwise propagation of twins invoke interactions between different twin variants. The twin–twin interactions are correlated with material hardening, so the incorporation of these interactions is warranted.

4.

The rapid lengthwise propagation of twins across a grain occurs after the basal–prismatic facet reaches a critical size that leads to its dissociation into twinning partials. These shoot forthwith across the grain or across several grains in case the stress value and state ahead of the twin at the neighboring grain is adequate for another twin nucleation by virtue of the autocatalysis behavior.

5.

Recent digital image correlation by the present authors emphasized substantial pseudoelasticity upon stress removal during three–point bending, which resulted from significant untwinning behavior. This extensive untwinning, not seen in simple compression or tension, is associated with the significant fraction of stopped elastic–twins in Garber’s notation as autocatalysis phenomena are arrested at the neutral axis, whereupon accommodation effects are not complete. Not capturing this phenomenon in crystal plasticity would result in significant errors in spring–back predictions for sheet-forming operations.

6.

Twin nucleation mechanisms exhibit the necessity for current crystal plasticity models to incorporate the combined effects of shear and shuffle during twin nucleation events. Few advances in site–specific nucleation criterion based on probabilistic and dissociation-based methods are discussed along with the preliminary results for twinning site–selection based on hydrostatic stress gradients.