A new constitutive model to describe evolving elastoplastic behaviours of hexagonal close-packed sheet metals

This study develops a new phenomenological constitutive model to capture the unique evolving cyclic elastoplastic behaviours of hexagonal close-packed (HCP) sheet metals. This new constitutive model is developed by adopting the concepts of multiple-yield surface approaches. Four phenomenological deformation modes, including Monotonic Compression (MC), Monotonic Tension (MT), Reverse Compression (RC), and Reverse Tension (RT), are considered to represent the hardening evolution of the materials, including the twining/untwining behaviours. Reference flow stress equations are introduced, and a Cazacu-Barlat 2004 (CB2004) type yield surface is employed to each deformation mode. In addition, the RT hardening parameters are defined as functions of plastic pre-strains to mitigate the interpolation error caused by parameter determination processes of existing models. For validation, the calculated stress–strain curves of AZ31B magnesium alloy are compared with experimental curves available from literature. Moreover, to show the accuracy of the proposed analytical model, the reproduced stress–strain curves are compared with those of an existing model—the modified homogeneous anisotropic hardening (HAH) model. The obtained results show that the new constitutive model can successfully reproduce experimental Tension–Compression-Tension (TCT) and Compression-Tension–Compression (CTC) stress–strain curves of HCP sheet metals with considerably less percentage errors.


Introduction
Over the last decade, attention on microscopic and macroscopic behaviours of HCP metals has been paid increasingly due to their unique features and industrial application growths. These metals have a wide range of applications in industry, including aerospace, automotive, electronics, etc. For example, magnesium alloys are the lightest structural components (about 1.71 g/cm 3 ) with the HCP crystal structure and become increasingly crucial for automotive industry [1]. However, in sheet metal forming, due to their crystalline structure, the number of slip systems is limited, and consequently, they exhibit poor formability [2].
Magnesium alloys show unique mechanical behaviours compared to conventional cubic crystalline metals such as aluminium alloys and steels [3]. Magnesium alloys' material properties and mechanical behaviours can differ based on the loading path and initial texture. For example, Styczynski et al. [4] showed a strong basal texture can be induced to wrought magnesium alloys in a rolling process. For most HCP metals, these unique material behaviours include the tension-compression asymmetry or eccentricity in the initial yield points, which is usually called "strength differential (SD)" or "yield asymmetry." Härtel et al. [5] reported that the compressive strength can be treated to be much lower than the tensile strength. Moreover, Lee et al. [3] study showed that the asymmetry can be observed in the subsequent plastic hardening, which is referred to as a "hardening differential" or "flow asymmetry." Hadadzadeh et al. [6] showed that the strong basal or nearly basal texture of coldrolled HCP metals is the main reason for the poor ductility and the strong tension-compression asymmetry.
Due to the basal crystallographic texture of HCP metals, main deformation modes are restricted to either prismatic {101̅ 0} < 112̅ 0 > or basal {0001} < 112̅ 0 > slips along the basal direction < a > in-plane tension. Because the < a > type slips cannot cause any deformations along the c-axis, the tensile twining {101̅ 2} < 101̅ 1 > is activated and results in an extension along the c-axis [7]. Tensile twinning is the predominant deformation mechanism during the in-plane compression, which causes a strong tension-compression asymmetry in the plastic behaviour of the materials [8]. Although the critical resolved shear stresses (CRSS) for slip on the basal plane are remarkably lower than for the prismatic and pyramidal planes, basal slip cannot shorten or elongate the c-axis . Consequently, it cannot satisfy the von-Mises criterion of five independent slip systems for homogeneous deformations [9].
Previous studies were carried out to model the unusual mechanical behaviours of magnesium alloys at room temperature. Usually, the constitutive relation was used to describe the plastic flow behaviours of materials. Generally, the constitutive models can be classified into two categories, including phenomenological and physically based models. Constitutive models based on crystal plasticity are classified as physically based models. In crystal plasticity, the constitutive equations are evaluated based on the deformation mechanisms of crystals, such as the slip mechanism of HCP metals. Because crystal plasticity is a theory for materials at a grain scale, a tremendous amount of calculation is necessary, and consequently, the practical application of this model is for FE-based numerical modelling. However, a large number of calculations at grain-scale of materials can be time-consuming and computationally insufficient.
Usually, the continuum plasticity is applied to phenomenological models. In the phenomenological continuum plasticity, the yield surface completes the plastic response of the material, which is the most critical constitutive equation. HCP metals, including magnesium alloys, have asymmetrical yield surfaces due to their strong eccentricity in yielding [8]. Several yield functions were proposed to model the anisotropy or the tension-compression asymmetry (TCA) of the materials. Cubic materials such as body-centred cubic (BCC) and face-centred cubic (FCC) metals show considerable anisotropy, while most do not exhibit tensioncompression asymmetry. Several studies were performed on the cyclic behaviours of such metals. Shariati and Mehrabi [10] studied the low cycle fatigue life of CK45 steel and SS316 stainless steel under strain-controlled loading. Moreover, they experimentally investigated the impact of mean strain and strain amplitude on fatigue life. Shariati and Mehrabi [11] performed an experimental study on ratcheting-fatigue interaction and ratcheting behaviour of Ck45 steels under stress-controlled uniaxial cyclic loading. Also, Shariati and Mehrabi [12] studied the influence of pre-fatigue on the fatigue behaviour of CK45 steel at room temperature. They recorded the fatigue life for various stress ratios and applied different mean stress models to predict the fatigue life. However, some cubic metals show considerable TCA in their plastic behaviours, which makes their plastic behaviour even more complicated. For example, Maeda et al. [13] showed that the DP980 steel sheets show considerable strength differential effect (SDE) under the in-plane uniaxial tension-compression test. Recently, Mehrabi et al. [14] proposed new analytical and numerical methods to model the tension-compression asymmetry in plastic behaviours of steels under bending. There are few continuum-based models reported in literature to consider the anisotropy and asymmetry of HCP metals. One of the noticeable and wellknown models is the Cazacu-Barlat family of yield surfaces developed since 2001, which can precisely model HCP metals' asymmetry and anisotropy. Cazacu and Barlat [15] proposed a generalisation of the invariants of the stress deviator and used this approach to extend Drucker's isotropic yield criterion. Cazacu and Barlat [16] applied this model to study the anisotropic plastic response of aluminium alloys. Moreover, Cazacu and Barlat [17] modified the proposed yield function to consider the tension-compression asymmetry. They modified the isotropic Drucker yield function [18] by applying the second and third invariants of stresses into the yield function. Mehrabi and Yang [19] developed a novel analytical method based on the Cazacu-Barlat 2004 asymmetric yield function and isotropic plastic hardening rule for modelling bending and springback behaviours of HCP metals. Furthermore, Mehrabi et al. [20] improved the Mehrabi and Yang [19] model by including the effects of the neutral surface shit and the variable elastic modulus into the model.
For the stress invariants-based yield criteria, not only the third stress invariant contributes the SD behaviours, but the first stress invariant can also play an essential role in such material behaviours. Cazacu et al. [21] developed a physically-based non-quadratic yield surface based on single and polycrystalline viscoplasticity to describe the material's anisotropy and initial yield asymmetry. In their research, the anisotropy was considered with the linear transformation of the deviatoric stress tensor, and the asymmetry parameter is controlled by the ratio of strength in tension and compression.
While the twining occurs under compression, the untwining happens under the subsequent tension. Lou et al. [22] showed that it initiates because of the contraction of the twined region without nucleation. Accurate modelling of twining/untwining behaviours is important since reverse loading usually occurs in sheet metal forming processes. Some phenomenological and physically-based models based on the continuum plasticity and crystal plasticity were presented to model magnesium alloys' twining/untwining behaviours. Although the proposed models can capture the twining/untwining behaviours of the material, further study is needed to improve the application and accuracy of these models. The existing models based on crystal plasticity such as Zhang et al. [23] have some complexities, which limit their applications. In addition, the available phenomenological models have some drawbacks. Li et al. [24] proposed a model based on von-Mises yield function and concepts of crystallographic models such as c-axis orientation of basal texture for each phenomenological deformation mode, which has the complexity of crystal plasticity concepts. Some models were developed based on applying the multiple-yield surface method, for example, the Lee et al. [8] model. However, since the material parameters in their models are related to pre-strains, the lack of independent functions to define these parameters as a function of compressive pre-strains during untwining is troublesome. In some cases, interpolating data from multiple experimental tests are required to determine material parameters for different pre-strains, which can cause interpolation errors. A recent model based on the HAH hardening model was presented by Lee et al. [3] for AZ31B magnesium alloy sheets at room temperature. They modified the HAH hardening model to capture the plastic behaviours of magnesium alloys under non-proportional strain path changes with constant material parameters. Overall, their model can reproduce the main characteristics of stress-strain curves. However, the reproducibility of their model was only shown acceptable results for a limited range of pre-strains and their model provided a limited accuracy as well. Recently, Lee et al. [25] modified the Lee et al. [3] model for AZ31B magnesium alloy sheets at elevated temperatures, which inherits the same drawback as their previous model. Therefore, there is still a technical need to develop a new constitutive model for HCP metals, capturing their complex behaviours precisely under dynamic loading conditions.
In this research, a new elastoplastic constitutive model is developed to capture the twining/untwining behaviours of HCP sheet metals during TCT and CTC cyclic loading. The developed model is based on the concepts of multiple yield surface approaches and isotropic hardening rule. The hardening response consists of four phenomenological modes, including monotonic compression (MC), monotonic tension (MT), reverse compression (RC), and reverse tension (RT), and a Cazacu-Barlat 2004 (CB2004) type yield surface is adopted for each phenomenological mode. Furthermore, a new parameter determination process is developed as the RT mode's hardening material parameters are defined as compressive pre-strains functions. Not only this method can make the parameter determination process convenient and straightforward, but it can mitigate the interpolation error of existing models. As a case study to validate the proposed analytical model, material behaviours of AZ31B magnesium alloy under tension-compression-tension (TCT) and compression-tension-compression (CTC) loading paths are calculated and compared with experimental results available in literature [3]. Moreover, the calculated results obtained via using the developed model are compared with those from the recent modified HAH model [3]. It is demonstrated that the developed analytic constitutive model can accurately capture the evolving asymmetric hardening behaviours of the magnesium alloys. The significance of the current constitutive model is that it can successfully model HCP sheet metal behaviours under cyclic loading with better accuracy in a more comprehensive pre-strain range compared to the existing models available in literature. Moreover, the hardening material parameters during all three deformation modes can be defined as pre-strain functions. Therefore, in contrast to those existing models in literature, a systematic parameter determination process can be applied to determine hardening parameters, which can mitigate interpolation errors. Finally, the developed model can well capture elastoplastic behaviours of the materials within any range of pre-strains, which verifies its versatility.

Modelling approach
In this study, phenomenological deformation modes of MC, MT, RC, and RT are defined to model evolving plastic behaviours of HCP sheet metals under cyclic loadings. Figure 1 shows the MT and MC modes including the initial tension and compression of the TCT and CTC curves, respectively. In addition, the RC mode consists of the reverse compression following tension, while the RT mode refers to the reverse tension following compression of the TCT and CTC curves. Young's moduli in MC and RT modes are considered differently, which aligns with experimental observations from literature, for example, experimental observations from Lou et al. [22]. Four CB2004 type yield surfaces are assigned to the deformation modes to develop the constitutive equations based on magnesium alloys' twining/untwining behaviour. Each yield surface is only active when its deformation mode is activated. A general reference fellow stress equation is introduced and developed in Sect. 2.4, representing the yield surfaces' hardening evolution at each deformation mode.
To indicate the direction of loading, an indicator function can be defined by summing the in-plane principal strain increments, d 1 and d 2 as where sgn represent the sign function. By assuming that the sheet metals are at least under a nonzero in-plane normal stress during deformation, the indicator function gives the value of 1 for d 1 + d 2 ≥ 0 (tension) and 0 for d 1 + d 2 < 0 (compression). Similar approaches to indicate the loading directions were also adopted by several studies in literature. For example see the Lee et al. [25] model. Moreover, to detect the reverse loading state, especially for numerical implementation purposes, a reverse loading criterion is adopted from literature [26]. As illustrated in Fig. 2, this criterion is defined using on the relative angle between two stresses, θ rel representing the previous, d old and the current loading directions, d new . In addition, θ ref is a prescribed reference angle with a typical value of π/2. For reverse loading to happen θ rel > θ ref , the relative angle can be defined as In Eq. (2), (m) and (m+1) are stress tensors of the previous and current time steps, respectively.

Constitutive equations
In order to model the asymmetric evolution of hardening responses of magnesium alloys, including untwining behaviours in reverse tension, a general constitutive equation is proposed as In Eq. (3), f ( ) is the yield function, yt and yut are the initial tensile and untwining yield stresses, and R MC iso , R MT iso and R RT iso are isotropic hardening functions in MC, MT and RT modes, respectively. Also, c is a material parameter, and c n 0 is the initial strength difference parameter of CB2004 yield criterion in the current mode and G RT and G RC are the initial conditions of isotropic hardening functions in RT and RC modes, respectively. The initial yield surface's size of the active mode equals the previous yield surface's size, which is updated using the initial conditions. In other words, the initial size can be determined using accumulated plastic strains in the previous mode. However, the yield surface expansion at the active mode is controlled by its local equivalent accumulated plastic strains. Moreover, v is an indicator function, indicating the untwining deformation during the RT mode and can be expressed as , σ ct is the initial compressive yield stresses, and n eq and n−1 eq are the current and previous modes' local equivalent accumulated plastic strains, respectively.

Modified yield function
The isotropic CB2004 yield function can be expressed as where τ y is the yield stress in pure shear, c is the strength difference parameter, J 2 and J 3 are the second and third invariants of the stress deviator tensor, respectively. The material � , the yield surface is convex [17]. While the asymmetric isotropic yield surface can represent individual yield surfaces for fixed accumulated plastic deformation levels, it cannot account for the continuous evolving TCA caused by the texture evolution during twinning deformation. To take this effect into account, the parameter c is defined as a function of local equivalent accumulated plastic strain and can be defined as In Eq. (7), s is a function of local equivalent accumulated plastic strain in the previous mode and can be expressed as Hill [27] showed that the local equivalent accumulated plastic strain associated with the asymmetric isotropic yield function by applying the plastic work-equivalence principle. In this study, the tensor components of Cauchy stress is denoted collectively by the (3 × 3) matrix σ ij or by the (9 × 1) column σ. Similarly, the strain tensor components are denoted by ε ij or ε. Based on this principle and considering any f (σ) can be made homogenous of degree one with respect to positive multipliers, the following expression can be written where the symmetric tensor β denotes a generic outward normal to the yield surface at a yield point. Then, by considering a unique β, it can be defined as By applying the associated flow rule and considering equivalent plastic strain increment as the plastic multiplier, the plastic strain can be defined as Note that, under the uniaxial loading condition, β is reduced to a scaler. More details can be found in Mehrabi = n eq et al. [20]. In this study, β refers to the first tensor component β 11 .

Phenomenological deformation modes
In this section, isotropic hardening responses in each deformation mode are presented. For this purpose, based on the hardening law proposed for each mode, a reference flow stress equation is given. The hardening laws in each deformation mode control the expansion of their yield surfaces based on the local equivalent accumulated plastic strain.

MC mode
The MC mode is active when the monotonic compressive loading is applied. Accordingly, the MC yield surface is active during this deformation mode. In this mode, according to Eqs.
(1) and (4), w = 0 and v = 1 and the reference flow stress equation considering the initial condition of G RC = 0 can be written as where, R MC iso is the isotropic hardening in monotonic compression and can be defined with a Boltzmann sigmoidaltype function as In Eq. (13), K MC , M MC and N MC are material parameters, respectively, which can be obtained from monotonic compressive experimental test data. Also, n 0 and MC 0 are the initial in current and MC modes and MC can be determined by applying the initial condition of R MC iso (0) = 0 . Moreover, according to Eq. (8), s = 0 and the strength difference parameter c is a function of local equivalent accumulated plastic strain in the current mode.

MT mode
Monotonic tension deformation mode and its yield surface are active when monotonic tensile loading is applied. During this mode, according to Eqs. (1) and (4), w = v = 1. Therefore, the general reference flow stress equation with no plastic prestrains can be written as iso is the isotropic hardening in monotonic tension and it can be expressed with a Voce-type hardening rule as In Eq. (15), K MT and N MT are material parameters, which can be obtained by curve fitting of monotonic tensile test data. Similar to the MC mode, parameter c is a function of local equivalent accumulated plastic strain in the current mode.

RC mode
This mode refers to the reverse compression following the previous tension. In this mode, the RC yield surface is the only active yield surface, and according to Eqs. (1) and (4), w = 0, and v = 1. From Eq. (3), the general reference flow for the RC yield surface by considering a non-zero initial condition ( G RC ≠ 0 ) can be written as The RC modes' initial condition, with the initial value of v 0 = 1 can be defined as where subscripts n and n-1 denote the current and the previous modes, respectively, and β 0 can be found from the initial yield surface ( n eq = 0 ). Moreover, according to Eq. (8), s = 1 and consequently the yield criterion strength difference parameter is zero (c = 0).

RT mode
The RT deformation mode and its yield surface are active during the reverse tension following the previous compression. During this deformation mode, w = 1 and parameter v can be 0 or 1 depending on the local equivalent accumulated plastic strain. Hence, the general reference flow stress equation can be expressed as The RT modes' yield surface's size updates with the initial condition G RT , which can be expressed as where c 0 0 = c MC 0 is the initial value. Moreover, similar to the RC mode, according to Eq. (8), s = 1 and, consequently, the yield criterion strength difference parameter is zero (c = 0). According to Eq. (18), the RT mode's isotropic hardening function R RT iso , can be active in this mode, which represents the hardening response of the material associated with untwining deformation. Note that the deformation mechanisms for magnesium alloys can be not as simple as pure twining/untwining or slip. However, for simplicity, the four phenomenological deformation modes were assumed with the dominant deformation mechanisms during the plastic evolution of magnesium alloys. The twining deformation mechanism is predominant during the in-plane compression of magnesium alloy sheets and exhausts with continuous compression. Abrupt grain reorientation, creation and disappearance of twining boundaries lead to the initiation of slip mechanism again, and consequently, the flow stress rises. Due to the untwining, the reoriented crystals rotate back to their original orientations during the subsequent inplane tension, and the basal texture appears again. Since the untwining does not need twin nucleation, the untwining stress is initially less than the twinning stress. However, the disappearance of existing twins exhausts untwining deformation, and consequently, the slip deformation mechanism becomes dominant, which leads to higher stress and results in an inflected sigmoidal curve. In fact, untwining eliminates the original twins by producing the reverse of the previous strain caused in twinning deformation [22]. In the current model, the amount of the equivalent plastic strain, which takes untwining to be exhausted, is assumed to be a fraction of the compressive equivalent pre-strain. Figure 3 shows the schematic CT stress-strain curve and its critical points in the RT mode. For the range of RT eq,1 Fig. 3 Critical points on the stress-strain curve during RT mode to RT eq,3 the untwining deformation is happening and v = 0. At P 1 , RT eq,1 = yut + G RT , where yut is the initial untwining yield stress and RT eq,1 = yut is its corresponding strain. Moreover, the at P 2 , RT eq,2 = yuth + G RT , where yuth is the stress, which the strain is halfway between the untwining yield strain RT eq,1 , and the strain at the upper plateau region RT eq,3 . At P 3 , untwining exhausts and the slip-dominant deformation initiates. From this point onward v = 1, however, at this point the reference fellow stress in Eq. (18) should be valid for both values of v = 1 and 0. Therefore, the RT mode's hardening function at this point can be expressed as As shown in this figure, MC eq is the equivalent pre-strain or the local equivalent plastic strain in the MC mode, and yuth is a material parameter, which can be determined from a CT uniaxial test.
A modified Boltzmann sigmoid-Voce function is proposed for the isotropic hardening response of the material in the RT mode as In Eq. (21), the first and second terms represent a modified Boltzmann sigmoid and an exponential Voce type hardening function, respectively. Moreover, in this equation, A is a function of the local equivalent accumulated plastic strain, which can be defined as In the abovementioned equations, K RT , N RT , M RT and L RT are hardening parameters, which are functions of the equivalent pre-strain. The parameter determination procedure is fully explained in the next section.

Determination of material parameters in the RT mode
In this section, the material parameter determination process in RT mode is explained. According to Fig. 3, A is defined as the stress difference between the top and bottom regions of the RT curve A = RT eq,3 − RT eq,1 . Moreover, as mentioned in the previous section, stress at P 3 can be found from the RT mode's reference flow stress with v = 1. Therefore, by applying reference flow stress and Eq. (22), the hardening parameter K RT can be expressed as Hardening parameter N RT is defined as a function of slope ratio of uniaxial compressive and tensile hardening curves at the compressive plastic pre-strain and can be expressed as where |.| denotes absolute value. Parameters M RT and L RT in Eq. (21) can be determined by considering the critical points on the RT hardening curve, as shown in Fig. 3. At P 1 , the local equivalent accumulated plastic strain is zero, and the first principal stress equals the initial RT yield surface size. The local equivalent accumulated plastic strain at P 2 is RT eq2 = MC 0,un MC eq ∕2 . Therefore, by substituting the hardening equations (Eqs. (15) and (21)) into Eq. (18) and considering the stress at P 2 , RT eq,2 the following relationship can be written By manipulating the abovementioned equation and rearrange it for parameter L RT , it can be written as By applying the parameter L RT from Eq. (26) into Eq. (20) and considering RT eq,3 ∕ RT eq,2 = 2, the following relationship can be obtained Equation (27) is a nonlinear equation with two real roots for M RT , which can be found by applying the bisection method via programming a simple MATLAB code. By solving Eq. (27) and using Eq. (26), two sets of parameters with the following relationship can be found where subscripts 1 and 2 denote the first and second roots, respectively. Two curves based on the two sets of parameters are illustrated in Fig. 3. As it can be found, Curve 1, which is related to the first roots, can successfully reproduce the shape of the hardening response in the RT mode.

Parametric study
This section performs a parametric study to evaluate the effect and sensitivity of hardening parameters during reverse tensile loading. For this purpose, a CT curve based on the model material AZ31B magnesium alloy [3] with a prestrain of − 0.06 was considered as a baseline, as shown in Fig. 4. The hardening parameters used for the baseline are listed in Table 1. Figure 5 shows the effects of the hardening parameters on the shape of the CT curve. As shown in Fig. 5a, parameter K RT controls the top plateau region of the reverse tensile curves. This parameter controls the stress difference between P 2 and P 3 by shifting the top region. That is to say, lower values of this parameter result in lower values of the stress for the top region. The effects of parameter L RT are shown in Fig. 5b. This parameter controls the hardening rate of the curve through the entire untwining deformation. By applying higher values of this parameter, the rate of the curve increases. Finally, the effects of parameters N RT and M RT are depicted in Fig. 5c, d, respectively. These two parameters control the hardening rate between P 2 and P 3 . As it can be found, higher values of N RT slightly decreases the hardening rate, while higher values of M RT significantly increases the hardening rate. It is worth mentioning that the parameter M RT has the most influence on the reverse tensile curve than the other parameters.

Elastic deformation behaviours
To the best of the authors' knowledge, no clear practical method to date was developed to model the evolving elastic deformations of magnesium alloys during cyclic TCT and CTC loadings. In the current model, the chord method, which was discussed thoroughly in Mehrabi et al. [20], was applied to take the variation of elastic modulus with plastic strain into account for the MT and MC modes. However, as shown in Fig. 1, the elastic modulus in the RT mode ( E RT eff ) is different from the elastic modulus in the MC mode ( E MC eff ). In this section, a new analytical model is developed to capture the elastic behaviours of magnesium alloys in the RT mode. A nonlinear equation can define the elastic modulus as where a RT E , b RT E and c RT E are material constants. In the proposed model, the maximum value of elastic modulus in RT mode is assumed to be equal to the saturated elastic modulus in monotonic compression E MC a . In fact, the above equation is only valid for strains between 0 and the corresponding equivalent strain of saturated elastic modulus MC eq a , and for higher strains, the effective elastic modulus is E MC a . By applying the mentioned boundary condition and considering the initial condition of E RT eff (0) = E 0 and one arbitrary data point n eq * , E * from experimental data into Eq. (29), the parameter c RT E can be expressed as By applying c ERT and a ERT into Eq. (29), the parameter b ERT can be defined by solving the following equation: Equation (31) is a nonlinear equation, which can be also solved using the bisection method. The elastic material Except for E * and n eq * , all chord method's parameters are discussed in Mehrabi et al. [20].  Table 2, the variations of elastic modulus in the RT mode versus compressive plastic pre-strain can be shown in Fig. 6. As it is shown, the proposed model can successfully capture the elastic behaviours of AZ31B magnesium alloy during reverse loading.

Results and discussion
In this section, material behaviours of the model material, AZ31B magnesium alloys during TCT and CTC loadings, are calculated based on the developed model, and the extracted results are compared with those experimental data from literature. Note that the current model is capable of representing HCP metals' twining/untwining and its continuous evolving asymmetry under cyclic loading. According to Yoon et al. [28], the dominant factor in the plastic behaviour of magnesium alloy sheets is the initial asymmetry in tension and compression. Hence, considering the aim of this study, the effects of initial anisotropy of the material were neglected. Figure 7 shows the yield surface evolution during MC and MT modes for eleven fixed levels of the local equivalent plastic strains ( n eq = 0, 0.01, 0.02, … , 0.1 ). The loading during MT and MC modes occurs within the third and first stress quadrants and are shown with solid lines in Fig. 7a, b, respectively. As it can be seen, the developed model with the modified CB2004 yield criterion can successfully capture the continuous evolving TCA of the material's plastic response due to the induced twining-slip texture change. The modelled plastic responses for AZ31B magnesium alloy during the MC and MT loadings algin with the experimental observations from Muhammad et al. [29] and theoretical results provided in Lee et al. [3]'s study. Figure 8 shows the yield surface evolution during the RC mode with different tensile equivalent pre-strains of MT eq = 0.016and0.036 . In the present model, the loading during the RC mode occurs within the third stress quadrant, which is shown by solid lines. The yield loci are plotted for eleven fixed levels of the local equivalent plastic strains ( n eq = 0, 0.01, 0.02, … , 0.1). Figure 9 shows the yield surface evolution during the RT mode with different compressive equivalent pre-strains of MC eq = 0.0125and0.0264 , respectively. The loading during the RT mode occurs within the first stress quadrant, which is shown by solid lines. The yield loci are plotted for eleven fixed levels of the local equivalent plastic strains ( n eq = 0, 0.01, 0.02, … , 0.1).

Yield surface prediction
The modelled plastic responses for AZ31B magnesium alloy during the RC and RT loadings show a similar trend as the experimental results provided in the literature by Muhammad et al. [29]. These yield loci can be further calibrated by the experimental results via applying anisotropy coefficients into the yield criterion, which will be considered for future studies.

Analytical results for cyclic behaviours of magnesium alloys
Additional to those RT mode hardening and elastic material parameters listed in Tables 1 and 2, and the rest of the plastic material parameters for the AZ31B magnesium alloy, including the hardening parameters during the MC and MT modes, are listed in Table 3. Superscripts T and C denote tension and compression modes, respectively. The experimental stress-strain curves are compared with analytical results based on the new and modified HAH models in Fig. 10. To determine the accuracy of the proposed model, for each case, the normalised root mean square error (NRMSE) of the two models are calculated and listed in Table 4. Moreover, number of check points and experimental repetitions for each curve is listed in this table. In this study, the NRMSE is defined as     where RMSE and SD are root mean squared error and standard deviation, respectively, and can be defined as where i is the measured stress from experimental results for the ith point, ̂ i is the predicted stress for the i th point, ̃ is the mean value of the measured stress from experiments, and n is the number of points or sample size, respectively. The model responses to CTC and TCT loading paths of AZ31B with different pre-strain values are shown in Fig. 10. In this figure, the numbers in front of T or C show the prestrain in each loading direction, i.e., T2C2T6 in Fig. 10a means pre-strains of 2% in both tension and reverse compression and pre-strain of 6% in the final tension. As it can be seen, the proposed elastic-plastic material model can successfully reproduce the experimental cyclic stress-strain curves under different loading paths.
For the current model, most of the discrepancies between calculated and experimental results are at the upper plateau region of the RT curves. These discrepancies are primarily due to the abrupt transition between untwining and slip dominant deformations. However, the results based on the Lee et al. [3] model show discrepancies in the elastic regions, reverse-compression under CTC loading path and the upper plateau region of the reverse tensile curves. In order to investigate the accuracy of the models, NRMSEs of the reproduced stress-strain curves with different loading paths are calculated and listed in Table 4 for the two models. Based on the given data in Table 4 It is worth mentioning that the highest and lowest errors for both models belong to C2T2C6 and T4C4T6 loading paths, respectively.
While the prediction accuracy of the developed model for paths with higher pre-strain values (C4T4C6 and T4C4T6) may be a little limited, the prediction accuracy of those with lower pre-strains (C2T2C6 and T2C2T6) are much high. This clearly shows that the developed model can have a high accuracy at a more comprehensive pre-strain range, which significantly enhances the accuracy and versatility of the model in general. The high prediction accuracy of the model in any pre-strain range is achieved with the application of variable RT hardening coefficients and yield function parameter c, which are defined as functions of plastic pre-strains.
Thus, the developed model shows excellent accuracy with an average error of approximately 12%, which is a significant improvement compared to the existing models in literature.

Concluding remarks
In this research, based on the concepts of multiple yield surface approaches and isotropic hardening, a new phenomenological elastoplastic model has been developed to capture twining/untwining behaviours of HCP metals under TCT and CTC cyclic loadings. This proposed constitutive model can be used to describe their elastoplastic behaviours by assigning a CB2004 type yield surface to four phenomenological deformation modes, including monotonic compression (MC), monotonic tension (MT), reverse compression (RC), and reverse tension (RT). In addition, the RT hardening parameters are defined as the functions of plastic pre-strains to mitigate interpolation errors caused by parameter determination processes of those existing models. In the showcase study, the proposed continuum-based constitutive model is able to capture the asymmetric hardening response of AZ31B magnesium alloy in monotonic and reverse loading paths. In addition, the results have been validated via comparing them with experimental and calculated results in literature-the modified homogeneous anisotropic hardening (HAH) model. The obtained results show that the new constitutive model can successfully reproduce experimental tension-compression-tension (TCT) and compressiontension-compression (CTC) stress-strain curves of HCP sheet metals with considerably less percentage errors.
Based on the obtained analytical results, the following conclusions can be drawn: • The proposed model can successfully capture the main characteristics of the stress-strain curves of AZ31B magnesium alloy in TCT and CTC loading paths. • The new material parameter identification method defines the hardening material parameters as functions of equivalent pre-strain. This straightforward parameter identification process leads to excellent reproduction accuracy of the model, which can make this model a powerful tool in sheet metal forming simulations. Due to the high accuracy of the developed model in predicting the material behaviours at any pre-strain range, it can have a great application in FEA studies of sheet metal forming, especially the bending process such as, pure bending, bending under tension, deep draw bending, etc.
• Compared to the recent modified HAH model of Lee et al. [3], the new constitutive model can reproduce the stress-strain curves of AZ31B magnesium alloys with more accuracy. Not only the fitting errors (NRMSE) of the new model for all the loading paths are less than those of the modified HAH model, but also the average fitting error (average NRMSE) of the new model is 4.35% less than that of the modified HAH model. • The modified CB2004 yield criterion and the defined hardening rules can successfully capture the continuous evolving TCA caused by the texture evolution during twinning deformation in the MC and MT modes. Moreover, the yield loci for the four loading modes show an acceptable trend for the plastic behaviours of HCP metals. • The proposed elastic model can well capture the elastic behaviour of the material in monotonic unloading and reverse loadings.
As for the future work to improve the current model, firstly, by applying the general anisotropic form of CB2004 yield surface in the constitutive model, the anisotropic effects can be included to improve the current model further. Secondly, as mentioned earlier, most of the discrepancies between calculated and experimental results are at the upper plateau region of the RT curves since the untwining to slip deformation transition happens abruptly. Therefore, it can be further developed to improve the accuracy, explicitly focusing on these abrupt transitions. Funding Open Access funding enabled and organized by CAUL and its Member Institutions. The first author was supported during his PhD study by the Postgraduate Scholarship Awards provided by the Western Sydney University, Australia. This research received no external funding.

Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.