Tests and Finite Element Simulation of Yield Anisotropy and Tension-Compression Strength Difference of an Extruded ZK60 Mg Alloy

: Multiroute yield experiments were carried out on extruded ZK60 magnesium alloy samples under 200 ◦ C and strain rate of 10 − 2 s − 1 by a thermomechanical simulator. A hypoelastoplastic large deformation constitutive relationship was employed to simulate large deformation yield subjected to uniaxial loading, biaxial loading, loading-reverse loading, and orthogonal loading in the ﬁnite element (FE) software with user subroutine. The results showed the following: (1) As the accumulative plastic strain increased, the convex yield surface in the 2D stress space gradually expanded or shrank, rotated, and distorted from the approximate ellipse at low accumulative plastic strain. At the same accumulative plastic strain, uniaxial tensile or compressive yield values along different directions were not the same. In addition, the tensile yield value differed considerably from the compressive one. The large deformation yield of ZK60 magnesium alloy showed signiﬁcant anisotropy, tension-compression strength difference, and evolution effect. (2) FE simulations based on the hypoelastoplastic constitutive relationship could accurately capture the strongly evolving asymmetric yield behavior under complex loading routes. The stress-strain relationships and hardening rates were in better accordance with the experimental results and reﬂected the yield behavior more realistically compared to simulations without the evolution effect or with traditional linear interpolation. Deformation at low temperature and high strain rate is important to save process costs and improve processing efﬁciency, but it may cause evolving asymmetric yield during large deformation. It is believed that the simulation approach used herein is reliable for the prediction and optimization of severe plastic deformation processes of hexagonal close-packed (HCP) alloys.


Introduction
Due to low density, high specific strength, good biocompatibility, and abundant reserves [1][2][3], magnesium and its alloys have attracted widespread attention for use in vehicles, aerospace, military equipment, intelligent technology, and medical devices [4][5][6][7]. As a typical lightweight hexagonal close-packed (HCP) material, magnesium alloys are praised as the most promising metal material in the 21st century [8]. However, the drawbacks, such as poor low-temperature formability, cumbersome processing procedure, and strong dependence on rare and precious elements, substantially restrict their extensive application. Developing low-cost and high-efficiency severe plastic deformation processes for plainified alloys is an effective way to better promote the application of magnesium in various industrial fields. Some severe plastic forming processes, such as high-ratio differential speed rolling [9][10][11][12], equal-channel angular pressure [13,14], friction stir processing [15,16], and high-pressure torsion [17,18], can achieve uniform and fine equiaxed grains, improving the forming and service performance. However, owing to the high cost (high temperature) and limited sample dimensions (small samples), these processes are only utilized in laboratories and have not been popularized in industrial production. 2 of 19 To effectively reduce the cost and improve processing efficiency, low-temperature and high-speed processes are needed. It is worth noting that at high temperature [19] and low strain rate, more slip systems, such as pyramidal and secondary pyramidal slips, are prone to be activated to coordinate deformation. However, at low temperature, high strain rate, and large deformation, magnesium alloys may exhibit strongly evolving yield anisotropy and tension-compression strength difference effect (SDE), making it difficult to predict and control the product profile and texture.
At low temperature, basal slip is almost the only active slip system of Mg and its alloys [20]. As the deformation accumulates, a strong basal texture is gradually formed, leading to significant yield anisotropy [21] and even earing defects [22][23][24]. While twinning initiates and coordinates deformation, and its directionality causes SDE [23,25]. Li et al. [26] indicated that during large deformation, stress-induced microstructure activities, including slip, twinning, dynamic recrystallization, and reorientation, lead to clear crystal rotation and continuous texture evolution with internal variables. The evolution of texture can cause continuous rotation and distortion of yield surfaces. Therefore, magnesium and its alloys exhibit strongly evolving yield anisotropy and SDE, which is unfavorable for material forming prediction and control. In continuous large deformation simulations, the strong evolution of yield surfaces cannot be ignored. To accurately capture the strongly evolving yield anisotropy and SDE at low temperature, high strain rate, and large deformation, it is necessary to develop a reasonable asymmetric yield function, large deformation constitutive relationship, and yield surface evolution description.
To simulate yield stress and ratio of strain rate in different directions (R-value), Barlat et al. [27], Banabic et al. [28], Barlat et al. [29], and Lou et al. [30] developed a variety of anisotropic yield functions that are mainly used for steel and aluminum sheets. With respect to HCP metals and their alloys, Plunkett et al. [31], Nixon et al. [32], and Yoon et al. [33] simulated the asymmetric yield with anisotropy and SDE of Zr, Ti, and Mg, respectively. Based on the above works, Yoon et al. [34] also developed an asymmetric yield function that can accurately describe hydrostatic pressure sensitivity and applied it to describe the yield of Zr and Ti. However, the evolution of asymmetric yield with internal variables caused by a strong microstructure and texture evolution has rarely been studied. The evolution, which may lead to great errors, should be theoretically and numerically described in large deformation simulations.
The phenomenological yield function and its constitutive description can only consider the yield function evolving with internal variables to describe the evolution effect. Mekonen et al. [19] and Tari et al. [35] set the yield model parameters as a monotonic function of accumulative plastic strain to model the evolution, which may lead to accumulated simulation error [19] at a large strain. Therefore, this approach is not suitable for accurate large deformation finite element (FE) simulation of HCP metals and their alloys with strong and nonmonotonic evolution effects. For example, at a large plastic strain, in-plane compression Zr presents a strong and nonmonotonic yield surface evolution due to the clear strain softening on the X-Y biaxial tension route [34]. Plunkett et al. [31] used a linear interpolation approach to describe yield surface evolution. However, the accuracy of the linear interpolation approach is strongly dependent on the gap between the interpolation points and whether the evolution is linear. Li et al. [27] innovatively proposed a standard continuum-based discontinuous (CBD) framework that can discontinuously describe the torsion and evolution of plasticity. In the FE simulation, the constitutive relationship of small deformation, the associated flow law, and the implicit integration algorithm were used to verify the validity of the CBD framework. However, it was also pointed out that the current CBD constitutive framework requires further extensive evaluation, development, and application in terms of describing complex forming conditions with space, abrupt strain route, or thermal-mechanical coupling load and greater nonlinear plastic behavior. Xie et al. [36] developed a hypoelastoplastic large deformation constitutive relationship that combines the evolution description based on stress states and internal variables, the asymmetric yield function, the nonassociated flow law, and the explicit incremental inte-gration algorithm, thereby accurately simulating the evolving asymmetric yield behaviors of Zr under complex loading routes.
The aims of the current study were to (1) reveal the significant yield anisotropy, SDE, and yield surface evolution behavior of plain ZK60 Mg alloy under low temperature and high strain rate and (2) validate the reliability of the approach used in the current study for the simulation of evolving large deformation asymmetric yield of Mg alloys.

Materials
The plain ZK60 magnesium alloy used in this study was received in the form of as-extruded cuboids. Its elasticity was assumed to be isotropic in the FE simulation. The Young's modulus was 43 GPa, and the Poisson's ratio was 0.35. Its chemical composition is shown in Table 1.

Electron Backscatter Diffraction Characterization
The electron backscatter diffraction (EBSD) characterization surface is the extrusion direction-normal direction (ED-ND) surface of the ZK60 magnesium alloy cuboid. The surface was first ground up to 3000 grit using SiC papers and then polished with diamond paste with the sequence of W3.5, W2.5, W1, W0.5, and W0.25. After being polished, the surface was etched using a mixture solution containing 10 mL HNO 3 , 30 mL acetic acid, 40 mL H 2 O, and 120 mL alcohol. Finally, the EBSD data were measured by a SU5000 Hitachi scanning electron microscope (SEM) (Hitachi High-Technologies GLOBAL, Tokyo, Japan) equipped with the Oxford Aztec Nordlys Max data acquisition system, and the microstructure and texture were analyzed using the HKL Channel 5 system (5.11.20405.0, Oxford Instruments, Abingdon, UK). The SEM and EBSD operation voltage, spot intensity, magnification, and step size were 20 KV, 70 pA, 150, and 0.5 µm, respectively. An inverse pole figure (IPF) map for the initial microstructure of the extruded ZK60 magnesium alloy is shown in Figure 1, and the pole figure maps are shown in Figure 2. The results showed that the average grain size of the material was 4.15 µm, the standard deviation was 3.25, the majority of c-axes were slightly inclined at around 3.35 • with respect to the ND, and the maximum intensity was 18.

Compression and Tension Tests
For convenience, the transverse direction (TD), the ED, and the 45 • directions between the TD and ED are hereafter written as X direction, Y direction, and 45 • direction, respectively. Wire cutting was performed on the original cuboid to obtain dog-bone-shaped tensile samples (1, 2, and 3 in Figure 3) and cylinder-shaped compressive samples (4, 5, and 6 in Figure 3). A temperature of 200 • C and a strain rate of 10 −2 s −1 [37] were used on the samples.

Compression and Tension Tests
For convenience, the transverse direction (TD), the ED, and the 45° directions between the TD and ED are hereafter written as X direction, Y direction, and 45° direction, respectively. Wire cutting was performed on the original cuboid to obtain dog-bone-shaped tensile samples (1, 2, and 3 in Figure 3) and cylinder-shaped compressive samples (4, 5, and 6 in Figure 3). A temperature of 200 °C and a strain rate of 2 1 10 s − − [37] were used on the samples.

Compression and Tension Tests
For convenience, the transverse direction (TD), the ED, and the 45° directions between the TD and ED are hereafter written as X direction, Y direction, and 45° direction, respectively. Wire cutting was performed on the original cuboid to obtain dog-bone-shaped tensile samples (1, 2, and 3 in Figure 3) and cylinder-shaped compressive samples (4, 5, and 6 in Figure 3). A temperature of 200 °C and a strain rate of 2 1 10 s − − [37] were used on the samples.   was reduced from 12 to 4 mm, and the tensile tests were performed until the samples were completely cracked. To minimize experimental error, the test for each direction was repeated in triplicate. The dog-bone-shaped tensile and cylinder-shaped compressive samples conformed to the Chinese national standard. The compressive and tensile samples did not slide during testing. The friction [38][39][40] that occurred during the testing was negligible and was not relevant for the purpose of FE simulation in this study. The primary force F versus displacement ΔL curves were calculated to the true tress σ and true strain ε curves. The equations are as follows: where 0 L represents the original length of the sample gauge zone, and 0 A represents the original cross-sectional area of the sample gauge zone. Figure 4 shows the experimental compressive and tensile equivalent stress versus accumulative plastic strain curves. The compressive sample was a cylinder with a diameter of 8 mm and a height of 12 mm. The tensile sample had a gauge size of 30 mm × 4 mm × 2 mm (GB/T 228-2002 (Chinese National Standard/Recommended)). Before the tests, the samples were heated to 200 • C using the thermocouple and kept for two minutes. The tests were performed on a Gleeble 3500 thermomechanical simulator (DATASCIENCESINTERNATIONAL, INC, Minnesota, America) with a true strain rate of 10 −2 s −1 . The compressive sample height was reduced from 12 to 4 mm, and the tensile tests were performed until the samples were completely cracked. To minimize experimental error, the test for each direction was repeated in triplicate. The dog-bone-shaped tensile and cylinder-shaped compressive samples conformed to the Chinese national standard. The compressive and tensile samples did not slide during testing. The friction [38][39][40] that occurred during the testing was negligible and was not relevant for the purpose of FE simulation in this study. The primary force F versus displacement ∆L curves were calculated to the true tress σ and true strain ε curves. The equations are as follows: where L 0 represents the original length of the sample gauge zone, and A 0 represents the original cross-sectional area of the sample gauge zone. Figure 4 shows the experimental compressive and tensile equivalent stress versus accumulative plastic strain curves.

Large Deformation Constitutive Relationship
There is no transformation between phases of different density and a massive amount of void nucleation in magnesium alloy [41,42], so there is no plastic dilatancy or only negligible plastic dilatancy [43]. Therefore, the following asymmetric yield function does not contain the first invariant of stress that causes plastic dilatancy. Applying the analytical interpolation of the stress invariant asymmetric yield function * f , the yield function ** f evolving with the internal variables 1 n ξ − can be derived:

Large Deformation Constitutive Relationship
There is no transformation between phases of different density and a massive amount of void nucleation in magnesium alloy [41,42], so there is no plastic dilatancy or only negligible plastic dilatancy [43]. Therefore, the following asymmetric yield function does not contain the first invariant of stress that causes plastic dilatancy. Applying the analytical interpolation of the stress invariant asymmetric yield functionf * , the yield functionf * * evolving with the internal variablesξ 1−n can be derived: whereĴ 2 = 1 2ŝ :ŝ is the second stress invariant of the isotropic plastic equivalent transformed stress tensorŝ ;Ĵ 3 = detŝ is the third stress invariant of another isotropic plastic equivalent transformed stress tensorŝ ; the subscripts (j) and (j + 1) represent the starting and terminal points of the interpolation based on the given experimental data, respectively; n represents the hardening exponent; ) n ; the intersections between the proportional stress route and the two adjacent surfaces are the stress statesσ app(1) andσ app (2) ; and dε p is the accumulative plastic strain.
To introduce the anisotropic coefficients, two isotropic plastic equivalent transformed stress tensorsŝ andŝ are transformed from the stress tensorσ by two fourth-order linear transformation tensors of L and L , respectively: where the two fourth-order linear transformation tensors introduce a total of 12 anisotropic parameters and are written in matrix form as follows: In addition, the stress tensor can be written in matrix form as follows: For a two-dimensional (2D) principal stress plane, only six parameters c 1 , c 2 , c 3 , c 1 , c 2 , and c 3 are relevant [34]. Based on uniaxial tensile and compressive tests in the X and Y directions together with equibiaxial tension and compression, a nonlinear error function can be built and optimized to evaluate the six anisotropic parameters by the downhill simplex approach [44].
In the process of anisotropic plastic deformation of HCP alloys, the nonassociated flow criterion is required. According to the work conjugate pair of the rotating Kirchhoff stress tensorσ and the rotating deformation rate tensord, the plastic potentialĝ, the plastic part of the rotating deformation rate tensor d p , the equivalent plastic rotating deformation rated p , and elastic constitutive relationship are introduced as follows: where C represents the isotropic stiffness tensor.
The accumulative plastic strain dε p is considered as the sole internal variable. The plastic modulus h will be obtained by consideringK =f * =σ T 11 and the relationship between the accumulative plastic strain rate and the equivalent plastic rotational deformation rate for the small elastic strain or rigid body rotation-free conditions.
where dσ T 11 = dσ T 11 = Jdt T 11 represents the stress increment of X uniaxial tension that is a rigid body rotation-free process, and dε p represents the accumulative plastic strain increment corresponding to the stress increment. The expression of the partial derivative of the asymmetric yield function obtained by continuous interpolation to the accumulative plastic strain is as follows [36]:

FE Algorithm
In the current study, the forward Euler explicit integration scheme was employed for the hypoelastoplastic constitutive relationship in the user material subroutine, and the asymmetric yield function based on the invariants of isotropic plastic equivalent transformed stress tensors was used to describe anisotropic and tension-compression SDE. For the problems with a single internal variable [34,45], based on the hardening effect with the function of the internal variable and the stress state, an analytical interpolation approach was developed to build the continuous evolution between two adjacent yield surfaces [36].
In this study, we conducted a multiroute loading simulation of extruded ZK60 magnesium alloy using the commercial FE software Abaqus with user subroutine defining the evolving asymmetric yield [36]. The routes included X uniaxial tension, X uniaxial compression, Y uniaxial tension, Y uniaxial compression, X-Y biaxial tension, X tension-Y compression, X compression-Y tension, X-Y biaxial compression, loading-reverse loading (X uniaxial tension-X uniaxial compression and Y uniaxial tension-Y uniaxial compression), and orthogonal loading (X uniaxial tension-Y uniaxial tension and Y uniaxial tension-X uniaxial tension).
In this study, 1D/2D stress problems were employed to validate the large deformation constitutive relationship considering asymmetric yield evolution. Uniaxial and biaxial loading were first simulated by the finite element method (FEM). For uniaxial loading, a single element cube with dimensions of 10 mm × 10mm × 10mm was built, and the displacement boundary conditions ( Table 2) were imposed on the end faces of the cube. For biaxial loading, the same cubic model was applied, and the displacement loads (Table 2) along two directions possessed the same loading increment at each time step. In the above simulations, the 3D stress hexahedral element (C3D8) and the static analysis step were employed. The displacement contours of deformed specimens of extruded ZK60 magnesium alloy are shown in Table 3.   Table 3. Displacement contours of deformed specimens under uniaxial and biaxial loading (unit: mm).  Table 3. Displacement contours of deformed specimens under uniaxial and biaxial loading (unit: mm).

Yield Surfaces
Experimental data in Table 4 shows the differences in tension and compression strength and the anisotropy. At an accumulative strain of 0.005, the X uniaxial tension-compression yield difference was 15.5 MPa and the relative difference reached 47.69%; the Y uniaxial tension-compression yield difference was 4.0 MPa and the relative difference reached 11.76%; the X uniaxial tension and Y uniaxial tension yield difference was 10.0 MPa and the relative difference reached 26.32%; and the X uniaxial compression and Y uniaxial compression yield difference was 1.5 MPa and the relative difference reached 4.62%. At an accumulative strain of 0.02, the X uniaxial tension-compression yield difference was 9.0 MPa and the relative difference reached 8.61%; the Y uniaxial tension-compression yield difference was 4.5 MPa and the relative difference reached  Table 3. Displacement contours of deformed specimens under uniaxial and biaxial loading (unit: mm).

Yield Surfaces
Experimental data in Table 4 shows the differences in tension and compression strength and the anisotropy. At an accumulative strain of 0.005, the X uniaxial tension-compression yield difference was 15.5 MPa and the relative difference reached 47.69%; the Y uniaxial tension-compression yield difference was 4.0 MPa and the relative difference reached 11.76%; the X uniaxial tension and Y uniaxial tension yield difference was 10.0 MPa and the relative difference reached 26.32%; and the X uniaxial compression and Y uniaxial compression yield difference was 1.5 MPa and the relative difference reached 4.62%. At an accumulative strain of 0.02, the X uniaxial tension-compression yield difference was 9.0 MPa and the relative difference reached 8.61%; the Y uniaxial tension-compression yield difference was 4.5 MPa and the relative difference reached

Yield Surfaces
Experimental data in Table 4 shows the differences in tension and compression strength and the anisotropy. At an accumulative strain of 0.005, the X uniaxial tension-compression yield difference was 15.5 MPa and the relative difference reached 47.69%; the Y uniaxial tension-compression yield difference was 4.0 MPa and the relative difference reached 11.76%; the X uniaxial tension and Y uniaxial tension yield difference was 10.0 MPa and the relative difference reached 26.32%; and the X uniaxial compression and Y uniaxial compression yield difference was 1.5 MPa and the relative difference reached 4.62%. At an accumulative strain of 0.02, the X uniaxial tension-compression yield difference was 9.0 MPa and the relative difference reached 8.61%; the Y uniaxial tension-compression yield difference was 4.5 MPa and the relative difference reached

Yield Surfaces
Experimental data in Table 4 shows the differences in tension and compression strength and the anisotropy. At an accumulative strain of 0.005, the X uniaxial tension-compression yield difference was 15.5 MPa and the relative difference reached 47.69%; the Y uniaxial tension-compression yield difference was 4.0 MPa and the relative difference reached 11.76%; the X uniaxial tension and Y uniaxial tension yield difference was 10.0 MPa and the relative difference reached 26.32%; and the X uniaxial compression and Y uniaxial compression yield difference was 1.5 MPa and the relative difference reached 4.62%. At an accumulative strain of 0.02, the X uniaxial tension-compression yield difference was 9.0 MPa and the relative difference reached 8.61%; the Y uniaxial tension-compression yield difference was 4.5 MPa and the relative difference reached

Yield Surfaces
Experimental data in Table 4 shows the differences in tension and compression strength and the anisotropy. At an accumulative strain of 0.005, the X uniaxial tension-compression yield difference was 15.5 MPa and the relative difference reached 47.69%; the Y uniaxial tension-compression yield difference was 4.0 MPa and the relative difference reached 11.76%; the X uniaxial tension and Y uniaxial tension yield difference was 10.0 MPa and the relative difference reached 26.32%; and the X uniaxial compression and Y uniaxial compression yield difference was 1.5 MPa and the relative difference reached 4.62%. At an accumulative strain of 0.02, the X uniaxial tension-compression yield difference was 9.0 MPa and the relative difference reached 8.61%; the Y uniaxial tension-compression yield difference was 4.5 MPa and the relative difference reached

Yield Surfaces
Experimental data in Table 4 shows the differences in tension and compression strength and the anisotropy. At an accumulative strain of 0.005, the X uniaxial tension-compression yield difference was 15.5 MPa and the relative difference reached 47.69%; the Y uniaxial tension-compression yield difference was 4.0 MPa and the relative difference reached 11.76%; the X uniaxial tension and Y uniaxial tension yield difference was 10.0 MPa and the relative difference reached 26.32%; and the X uniaxial compression and Y uniaxial compression yield difference was 1.5 MPa and the relative difference reached 4.62%. At an accumulative strain of 0.02, the X uniaxial tension-compression yield difference was 9.0 MPa and the relative difference reached 8.61%; the Y uniaxial tension-compression yield difference was 4.5 MPa and the relative difference reached  Table 3. Displacement contours of deformed specimens under uniaxial and biaxial loading (unit: mm).

Yield Surfaces
Experimental data in Table 4 shows the differences in tension and compression strength and the anisotropy. At an accumulative strain of 0.005, the X uniaxial tension-compression yield difference was 15.5 MPa and the relative difference reached 47.69%; the Y uniaxial tension-compression yield difference was 4.0 MPa and the relative difference reached 11.76%; the X uniaxial tension and Y uniaxial tension yield difference was 10.0 MPa and the relative difference reached 26.32%; and the X uniaxial compression and Y uniaxial compression yield difference was 1.5 MPa and the relative difference reached 4.62%. At an accumulative strain of 0.02, the X uniaxial tension-compression yield difference was 9.0 MPa and the relative difference reached 8.61%; the Y uniaxial tension-compression yield difference was 4.5 MPa and the relative difference reached  Table 3. Displacement contours of deformed specimens under uniaxial and biaxial loading (unit: mm).

Yield Surfaces
Experimental data in Table 4 shows the differences in tension and compression strength and the anisotropy. At an accumulative strain of 0.005, the X uniaxial tension-compression yield difference was 15.5 MPa and the relative difference reached 47.69%; the Y uniaxial tension-compression yield difference was 4.0 MPa and the relative difference reached 11.76%; the X uniaxial tension and Y uniaxial tension yield difference was 10.0 MPa and the relative difference reached 26.32%; and the X uniaxial compression and Y uniaxial compression yield difference was 1.5 MPa and the relative difference reached 4.62%. At an accumulative strain of 0.02, the X uniaxial tension-compression yield difference was 9.0 MPa and the relative difference reached 8.61%; the Y uniaxial tension-compression yield difference was 4.5 MPa and the relative difference reached

Yield Surfaces
Experimental data in Table 4 shows the differences in tension and compression strength and the anisotropy. At an accumulative strain of 0.005, the X uniaxial tensioncompression yield difference was 15.5 MPa and the relative difference reached 47.69%; the Y uniaxial tension-compression yield difference was 4.0 MPa and the relative difference reached 11.76%; the X uniaxial tension and Y uniaxial tension yield difference was 10.0 MPa and the relative difference reached 26.32%; and the X uniaxial compression and Y uniaxial compression yield difference was 1.5 MPa and the relative difference reached 4.62%. At an accumulative strain of 0.02, the X uniaxial tension-compression yield difference was 9.0 MPa and the relative difference reached 8.61%; the Y uniaxial tension-compression yield difference was 4.5 MPa and the relative difference reached 5.59%; the X uniaxial tension and Y uniaxial tension yield difference was 28.5 MPa and the relative difference reached 33.53%; The X uniaxial compression and Y uniaxial compression yield difference was 24.0 MPa and the relative difference reached 29.81%. The extruded ZK60 magnesium alloy presented clear SDE and anisotropy. Based on the experimental data and the optimization approach, nine groups of anisotropic parameters (Appendix A and Table A1) corresponding to different accumulative plastic strains were calculated according to the method in Section 4.1 of [36], and the convex yield surfaces (the convexity can be validated according to the Appendix D of [36]) on the 2D principal stress plane are shown in Figure 5. As can be seen, the convex yield surfaces could accurately capture the experimental yield points on the six loading routes. In addition, the dashed line, which builds an ellipse with its long axis inclined by an angle of 45 • with respect to the x direction, is the isotropic yield surface. Its yield stress values of X uniaxial tension, Y uniaxial tension, X uniaxial compression, and Y uniaxial compression were completely equal. At an accumulative strain of 0.002, the calculated yield surface approximated the ellipse. However, it gradually expanded or shrank, rotated, and distorted with increasing accumulative plastic strain, indicating the evolution effect with complicated strain hardening and softening along different routes.

Uniaxial Loading
In this study, the large deformation uniaxial loading of the extruded ZK60 magnesium alloy was first simulated. The experimental results (yield points as shown in Figure 5), the FE stress-strain relationship, and the hardening rate calculated using nonevolution, linear interpolation, and the analytical interpolation approach [36] are shown in Figure 6. dition, the dashed line, which builds an ellipse with its long axis inclined by an angle of 45° with respect to the x direction, is the isotropic yield surface. Its yield stress values of X uniaxial tension, Y uniaxial tension, X uniaxial compression, and Y uniaxial compression were completely equal. At an accumulative strain of 0.002, the calculated yield surface approximated the ellipse. However, it gradually expanded or shrank, rotated, and distorted with increasing accumulative plastic strain, indicating the evolution effect with complicated strain hardening and softening along different routes. As can be seen from the stress-strain relationships and strain hardening rates of continuous large deformation X uniaxial tension in Figure 6a, across the entire strain range of 0-0.1, the FEM simulation results with and without considering evolution were exactly the same and in good agreement with the experimental result because of the absence of yield surface evolution along this loading route [36].
For the X uniaxial compression (Figure 6b), the simulated stress-strain relationships considering evolution conformed to the experimental result significantly better than the one without considering evolution across the entire accumulative strain range of 0-0.1. At an accumulative strain of 0.005, the stress was corrected from 47.3 to 32.9 MPa, which was approximately equal to the experiment's 32.5 Mpa. At an accumulative strain of 0.03, the stress was effectively corrected from 94.0 to 109.8 Mpa, which was approximately equal to the experiment's 109.2 Mpa. At an accumulative strain of 0.06, the stress was corrected from 94.9 to 119.6 Mpa, which was approximately equal to the experiment's 119.0 Mpa. At an accumulative strain of 0.08, the stress was effectively corrected from 91.9 to 124.7 Mpa, which was approximately equal to the experiment's 124.0 Mpa. It is worth noting that in the strain range of 0.005-0.02, the traditional linear interpolation approach [31] showed an increasing hardening rate, while the analytical interpolation approach [36] showed a decreasing hardening rate. The experimental hardening exponent was in the range of 0-1, and therefore the decreasing hardening rate simulated by the analytical interpolation approach was more reasonable. In addition, in the same strain range, the analytical interpolation approach still captured the upward convex stress-strain relationship well, even when the numerical point gap was big; however, the traditional linear interpolation approach had a comparatively large error. It can be concluded that the traditional linear interpolation approach depends on whether the hardening is linear and how big the interpolation gap is and that the analytical interpolation approach requires fewer data and is more accurate and reasonable.

Uniaxial Loading
In this study, the large deformation uniaxial loading of the extruded ZK60 magnesium alloy was first simulated. The experimental results (yield points as shown in Figure  5), the FE stress-strain relationship, and the hardening rate calculated using nonevolution, linear interpolation, and the analytical interpolation approach [36] are shown in Figure 6. As can be seen from the stress-strain relationships and strain hardening rates of continuous large deformation X uniaxial tension in Figure 6a, across the entire strain range of 0-0.1, the FEM simulation results with and without considering evolution were exactly the same and in good agreement with the experimental result because of the absence of yield surface evolution along this loading route [36].
For the X uniaxial compression (Figure 6b), the simulated stress-strain relationships considering evolution conformed to the experimental result significantly better than the For Y uniaxial compression (Figure 6d), the simulated stress-strain relationship considering evolution was in significantly better accordance with the experimental result than the one without corrections across the entire accumulative strain range of 0-0.1. At an accumulative strain of 0.05, the stress was corrected from 85.0 to 105.5 Mpa, which was approximately equal to the experiment's 105.2 Mpa. At an accumulative strain of 0.06, the stress was corrected from 86.9 to 108.9 Mpa, which was approximately equal to the experiment's 108.5 Mpa. At an accumulative strain of 0.08, the stress was effectively corrected from 84.0 to 112.9 Mpa, which was approximately equal to the experiment's 112.5 Mpa. At an accumulative strain of 0.1, the stress was effectively corrected from 80.1 to 110.4 Mpa, which was approximately equal to the experiment's 110.0 Mpa. The correction term clearly corrected the strain hardening rate in the accumulative strain range of 0.02-0.1. In addition, in the accumulative strain ranges of 0-0.03 and 0.05-0.1, the simulation results of the traditional linear interpolation and analytical interpolation were substantially identical to each other. It is worth noting that in the accumulative strain range of 0.03-0.05, the traditional linear interpolation approach showed an increasing hardening rate and a clear error in the stress-strain relationship, while the analytical interpolation approach showed a decreasing hardening rate and captured the upward convex stress-strain relationship much better. This indicated that the analytical interpolation approach was more reasonable and accurate.
Clearly, for continuous large deformation uniaxial loading of extruded ZK60 magnesium alloy, the simulation results without evolution showed clear errors, and the simulation results considering evolution could better reflect the stress-strain relationships and hardening rates. In addition, for X uniaxial compression and Y uniaxial compression, the analytical interpolation approach was more reliable and efficient than the traditional linear interpolation approach.

Biaxial Loading
To further validate the evolution effect described by the proposed asymmetric yield function and constitutive relationship for continuous large deformations and to determine the accuracy of the evaluation with corrections, simulations of the displacement loading routes of equibiaxial tension, equibiaxial compression, X tension-Y compression, and X compression-Y tension of extruded ZK60 magnesium alloy were performed. The simulation data of the uniaxial and biaxial loading and experimental yield surfaces on the 2D principal stress plane are shown in Figure 7a the analytical interpolation approach was more reliable and efficient than the traditional linear interpolation approach.

Biaxial Loading
To further validate the evolution effect described by the proposed asymmetric yield function and constitutive relationship for continuous large deformations and to determine the accuracy of the evaluation with corrections, simulations of the displacement loading routes of equibiaxial tension, equibiaxial compression, X tension-Y compression, and X compression-Y tension of extruded ZK60 magnesium alloy were performed. The simulation data of the uniaxial and biaxial loading and experimental yield surfaces on the 2D principal stress plane are shown in Figure 7a,b. As shown in Figure 7, the simulated yield points with and without considering evolution were exactly the same along the nonevolving X uniaxial tension route and conformed to the experimental yield surfaces very well. However, in the other three uniaxial loading routes, the simulated yield results without considering evolution had clear errors, and those considering evolution effect better conformed to the yield surfaces.
For X-Y biaxial tension and X-Y biaxial compression, at a low accumulative strain of 0.002, the simulated yield points considering evolution coincided with those without considering evolution because the evolution effect was insignificant at the small deformation stage. However, with increasing strain, the evolution effect gradually became As shown in Figure 7, the simulated yield points with and without considering evolution were exactly the same along the nonevolving X uniaxial tension route and conformed to the experimental yield surfaces very well. However, in the other three uniaxial loading routes, the simulated yield results without considering evolution had clear errors, and those considering evolution effect better conformed to the yield surfaces.
For X-Y biaxial tension and X-Y biaxial compression, at a low accumulative strain of 0.002, the simulated yield points considering evolution coincided with those without considering evolution because the evolution effect was insignificant at the small deformation stage. However, with increasing strain, the evolution effect gradually became clear. Substituting the simulated stress state into the yield surface equation [36] corresponding to the accumulative plastic strain, the intersection of the proportional yield surface and the X-axis, which is hereafter defined as simulated results, was obtained.
For X-Y biaxial tension, in the accumulative strain range of 0.005-0.1, the stress was effectively corrected. At an accumulative strain of 0.03, the simulated stress value was greatly corrected from 129.8 to 116.1 Mpa, which was approximately equal to the experiment's 116.0 Mpa.
For X-Y biaxial compression, in the accumulative strain range of 0.005-0.1, the stress was effectively corrected. At an accumulative strain of 0.03, the simulated stress was greatly corrected from 130.0 to 116.6 Mpa, which was approximately equal to the experiment's 116.0 Mpa.
Along the routes of X tension-Y compression and X compression-Y tension, at a low accumulative strain of 0.002, the simulated yield points considering evolution coincided with those without considering evolution because the evolution effect was insignificant at the small deformation stage. However, with increasing strain (greater than or equal to 0.005), the evolution effect gradually became clear, and the effect of corrections could be observed. For X tension-Y compression, in the accumulative strain range of 0.005-0.1, the stress was effectively corrected. At an accumulative strain of 0.05, the stress was greatly corrected from 110.7 to 115.3 Mpa, which was approximately equal to the experiment's 115.0 Mpa, and the correction effect was the most obvious. For X compression-Y tension, in the accumulative strain range of 0.005-0.1, the stress was effectively corrected. At an accumulative strain of 0.1, the stress was greatly corrected from 84.1 to 110.5 Mpa, which was approximately equal to the experiment's 110.0 Mpa. Because the true strain of tension was clearly different from that of compression at large deformation stages, the biaxial stresses gradually exhibited significant asymmetry.
By considering the partial derivative of the continuously interpolated yield function ∂f * * ∂( dε p ) [36] with respect to the accumulative plastic strain, the FEM simulations accurately captured the evolving asymmetric yield of the extruded ZK60 magnesium alloy during continuous large deformations. To further validate that the large deformation constitutive relationship and the analytical interpolation approach had a wide range of applicability, continuous large deformation loading-reverse loading and orthogonal loading FE simulations of the extruded ZK60 magnesium alloy were performed in this study.

Loading-Reverse Loading
The FEM simulation results of loading-reverse loading are shown in Figure 8a,b. The displacement loading route was first along X uniaxial tension (Figure 8a) or Y uniaxial tension (Figure 8b) from the accumulative plastic strain of 0 to 0.03 (points: 0 → 1 → 2 → 3 → 4 ). The unloading and reverse loading were then conducted, and the accumulative plastic strain increment kept to zero until point 5 (this process is an elastic deformation, and the yield stress at point 5 was different from that at point 4 due to the tension-compression SDE). The subsequent yield was from the accumulative plastic strain of 0.03 to 0.1 (points: 5 → 6 → 7 → 8 → 9 → 10 ) with increasing reverse loading. For loading-reverse loading, the simulations considering evolution were in significantly better agreement with the experimental yield surfaces than that without considering evolution. In addition, there was no clear difference between the traditional linear interpolation and the analytical interpolation in the current simulations.

Orthogonal Loading
The FEM simulation results of orthogonal loading are shown in Figure 9a,b. The displacement loading route was first along X uniaxial tension (Figure 9a) or Y uniaxial tension (Figure 9b) from the accumulative plastic strain of 0 to 0.03 (points: 0 1 2 3 4 → → → → ), and then the unloading was conducted until the stress decreased to zero (points: 4 0 → , which is an elastic deformation). Subsequently, the Y uniaxial tension ( Figure 9a) and X uniaxial tension (Figure 9b) are conducted, respectively, and the accumulative plastic strain increment keeps to be zero until point 5 (this process is also an elastic deformation, and the yield stress at point 5 was different from that at point 4 due to the anisotropy). The subsequent yield was from the accumulative plastic strain of 0.03 to 0.1 (points: 5 6 7 8 9 10 → → → → → ) with increasing Y uniaxial tension. For orthogonal loading, the simulations considering evolution were in significantly better agreement with the experimental yield surfaces than that without considering evolution. In addition, there was no clear difference between the traditional linear interpolation and the analytical interpolation in the current simulations.

Orthogonal Loading
The FEM simulation results of orthogonal loading are shown in Figure 9a,b. The displacement loading route was first along X uniaxial tension (Figure 9a) or Y uniaxial tension (Figure 9b) from the accumulative plastic strain of 0 to 0.03 (points: 0 → 1 → 2 → 3 → 4 ), and then the unloading was conducted until the stress decreased to zero (points: 4 → 0 , which is an elastic deformation). Subsequently, the Y uniaxial tension (Figure 9a) and X uniaxial tension (Figure 9b) are conducted, respectively, and the accumulative plastic strain increment keeps to be zero until point 5 (this process is also an elastic deformation, and the yield stress at point 5 was different from that at point 4 due to the anisotropy). The subsequent yield was from the accumulative plastic strain of 0.03 to 0.1 (points: 5 → 6 → 7 → 8 → 9 → 10 ) with increasing Y uniaxial tension. For orthogonal loading, the simulations considering evolution were in significantly better agreement with the experimental yield surfaces than that without considering evolution. In addition, there was no clear difference between the traditional linear interpolation and the analytical interpolation in the current simulations.
Although the results of the traditional linear interpolation at some loading routes were similar to the results of the analytical interpolation approach, the linear assumption is still rough, and the accuracy is strongly dependent on whether the hardening is linear and how big the interpolation gap is. Therefore, the accuracy of the results between the traditional linear interpolation starting and terminal points cannot be guaranteed, and the linear assumption also cannot well describe the strong evolution (see the X uniaxial compression and Y uniaxial compression of the ZK60 magnesium alloy in Figure 6b,d). As opposed to this assumption, the analytical interpolation approach considering the physical hardening effect and the Lode parameters is independent of numerous experimental data and model parameters and presents higher accuracy and reasonability. Although the results of the traditional linear interpolation at some loading routes were similar to the results of the analytical interpolation approach, the linear assumption is still rough, and the accuracy is strongly dependent on whether the hardening is linear and how big the interpolation gap is. Therefore, the accuracy of the results between the traditional linear interpolation starting and terminal points cannot be guaranteed, and the linear assumption also cannot well describe the strong evolution (see the X uniaxial compression and Y uniaxial compression of the ZK60 magnesium alloy in Figure 6b,d). As opposed to this assumption, the analytical interpolation approach considering the physical hardening effect and the Lode parameters is independent of numerous experimental data and model parameters and presents higher accuracy and reasonability.

Conclusions
Under the deformation conditions of low temperature (200 °C) and high processing strain rate ( wire-cutting samples along three directions of 0, 45, and 90° in the transverse extrusion plane of extruded ZK60 magnesium alloy. The yield surface anisotropy parameters were obtained using the true stress-accumulative plastic strain data and the optimization algorithm, and the yield surfaces at different accumulative plastic strains were determined. In addition, the hypoelastoplastic large deformation constitutive relationship based on a stress invariant asymmetric yield function considering evolution, an analytical interpolation approach related to the accumulative plastic strain and stress state, and a conjugate pair of rotational Kirchhoff stress-rotational deformation rate was used to simulate evolving large deformation yield subjected to uniaxial, biaxial, loading-reverse loading, and orthogonal loading in the commercial FE software Abaqus with user subroutine. Some conclusions can be drawn as follows: (1) As a typical plain Mg alloy, extruded ZK60 alloy exhibited significant anisotropy, tension-compression SDE, and evolution effect during low-temperature continuous large deformation yield. At an accumulative strain of 0.002, the shape of the yield surface in 2D stress space approximated an ellipse. As the accumulative plastic strain increased, the yield surface not only expanded or shrank in proportion to the ratio of the X uniaxial tensile stresses but also simultaneously rotated and distorted,

Conclusions
Under the deformation conditions of low temperature (200 • C) and high processing strain rate (10 −2 s −1 ), uniaxial tension and compression tests were conducted on wirecutting samples along three directions of 0, 45, and 90 • in the transverse extrusion plane of extruded ZK60 magnesium alloy. The yield surface anisotropy parameters were obtained using the true stress-accumulative plastic strain data and the optimization algorithm, and the yield surfaces at different accumulative plastic strains were determined. In addition, the hypoelastoplastic large deformation constitutive relationship based on a stress invariant asymmetric yield function considering evolution, an analytical interpolation approach related to the accumulative plastic strain and stress state, and a conjugate pair of rotational Kirchhoff stress-rotational deformation rate was used to simulate evolving large deformation yield subjected to uniaxial, biaxial, loading-reverse loading, and orthogonal loading in the commercial FE software Abaqus with user subroutine. Some conclusions can be drawn as follows: (1) As a typical plain Mg alloy, extruded ZK60 alloy exhibited significant anisotropy, tension-compression SDE, and evolution effect during low-temperature continuous large deformation yield. At an accumulative strain of 0.002, the shape of the yield surface in 2D stress space approximated an ellipse. As the accumulative plastic strain increased, the yield surface not only expanded or shrank in proportion to the ratio of the X uniaxial tensile stresses but also simultaneously rotated and distorted, evolving into convex and smooth closed curves with various shapes. At an accumulative strain of 0.02, the yield difference between X uniaxial compression and Y uniaxial compression was 24.0 MPa and the relative difference was as high as 29.81%. At an accumulative strain of 0.03, the yield difference between X uniaxial tension and Y uniaxial tension was 29.5 MPa and the relative difference reached 34.1%. At an accumulative strain of 0.005, the X tension-compression yield difference was 15.5 MPa and the relative difference reached 47.69%. At a strain of 0.08, the Y tensioncompression yield difference was 15.5 MPa and the relative difference reached 15.98%. (2) For the stress-strain relationships and hardening rates of all loading routes, the simulation results considering the yield surface evolution effect were approximately the same as the experiments, and the analytical interpolation approach was more reasonable and accurate than the linear interpolation approach. When the accumulative plastic strain was large, the simulation results without considering the evolution effect, excluding X uniaxial tension, clearly deviated from the experiments (X uniaxial tension loading route had no evolution effect).
Low temperature and high strain rate can effectively save process cost and improve processing efficiency. However, asymmetric yield evolution behavior during the processes may significantly influence the forming prediction and mechanical performance of magnesium alloys. The simulation approach presented herein, which can accurately simulate large deformation yield of magnesium alloys under complex loading routes, is promising for the prediction and optimization of severe plastic processes of Mg alloys.

Data Availability Statement:
The raw/processed data will be provided as request.

Acknowledgments:
The authors are grateful for the support from the Natural Science Foundation of Zhejiang Province (No. LY18A020003), the Natural Science Foundation of Ningbo City (No. 202003N4083), and the K. C. Wong Magna Fund from Ningbo University.

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