Crystal Plasticity Modeling of Anisotropic Hardening and Texture Due to Dislocation Transmutation in Twinning

In crystalline materials, dislocations are three-dimensional lattice distortions that systematically distort twin interfaces that they encounter. This results in dislocation dissociation events and changes in the atomic structure of the interface. The manner in which the interface distorts drive the product of the dissociation event, and consequently, the incident dislocation core and the magnitude and relative direction of the Burgers vector govern these slip-twin interaction phenomena. Recent characterization studies using transmission electron microscopy as well as advanced molecular dynamic simulations have shown that slip dislocations, whether striking or struck by a {101¯2} twin boundary, dissociate into a combination of twinning disconnections, interfacial disclinations (facets), jogs, and other types of dislocations engulfed inside the twin domains, called transmuted dislocations. While twinning disconnections were found to promote twin propagation, the dislocations incorporated inside the twin are of considerable importance to hardening and damage initiation as they more significantly obstruct slip dislocations accommodating plasticity of the twins. In this work, the dislocation transmutation event and its effect on hardening is captured using a dislocation density based hardening model contained in a visco-plastic self-consistent mean-field model. This is done by allowing the twins to increase their dislocation densities, not only by virtue of slip inside the twin, but also through dislocations that transmute from the parents as the twin volume fraction increases. A correspondence matrix rule is used to determine the type of converted dislocations while tracking and parameterizing their evolution. This hypothesis provides a modeling framework for capturing slip-twin interactions. The model is used to simulate the mechanical response of pure Mg and provides a more physically based approach for modeling stress-strain behavior.


Introduction
The utilization of hexagonal-close-packed (HCP) magnesium (Mg) and magnesium based alloys is of interest to industry because of their low density and high specific strength. Their low crystal symmetry leads to mechanically anisotropic behavior, especially in highly textured polycrystals, in which significant nucleation and growth of {1 0 1 2} tensile twins may occur. {1 0 1 1} compression The findings of Oppedal et al. [14] were bolstered by the more recent studies of slip-twin interactions in the works of El Kadiri et al. [20], Barrett et al. [21], and Wang et al. [22], in which the behavior of basal dislocations encountering {1 0 1 2} twin boundaries were investigated by means of molecular dynamic simulations and transmission electron microscopy. These studies showed that such dislocations readily disassociate into twinning disconnections bolstering twin propagation, interfacial disclination dipoles (facets), and other dislocations that transmute into the interior of the twin volume, called transmuted dislocations. Though the complete dislocation reaction at the twin interface was only recently identified by means of atomistic simulations and advanced interfacial defect theory, the part of the dissociation leading to dislocation transmutation has been hypothesized as early as the nineteen sixties by Bilby and Saxl and was later known as the Basinski mechanism [23][24][25]. Similar to the Basinski mechanism [25], the capacity of dislocations to transmute from parent grains across twin boundaries represents a potential vehicle by which the dislocation density of twin grains might be increased in excess of statistically and geometrically necessary dislocations typically formed to accommodate slip within the twin. Transmutation could then account for the increased dislocation density enforced by the twin storage factor of Oppedal et al. [14]. While previous studies have conceded that a Basinski type transmutation effect may be occuring in twinning Mg, experimental evidence has been sparse thus far [2,15,16,18]. The work of El Kadiri and Oppedal [26] provides a crystal plasticity framework by which such transmutation effects might be modeled in VPSC simulations. They developed a transmutation matrix α where each element α ij represents the proportion of dislocation density that is transmuted to the ith slip mode inside the twin volume from the jth slip mode in the parent volume.
In order to investigate the potential role of dislocation transmutation effects in the hardening of Mg, the work presented herein presents a modified version of the crystal plasticity based dislocation transmutation theory of El Kadiri and Oppedal [26], further adapted to capture the system-wise nature slip dislocation interaction with twin boundaries and allow for the inclusion of both transmutation and dissociation behavior of these dislocations. In this way, a physically based model for the complex interactions of slip and twinning in polycrystal Mg is developed. This model is implemented in the LANL code VPSC-7d, and is calibrated against data from rolled Mg under multiple compression load paths. Finally, it is used to simulate the anisotropic hardening behavior of the material.

Model
The hardening model for this work is based on the approach of Beyerlein and Tomé [19] in which the hardening behavior of HCP polycrystals is simulated by the evolution of the critical resolved shear stress (CRSS) used in Hutchinson's strain rate sensitive constitutive equation for slip. In this approach, the CRSS for slip for a given slip system s as part of the slip mode, i, is defined as a function of the dislocation density on that system. The CRSS of each system is given by the additive decomposition Here, τ i 0 , τ i f orest , τ i deb , and τ s HP are the initial critical stress values, hardening from forest dislocations, hardening from debris formation, and hardening from Hall-Petch effects, respectively. They can be defined as χ and µ are the dislocation interaction coefficient and shear modulus, respectively. Equation (2b) also contains the debris formation parameter k deb and debris density ρ deb . The Hall-Petch contributions in Equation (2c) depend on the presence of twins. In grains where there is no active twinning, these contributions are the same for all s ∈ i, with a bulk interaction term HP i and grain size d g .
For grains containing twins, the values of τ s HP cannot be assumed to be the same for all s ∈ i. In this case, the interaction term HP ij describes the interaction between twin mode j, which contains the predominate twin system, and slip mode i. Here, d s,PTS m f p is the mean free path between twin lamella for dislocations on slip system s as described in Proust et al. [17].
In the simulations presented in this work, the forest and debris contributions remain unchanged from the forms presented in Beyerlein and Tomé [19]. However, the constitutive model of Beyerlein and Tomé [19] is based on the Composite Grain (CG) model of Proust et al. [17] and treats twin boundaries as full-stop barriers to dislocation glide. As twin boundaries are not considered to act as barriers in this work by hypothesis, the Hall-Petch contributions to CRSS are treated as being identically zero, hence In the model of Beyerlein and Tomé [19], the change in dislocation density for a given slip mode α with change in shear on the same system is calculated by a Mecking-Kocks [27] type equation and adopted in this work to calculate the change in dislocation density with shear as Equation (4) introduces the phenomenological generation and rate and temperature dependent parameters k i 1 and k i 2 . In Beyerlein and Tomé [19], this leads to a relationship between these two parameters, Above, g i , D i are the normalized activation energy and drag stress, respectively, and˙ and˙ 0 are the strain rate and reference strain rate.
In the present model, the dislocation densities were modified based on changes in twin volume fraction of the composite grain, which were added: In this way, the hardening contributions from dislocation density to the CRSS of slip systems inside the twin volume fraction of each grain are differentiated from those systems inside the parent volume fraction.
The function ρ iTF (V) from the work of El Kadiri and Oppedal [26] defining the increase in dislocation density of the ith slip mode type inside the jth twin volume fraction is expressed as In the present work, the matrix α includes an additional row in order to account for dislocation transmutation into sessile, higher order dislocations and Equation (7) is modified in order to account for the dissociation of parent dislocations at twin grain boundaries. A normalized dissociation parameter η is introduced, representing the proportion of dislocation density that fails to transmute from the jth slip system in the parent volume to the twin volume and does not contribute to dislocation density in the twin. This leads to the final, modified constitutive law The model of Beyerlein and Tomé [19] also contains constitutive laws for twin activation, propagation, and interaction. These remained unchanged in the current work and are summarized in Appendix A.

Implementation and Calibration
The above described changes to the dislocation density model were implemented in the LANL code VPSC-7b by making changes directly to the subroutines UPDATE_CRSS_CG_disl and DATA_CRYS and were compiled using the GFortran GNU compiler. The model was calibrated using a three stage process. First, the values for α were assigned based on the correspondence method introduced by Bilby and Crocker [23]. After this, values for the elements η j were assigned, and finally, the slip and twinning parameters respectively associated with dislocation generation, twin nucleation, and twin propagation were adopted from Oppedal et al. [14].

Correspondence Method for Transmutation for the Construction of the α matrix
A generalized correspondence method for mapping a vector on a slip plane in a parent grain on to a corresponding slip plane inside a twin grain was first presented in the work of Bilby and Crocker [23]. This method was applied to twin modes in FCC and HCP crystals in the works of Niewczas [28,29]. Summarized in Appendix B, the precise mathematical expression of the general theory of the correspondence method shown in Niewczas [29] in the case of Mg was adapted for the purpose of assigning values to the transmutation matrix α.
The onto mapping of slip systems in the parent volume to corresponding slip systems in the twin volume implemented by Niewczas [29] provides explicit calculations for transmuting the plane normal and directional vectors for each slip system in a parent grain to its corresponding slip system inside a twin. Rather than using the method to calculate values of α at each deformation step, the transmutation matrix for {1 0 1 2} twins was constructed in pre-processing. {1 0 1 1} twins were not assumed to contribute to transmutation in this work as the onset of compression twinning in Mg is taking place at higher strains, and very closely associated with the nucleation of damage in the material, leading quickly to brittle failure. As such, the contributions of compression twins to dislocation transmutation were assumed to be negligible. Strictly for the purposes of the element wise construction of α, the following assumptions were made: 1. The entirety of the dislocation density from a slip system in the volume fraction of the parent grain overtaken by a twin mode is considered to transmute to its corresponding slip system inside the twin grain. This pairing of slip systems is defined by the mappings defined by the correspondence method used by Niewczas [29] for HCP materials. 2. The dislocation density of slip mode in the volume fraction of the parent grain overtaken by a twin volume is considered to be evenly distributed across the systems of that mode. 3. In this simulation, the dislocation density of any parent slip system corresponding to a slip system inside the twin grain that was not part of the prismatic, basal, or 2nd order pyramidal < c + a > slip modes was assumed to contribute to debris formation inside the twin as part of ρ deb .
Working from these assumptions, the value of each element α ij can be calculated as the proportion of systems from the jth slip mode in the parent volume fraction which contribute to dislocation density in the ith slip mode in the twinned volume. This can be expressed as where n C j is the number of slip systems in the jth slip mode of parent volume mapped onto slip systems in the ith slip mode inside the twin volume and n tot j is the total number of slip systems in the jth slip mode of the parent volume. The values of α ij for Mg used in this work are summarized in Appendix C.

Parameters for Dissociation
The molecular dynamics simulations in the work of El Kadiri et al. [20] suggest that the type of dislocation and its orientation relative to the advancing twin boundary play a significant role in determining whether a given dislocation will transmute across the boundary or dissociate upon contact. In the simulations presented in this work, screw type basal dislocations were allowed to transmute perfectly across tensile twin boundaries. Edge type dislocations would either transmute across these twin boundaries or dissociate into twinning disconnections according to their orientation relative to the twin boundary, with positively oriented dislocations transmuting across the boundary and negatively oriented ones dissociating.
The simulations presented in this work were conducted at relatively low strain regimes, with ≤ 0.25 for all simulations. As such, the following assumptions were made: 1. Dislocations in the prismatic, basal, and 2nd order pyramidal < c + a > modes are assumed to transmute. 2. The incidence of screw type dislocations at low strain regimes is quite low. As such, for the purposes of this work, it was assumed that the contributions to dislocation density inside the twin volume fractions made by these types of dislocations were negligible. 3. Relative to twin boundaries, it was assumed that dislocations with positive and negative Burgers vector occur in equal measure.
The proportion of dislocation density dissociated at twin boundaries rather than transmuted is set by the η parameter at and summarized again in Appendix C.

Parameters for Dislocation Generation and Twin Nucleation and Propagation
Simulations from the work of Oppedal et al. [14] were taken as a baseline for the presented work, and parameters for the constitutive laws governing dislocation generation, drag stress, normalized activation energy, critical stress values, and deactivated Hall-Petch effects were taken from this work. Transmutation effects were used in place of twin storage factor, reducing the latter value to zero. The values of these parameters are summarized in Appendix A.

Simulation
Experimental data from the work of Oppedal et al. [14] was utilized to establish a control data set, where possible. In the case where such data were not available, comparisons were made to simulated results from the same work. In this work, rolled, pure magnesium with a highly basal texture in the normal direction was subjected to uniaxial compression along the normal and transverse directions, hereafter referred to as through thickness compression (TTC) and in-plane compression (IPC), respectively. This was done in order to capture the anisotropic mechanical behavior of textured magnesium, and, in particular, to motivate both scant and profuse twin volume growth. The chemical composition of the pure magnesium used to generate the experimental data in Oppedal et al. [14] is summarized in Table 1. The initial texture data was obtained via neutron diffraction at LANL and pole figures generated from this data are shown in Figure 1. Shown in Figure 2, the load paths for the simulations presented in this work were performed along both TTC and IPC directions on polycrystals consisting of 1944 grains with a highly basal texture in order to recreate the conditions under which the literature data were collected. The loading was conducted at room temperature conditions, with a strain increment of ∆ = 0.001, up to a total compressive true strain of = 0.25. The initial texture of the compressed Mg specimen is shown in Figure 1.

TTC Load Path
As highlighted in Figure 3, the TTC simulations reproduced both the texture and hardening curves from the experimental literature data with good accuracy. Although the strain softening observed in the experimental data at approximately = 0.1 was not reproduced in the simulation, capturing such behavior lies beyond the scope of the current work. The simulated modal activity, being defined as the proportional contribution to strain from each mode at each deformation increment, and evolution of twin volume fraction shown in Figure 4, both followed trends seen in simulations utilizing twin storage factor effects. Good agreement between the two simulation approaches was observed, however, tension twin volume was somewhat under predicted relative to simulations utilizing the twin storage factor approach.
From the literature, the texture in TTC simulations was not expected to change much over the course of loading. This was evident in Figure 5 showing results of the simulations utilizing the transmutation scheme, with the basal texture remaining consistent throughout the compression loading.

IPC Load Path
The sigmoidal stress curve associated with a high incidence of twinning was successfully captured in the IPC simulations. The associated simulated hardening rate was reasonably similar to experimental data from the literature. These results are shown in Figure 6. Shown in Figures 7 and 8, modal contributions to strain were nearly identical to simulation results from previous simulation literature as well. The growth of secondary compression twin volume fraction did exceed that seen in simulations utilizing the twin storage factor approach, but it is assumed that this can be neglected as the onset of compression twinning is, at higher strains, associated with void nucleation along compression twin grain boundaries leading to brittle failure. This fact was reinforced by the failure of the experimental control specimen coinciding with the formation of significant (>10%) volume fraction of secondary compressive twins within the primary twin volume.
In Figure 9, the textural evolution is shown to be consistent with experimental data from the literature, with high concentrations of orientations in the neighborhood of 90°from the normal specimen axis. This type of reorientation is consistent with the high degree of tensile twinning.

Discussion
As the evolution of yield stress and twin volume with strain of pure magnesium is heavily anisotropic, it becomes necessary to consider a measure for "goodness of fit" of simulated data along multiple load paths when evaluating the degree to which the transmutation model and the method for calibrating its parameters are validated. This measure must also be extended to simulated data obtained from simulations utilizing the TSF model in order to make truly meaningful conclusions regarding connections between the mechanical behavior resulting from artificially increased dislocation density in twin grains and the increased dislocation density resulting from transmutation effects. To this end the normalized mean squared error (NMSE) was calculated across the simulated strain sub-domains for which experimental data was available ( = 0 to 0.2 for TTC and = 0 to 0.125 for IPC). Summarized in Table 2, these values were calculated for both load paths and for both modeling approaches. For simulations of the TTC load path, the normalized mean squared error of the transmutation model is approximately 10% lower than that of the TSF model. In the IPC load path simulations, the normalized mean squared error is roughly 5% higher. Taking both of these comparisons into consideration, it can be concluded that while the transmutation model does suffer from some loss of accuracy in IPC simulations compared to the TSF approach, this loss is marginal at worst and more than offset by gains in accuracy in TTC simulations. This conclusion is further reinforced when one considers that the TSF model is a more empirical approach and, therefore, more limited than transmutation model by the availability of sound experimental data to insure its predictive capabilities.
η Sensitivity Although stress-strain curves, texture evolution, and modal activities were consistent both with previous experimental and simulation work, further validation of the modeling approach is necessary to conclude that predicted results are indeed the product of increased dislocation density inside the twinned volume fraction of the simulated polycrystal.
In experiments and in simulations using both the TSF method and the method from the current work, the saturation stress observed in the TTC and IPC load paths differs significantly by approximately 60 MPa. In order to confirm that the "gulf" observed between the saturation stress levels of TTC and IPC simulations was the result of increased hardness due to higher dislocation density brought on by dislocation transmutation, additional simulations in which transmutation effects were deactivated were performed for both load paths. This corresponds to simulations in which all of the dislocations dissociate upon reaching the twin boundaries, leading to dissociation parameter values of η = 1.0. In these cases, the saturation stress levels for both TTC and IPC simulations were shown to be nearly equal, at roughly σ saturation = 180 MPa. These results are summarized in Figure 10.  Figure 11a shows the evolution of the dislocation density of each slip as a function of IPC compression. Dislocation activity correlates well with the modal activity trends shown in Figure 7a, showing increasing dislocation density across strain sub-domains in modes with more modal activity. As shown in Figure 11b, the dislocation density evolution in IPC simulations is consistent with the TSF parameters of Oppedal et al. [14]. At the point in which twin volume fraction was fully realized (roughly 80% total volume fraction) the dislocation density of the twin volume fraction lies in the range of twice that of the the dislocation density of the parent volume fraction. The assumption that dislocations exist on the prismatic a and 2nd order pyramidal c + a was also put to the test. IPC compression was simulated under the assumption that only basal dislocations would be allowed to transmute across twin grain boundaries. This corresponds to dissociation input parameters of η basal = 0.5 and η prismatic = η 2 nd order pyramidal = 1.0. Shown in Figure 12, the stress levels after the onset of twinning were noticeably under predicted. While stress levels for simulations with only basal dislocations actively transmuting were closer to approximating experimental stress levels than simulations in which no transmutation was activated, the relative similarities can be accounted for by acknowledging that deformation by basal slip, and therefore basal type dislocation density, is significantly higher than both prismatic a and 2nd order pyramidal c + a in general. This is supported by observations of simulated dislocation density by mode, summarized in Figure 11a.

Conclusions
In summary, methods for simulating the effects of twin transmutation and dissociation at twin grain boundaries were implemented in a visco-plastic self-consistent framework in place of more conventional approaches that utilize Hall-Petch effects, and used to simulate uniaxial compression along multiple load paths in order to capture the behavior of magnesium under conditions which twin volume growth was both profuse and sparse. The results of these simulations were compared to both experimental data and simulated data from previous approaches. The results permit highlighting the following points:

1.
In place of Twin Storage Factor and Hall-Petch effects, a model for dislocation and twin boundary interactions was implemented. These methods and parameter selections were used to simulate the behavior of pure, basal textured, rolled magnesium subjected to uniaxial compression along (TTC) and perpendicular to (IPC) the dominant c -axis of the texture. The simulation results for stress, hardening, and texture development were consistent with observed experimental results. Modal activity and twin volume growth were largely similar to simulation work performed using the TSF approach.

2.
The large difference between the saturation stress of the simulated IPC and TTC load paths was confirmed to be the result of transmutation effects included in the model by disabling the contributions of transmuted dislocations to the dislocation density of twin volume fractions.

3.
It can be stated that Hall-Petch effects cannot be assumed to be the sole or even primary source of twinning induced anisotropy in the mechanical behavior of pure magnesium. Government. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: Hardening and dislocation generation parameters adopted from Oppedal et al. [14]: Table A1. Parameters for hardening evolution of slip.

Appendix B. Correspondence Method; Summary of Equations
Presented here is a summary of the general Correspondence Method equations from Niewczas [28], Niewczas [29] and Bilby and Crocker [23]. In this, a vector u M on in the parent lattice is mapped onto the vector v T in the twin vector using the matrices S, X, and U. These matrices are constructed from elements of m and n describing the twinning direction and normal to the twin plane.