Application of an Evolving Non-Associative Anisotropic-Asymmetric Plasticity Model for a Rare-Earth Magnesium Alloy

: Magnesium sheet metal alloys have a hexagonal close packed (hcp) crystal structure that leads to severe evolving anisotropy and tension-compression asymmetry as a result of the activation of different deformation mechanisms (slip and twinning) that are extremely challenging to model numerically. The low density of magnesium alloys and their high speciﬁc strength relative to steel and aluminum alloys make them promising candidates for automotive light-weighting but standard phenomenological plasticity models cannot adequately capture the complex plastic response of these materials. In this study, the constitutive plastic behavior of a rare-earth magnesium alloy sheet, ZEK100 (O-temper), was considered at room temperature, under quasi-static conditions. The CPB06 yield criterion for hcp materials was employed along with a non-associative ﬂow rule in which the yield function and plastic potential were calibrated for a range of plastic deformation levels to account for evolving anisotropy under proportional loading. The non-associative ﬂow rule has not previously been applied to magnesium alloys which require the use of ﬂexible constitutive models to capture the severe anisotropy and its evolution with plastic deformation. The non-associative ﬂow rule can provide the required ﬂexibility by decoupling the yield function and plastic potential. For the associative ﬂow rule, such ﬂexibility can only be achieved by multiple linear transformations of the stress tensor resulting in expensive models for calibration and simulations. The constitutive model was implemented as a user material subroutine (UMAT) within the commercial ﬁnite element software, LS-DYNA, for general 3-D stress states along with an interpolation technique to consider the evolution of anisotropy based upon the plastic work. To evaluate the accuracy of the implemented model, predictions of a single-element model were compared with the experimental results in terms of ﬂow stresses and plastic ﬂow directions under various proportional loading conditions and along different test directions. Finally, to assess the predictive capabilities of the model, full-scale simulations of coupon-level formability experiments were performed and compared with experimental results in terms of far-ﬁeld load-displacement and local strain paths. Using these experiments, the constitutive model was evaluated across the full range of representative stress states for sheet metal forming operations. It was shown that the predictions of the model were in very good agreement with experimental data.


Introduction
The growing demand for vehicle weight reduction to decrease fuel consumption and protect the environment has motivated research on light-weight alloys and among these, magnesium appears to be a promising candidate for significant weight reduction. In particular, the development of rare-earth magnesium alloys, such as ZEK100, have the potential to be used as light-weight structural alloys since they possess a weakened crystallographic texture that leads to superior ductility over conventional magnesium sheets. However, there are some complications for commercial applications of these alloys that are mainly attributed to low corrosion resistance and significant plastic anisotropy [1]. Challenges with corrosion can be overcome with emerging coating techniques (e.g., [2]). In terms of anisotropic constitutive response, these alloys require extensive material characterization to develop and calibrate advanced constitutive plasticity models for design of metal forming processes of automotive components and to predict crash performance. Most commercial alloys used in sheet metal forming applications exhibit an orientation-dependent plastic response although some of these materials can be reasonably idealized as isotropic for metal forming and vehicle crash simulations (e.g., [3]). However, for light metals such as aluminum and magnesium, the plastic anisotropy can be severe and non-trivial. To consider the plastic response of such materials, many anisotropic yield functions have been proposed in the literature including the Barlat family of yield criteria [4][5][6][7][8] in which linear transformations are applied on the stress tensor to account for anisotropy. However, these models were intended for bcc and fcc cubic materials with slip-dominated deformation mechanisms, while magnesium alloys have an hexagonal close packed (hcp) crystal structure. In hcp materials, plastic deformation occurs by slip and twinning mechanisms and due to directional sensitivity of twining mechanisms, strong tension-compression asymmetry has been observed in the yielding behavior (e.g., [9][10][11][12][13]). In order to account for anisotropy and strength differential effects for hcp materials, Cazacu and Barlat [14] generalized Drucker's isotropic yield criterion [15] followed by an anisotropic extension developed by Cazacu et al. [16] known as the CPB06 yield criterion based on the linear transformation approach. Plunkett et al. [17] demonstrated that the accuracy of the CPB06 model could be significantly improved if more than one linear stress transformation was applied. The designation "CPB06ex2" indicates that two transformations have been applied to the stress deviator. Despite the improvements in the CPB06 framework, the initial anisotropy at the onset of yielding can be well described but unfortunately, anisotropy significantly evolves with plastic deformation in hcp materials.
Due to texture evolution during plastic deformation, the shape of yield surface of hcp alloys does not remain constant, thus, isotropic hardening models cannot capture the materials behavior accurately. To consider evolving anisotropy, Plunkett et al. [17] proposed a piece-wise linear interpolation of the CPB06 yield function calibrated at different levels of plastic deformation. Similar approaches were adopted in Reference [11,[18][19][20][21] for magnesium alloys and in Reference [22][23][24][25] for titanium alloys. It should be noted that evolving plastic anisotropy and distortional hardening in hcp alloys occurs even during monotonic proportional loading. The role of non-proportionality on the yielding behavior of materials is well-established with concepts of kinematic hardening [26] and distortional hardening in the homogeneous anisotropic hardening (HAH) models of Barlat et al. [27,28]. These concepts were also utilized in the literature to capture the evolution of anisotropy in hcp materials under proportional loadings. For instance, Steglich et al. [29] employed an associative flow rule and combined two dissimilar yield functions of CPB06 and Yld91 [5] in different regions of the yield loci along with a kinematic hardening function to account for evolving anisotropy of magnesium alloys. More recently, Kondori et al. [30] used an approach that was similar to that of Steglich et al. [29] that is, utilizing two yield surfaces of CPB06 and Yld91. In the study by Kondori et al. [30] the model was calibrated by accounting for 3-D stress states for an AZ31B plate with promising results for in-plane and out-of-plane uniaxial tension and compression stress states.
In addition, to address evolving anisotropy of magnesium alloys, Mekonen et al. [18] and Ghaffari Tari et al. [11] employed saturating exponential-type functions to evolve the yield function parameters for AZ31B in which it was assumed that anisotropy will stop evolving at some saturation strain level. Using this approach, there is a possible bias introduced by the selection of the functional form and the number of parameters in these interpolations can restrict calibration accuracy, specifically for ZEK100 which displays severe anisotropy, since the parameters are forced to follow a pre-assumed trend. Also there is no guarantee that the interpolated coefficients produce loci with interpolated shapes. Thus, in the present paper a different methodology is utilized to directly interpolate the yield surface and plastic potential rather than their coefficients. This modelling framework is developed based on the experimental data of Abedini et al. [12] for ZEK100-O. Note that the work by Abedini et al. [12] was mainly focused on experimental constitutive characterization of this material under different stress states and discussion of its asymmetry in different quadrants of yield locus while the present paper aims at developing a modelling framework to simulate the response of the material under a range of mechanical loading conditions. The studies mentioned above were all based upon an assumption of the associative flow rule (AFR) in which the yield function is also the plastic potential for calculating the plastic strain increments (normality rule). The classical work of Bishop and Hill [31] demonstrated that the AFR was valid for metals based on a crystal plasticity model. The use of plasticity models that employ a non-associative flow rule (non-AFR) has been on the rise in recent years after the pioneering work of Stoughton [32] and Stoughton and Yoon [33] who developed a series of necessary constraints upon the non-AFR models to ensure that Drucker's postulate [34] is satisfied. The flexibility of non-AFR models to efficiently and accurately characterize anisotropy in both the stress and strain responses has led to successful applications such as in References [35][36][37]. Considering materials with severe anisotropy, such as magnesium alloys, the non-AFR provides a high degree of flexibility for calibrating yield stresses and r-values. Furthermore, models with the non-AFR enable the possibility that the plastic potential and yield function can evolve independently, a feature that is not possible in models based on the AFR. These characteristics make the non-AFR of potential interest to model materials such as ZEK100-O magnesium which exhibit strong evolving anisotropic and asymmetric behavior [9,12]. Despite these advantages, to the best of authors' knowledge, models based on the non-AFR have not been evaluated in the literature for hcp alloys.
The objectives of the present work are to investigate and model the anisotropic plastic response of a rare-earth magnesium alloy sheet (ZEK100-O) under quasi-static conditions. There is a clear and present need to investigate a 3-D CPB06 framework using the non-AFR that can be efficiently calibrated to high accuracy across a range of loading conditions. Such a modeling framework must be able to capture the evolving anisotropic-asymmetric response of the material. The present paper addresses this need. To this end, the non-AFR is employed by calibrating the yield function and plastic potential at different plastic deformation levels using the CPB06 formulation with two linear stress transformations, that is, CPB06ex2. The full 3-D model is implemented into the commercial finite element package, LS-DYNA, along with an interpolation technique to consider the evolving anisotropy of the material. Predictions of the model were evaluated using single-element simulations in which the finite element model comprises a 3-D element that is subjected to various stress states in different test orientations with respect to the rolling direction. The outcomes of each single-element simulation in terms of flow stresses and r-values are compared to experimental results of Abedini et al. [12]. Furthermore, to assess the predictive capabilities of the model across the wide range of stress states experienced in practical sheet metal forming operations, full-scale simulations of coupon-level formability experiments were performed and compared with results of experiments conducted in the present study with the aid of full-field stereoscopic digital image correlation (DIC) techniques. Finally, to highlight the importance of considering anisotropy evolutions, predictions of evolving and non-evolving models are compared in terms of their accuracy to capture the experimental data.

Material
A rare-earth magnesium alloy ZEK100-O rolled sheet with a nominal thickness of 1.55 mm was used in the present study. Recently, Abedini et al. [12] performed X-Ray Diffraction (XRD) analysis on this material and reported that the texture of ZEK100-O exhibits a spread of the basal poles along the transverse direction which is typical of rare-earth magnesium rolled sheets ( Figure 1). Compared to the more common commercial AZ31B magnesium alloy, ZEK100-O possesses a relatively weaker basal texture that translates to higher ductility; however, a strong evolving anisotropy persists which is significantly more severe than that of AZ31B [9,12,13].

Plastic Response
An extensive experimental investigation into the plastic response of this same lot of material was performed by Abedini et al. [12] and this test data will be utilized in the present paper to develop the constitutive plasticity model. The experiments reported by Abedini et al. [12] were performed at room temperature, under a quasi-static strain rate of 0.001 s −1 for a range of stress states. The constitutive behavior of the material is presented in Figure 2 where it can be seen that ZEK100-O exhibits significant anisotropy. It is apparent from Figure 2a-c that the material displays a tensioncompression asymmetric response that is indicated by the sigmoidal shape of the compressive stressstrain curves and is due to the twinning mechanisms that are more dominant under the compression mode [9]. Moreover, the evolving tensile r-values in Figure 2d highlight the significant evolution of the texture with deformation. For instance, the r-value in the transverse direction is about 0.20 at the onset of yielding and evolves to be greater than unity by the onset of necking in the tensile test.
The response of the material under an equal-biaxial tension state was also determined using through-thickness compression experiments under the assumption that yielding is independent of the pressure [11,38]. Note that the modelling approach used in the subsequent sections is based on the CPB06 framework that ignores pressure sensitivity in the plastic response of the material. The result of the equal-biaxial tensile response is shown in Figure 2e in which the uniaxial tensile response of the material in the rolling direction (RD) and transverse direction (TD) are also presented for comparison purposes. It can be seen that the equal-biaxial response of the material lies between that of the uniaxial tensile states in the RD and TD. In addition, Figure 2f presents the behavior of the material undergoing shear deformation in which anisotropy is apparent in terms of the initial yield strength and hardening rates. The experimental results shown in Figure 2 show the challenging nature of ZEK100-O in terms of modelling and characterization and highlight the need for accurate constitutive plasticity models that can capture the evolving anisotropy of the material.

Plastic Response
An extensive experimental investigation into the plastic response of this same lot of material was performed by Abedini et al. [12] and this test data will be utilized in the present paper to develop the constitutive plasticity model. The experiments reported by Abedini et al. [12] were performed at room temperature, under a quasi-static strain rate of 0.001 s −1 for a range of stress states. The constitutive behavior of the material is presented in Figure 2 where it can be seen that ZEK100-O exhibits significant anisotropy. It is apparent from Figure 2a-c that the material displays a tension-compression asymmetric response that is indicated by the sigmoidal shape of the compressive stress-strain curves and is due to the twinning mechanisms that are more dominant under the compression mode [9]. Moreover, the evolving tensile r-values in Figure 2d highlight the significant evolution of the texture with deformation. For instance, the r-value in the transverse direction is about 0.20 at the onset of yielding and evolves to be greater than unity by the onset of necking in the tensile test.
The response of the material under an equal-biaxial tension state was also determined using through-thickness compression experiments under the assumption that yielding is independent of the pressure [11,38]. Note that the modelling approach used in the subsequent sections is based on the CPB06 framework that ignores pressure sensitivity in the plastic response of the material. The result of the equal-biaxial tensile response is shown in Figure 2e in which the uniaxial tensile response of the material in the rolling direction (RD) and transverse direction (TD) are also presented for comparison purposes. It can be seen that the equal-biaxial response of the material lies between that of the uniaxial tensile states in the RD and TD. In addition, Figure 2f presents the behavior of the material undergoing shear deformation in which anisotropy is apparent in terms of the initial yield strength and hardening rates. The experimental results shown in Figure 2 show the challenging nature of ZEK100-O in terms of modelling and characterization and highlight the need for accurate constitutive plasticity models that can capture the evolving anisotropy of the material.

Plasticity Model
Sheet metal forming and vehicle crash simulations rely on an accurate description of the plastic behavior of materials. In general, there are two major approaches to describe plasticity or yielding behavior of polycrystalline materials. The first is the crystal plasticity approach which is based on physical aspects of plastic deformation and on averaging the response over a large number of grains. The crystallographic texture is the main input for these models. Crystal plasticity models have been implemented into finite element codes and have successfully predicted the plastic response of magnesium alloys [39][40][41][42]. Unfortunately, such finite element calculations are computationally expensive, thus limiting their applicability for industrial large-scale simulations. The second approach uses phenomenological continuum-based models that have shorter computational times compared to the crystal plasticity approaches. In addition, phenomenological models are easier to implement into finite element codes. In the present study, a continuum-based phenomenological approach has been selected. This section is devoted to describing the plasticity model, development of the interpolation method for evolving anisotropy and its implementation into a finite element software package.

Yield Function and Plastic Potential
The phenomenological yield criterion CPB06 proposed by Cazacu et al. [16] was adopted in the present study to investigate the yielding behavior of ZEK100-O. In analogy to Plunkett et al. [17], two linear stress transformations were performed to increase the degree of flexibility of the model (denoted as CPB06ex2) which is defined as: where k and k are material parameters that account for strength differential effects and a is the exponent of the yield function. Also Σ 1-3 and Σ 1-3 are the principal values of the transformed stress deviators written as: in which s is the deviatoric stress tensor and the fourth-order transformation tensors C and C are given by: Note that the symbol ":" denotes the doubled contracted product between two tensors and boldface letters refer to tensor quantities. Second-order tensors are italicized such as s whereas fourth-order tensors are non-italicized such as C. In Equation (1), B and B are given by: where It is important to mention that the CPB06 criterion can be generalized to n-transformations as done in References [11,17] who considered up to four stress transformations. As each transformation introduces additional parameters, the flexibility of the yield function increases but issues of calibration arise since it can be challenging to obtain the required number of experimental data points. Alternatively, the non-AFR formulation of CPB06ex2 was selected in the present work with two transformations for both the yield function and the plastic potential. Effectively, this leads to the same number of parameters as for the AFR model with four transformations but with the advantage of decoupling the stresses and r-values that greatly simplifies calibration of the parameters and can improve accuracy over an equivalent four-transformation model as demonstrated in Section 6. Thus, the non-associative flow rule was employed with the same functional form as Equation (1) to define the plastic potential, Ψ. The plastic strain increments are normal to Ψ and are: where dε p is the incremental plastic strain tensor and dλ is the plastic multiplier. Using the principle of plastic work equivalence and assuming that the plastic potential is a first-order homogeneous function, the following relation can be used to calculate the equivalent plastic strain Reference [35]: in the case of the associative flow rule (Φ = Ψ), Equation (10) reduces to: In terms of stability, Drucker [34] showed that models based on the AFR are always stable; however, Stoughton and Yoon [33] discussed that the AFR is not the only postulate for stability. These authors proved that by imposing the following constraint, the rate of change of stress and plastic strain are uniquely defined which is a necessary condition for stability: where σ is the equivalent stress and L is the fourth-order stiffness tensor. There are additional constraints that have been mentioned by Stoughton [32] and Stoughton and Yoon [33] for stability of non-AFR models but the constraint of Equation (12) is the most important for numerical implementation as discussed in Reference [43]. The reader is referred to Stoughton [32] and Stoughton and Yoon [33] for details on the proof of the uniqueness of stress and strain states as well as the proof of stability of the non-AFR and consistency in terms of plastic dissipation. The values of the coefficients of transformation tensors, C and C and the strength differential parameters, k and k , can be determined from an optimization approach to minimize the difference between the experimental data and the values predicted by the yield function and plastic potential. In the present study, the genetic algorithm (GA) which is a global optimizer available in Matlab ® (R2013b, MathWorks, Natick, MA, USA) was used to minimize "Error" given by: (13) where j, k and l are the number of tensile, compressive and shear stresses, respectively. Note that the superscripts t and c show the tension and compression values. In addition, the subscripts exp and CPB indicate whether the values are from experiments or from the model. To calibrate the plasticity model, seven uniaxial tensile yield stresses and r-values (in 15 • increments from the RD, that is, 0 • , 15 • , 30 • , 45 • , 60 • , 75 • and 90 • ) were used along with an equal-biaxial yield stress and r-value (j = 8). Moreover, three uniaxial compression yield stresses (k = 3) and three shear stresses (l = 3) were utilized. Note that the first three terms in Equation (13) are considered for calibrating the yield function while only the last term in Equation (13) is needed to calibrate the plastic potential.
Since the magnesium alloy exhibits an evolving anisotropic response, in the present study the experimental values corresponding to nine levels of plastic work per unit volume, w p , associated with equivalent plastic strains of 0.01 to 0.09 for uniaxial tension in reference direction of the RD were considered to calibrate the yield function and plastic potential. The values of the plastic work are given in Table 1 and the yield function and plastic potential were calibrated at each of these individual plastic work levels. It is worth mentioning that in the absence of through-thickness shear data to calibrate C 55 , C 55 , C 66 and C 66 , the yield function and plastic potential each include 17 coefficients for each calibration level to be determined from experimental data. The upper bound value of the plastic work level for calibration of the anisotropic yield function and plastic potential is limited by the onset of plastic instability (necking) for the tensile test in the RD direction which occurs at a strain of approximately 0.09. Due to this limitation, it was assumed that the shape of the yield loci and plastic potential remain constant after this level of plastic deformation.
It is important to note that stacked compressive tests were performed to extract the yielding response of the material. It was observed that the compressive r-values had a large variation and were deemed unreliable. Therefore, although the yield function does capture the asymmetric yielding behavior of the material in terms of compressive and tensile strengths, in the absence of reliable compressive r-values, the plastic potential was assumed to be symmetric. This leads to the simplification that the r-values in tension and compression are assumed to be equal and is a limitation of the present work. Note that the model is able to capture asymmetric yield stresses in tension versus compression that is critical from the point of view of forming and crash which are the intended applications of the present work. The compressive r-values might be better measured using an anti-buckling device that will be examined in future work.
Moreover, Swift effects (the occurrence of axial plastic deformation during monotonic free-end torsion of solid rods or tubes which was first observed by Swift [44]) were neglected for the shear state of the yield function. In other words, it was assumed that the principal strains are equal and opposite under a shear state (see also [45]). This assumption is in agreement with experimental investigations conducted by Abedini et al. [46] where equal and opposite principal strains were reported under shear state for small plastic deformation.
It is worth mentioning that the measured instantaneous r-values were directly used in Equation (13) to calibrate the model. The measured strains could also be directly used for calibration as advocated by Kondori et al. [30]. However, it is cautioned that general agreement of the predicted and measured total strains does not guarantee accuracy in terms of the instantaneous r-value which is a ratio of the normal vectors of the plastic potential that evolves with deformation for ZEK100-O.

Implementation into a Finite Element Code
The constitutive plasticity model was implemented as an explicit user material subroutine (UMAT) within the commercial finite element software package, LS-DYNA (R8.1.0), for solid elements. At the start of each time step, the current strain increments (time = t n+1 ) and the stress and plastic strain from the previous time step (t n ) are provided. The role of the UMAT is to update the stress and plastic strain for the current time step and return these values to LS-DYNA. To do so, the following stress integration algorithm was used in the present study which is based on the semi-explicit convex cutting plane (CCP) algorithm of Ortiz and Simo [47]. In this algorithm, all derivatives are calculated with variables from the previous time step, thus small time steps are required to ensure that this approximation is valid and that the hardening rate and directions of plastic flow are taken as constant over the increment. The algorithm is semi-explicit in that an iterative loop is used to determine the plastic multiplier at the end of the step. Note that all the variables described in this section correspond to the time increment of t n+1 unless stated otherwise.
At the start of a time step, the strain increment, dε, is assumed to be totally elastic and a trial stress is calculated as: It should be mentioned that the elastic response of the material was assumed to be isotropic with a Young's modulus of 39 GPa and Poisson's ratio of 0.3. If the trial stress lies inside or on the yield surface, the condition of Equation (15) given below holds and the trial stress is accepted as the current stress: If the condition of Equation (15) does not hold, the trial stress is outside of the yield locus hence plastic deformation occurs and the stress and plastic strain need to be updated such that: where TOL is a tolerance that was set to 10 −6 MPa in the present study. The stress and plastic strain increments are related by: In which dε e is the elastic strain increment tensor. Inserting Equation (9) into Equation (17) and using the consistency condition, the plastic multiplier can be obtained by: To determine the current yield function and plastic potential, an interpolation-based methodology similar to that presented by Gilles et al. [22] and Li et al. [24] was adopted. However, as opposed to these studies in which the AFR was used, in the present study the evolution of the yield function and plastic potential should be captured independently using a piece-wise linear interpolation technique: where Φ m (Ψ m ) and Φ m+1 (Ψ m+1 ) are the calibrated yield function (plastic potential). The subscripts "m" and "m + 1" are used to denote the plastic work levels of w p m and w p m+1 for which calibration of the yield function and plastic potential have been performed (see Table 1) and w Similarly, the interpolation technique is used to determine derivatives of the yield function and plastic potential: It should be mentioned that the derivatives of the CPB06 model for plane-stress conditions are given in References [16,17]. Extension from plane-stress conditions to general 3-D states can be obtained through the Cayley-Hamilton solution given by Barlat et al. [8] and is not provided here for brevity. In the UMAT constitutive subroutine, using Equations (17)-(23) the plastic strains are updated and the procedure is iterated until Equation (16) is satisfied, at which point the updated stress and plastic strain tensors are returned to LS-DYNA.
It is noted that the interpolation approach used in the present study differs from that adopted in Reference [11,18,20,25]. In those studies, the coefficients of the CPB06 yield function were interpolated as a function of equivalent strain. It should be noted that interpolating the values of two sets of yield function coefficients calibrated at two specific levels of plastic deformation, does not guarantee that the shape of the resulting yield locus is an interpolation of the shapes of the two calibrated yield loci. Therefore, in the present study, the methodology given in Equations (19)-(23) was utilized to directly interpolate Φ and Ψ rather than their coefficients.
Also, it is emphasized that in the most published studies mentioned above, the CPB06 yield function was reduced to 2-D stress states with shell elements (e.g., [11,21,48]), whereas the current implementation accounts for a general 3-D stress state. Shell elements are not suitable for the evaluation of local fields after the onset of necking while solids elements can provide more accurate predictions of localized deformation and through-thickness gradients of stress and strain fields. Solid elements are specifically required for fracture characterization of materials where resolving local fields is critical [49].

Hardening Function
In terms of the hardening response of the material, a Hockett-Sherby [50] function with the relation given in Equation (24) was adopted and calibrated with the experimental uniaxial tensile data in reference direction of the RD: where σ y is the initial yield stress and b, c and d are constant values that are presented in Table 2.
Predictions of the hardening function are compared with the experimental data in Figure 3a which shows good correlation. Note that the uniaxial tensile experimental data is limited by the onset of necking at a plastic strain of approximately 9%. For larger strains, the hardening behavior was approximated using an area-reduction method in which the flow stress of the material is obtained from the deformed cross-sectional area determined from DIC measurements at the localization zone (see [51]). This procedure was followed by iteratively adjusting the flow curve in a finite element model of the tensile specimen until the measured and predicted load-displacement curves were matched [52]. The load-displacement responses of the model and experiments are shown in Figure 3b. Note that the modelling approach for this analysis was consistent with that described in Section 5. This hardening function was used along with the plasticity model to simulate the behavior of the material under different loading conditions. (see [51]). This procedure was followed by iteratively adjusting the flow curve in a finite element model of the tensile specimen until the measured and predicted load-displacement curves were matched [52]. The load-displacement responses of the model and experiments are shown in Figure  3b. Note that the modelling approach for this analysis was consistent with that described in Section 5. This hardening function was used along with the plasticity model to simulate the behavior of the material under different loading conditions.

Coupon-Level Experiments
To evaluate the accuracy of the model, a series of experiments typical for sheet metal forming operations were performed on several specimen geometries. These tests feature multi-axial loading states; thus they can serve to assess the predictive capabilities of the developed plasticity model. Depending on the nature of the applied load, these tests are categorized into: in-plane and out-ofplane tests. The applied loads are within the plane of the sheet for the in-plane tests while the Nakazima dome tests utilize a punch that applies out-of-plane loads. The experimental procedure and the method of extracting data are presented in this section along with the characteristics of fullfield DIC strain measurements.

In-Plane Tensile Tests
Two different types of tensile tests were considered in the present paper using the geometries depicted in Figure 4a,b. These tests were performed using an MTS Criterion Model 45 servo-electric tensile frame and the results of the experiments were used to assess global and local predictions of the plasticity model described in Section 3.

Coupon-Level Experiments
To evaluate the accuracy of the model, a series of experiments typical for sheet metal forming operations were performed on several specimen geometries. These tests feature multi-axial loading states; thus they can serve to assess the predictive capabilities of the developed plasticity model. Depending on the nature of the applied load, these tests are categorized into: in-plane and out-of-plane tests. The applied loads are within the plane of the sheet for the in-plane tests while the Nakazima dome tests utilize a punch that applies out-of-plane loads. The experimental procedure and the method of extracting data are presented in this section along with the characteristics of full-field DIC strain measurements.

In-Plane Tensile Tests
Two different types of tensile tests were considered in the present paper using the geometries depicted in Figure 4a,b. These tests were performed using an MTS Criterion Model 45 servo-electric tensile frame and the results of the experiments were used to assess global and local predictions of the plasticity model described in Section 3.
The hole tension specimen shown in Figure 4a has become commonly adopted to characterize the plasticity and fracture of sheet materials [49,[53][54][55]. This test is of particular interest to the fracture community since the stress state at the edge of the hole remains close to uniaxial tension up to large values of deformation. However, the gauge area or ligament of the specimen experiences a range of stress states from uniaxial tension at the free boundaries towards the plane-strain state at some location behind the hole edge. Therefore, given the wide range of stress states, this test can be considered as a good evaluation test for the predictions of the plasticity model. The hole tension tests were performed in the three orientations of the rolling (RD), diagonal (DD) and transverse (TD) directions with a cross-head velocity of 0.005 mm/s resulting in a von Mises equivalent strain rate of approximately 0.001 s −1 at the gauge area. The hole tension specimen shown in Figure 4a has become commonly adopted to characterize the plasticity and fracture of sheet materials [49,[53][54][55]. This test is of particular interest to the fracture community since the stress state at the edge of the hole remains close to uniaxial tension up to large values of deformation. However, the gauge area or ligament of the specimen experiences a range of stress states from uniaxial tension at the free boundaries towards the plane-strain state at some location behind the hole edge. Therefore, given the wide range of stress states, this test can be considered as a good evaluation test for the predictions of the plasticity model. The hole tension tests were performed in the three orientations of the rolling (RD), diagonal (DD) and transverse (TD) directions with a cross-head velocity of 0.005 mm/s resulting in a von Mises equivalent strain rate of approximately 0.001 s −1 at the gauge area.
Furthermore, using notch tensile specimens, such as the one depicted in Figure 4b, is a classical way to control the stress triaxiality for plasticity and fracture characterization of materials. Due to the geometry of the notch tensile specimen, the stress triaxiality continuously increases with deformation leading to a severe triaxial stress state in the gauge area of the specimen, thus the notch test can also be considered as a good assessment test to analyze the predictive ability of the developed plasticity model. In a similar manner to the hole tension test, the notch tests were performed in the three orientations of RD, DD and TD with a cross-head velocity of 0.008 mm/s leading to an equivalent strain rate of approximately 0.001 s −1 at the gauge area.

Nakazima Dome Tests
Limiting dome height tests, proposed by Nakazima et al. [56], are standard experiments used for formability characterization of materials. In these tests a blank is clamped between a die and a binder and a hemispherical punch forms the specimen until failure. The blank geometries, shown in Figure 4c,d, conformed to the ISO12004-2 standard and were formed within an MTS hydraulic press with the tooling depicted in Figure 5. Dome tests with the specimen geometry shown in Figure 4c lead to a nearly plane-strain state at the center of the specimen. Furthermore, equal-biaxial tests were performed with the square blank depicted in Figure 4d. Dome tests involve tensile and compressive stress states, an advantage that assists with better evaluations of the tension-compression asymmetric predictions of the developed plasticity model. Moreover, compared to the in-plane tensile tests, the stress states in the dome tests are expected to remain closer to plane-stress conditions, providing a good assessment of the model in these states of stress. Furthermore, using notch tensile specimens, such as the one depicted in Figure 4b, is a classical way to control the stress triaxiality for plasticity and fracture characterization of materials. Due to the geometry of the notch tensile specimen, the stress triaxiality continuously increases with deformation leading to a severe triaxial stress state in the gauge area of the specimen, thus the notch test can also be considered as a good assessment test to analyze the predictive ability of the developed plasticity model. In a similar manner to the hole tension test, the notch tests were performed in the three orientations of RD, DD and TD with a cross-head velocity of 0.008 mm/s leading to an equivalent strain rate of approximately 0.001 s −1 at the gauge area.

Nakazima Dome Tests
Limiting dome height tests, proposed by Nakazima et al. [56], are standard experiments used for formability characterization of materials. In these tests a blank is clamped between a die and a binder and a hemispherical punch forms the specimen until failure. The blank geometries, shown in Figure 4c,d, conformed to the ISO12004-2 standard and were formed within an MTS hydraulic press with the tooling depicted in Figure 5. Dome tests with the specimen geometry shown in Figure 4c lead to a nearly plane-strain state at the center of the specimen. Furthermore, equal-biaxial tests were performed with the square blank depicted in Figure 4d. Dome tests involve tensile and compressive stress states, an advantage that assists with better evaluations of the tension-compression asymmetric predictions of the developed plasticity model. Moreover, compared to the in-plane tensile tests, the stress states in the dome tests are expected to remain closer to plane-stress conditions, providing a good assessment of the model in these states of stress.
To reduce friction, four layers of Teflon ® film lubricated with Vaseline ® were utilized between the punch and the blank. Initially, a die set with a lockbead was employed but cracking initiated at the lockbead during die closure. Instead, a flat die and binder were used with a clamping force of 660 kN that was sufficient to avoid draw-in during the test. A punch velocity of 0.25 mm/s was used for all Nakazima dome tests and this velocity corresponds to a von Mises equivalent strain rate of approximately 0.001 s −1 at the center of the specimens. The plane-strain dome tests were performed in the RD and TD directions. To reduce friction, four layers of Teflon ® film lubricated with Vaseline ® were utilized between the punch and the blank. Initially, a die set with a lockbead was employed but cracking initiated at the lockbead during die closure. Instead, a flat die and binder were used with a clamping force of 660 kN that was sufficient to avoid draw-in during the test. A punch velocity of 0.25 mm/s was used for all Nakazima dome tests and this velocity corresponds to a von Mises equivalent strain rate of approximately 0.001 s −1 at the center of the specimens. The plane-strain dome tests were performed in the RD and TD directions.

Strain Measurements
The experiments described above were recorded using two digital cameras to obtain full-field stereoscopic logarithmic strain measurements using DIC techniques. The DIC images were recorded using two Point Grey 4.1 MP cameras. Depending on the size of the specimens, different types of lenses were utilized, resulting in various image resolutions while a virtual strain gauge length (VSGL) of 0.3 mm was used for all the tests (Table 3). More information regarding the VSGL size and its influence on the DIC strain measurements can be found in Reu [57] and Rahmaan et al. [58]. The commercial VIC3D DIC (7.2.4) software from Correlated Solutions, Inc. was used for the DIC analysis. Images in each experiment were recorded with a prescribed frequency to provide at least 300 images from the initial deformation up to fracture. For the in-plane tensile tests, virtual extensometers with initial lengths of 20 mm and 11 mm were used to report far-field displacements of the hole tension and notch specimens, respectively. Note that for all the experiments reported in this paper, at least four specimens were tested to ensure the repeatability of the tests.

Strain Measurements
The experiments described above were recorded using two digital cameras to obtain full-field stereoscopic logarithmic strain measurements using DIC techniques. The DIC images were recorded using two Point Grey 4.1 MP cameras. Depending on the size of the specimens, different types of lenses were utilized, resulting in various image resolutions while a virtual strain gauge length (VSGL) of 0.3 mm was used for all the tests (Table 3). More information regarding the VSGL size and its influence on the DIC strain measurements can be found in Reu [57] and Rahmaan et al. [58]. The commercial VIC3D DIC (7.2.4) software from Correlated Solutions, Inc. was used for the DIC analysis. Images in each experiment were recorded with a prescribed frequency to provide at least 300 images from the initial deformation up to fracture. For the in-plane tensile tests, virtual extensometers with initial lengths of 20 mm and 11 mm were used to report far-field displacements of the hole tension and notch specimens, respectively. Note that for all the experiments reported in this paper, at least four specimens were tested to ensure the repeatability of the tests.

Finite Element Models
To assess predictability of the plasticity approach, finite element models of each of the test geometries were created, as depicted in Figure 6 and were simulated in LS-DYNA. Eight-node linear brick elements were used for discretization of only one-quarter of the geometries due to symmetry and appropriate symmetry boundary conditions were applied to simulate the behavior of the full specimens. For simplicity, the regions clamped inside the grips were not considered for the finite element model of the in-plane tensile specimens. To capture through-thickness gradients of stress and strain fields, five elements were used through the half thickness of the models resulting in an average element size of approximately 0.15 mm. This modelling methodology agrees with the numerical guidelines of Dunand and Mohr [49] for simulations of tensile-based test specimens. Smaller element sizes were also examined that resulted in remarkable computational time increase with insignificant changes in the global response of the models. brick elements were used for discretization of only one-quarter of the geometries due to symmetry and appropriate symmetry boundary conditions were applied to simulate the behavior of the full specimens. For simplicity, the regions clamped inside the grips were not considered for the finite element model of the in-plane tensile specimens. To capture through-thickness gradients of stress and strain fields, five elements were used through the half thickness of the models resulting in an average element size of approximately 0.15 mm. This modelling methodology agrees with the numerical guidelines of Dunand and Mohr [49] for simulations of tensile-based test specimens. Smaller element sizes were also examined that resulted in remarkable computational time increase with insignificant changes in the global response of the models.  For the Nakazima dome test simulations (see Figure 7), a similar modelling approach was employed with a solid element size of 0.15 mm within the blank. The tooling (punch, blank holder and die) were modelled as rigid bodies and a penalty-based contact algorithm with a constant Coulomb friction coefficient of 0.05 was used between the punch and the blank. This friction coefficient was taken from [59] where twist-compression tests were used to experimentally determine the friction coefficient of ZEK100-O using Teflon ® film lubrication at room temperature. Moreover, a friction coefficient of 0.4 with the same penalty-based contact algorithm was enforced between the blank and the die and blank holder (unlubricated). A constant blank holder force of 165 kN (quarter-symmetry) was applied to the binder to simulate the experimental clamping condition. In Section 6, the results of the finite element simulations are compared with the experimental results in terms of local strain paths, as well as the far-field load-displacements. the friction coefficient of ZEK100-O using Teflon film lubrication at room temperature. Moreover, a friction coefficient of 0.4 with the same penalty-based contact algorithm was enforced between the blank and the die and blank holder (unlubricated). A constant blank holder force of 165 kN (quartersymmetry) was applied to the binder to simulate the experimental clamping condition. In Section 6, the results of the finite element simulations are compared with the experimental results in terms of local strain paths, as well as the far-field load-displacements.

Yield Function and Plastic Potential
The coefficients of the yield function and plastic potential calibrated for ZEK100-O at different levels of plastic work are given in Tables 4 and 5 respectively. These coefficients were obtained with the calibration approach described in Section 3.1. Note that for both the yield function and plastic potential, C11 and C′11 coefficients are fixed and equal to 1.0 as suggested by Cazacu et al. [16]. Also a fixed exponent of a = 8.0 was utilized in the models. The other fixed coefficients are K and Kʹ for the "plastic potential" that are equal to zero due to the assumed symmetry explained in Section 3.1.
Due to brevity, only three levels of plastic deformation were chosen (plastic work levels of 2.

Yield Function and Plastic Potential
The coefficients of the yield function and plastic potential calibrated for ZEK100-O at different levels of plastic work are given in Tables 4 and 5 respectively. These coefficients were obtained with the calibration approach described in Section 3.1. Note that for both the yield function and plastic potential, C 11 and C 11 coefficients are fixed and equal to 1.0 as suggested by Cazacu et al. [16]. Also a fixed exponent of a = 8.0 was utilized in the models. The other fixed coefficients are K and K' for the "plastic potential" that are equal to zero due to the assumed symmetry explained in Section 3.1.   It can be seen from Figures 8-10 that the CPB06ex2 yield function and plastic potential correlate with the measured data with good accuracy. For the smallest plastic work level of 2.24 MPa, which is close to the initial yielding of the material, it can be seen from Figure 8 that the material has a clear tension-compression asymmetry with the tension region having larger yield stresses than the compression region. In the literature, this behavior has been attributed to plastic deformation mechanisms where slip-dominated mechanisms are responsible for plastic deformation under tension while compression leads to twinning-dominated mechanisms with lower critical resolved shear stress (CRSS) [9,10,60]. The tension-compression asymmetry of the yield loci reduces with deformation due to the high hardening rate in compression, offsetting the initially lower compressive yield strength with respect to tension as shown in Figure 9 for the plastic work level of 14.61 MPa. At the plastic work level of 22.46 MPa (Figure 10), the yield stresses in compression have grown larger than in tension which is opposite to that observed for the material at the onset of yielding. Note that this material exhibits asymmetry not only in the tension-compression regions but also in the shear regions represented by the second and fourth quadrants of yield loci which is correlated to the sign of the principal stresses acting under the shear state, activating different types of slip or twinning deformation mechanisms [12]. Thus, Figures 8-10 highlight the complexities associated with the material in terms of modelling and characterization with its significant anisotropy and asymmetry that evolve with plastic deformation  It can be seen from Figures 8-10 that the CPB06ex2 yield function and plastic potential correlate with the measured data with good accuracy. For the smallest plastic work level of 2.24 MPa, which is close to the initial yielding of the material, it can be seen from Figure 8 that the material has a clear tension-compression asymmetry with the tension region having larger yield stresses than the compression region. In the literature, this behavior has been attributed to plastic deformation mechanisms where slip-dominated mechanisms are responsible for plastic deformation under tension while compression leads to twinning-dominated mechanisms with lower critical resolved shear stress (CRSS) [9,10,60]. The tension-compression asymmetry of the yield loci reduces with deformation due to the high hardening rate in compression, offsetting the initially lower compressive yield strength with respect to tension as shown in Figure 9 for the plastic work level of 14.61 MPa. At the plastic work level of 22.46 MPa (Figure 10), the yield stresses in compression have grown larger than in tension which is opposite to that observed for the material at the onset of yielding. Note that this material exhibits asymmetry not only in the tension-compression regions but also in the shear regions represented by the second and fourth quadrants of yield loci which is correlated to the sign of the principal stresses acting under the shear state, activating different types of slip or twinning deformation mechanisms [12]. Thus, Figures 8-10 highlight the complexities associated with the material in terms of modelling and characterization with its significant anisotropy and asymmetry that evolve with plastic deformation.

Evaluation of the Model Using Single-Element Simulations
In order to evaluate the accuracy of the implemented plasticity model, simulations utilizing a single 3-D brick element under different proportional loading conditions were used to predict the stress-strain response and evolution of r-values of the material. An important advantage of singleelement simulations is that any desired stress state can be achieved without structural complexities such as buckling, plastic instability and necking. Furthermore, single-element simulations can better highlight deficiencies in the implementation of the models as opposed to full-scale simulations where such potential flaws might be smeared out. Therefore, single-element simulations can provide firstorder evaluations of the implemented models. Note that full-scale assessments of the model are also considered in Section 6.3. Figure 11 includes comparisons between predictions of the single-element model with experimental results under various stress states with the evolving non-AFR model providing excellent agreement with the experimental trends. This fidelity is expected due to the calibration of the model in these stress states. It can be observed from Figure 11a-c that the tension-compression asymmetric responses of the material are reproduced to high accuracy. Comparing the model

Evaluation of the Model Using Single-Element Simulations
In order to evaluate the accuracy of the implemented plasticity model, simulations utilizing a single 3-D brick element under different proportional loading conditions were used to predict the stress-strain response and evolution of r-values of the material. An important advantage of single-element simulations is that any desired stress state can be achieved without structural complexities such as buckling, plastic instability and necking. Furthermore, single-element simulations can better highlight deficiencies in the implementation of the models as opposed to full-scale simulations where such potential flaws might be smeared out. Therefore, single-element simulations can provide first-order evaluations of the implemented models. Note that full-scale assessments of the model are also considered in Section 6.3. Figure 11 includes comparisons between predictions of the single-element model with experimental results under various stress states with the evolving non-AFR model providing excellent agreement with the experimental trends. This fidelity is expected due to the calibration of the model in these stress states. It can be observed from Figure 11a-c that the tension-compression asymmetric responses of the material are reproduced to high accuracy. Comparing the model predictions with the experimental data, equal-biaxial tension and shear states are in good agreement as shown in Figure 11d,e, respectively. Moreover, the evolutions of the r-values are well described to reasonable accuracy in Figure 11f.  Figures 12 and 13 compare the measured and predicted load-displacement responses of the inplane tensile and Nakazima dome tests, respectively. The experimental data is presented up to appearance of the first visible surface crack on the DIC image which is followed by a sudden load drop. Furthermore, the predicted and measured (DIC) contours of major and minor strains at the last image before the onset of fracture are presented in Figure 14. For brevity, only the strain contours of the tests and simulations in the rolling direction are shown in Figure 14 while the strain paths to fracture for all test orientations are compared in Figure 15 in terms of major versus minor strains. In general, it can be seen that the agreement between the model and experiments are reasonably good

Evaluation of the Model for Full-Scale Simulation
Full-scale finite element simulations of the experiments described in Section 4 were also performed to evaluate the predictive capabilities of the current plasticity model. Note that these tests were not utilized for calibration of the model and can be considered as evaluation tests. In these tests, the material undergoes different stress states that are commonly experienced in practical sheet metal forming operations and crash events. The experimental and simulation data are compared in terms of global load-displacements as well as local principal strain paths at the critical locations where fracture is expected to initiate. The fracture locations are at the center of the gauge area of all specimens except the hole tension test where fracture initiates from edges of the hole. To report displacements of the in-plane tensile tests, the gauge lengths defined in Section 4.3 (20 mm and 11 mm for the hole tension and notch tension tests, respectively) were used consistently in the experiments and simulations. Furthermore, the local strains measured in the DIC software at the fracture initiation locations were obtained using a "box inspector" tool in VIC3D with a dimension of 0.3 mm × 0.3 mm. The DIC software reports the average of strains measured within the box and the same approach was used in the finite element models, that is, strains of surface elements within a 0.3 mm × 0.3 mm square were averaged and reported. Figures 12 and 13 compare the measured and predicted load-displacement responses of the in-plane tensile and Nakazima dome tests, respectively. The experimental data is presented up to appearance of the first visible surface crack on the DIC image which is followed by a sudden load drop. Furthermore, the predicted and measured (DIC) contours of major and minor strains at the last image before the onset of fracture are presented in Figure 14. For brevity, only the strain contours of the tests and simulations in the rolling direction are shown in Figure 14 while the strain paths to fracture for all test orientations are compared in Figure 15 in terms of major versus minor strains. In general, it can be seen that the agreement between the model and experiments are reasonably good for different test geometries and orientations.
It can be seen from Figure 12 that the correlation between the model and experiments are good for the in-plane tensile tests up to the points of maximum force. After this point, the accuracy of the model deteriorates with further deformation which leads to an over prediction of the applied forces. Given the complexity in the plastic response of the material, this discrepancy can be attributed to the hardening function or to the extent of evolution of the plasticity model. As mentioned earlier, the developed plasticity model evolves up to the limiting plastic work level of 22.46 MPa (9% plastic strain for uniaxial tension test in the RD) and then the shapes of yield function and plastic potential remain constant with further plastic deformation. It is interesting to note that the evolution of the r-values exhibits a tendency to saturate with deformation (see Figure 1d) that is believed to be the underlying reason for the good agreement between the predicted and experimental strain paths in Figure 15. Furthermore, it can be seen from Figure 14 that the contours of predicted strains correlate well with the DIC results. Nevertheless, it is expected that the anisotropy evolution and distortional hardening of the material persists in terms of the yield function after the onset of necking in the tensile tests but this could not be characterized experimentally. The evolution of the plastic potential is anticipated to be less significant. Moreover, it should be mentioned that in the tensile tests, upon localization the stress state moves towards a plane-strain condition and unfortunately the experimental stress response for this state is unknown since plane-strain tension tests were not performed in this study. Extracting the plane-strain behavior of materials is not straightforward because achieving a uniform plane-strain condition is experimentally difficult (see for instance [61]). It is expected that calibrating the model with accurate plane-strain data could improve post-localization prediction of the developed plasticity model. experimental stress response for this state is unknown since plane-strain tension tests were not performed in this study. Extracting the plane-strain behavior of materials is not straightforward because achieving a uniform plane-strain condition is experimentally difficult (see for instance [61]). It is expected that calibrating the model with accurate plane-strain data could improve postlocalization prediction of the developed plasticity model.   performed in this study. Extracting the plane-strain behavior of materials is not straightforward because achieving a uniform plane-strain condition is experimentally difficult (see for instance [61]). It is expected that calibrating the model with accurate plane-strain data could improve postlocalization prediction of the developed plasticity model.  It can be observed from Figure 13a that the punch force-displacement predictions from the finite element model are in accordance with experimental results for the plane-strain dome tests in both orientations of the RD and TD. Furthermore, the model predicts the strain paths of these tests with good accuracy as shown in Figure 15. However, the model over predicts the punch force for the equal-biaxial dome test by approximately 12% as shown in Figure 13b. From the strain path of the equal-biaxial test in Figure 15, it appears that the model does not accurately reproduce the experimental data in terms of change of the strain path due to the localization. A large region of the dome specimen experiences the equal-biaxial tensile condition and in the plasticity model, this state of stress should be calibrated to high accuracy so that precise load-displacements can be obtained from the equal-biaxial dome test simulations. It can be seen from Figure 11d that the plasticity model is well-calibrated under this state of stress. Therefore, the discrepancy between the model and experiments might be related to the type of the characterization test (through-thickness compression test) that was used for the equal-biaxial condition. Obtaining the equal-biaxial tensile stress-strain response of sheet materials is difficult and the through-thickness compression test (e.g., [11,38]) is used as an alternative to cruciform (e.g., [62]) or hydraulic bulge (e.g., [63]) tests. It should be noted that the cruciform and bulge tests have their own challenges, for instance, the cruciform test is only suitable for small plastic deformation and it requires complex specimen preparation [64] and the stress state in bulge tests is only approximately equal-biaxial [65]. Keeping in mind these discrepancies, further investigations with other biaxial testing methods of cruciform and bulge tests are required to provide a clearer picture of the equal-biaxial response of the material and will be considered in future work. It can be observed from Figure 13a that the punch force-displacement predictions from the finite element model are in accordance with experimental results for the plane-strain dome tests in both orientations of the RD and TD. Furthermore, the model predicts the strain paths of these tests with good accuracy as shown in Figure 15. However, the model over predicts the punch force for the equalbiaxial dome test by approximately 12% as shown in Figure 13b. From the strain path of the equalbiaxial test in Figure 15, it appears that the model does not accurately reproduce the experimental data in terms of change of the strain path due to the localization. A large region of the dome specimen experiences the equal-biaxial tensile condition and in the plasticity model, this state of stress should be calibrated to high accuracy so that precise load-displacements can be obtained from the equalbiaxial dome test simulations. It can be seen from Figure 11d that the plasticity model is wellcalibrated under this state of stress. Therefore, the discrepancy between the model and experiments might be related to the type of the characterization test (through-thickness compression test) that was to cruciform (e.g., [62]) or hydraulic bulge (e.g., [63]) tests. It should be noted that the cruciform and bulge tests have their own challenges, for instance, the cruciform test is only suitable for small plastic deformation and it requires complex specimen preparation [64] and the stress state in bulge tests is only approximately equal-biaxial [65]. Keeping in mind these discrepancies, further investigations with other biaxial testing methods of cruciform and bulge tests are required to provide a clearer picture of the equal-biaxial response of the material and will be considered in future work.

Influence of Evolving Anisotropy on the Predicted Structural Response
One convenient approach to consider anisotropy is to assume that the shapes of yield function and plastic potential remain constant with deformation. This approach is typically sufficient for bcc and fcc cubic materials; however, this might not be the case for the hcp materials. To analyze the importance of considering the evolving anisotropy, evolving and non-evolving models are compared for the current magnesium alloy. Evolution of the plasticity model was disabled by considering anisotropy to be constant after the onset of yielding which corresponds to the lowest plastic work level of 2.24 MPa (see Figure 8). Simulations of the hole tension tests were performed assuming that the shape of the yield function and plastic potential remain constant with deformation. Figure 16a compares the load-displacement predictions of the non-evolving model with experimental results from which it can be seen that, as expected, only the initial yielding response is well-captured by the model for the TD and DD directions and the hardening rates are not in agreement with the experimental data. Furthermore, the predicted strain paths for all orientations shown in Figure 16b display significant discrepancies with experimental data.
In addition, it can be seen from Figure 16a that simulations with the non-evolving model in the DD and TD directions lead to significantly earlier peak loads than observed in the experiments due to premature localization as a result of lower hardening rates. This behavior is expected to have significant influence on fracture predictions, for instance, using popular phenomenological approaches (e.g., [66,67]) where the stress triaxiality is used to monitor the variations of stress states

Influence of Evolving Anisotropy on the Predicted Structural Response
One convenient approach to consider anisotropy is to assume that the shapes of yield function and plastic potential remain constant with deformation. This approach is typically sufficient for bcc and fcc cubic materials; however, this might not be the case for the hcp materials. To analyze the importance of considering the evolving anisotropy, evolving and non-evolving models are compared for the current magnesium alloy. Evolution of the plasticity model was disabled by considering anisotropy to be constant after the onset of yielding which corresponds to the lowest plastic work level of 2.24 MPa (see Figure 8). Simulations of the hole tension tests were performed assuming that the shape of the yield function and plastic potential remain constant with deformation. Figure 16a compares the load-displacement predictions of the non-evolving model with experimental results from which it can be seen that, as expected, only the initial yielding response is well-captured by the model for the TD and DD directions and the hardening rates are not in agreement with the experimental data. Furthermore, the predicted strain paths for all orientations shown in Figure 16b display significant discrepancies with experimental data.
In addition, it can be seen from Figure 16a that simulations with the non-evolving model in the DD and TD directions lead to significantly earlier peak loads than observed in the experiments due to premature localization as a result of lower hardening rates. This behavior is expected to have significant influence on fracture predictions, for instance, using popular phenomenological approaches (e.g., [66,67]) where the stress triaxiality is used to monitor the variations of stress states at critical locations to predict a triaxiality-dependent fracture strain and trigger element deletion. Early localization of the non-evolving model directly affects the stress triaxiality at the critical element, as shown in Figure 17 for the hole tension test in the TD direction, for example. This test direction was selected since the model exhibits the largest discrepancy from the experimental hardening rates. Compared to the evolving model, high stress triaxialities and plastic strains at low global deformation are apparent for the non-evolving model as observed in Figure 17. These large differences show that employing evolving plasticity models to capture the behavior of magnesium alloys is necessary and highlight its importance for simulations (at the cost of inevitable increased model complexity and computational cost).
at critical locations to predict a triaxiality-dependent fracture strain and trigger element deletion. Early localization of the non-evolving model directly affects the stress triaxiality at the critical element, as shown in Figure 17 for the hole tension test in the TD direction, for example. This test direction was selected since the model exhibits the largest discrepancy from the experimental hardening rates. Compared to the evolving model, high stress triaxialities and plastic strains at low global deformation are apparent for the non-evolving model as observed in Figure 17. These large differences show that employing evolving plasticity models to capture the behavior of magnesium alloys is necessary and highlight its importance for simulations (at the cost of inevitable increased model complexity and computational cost).  at critical locations to predict a triaxiality-dependent fracture strain and trigger element deletion. Early localization of the non-evolving model directly affects the stress triaxiality at the critical element, as shown in Figure 17 for the hole tension test in the TD direction, for example. This test direction was selected since the model exhibits the largest discrepancy from the experimental hardening rates. Compared to the evolving model, high stress triaxialities and plastic strains at low global deformation are apparent for the non-evolving model as observed in Figure 17. These large differences show that employing evolving plasticity models to capture the behavior of magnesium alloys is necessary and highlight its importance for simulations (at the cost of inevitable increased model complexity and computational cost).

Discussion
This research has demonstrated that the plastic anisotropy of ZEK100-O is extremely complex. Considering all the challenges associated with the asymmetry and anisotropy of the material, the current evolving non-associative model was able to capture the behavior of this alloy with good accuracy in practical stress states for forming operations. It should be mentioned that models with AFR have also been successfully used in the past to capture the response of magnesium alloys. For instance, Ghaffari Tari et al. [11] employed up to four linear transformations for the CPB06 model with an AFR to be able to reproduce behavior of AZ31B. Such a large number of transformations can significantly complicate implementation of models into FE codes and increase computational time as described in Reference [11]. The advantages of the framework developed in the present paper are simplified implementation, calibration and improved accuracy across a comprehensive range of stress states and loading conditions. Note that compared to the model developed in this paper, the recent two-surface model of Kondori et al. [30] has fewer parameters and was used successfully in uniaxial tension and compression.
Despite the capabilities of the present model discussed in Section 6, there are a few limitations that are remaining to be addressed in future work. This section summarizes the strengths and weaknesses of the model.

•
The proposed plasticity model was developed for proportional loading conditions and it was shown that the model performs well under these conditions; however, the performance of the model under severe non-proportional conditions has to be evaluated in future.

•
Consideration of the evolution of anisotropy was shown to be critical for global response and local strain/stress histories. A limitation of the approach is that the evolving anisotropy is restricted to the onset of necking at relatively low strain levels after which the anisotropy evolution was assumed to saturate. New experimental techniques that can suppress necking are required to characterize and expose the behavior of the material at larger strains under proportional conditions. • The work hardening rate under equal-biaxial tension was obtained using the through-thickness compression test and should be revisited by performing cruciform or bulge tests. Re-calibrating the model with the results of these alternative tests might improve its accuracy under biaxial-dominant states such as in the equal-biaxial Nakazima dome test.

•
This research was focused on the room temperature behavior of ZEK100-O under quasi-static conditions. Forming operations of magnesium alloys may also be performed at elevated temperatures and strain rates but the severity of anisotropy and asymmetry of the material is reduced by increasing the temperature [68,69], thus it is expected that the present modelling framework will perform well at elevated temperatures. However, this needs further investigations that will be considered in future work.

Conclusions
The room temperature constitutive plastic behavior of a rare-earth magnesium alloy sheet, ZEK100-O, was studied under different stress states. It was demonstrated that the material exhibits severe plastic anisotropy that evolves with deformation. It was shown that adopting the non-AFR along with an evolving CPB06 formulation with two linear stress transformations resulted in good agreement between the finite element predictions and experimental data. The strategy adopted to consider evolution was piece-wise linear interpolation between yield functions and plastic potentials that were calibrated at different plastic work levels. The model was implemented for general 3-D stress states to capture the through-thickness variation of stress and strain fields that is essential for accurate fracture prediction of materials. Full-scale simulations of different experiments including hole tension, notch tension and Nakazima dome tests showed that predictions of the model were in good agreement with the experimental results except for the equal-biaxial dome test simulations where the punch load was over predicted. Finally, it was shown that considering the evolution of yield function and plastic potential in finite element models is critical to precisely capture experimental load-displacements and local strains.