Numerical Investigation of Secondary Deformation Mechanisms on Plastic Deformation of AZ31 Magnesium Alloy Using Viscoplastic Self-Consistent Model

Uniaxial tension and compression of AZ31 magnesium alloy were numerically investigated via the viscoplastic self-consistent (VPSC) model to shed a light on the effect of secondary deformation mechanisms (prismatic  slip, pyramidal  slip, and { 10 1 ¯ 1 } contraction twinning) during plastic deformation. The method adopted in the present study used different combinations of deformation mechanisms in the VPSC modeling. In terms of the pyramidal  slip, it served as the first candidate for sustaining the extra plastic strain during the plastic deformation. The improvement of activity in the pyramidal  slip contributed to the increase in the mechanical response and the splitting of pole densities in { 0002 } pole figure during uniaxial tension. As for the prismatic  slip, its increasing activity was not only conducive to the improvement of flow stress in mechanical response, but also responsible for the splitting of pole densities in { 0002 } pole figure during uniaxial compression. With respect to the { 10 1 ¯ 1 } contraction twinning, it had a negligible influence on the plastic deformation of AZ31 magnesium alloy in terms of the mechanical response as well as the slip and the twinning activities. However, it is better to include the { 10 1 ¯ 1 } contraction twinning in the VPSC modeling to more accurately predict the texture evolution.


Introduction
The exploitation of magnesium and its alloys for an increasing diversity of applications in the fields of materials science and industry is derived from their remarkable attractiveness including low density, high specific stiffness and strength, good damping property, etc. [1,2]. These characteristics make magnesium and its alloys one of the light-weight materials with a huge potential in tackling the current energy crisis [3]. Among all magnesium alloys, the AZ31 magnesium alloy would be one of the most widely applied. However, owing to its anisotropy in mechanical response and poor formability at room temperature, many challenges with respect to metal forming and structural application still remain to be solved. The anisotropic plastic deformation and poor formability of the AZ31 magnesium alloy at room temperature resulted from the activated deformation mechanisms during the plastic deformation. In general, it is accepted that the ductility of metallic materials is mainly associated with the number of activated slip/twinning systems during the plastic deformation. As for the AZ31 magnesium alloy with the hexagonal close-packed (HCP) crystal structure, the known slip/twinning systems include {0001} < 1120 > (basal <a>) slip, 1010 < 1120 > (prismatic <a>) slip, 1122 < 1123 > (pyramidal <c+a>) slip, 1012 extension twinning, and 1011 contraction twinning [4]. Among these deformation mechanisms, the basal <a> slip and the 1012 extension twinning had been confirmed to be the main deformation mechanisms at room temperature to accommodate the plastic strain along and perpendicular to the crystal c-axis, respectively, as they possess a relatively lower critical resolved shear stress (CRSS) in comparison with the other deformation mechanisms [4,5]. However, the remaining deformation mechanisms including the non-basal slip (prismatic <a> slip and pyramidal <c+a>) slip) as well as the 1011 contraction twinning also play important roles in the plastic deformation of AZ31 magnesium alloy at room temperature as determined through the experimental and modeling investigation [6,7]. In the present study, due to their contributions to plastic deformation, these three deformation mechanisms are termed as secondary deformation mechanisms. In addition, our literature search has shown that there seems to be a lack of study focusing on the roles of these secondary deformation mechanisms during the plastic deformation of AZ31 magnesium alloy.
To date, there is still a lack of efficient and convenient experimental observation techniques in revealing the plastic deformation mechanisms of metallic materials. Consequently, the numerical investigation on the basis of crystal plasticity theory has become a perfect candidate for clarifying the influence of deformation mechanisms during plastic deformation. As for the AZ31 magnesium alloy, researchers have developed numerous models in consideration of the dominant deformation mechanisms to predict the anisotropic mechanical response and texture evolution [8][9][10][11][12][13]. Among these models, the viscoplastic self-consistent (VPSC) model can be termed as one of the most adopted models due to its efficiency and accuracy in modeling the plastic deformation of metallic materials. Therefore, in the present study, the VPSC model is used to shed a light on the roles of these secondary deformation mechanisms of AZ31 magnesium alloy during the plastic deformation, which could lay foundations for extending its application at room temperature. The method adopted in the present study is to systematically change the embedded deformation mechanisms during the VPSC modeling; however, the input data in terms of initial orientation and boundary conditions remain unchanged.

VPSC-Based Crystal Plasticity Modeling
In the present study, the viscoplastic self-consistent (VPSC) model developed by Lebensohn et al. [14] is utilized to predict the macroscopic stress-strain response, slip and twinning activities, and texture evolution of the AZ31 magnesium alloy during uniaxial tension and compression. The core of VPSC model is that it can account for the full anisotropy of properties and responses of the single crystals by treating the constituting individual grains as anisotropic viscoplastic ellipsoidal inclusions [15] inside a homogeneous effective medium (HEM), which represents the average constitutive behavior of the polycrystalline aggregate. The interaction between each grain and HEM can be described by different linearization assumptions embedded in the VPSC model [16]. To retain self-consistency and to ensure strain compatibility as well as stress equilibrium during the numerical modeling, the condition that the average strain rate of individual grains has to be consistent with the applied macroscopic magnitude is enforced in the VPSC model [17]. As the detailed procedure for calculation in the VPSC model has been well documented by Li et al. [18], only the key equations about the hardening model and twinning model are briefly presented here.
The VPSC model is based on the physical shear mechanisms of slip and twinning during the plastic deformation. The activation of the physical shear mechanisms is mainly associated with the corresponding CRSS, whose evolution is totally described by the hardening model. To date, several hardening models have been successfully developed and applied in the VPSC model including a modified Voce model, the mechanical threshold stress (MTS) model, and some dislocation density-based models [18][19][20]. In the present study, though it is in fact empirical, a modified Voce model is applied because of its simplification during parameter calibration and its accuracy in terms of predicting the mechanical response and the texture evolution during the plastic deformation [21,22]. The modified Voce model relates the evolution of CRSS with the accumulative shear strain of all the activated deformation mechanisms including the slip and the twinning systems.
where τ α c represents the instantaneous CRSS of each slip/twinning system α. The expressions τ α 0 and θ α 0 stand for initial CRSS and initial hardening rate, respectively, whereas τ α 1 is associated with the back-extrapolated CRSS and θ α 1 refers to the asymptotic hardening rate. The expression Γ stands for the total accumulated shear strain in all slip and twinning systems.
It is generally accepted that there exist interactions among the slip and the twinning systems, hence the propagation in a slip or a twinning system would be impeded by other slip and twinning systems. This phenomenon has been taken into consideration in the VPSC model by introducing the coupling coefficient h αα .
where ∆τ α c refers to the increment of CRSS in a slip or a twinning system, and it mainly depends on the increment of shear strain in each slip and twinning systems. The expression h αα in the case of (α = α ) represents the self-hardening coefficient, whereas h αα in the case of (α = α ) is known as the latent-hardening coefficient.
The embedding of twinning in modeling the plastic deformation of metallic materials has been successfully obtained by applying the predominant twin reorientation (PTR) scheme proposed by Tomé et al. [23]. In the PTR approach, both the shear strain γ α and the corresponding twin volume fraction f α = γ α /γ tw associated with each twin system α within a single grain are tracked at each incremental step (γ tw refers to the characteristic twinning shear). Since treating each twinned fraction as a newly generated grain is not numerically feasible in the PTR scheme a given grain is entirely reoriented only if the predominant twinning system exceeds a threshold value V thres as follows.
where A th1 and A th2 are the two parameters associated with the evolution of twin volume fraction during the plastic deformation, and they could have different values with respect to various twinning modes. In the present study, for the 1012 extension twinning, A th1 = 0.7 and A th2 = 0.1, whereas for the 1011 contraction twinning, A th1 = 0.05 and A th2 = 0.05 [24]. Here, V e f f refers to the volume fraction of the twin-terminated grains and V accum represents the volume fraction of the twinned regions.

Results and Discussion
In the present study, a hot-rolled AZ31 magnesium alloy plate with a strong basal texture is considered as shown in Figure 1, where the initial {0002} pole figure demonstrates that the c-axis of most grains in the sheet is close to the normal direction (ND). The specific contents about fabricating this AZ31 magnesium alloy sheet and the procedure of measuring its initial texture have been well documented in the literature [24]; hence they are not discussed here. The input data about crystallographic orientations that are needed in the VPSC modeling are then extracted from the measured texture. By using the method of discretizing the orientation distribution function (ODF) via the toolbox MTEX in MATLAB code, which is depicted in the literature [25], 1600 discrete orientations with equal volume fractions are finally used for the VPSC modeling. As for the boundary conditions, which is also key input information in the VPSC modeling, a combination of velocity gradient components and stress components are imposed in the present study to predict the uniaxial tension and compression of the AZ31 magnesium alloy in different loading paths. This operation allows for the ovalization of the deformed sample [26]. In all simulations, a true strain of 0.15 is enforced and a strain rate 0.001 s −1 is applied in each loading direction, which is equal to the macroscopic strain rate in uniaxial tension and compression experiments as depicted in the literature [24]. The mainly reported deformation mechanisms involving the plastic deformation of AZ31 magnesium alloy at room temperature are {0001} < 1120 > (basal <a>) slip, 1010 < 1120 > (prismatic <a>) slip, 1122 < 1123 > (pyramidal <c+a>) slip, 1012 extension twinning, and 1011 contraction twinning. These deformation mechanisms are considered in VPSC simulations. Moreover, 1011 − 1012 double twinning is also allowed in the present study. In addition, the affine linearization scheme is applied in the present study to relate the responses of single grains to the response of the polycrystal aggregate. Wang et al. [16] have confirmed that the affine scheme gives the best accordance between the experimental results and the corresponding simulated ones for the AZ31 magnesium alloy. The exponent of rate sensitivity n = 20 is applied in the present study, which is a recommended value in the VPSC manual to simulate an almost rate-insensitive deformation. Figure 1 shows the mechanical responses during through-thickness compression (TTC) along the normal direction (ND), in-plane compression (IPC) along the rolling direction (RD), and in-plane tension (IPT) along the transverse direction (TD). These experimental data are from the work of Chapuis et al. [24]. Using the corresponding material parameters identified by Chapuis et al. [24], as shown in Table 1, the predicted stress-strain curves are obtained and are also shown in Figure 1. Though during IPC along RD, the flow stresses between the true strain of 0.04 and 0.07 exhibit some discrepancy, the consistency between the global mechanical responses and the corresponding simulated ones confirms the validity of the adopted materials.
orientations with equal volume fractions are finally used for the VPSC modeling. As for the boundary conditions, which is also key input information in the VPSC modeling, a combination of velocity gradient components and stress components are imposed in the present study to predict the uniaxial tension and compression of the AZ31 magnesium alloy in different loading paths. This operation allows for the ovalization of the deformed sample [26]. In all simulations, a true strain of 0.15 is enforced and a strain rate 0.001 s −1 is applied in each loading direction, which is equal to the macroscopic strain rate in uniaxial tension and compression experiments as depicted in the literature [24]. The mainly reported deformation mechanisms involving the plastic deformation of AZ31 magnesium alloy at room temperature are {0001}<1120> (basal <a>) slip, {10 10}<1120> (prismatic <a>) slip, {1122}<1123> (pyramidal <c+a>) slip, {10 12} extension twinning, and {10 11} contraction twinning. These deformation mechanisms are considered in VPSC simulations.
Moreover, {10 11}-{10 12} double twinning is also allowed in the present study. In addition, the affine linearization scheme is applied in the present study to relate the responses of single grains to the response of the polycrystal aggregate. Wang et al. [16] have confirmed that the affine scheme gives the best accordance between the experimental results and the corresponding simulated ones for the AZ31 magnesium alloy. The exponent of rate sensitivity 20  n is applied in the present study, which is a recommended value in the VPSC manual to simulate an almost rate-insensitive deformation. Figure 1 shows the mechanical responses during through-thickness compression (TTC) along the normal direction (ND), in-plane compression (IPC) along the rolling direction (RD), and in-plane tension (IPT) along the transverse direction (TD). These experimental data are from the work of Chapuis et al. [24]. Using the corresponding material parameters identified by Chapuis et al. [24], as shown in Table 1, the predicted stress-strain curves are obtained and are also shown in Figure 1. Though during IPC along RD, the flow stresses between the true strain of 0.04 and 0.07 exhibit some discrepancy, the consistency between the global mechanical responses and the corresponding simulated ones confirms the validity of the adopted materials.     [5,27,28]. The focus of the present study is to numerically investigate the effects of these secondary deformation mechanisms during the plastic deformation of AZ31 magnesium alloy, hence four different combinations of deformation mechanisms are constructed and then applied in all subsequent VPSC simulations, as shown in Table 2. The first combination of deformation mechanisms is termed as combination A and it consists of basal <a> slip, prismatic <a> slip, pyramidal <c+a> slip, 1012 extension twinning, and 1011 contraction twinning. It serves as a benchmark in VPSC modeling. The second combination of deformation mechanisms is termed as combination B, which contains basal <a> slip, pyramidal <c+a> slip, 1012 extension twinning, and 1011 contraction twinning. The third combination of deformation mechanisms is termed as combination C, which includes basal <a> slip, prismatic <a> slip, 1012 extension twinning, and 1011 contraction twinning. The fourth combination of deformation mechanisms is termed as combination D and it is made up of basal <a> slip, prismatic <a> slip, pyramidal <c+a> slip, and 1012 extension twinning.  Figure 2 demonstrates the predicted stress-strain curves during IPC along RD and the corresponding slip and twinning activities in the cases of different combinations of deformation mechanisms. As shown in Figure 2a, there exists little discrepancy between the stress-strain curve in the case of combination A and the corresponding one in the case of combination D. On the contrary, both the stress-strain curves in the cases of combination B and combination C are higher than the corresponding one in the case of combination A. It can be seen from Figure 2b that when the prismatic <a> slip is eliminated from the adopted deformation mechanisms, there exist almost no changes with respect to the activities in the 1012 extension twinning and the 1011 contraction twinning, and only a minor change appears in the activity of basal <a> slip. However, the evolution of activity in the pyramidal <c+a> slip changes obviously and the increasing activity of the pyramidal <c+a> slip in the case of combination B is nearly equal to the total activity of the prismatic <a> slip in the case of combination A. In addition, it is important to note that the actual activity of the pyramidal <c+a> slip in the case of combination A and combination B are quite close to each other at the true stain of 0.15. As for the pyramidal <c+a> slip, when it is excluded during IPC along RD, as shown in Figure 2c, the evolution laws in the basal <a> slip and the prismatic <a> slip as well as in the 1012 extension twinning and the 1011 contraction twinning change obviously. Moreover, the 1011 contraction twinning possesses the biggest increase in the case of combination C. With respect to the 1011 contraction twinning, when it is not embedded in VPSC modeling, there seems to be little influence on the activities of other deformation mechanisms as shown in Figure 2d. Based on the aforementioned analysis with respect to the slip and twinning activities in VPSC modeling with different combinations of deformation mechanisms, it can be noted that the aforementioned discrepancy with respect to the stress-strain curves during IPC along RD is actually derived from the various roles of each deformation mechanism during the plastic deformation of AZ31 magnesium alloy, which has been reported by Chapuis et al. [24]. As for the pyramidal <c+a> slip, when it is excluded during IPC along RD, as shown in Figure 2c, the evolution laws in the basal <a> slip and the prismatic <a> slip as well as in the {10 12} extension twinning and the {10 11} contraction twinning change obviously. Moreover, the {10 11} contraction twinning possesses the biggest increase in the case of combination C. With respect to the {10 11} contraction twinning, when it is not embedded in VPSC modeling, there seems to be little influence on the activities of other deformation mechanisms as shown in Figure 2d. Based on the aforementioned analysis with respect to the slip and twinning activities in VPSC modeling with different combinations of deformation mechanisms, it can be noted that the aforementioned discrepancy with respect to the stress-strain curves during IPC along RD is actually derived from the various roles of each deformation mechanism during the plastic deformation of AZ31 magnesium alloy, which has been reported by Chapuis et al. [24].  Figure 3 demonstrates the mechanical responses as well as the slip and twinning activities during IPT along TD. Similar to the evolution laws of stress-strain curves as shown in Figure 2a, the mechanical responses in Figure 3a possess larger flow stresses in the cases of combination B and combination C than in the cases of combination A and combination D. Moreover, in the case of combination B, the growth with respect to the flow stress is bigger than the corresponding one in the case of combination C. The difference in mechanical responses can also be related to the slip and twinning activities. As shown in Figure 3b, in the case of combination A, the basal <a> slip and the prismatic <a> slip possess the dominant activities during IPT along TD. When the prismatic <a> slip is not considered in the VPSC modeling, namely in the case of combination B, the increment in the activity of pyramidal <c+a> slip is the highest in comparison with the corresponding increment with respect to the basal <a> slip. As for the pyramidal <c+a> slip, when it is not embedded in the adopted deformation mechanisms during the VPSC modeling, it only results in a slight increase of activity in Figure 2. (a) The predicted stress-strain curves of AZ31 magnesium alloy under IPC along RD; (b-d) the corresponding slip and twinning activities predicted by the VPSC modeling in the cases of combination B, combination C, and combination D, respectively. Figure 3 demonstrates the mechanical responses as well as the slip and twinning activities during IPT along TD. Similar to the evolution laws of stress-strain curves as shown in Figure 2a, the mechanical responses in Figure 3a possess larger flow stresses in the cases of combination B and combination C than in the cases of combination A and combination D. Moreover, in the case of combination B, the growth with respect to the flow stress is bigger than the corresponding one in the case of combination C. The difference in mechanical responses can also be related to the slip and twinning activities. As shown in Figure 3b, in the case of combination A, the basal <a> slip and the prismatic <a> slip possess the dominant activities during IPT along TD. When the prismatic <a> slip is not considered in the VPSC modeling, namely in the case of combination B, the increment in the activity of pyramidal <c+a> slip is the highest in comparison with the corresponding increment with respect to the basal <a> slip. As for the pyramidal <c+a> slip, when it is not embedded in the adopted deformation mechanisms during the VPSC modeling, it only results in a slight increase of activity in the prismatic <a> slip, as shown in Figure 3c. With respect to the 1011 contraction twinning, Figure 3d illustrates that whether it exists or not, the activities in other deformation mechanisms would not be influenced during IPT along TD. Hence, it has a negligible influence on the mechanical response of AZ31 magnesium alloy during IPT along TD. the prismatic <a> slip, as shown in Figure 3c. With respect to the {10 11} contraction twinning, Figure 3d illustrates that whether it exists or not, the activities in other deformation mechanisms would not be influenced during IPT along TD. Hence, it has a negligible influence on the mechanical response of AZ31 magnesium alloy during IPT along TD. It is generally accepted that the slip and twinning activities have a direct influence on the texture evolution of AZ31 magnesium alloy during the plastic deformation; hence, the predicted textures during IPC along RD and IPT along TD are also analyzed in the present study using VPSC modeling with different combinations of deformation mechanisms, as shown in Figure 4. As for the IPC along RD, Figure 4a The aforementioned difference with respect to the predicted texture can also be related to the slip and twinning activities during the plastic deformation. Based on the aforementioned analysis on the slip and twinning activities during IPC along RD and IPT along TD, an agreement can be made that the pyramidal <c+a> slip is the first candidate to sustain the extra plastic strain during plastic deformation, which is consistent with the conclusion made by Styczynski et al. [28] that the pyramidal <c+a> slip is a "compatibility first" deformation mechanism of the AZ31 magnesium alloy. As for the case of combination B, the activity of pyramidal <c+a> slip increases obviously during IPT along TD and further results in the splitting of pole densities around ND in the {0002} pole figure. However, It is generally accepted that the slip and twinning activities have a direct influence on the texture evolution of AZ31 magnesium alloy during the plastic deformation; hence, the predicted textures during IPC along RD and IPT along TD are also analyzed in the present study using VPSC modeling with different combinations of deformation mechanisms, as shown in Figure 4. As for the IPC along RD, Figure 4a,c,g demonstrate that the distribution of pole densities in {0002} and 1010 pole figures are nearly the same in the cases of combination A, combination B, and combination D. However, in the case of combination C, the splitting of pole densities is obvious around RD, as shown in Figure 4e, and there also exist some pole densities close to ND. In contrast, during IPT along TD, the {0002} and 1010 pole figures exhibit only minor changes in the case of combination A, combination C, and combination D, as shown in Figure 4b,f,h, respectively. However, in the case of combination B, there exists an obvious splitting of pole densities around ND in the {0002} pole figure.
The aforementioned difference with respect to the predicted texture can also be related to the slip and twinning activities during the plastic deformation. Based on the aforementioned analysis on the slip and twinning activities during IPC along RD and IPT along TD, an agreement can be made that the pyramidal <c+a> slip is the first candidate to sustain the extra plastic strain during plastic deformation, which is consistent with the conclusion made by Styczynski et al. [28] that the pyramidal <c+a> slip is a "compatibility first" deformation mechanism of the AZ31 magnesium alloy. As for the case of combination B, the activity of pyramidal <c+a> slip increases obviously during IPT along TD and further results in the splitting of pole densities around ND in the {0002} pole figure. However, during IPC along RD, the increment with respect to the activity of pyramidal <c+a> slip is small, hence its influence on the predicted texture is very small. In the case of combination C, during IPC along RD, the appearance of pole densities close to ND can be attributed to the increase of activity in the 1011 contraction twinning, which mainly results in the 56 • rotation about c-axis in the crystals of AZ31 magnesium alloy. Moreover, the splitting of pole densities around RD can be due to the increasing activity in the prismatic <a> slip, as Chapuis et al. [24] have found that the prismatic <a> slip can generate a rotation around its c-axis in HCP metals. In terms of combination D, it can be concluded that whether it exists or not, the 1011 contraction twinning would not affect the predicted textures during IPC along RD and IPT along TD.
during IPC along RD, the increment with respect to the activity of pyramidal <c+a> slip is small, hence its influence on the predicted texture is very small. In the case of combination C, during IPC along RD, the appearance of pole densities close to ND can be attributed to the increase of activity in the {10 11} contraction twinning, which mainly results in the 56° rotation about c-axis in the crystals of AZ31 magnesium alloy. Moreover, the splitting of pole densities around RD can be due to the increasing activity in the prismatic <a> slip, as Chapuis et al. [24] have found that the prismatic <a> slip can generate a rotation around its c-axis in HCP metals. In terms of combination D, it can be concluded that whether it exists or not, the {10 11} contraction twinning would not affect the predicted textures during IPC along RD and IPT along TD.

Conclusions
The effect of secondary deformation mechanisms (prismatic <a> slip, pyramidal <c+a> slip, and 1011 contraction twinning) on the plastic deformation of AZ31 magnesium alloy during in-plane compression (IPC) along the rolling direction (RD) and in-plane tension (IPT) along the transverse direction (TD) had been numerically investigated by VPSC modeling with different combinations of deformation mechanisms. The following conclusions can be drawn.
(1) As for the pyramidal <c+a> slip, it is the first candidate for sustaining the extra plastic strain during the plastic deformation. The increment in the activity of the pyramidal <c+a> slip contributes to the increase of flow stress in mechanical response during IPC along RD and IPT along TD. Moreover, during IPT along TD, its great growth contributes to the splitting of pole densities around ND in the {0002} pole figure.
(2) The increment in the activity of prismatic <a> slip is also conducive to the improvement of flow stress in mechanical response during IPT along TD. In addition, during IPC along RD, its great growth contributes to the splitting of pole densities around RD in the {0002} pole figure.
(3) The 1011 contraction twinning has a minor influence on the plastic deformation of AZ31 magnesium alloy during IPC along RD and IPT along TD as well as the corresponding mechanical responses. However, during IPC along RD, its great growth contributes to the appearance of pole densities close to ND. Consequently, it is better to include the 1011 contraction twinning in VPSC modeling to predict the texture more accurately during the plastic deformation of AZ31 magnesium alloy.