**1. Introduction**

Mg alloys can significantly impact many industrial applications such as automobile, aerospace, and transportation, not to mention the energy sector, by effectively reducing the weight of the designed system. This can lead to a significant improvement in fuel economy and reducing emissions [1]. However, the current understanding of commercial Mg alloys and their physical metallurgy is less developed compared to other established structural metals and alloys. To overcome many of the challenges in using Mg alloys in these applications, different aspects of these alloys, including strength, formability, and fatigue resistance, should be enhanced. To fulfil any of these objectives, the underlying deformation mechanisms of Mg alloys should be thoroughly understood.

Two important deformation mechanisms of Mg alloys are plastic slip and extension twinning. Mg alloys have Hexagonal Closest Packed (HCP) crystal structure. In the case of slip, different modes can be activated, including basal a-{0001}1120, prismatic a-10101120, pyramidal a-10111120, and pyramidal c + a-11221123. In the case of unalloyed Mg polycrystals, the basal mode *a* has a very low critical resolved shear stress (CRSS) and becomes the dominant plasticity deformation mechanism in most of the loading conditions. In the case of Mg alloys, although the basal mode CRSS increases, it still becomes the lowest among the slip modes. Although the basal mode CRSS is very low, its corresponding resolved shear stress disappears during the tensile deformation along the c-axis. The extension twinning 10121011 becomes an important deformation mechanism of Mg alloys in this loading condition [2–9]. Similarly, the extension twinning is the controlling deformation mechanism during the compression loading parallel to the

**Citation:** Yaghoobi, M.; Voyiadjis, G.Z.; Sundararaghavan, V. Crystal Plasticity Simulation of Magnesium and Its Alloys: A Review of Recent Advances. *Crystals* **2021**, *11*, 435. https://doi.org/10.3390/ cryst11040435

Academic Editor: Sergio Brutti

Received: 17 March 2021 Accepted: 15 April 2021 Published: 17 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

basal plane. The highly asymmetric response of Mg alloys for some special orientations can be attributed to the activation of extension twinning, which has a complex morphology and micromechanics [3,10–16]. Although consideration of slip and twinning mechanisms is enough to describe the response of Mg alloys in most of the loading types, which are monotonic, this is not the case for the loading types, which involve unloading, such as cyclic loading. These loadings are especially important in many applications, such as fatigue. In the case of cyclic loading parallel to the basal plane, twinning occurs during compression loading, which results in the 86.3◦ reorientation of the basal pole. During subsequent reversed cyclic loading, i.e., tensile loading parallel to the basal plane, the size of the twinned regions shrinks, and in some cases, the twin fully disappears, a phenomenon that is commonly known as detwinning [7,17–22]. This twinning and detwinning occur alternately during the cyclic loading of Mg alloys.

In-situ experiments have a key role in unraveling the underlying deformation mechanisms of Mg alloys. Three in-situ experimental schemes, that were used to investigate the Mg alloy response were Digital Image Correlation (DIC), Synchrotron X-ray techniques, and neutron diffraction. In-situ DIC is another experimental technique that is used to address the deformation mechanisms of Mg alloys. Different aspects of Mg alloys have been studied by in-situ DIC, including the effect of twin nucleation on strain field [23], twin morphology [24], and the effect of aging on the accumulation of microscale plasticity [25]. Two common in-situ Synchrotron X-ray techniques of three-dimensional X-ray diffraction (3DXRD) and high-energy diffraction microscopy (HEDM) were incorporated in the investigation of Mg alloy behavior. Aydıner et al. [26] investigated the response of the AZ31 Mg alloy using the 3DXRD technique and analyzed the evolution of stress during deformation. It was observed that the stress state inside the twinned region was different from the untwinned grain. Other studies of Mg alloys have incorporated the in-situ 3DXRD and HEDM to analyze the micromechanics of twinning [27], the effect of age hardening on the deformation behavior [28], evaluate the CRSS for deformation modes [29,30], twin growth [31], and detwinning [7,22,32]. Using in-situ neutron diffraction, it was shown that the twinning in extruded Mg-7.7 at.% Al alloy was manifested as the change in the intensity of diffraction peaks during loading [33]. The in-situ neutron diffraction experiment on the Mg alloy showed that the stress relaxation occurred in the twinned region compared to the untwinned region of grain [34]. In-situ neutron experiments were also used to address the deformation detwinning in Mg alloys [35–39]. Although these in-situ experiments can provide invaluable information regarding the underlying deformation mechanisms of Mg alloys, they are very expensive, time-consuming, and complex. Potential simulation frameworks benefit the research community to complement these experimental observations.

Namakian and Voyiadjis [40] incorporated atomistic simulations and proposed a new mechanism for 10121011 twinning in which it was postulated that partial stacking faults (PSFs) play a key role in the formation of 1012 twins. A detailed crystallographic study was conducted for different HCP metals. Accordingly, PSFs were created by 13 1010 displacement vectors on every other basal plane and created a faulting plane on 1012 plane. Figure 1 shows the geometrical properties associated with the transformation of the untwinned parent into the twinned child unit cell. The model defines the properties of HCP crystals in a way that atoms can undergo spontaneous cooperative movements within these crystals. Accordingly, considerable stress relaxation occurs due to twin formation. This can be described using a deformation gradient, which shows that the macroscopic effect of the current mechanism is a simple shear process. Namakian et al. [41] incorporated atomistic simulations to model both deformation mechanisms of slip and twinning in single crystal Mg along with the interrelationships on the loose π\_1L and dense π\_1D first-order crystallographic planes. The generalized stacking fault energy (GSFE) analysis was used to investigate the slip mechanisms to study the c+a dislocation's core structure, dissociation mechanism, and mobility. Namakian et al. [41] stated that the screw component of a dissociated pyramidal-I a dislocation had a key role in compression twinning (CTW) nucleation. They also showed that CTW can grow by activating pyramidal-I a slip on the preexisting twin boundaries. Figure 2 shows the minimum energy path (MEP) of four deformation mechanisms associated with π\_1 plane, including homogeneous CTW nucleation, heterogeneous nucleation mechanism of pyramidal-I c+a slip, i.e., 1/31123 -1011 slip to CTW 1012 -1011 (or zonal twinning mechanism), heterogeneous nucleation mechanism of pyramidal- I a slip 1/31210 -1011 to CTW 1012 -<sup>1011</sup>*π*1\_*D*, and pyramidal-I c+a slip 1/31123 -1011 using Liu et al. [42] embedded-atom method (EAM) and Wu et al. [43] modified embedded-atom method (MEAM) interatomic potentials. Namakian et al. [41] showed the importance of the interrelationships of non-basal slips and CTW on π\_1 plane. They concluded that these interrelationships should be incorporated in simulation frameworks with larger length scales, including crystal plasticity, to accurately capture the mechanical response of Mg.

**Figure 1.** The schematics of the transformation of the untwinned parent unit cell (blue segments and red and blue atoms) into the twinned children unit cell (cyan segments and magenta and cyan atoms for the first two layers close to the twin boundary) during -10121011 twinning (After Namakian and Voyiadjis [40]).

The deformation mechanisms and response of metallic systems and alloys have been studied using frameworks at different length scales spanning from density functional theory (DFT) [44,45], atomistic simulations [11,40,41,46–56], discrete dislocation dynamics [57,58], and continuum mechanics [59]. However, all these schemes except continuum mechanics suffer from length and time scale gaps to model a real-size sample or experiment. Crystal plasticity (CP) is a very powerful continuum mechanics framework to model the response of Mg alloys with both slip and twinning deformations. Various CP models have been developed to simulate slip and deformation in HCP polycrystals [60–69]. A combination of CP and finite elements, which is commonly known as CP finite element (CPFE), benefits the advantages of both frameworks. Two main categories of CPFE frameworks were developed to capture both slip and twinning deformations, which were developed by Staroselsky and Anand [63] using rate-independent formulation and Kalidindi [62]

using rate-dependent formulation. In the case of Mg and its alloys, the rate-independent framework developed by Staroselsky and Anand [63] was used to model the slip along with twinning [66,68,69]. Kalidindi [62] introduced a modified plastic velocity gradient tensor by incorporating twinning. Various researchers have included twinning into their CPFE framework following the formulations introduced by Kalidindi [62], such as Abdolvand and his coworkers [65,70–72] and Zhang and Joshi [67]. Qiao et al. [73] included the stress relaxation due to twinning into the CPFE formulation, a CPFE framework that can capture the stress relaxation effect for twinning. Hama et al. [74] incorporated twinning into the CPFE framework to model the anisotropy of Mg alloy sheets subjected to a two-step loading. Prasad et al. [75] simulated the ductile fracture in pure Mg by adding the twinning to plane strain CPFE.

**Figure 2.** Minimum energy paths of four deformation mechanisms associated with π\_1 plane, i.e., homogeneous CTW nucleation, heterogeneous nucleation mechanism of pyramidal-I c+a slip, i.e., 1/31123 -1011 slip to CTW 1012 -<sup>1011</sup>, heterogeneous nucleation mechanism of pyramidal-I a slip 1/31210 -1011 to CTW 1012 \_(π1\_D), and pyramidal-I c+a slip 1/31123 -1011 (**a**) Liu et al. [42] EAM potential (**b**)Wu et al. [43] MEAM potential (After Namakian et al. [41]).

The CP models were used in combination with other schemes to model twinning. One of these candidates was a viscoplastic self-consistent (VPSC), which was first used

by Lebensohn and Tomé [76] to model the twinning that occurred in zirconium alloys. Agnew et al. [77] used the VPSC to model the slip and twinning in Mg alloys. Beyerlein and Tomé [12] simulated the twinning in pure zirconium using VPSC by introducing a probabilistic model for the twin nucleation. Kumar et al. [78] used a full-field elastoviscoplastic model based on fast Fourier transformation (FFT), which included twinning to model the response of Mg. <sup>L</sup>évesque et al. [79] introduced twinning into a Taylor-type CP model to simulate the response of AM30 and AZ31B Mg alloys. The multiscale framework of VPSC-FE was incorporated to simulate the slip and twinning deformations [80–82]. Accordingly, the response of each FE material point was defined using an RVE, which was modeled by VPSC. A similar multiscale scheme has been used by Ardeljan et al. [83] and Feather et al. [84] to simulate the response of Mg alloys in which an RVE modeled by a Taylor-type CP defines the constitutive model of an FE material point. A CP model has also been coupled with a phase field method to better describe the twin morphology in HCP metals [85,86].

Most of the studies have focused on the response of the Mg alloys during monotonic loading, which only requires the deformation slip and twinning. In the case of more complex strain paths such as cyclic loading, however, detwinning should also be defined in the CP model. A meso-scale composite grain (CG) model was the first detwinning model, which was included in the VPSC framework [87,88]. An empirical model defines the twin nucleation and propagation in the CG model. Guillemer et al. [89] used a simple phenomenological detwinning model along a self-consistent model to capture the cyclic behavior of extruded Mg. The first physically-based model, which includes both mechanisms, was the Twinning-Detwinning (TDT) model developed by Wang and his coworkers [90–92], which was used along the EVPSC framework. Qiao et al. [93] modeled the cyclic response of the ZK60A Mg alloy using the TDT/EVPSC scheme. The twinning and detwinning can be effectively handled by EVPSC framework by introducing new grains as twin nucleation and removing that twinned grain for complete detwinning, while the change in the volume of the twinned grain manifests the twin growth or shrinkage. However, in the case of studying the local variable information in a polycrystal, one cannot use the EVPSC model. One way to obtain the local information is using the CPFE framework to capture twinning and detwinning mechanisms. However, most of the CPFE models introduced to capture the twinning and detwinning mechanisms are phenomenological and neglect some of the main elements of the deformation twinning and detwinning. The biggest challenge in CPFE is that the introduction of twin nucleation, growth, shrinkage, and detwinning is not as straightforward as EVPSC. In the case of Mg and its alloys, the CPFE method has been used recently to model the cyclic behavior, including twinning and detwinning [22,94–97]. A material point is treated in these models as two states of twinned or not twinned, and the reorientation criteria are commonly based on the most active twinning system, which is called the predominant twinning reorientation (PTR) scheme [61]. Before the reorientation, the stress inside the twinned region does not contribute to the stress state of the material point, while after reorientation, the stress state inside the parent grain is neglected. However, the stress state of the twinned region and untwinned region are not similar, which was shown using in-situ experiments [26,34]. Yaghoobi et al. [9] introduced a multiscale CPFE framework in which the deformation twinning and detwinning were captured using a physically-based TDT model proposed by Wang and his coworkers [90–92]. They used a Taylor-type homogenization scheme to model both twinned and untwinned regions, coexisting at a material point.

Traditional crystal plasticity models were developed largely without a connection to grain size and shape effects. In Magnesium alloys, it has been observed that the grain size strengthening follows the Hall–Petch effect [98]. The incorporation of grain size effect into constitutive models for single slip began in 1962 with Armstrong et al. [99], who modified the Hall–Petch equation to correspond to the flow stress on a slip system (the "micro-Hall–Petch relation"). The interrelationship between grain size and texture was not considered until 1983, when Weng [100] employed the mean grain size in the equation for slip system resistance through the micro-Hall–Petch relation. Sun and Sundararaghavan [101] presented a crystal plasticity model for grain size and shape effects where the slip length of each slip system at a material point was considered in a micro-Hall–Petch strengthening term for each slip system. Recent HREBSD experiments [102] in Mg-4Al alloy have revealed the interaction of slip systems and grain boundaries is dependent on not only the slip length but also strongly on the angle between the incoming and emerging slip planes at a grain boundary. In the future, there is a need to better incorporate such insights in crystal plasticity.

In the current work, a review of works available in the literature, which addresses different properties of Mg and its alloys using crystal plasticity models, is presented. In addition to slip and twinning, detwinning is another deformation mechanism, which is activated in Mg and its alloys. The different models that capture detwinning will also be addressed here. Finally, the recent experimental frameworks, such as in-situ neutron diffraction, 3D high energy synchrotron X-ray techniques, and digital image correlation under scanning electron microscopy (SEM-DIC), are incorporated along crystal plasticity models to investigate the properties of Mg and its alloys are addressed. Future research directions towards improving the deformation response of Mg and its alloys are identified, which can lead to increased deployment of the lightest structural metal in engineering applications.

#### **2. Crystal Plasticity Models**

*2.1. Overview*

The CP theory is developed according to the assumption that the plasticity in metals and alloys occurs as a result of slip-on prescribed slip systems. The finite deformation continuum mechanics is usually incorporated to describe the formulation. The deformation gradient tensor **F** is decomposed as below:

$$\mathbf{F} = \mathbf{F}^{\mathbf{e}} \mathbf{F}^{\mathbf{p}} \tag{1}$$

where **F**e and **F**p are the elastic and plastic deformation gradients, respectively. Figure 3 shows the deformation described in Equation (1). Elastic distortion of the crystal lattice and plastic slip are two deformation mechanisms, which sustain the applied deformation in this formulation. Accordingly, the macroscopic velocity gradient **L** can be decomposed as below:

$$\mathbf{L} = \mathbf{L}^{\mathbf{e}} + \mathbf{L}^{\mathbf{P}} \tag{2}$$

where **L**e and **L**p are the elastic and plastic velocity gradients, respectively. The novel relation of CP is to link the macroscopic velocity gradient **L**p to micro deformation . *γ*<sup>α</sup> as below:

$$\mathbf{L}^{\mathbf{P}} = \sum\_{\alpha=1}^{N\_s} \dot{\gamma}^{\alpha} \mathbf{S}^{\alpha} \tag{3}$$

where . *γα* is the shearing rate on slip system *α*, *Ns* is the number of slip systems, and **S**α is the Schmid tensor for the slip system *α*, which can be defined as follows:

$$\mathbf{S}^{\alpha} = \mathfrak{m}^{\alpha} \bigotimes \mathfrak{n}^{\alpha} \tag{4}$$

where unit vectors **m**<sup>α</sup> and **n**<sup>α</sup> denote the slip direction and slip plane normal, respectively, in the deformed configuration.

**Figure 3.** Different configurations of the body in finite strain framework (After Yaghoobi et al. [103]).

The resolved shear stress on the slip system *α* can be obtained as follows:

$$
\boldsymbol{\tau}^{\alpha} = \boldsymbol{\sigma} \cdot \mathbf{S}^{\alpha} \tag{5}
$$

where σ is the Cauchy stress tensor and · operator denotes the standard inner product of tensors. The shearing rate on the *α*th slip system *γ<sup>α</sup>* can be calculated depending on if the CP formulation is rate-independent or rate-dependent. In the case of a rate-independent model, the yield surface is defined for each system, i.e., *f α*, as follows:

$$f^{a} = \lfloor \tau^{a} \rfloor - s^{a} \tag{6}$$

where *s<sup>α</sup>* is the slip resistance for slip system α. The values of .*γα* is then obtained using the yield surface and considering the hardening model, which describes the evolution of slip resistance.

In the case of rate-dependent formulation, the evolution of shearing rate on the *α*th slip system *γ<sup>α</sup>* can be defined as follows:

$$\dot{\gamma}^a = \dot{\gamma}\_0 \left| \frac{\tau^a}{s^a} \right|^{1/m} \text{sign}[\tau^a] \tag{7}$$

where .*γ*0 and *m* are the material parameters denoting the reference shearing rate and rate sensitivity of the material. One should note that the kinematic hardening is neglected in both Equations (6) and (7).
