Investigation of Yield Surfaces Evolution for Polycrystalline Aluminum after Pre-Cyclic Loading by Experiment and Crystal Plasticity Simulation

This study aims at introducing the back stress of anisotropic strain-hardening into the crystal plasticity theory and demonstrating the rationality of this crystal plasticity model to describe the evolution of the subsequent yield surface of polycrystalline aluminum at the mesoscopic scale under complex pre-cyclic loading paths. By using two different scale finite element models, namely a global finite element model (GFEM) as the same size of the thin-walled tube specimen used in the experiments and a 3D cubic polycrystalline aggregate representative volume element (RVE) model, the evolution of the subsequent yield surface for different unloading cases after 30 pre-cycles is further performed by experiments and numerical simulations within a crystal plasticity finite element (CPFE) frame. Results show that the size and shape of the subsequent yield surfaces are extremely sensitive to the chosen offset strain and the pre-cyclic loading direction, which present pronounced anisotropic hardening through a translation and a distortion of the yield surface characterized by the obvious “sharp corner” in the pre-deformation direction and “flat” in the reverse direction by the definition of small offset strain, while the subsequent yield surface exhibits isotropic hardening reflected by the von Mises circle to be distorted into an ellipse by the definition of large offset strain. In addition, the heterogeneous properties of equivalent plastic strain increment are further discussed under different offset strain conditions. Modeling results from this study show that the heterogeneity of plastic deformation decreases as a law of fraction exponential function with the increasing offset strain. The above analysis indicates that anisotropic hardening of the yield surface is correlated with heterogeneous deformation caused by crystal microstructure and crystal slip. The crystal plasticity model based on the above microscopic mechanism can accurately capture the directional hardening features of the yield surface.


Introduction
Polycrystalline materials often exhibit more complicated plastic anisotropy during the forming and manufacturing processes of metal components [1,2]. The initial anisotropy induced by plastic deformations gives rise to the complex distorted evolution and translation of subsequent yield surfaces under complex pre-loading paths, which include cyclic loading histories. Under the monotonous proportional and non-proportional pre-loading paths, a series of studies has already demonstrated that the subsequent yield surface will change size, location, and shape [3][4][5][6][7][8], which are derived from Materials 2020, 13 the isotropic, kinematic, and distortional hardening mechanisms in plastic constitutive models [9]. Under the pre-cyclic loading paths, a few studies of the yield surfaces were performed [10,11] and the results showed that cyclic hardening, softening, ratcheting, Bauschinger effect, and cyclic strain amplitude resulted in a complex evolution of subsequent yield surfaces. For these complex subsequent yielding phenomena of materials, the study of yield surface evolution is the key to establish a reasonable constitutive model to describe metal plastic behavior accurately under complex pre-loading paths [10,12]. In general, the shape of the yield surface for most metals is sensitive to the yield definition. Numerous experimental results have shown the distorted evolution of the subsequent yield surfaces characterized by a triangle vertex in the loading direction but perfectly flattening in the reverse loading direction by yield definition of small offset strain [3,[5][6][7]13], while the subsequent yield surface based on large offset strain value is closer to von-Mises surface [14,15]. For instance, using the same experimental materials as stainless steel 304, the shape of subsequent yield surfaces becomes an ellipse with the definition of yield as 50 µε offset strain [16,17], but exhibits directional distortion with the definition of yield as 5 µε offset strain [3]. The previous experimental studies about the yield surface were performed by the smaller proof strain in the range of 5-50 µε using a single specimen method, which all yield points on a subsequent yield surface were determined by a thin-walled tube specimen. Because the accumulated plastic strain caused by the previous yield points probed and the accuracy of measuring instrument is negligible, some yield points on the subsequent yield surface by a single specimen methods may be have large deviation, but numerical simulations and multi-specimen test method can almost eliminate the effect of these factors using this inaccurate measurement methods. The experimental studies of the yield surfaces using multiple sample methods are rarely reported with various offset strain ranging from 20 µε to 1000 µε.
Based on experimental observations, an early documented attempt to describe the macroscopic anisotropy yielding of metals introduced the appropriate kinematic and distortional hardening rules, based on the classical isotropic hardening, into plasticity theories [18][19][20][21][22][23]. For instance, a non-quadratic yield surface had been adopted for modeling the anisotropic responses of aluminium alloys [24,25]. Some researchers introduced a distortional stress tensor into the von Mises yield criterion, and then the yield surface described as an "egg shape" [21]. In this study, the hardening rules were derived from thermodynamics for both the back-stress and the fourth order tensor. These macroscopic continuum models make no attempt to incorporate the microscopic behavior of the material, and the yield functions are defined by the stress field and a set of internal variables by the choice of distortion parameters. Most of the above models predict the von Mises circle to be distorted into an ellipse or a "round nose" rather than to become high curvature of yield surface distortion as experimental observations. The evolution of the yield surface depends markedly on the plastic slip mechanisms and the crystallographic texture for polycrystal material [25]. The crystal plasticity model can be employed to accurately describe the anisotropic evolution of the subsequent yield surfaces [26,27].
In this work, the evolution of the subsequent yield surfaces for different unloading cases with various yield definitions after different pre-cyclic loading is studied in detail by both experiments and crystal plasticity simulations using multiple identical specimens. The work is highlighted within a crystal plasticity finite element (CPFE) scheme using two different scale models: a thin-walled circular tube finite element model of the same size as the test specimen and a representative 3D cubic polycrystalline aggregate representative volume element RVE model. Furthermore, by comparing the experimental data with the numerical simulations, the ability of the crystal plasticity model considering a back stress for anisotropic strain-hardening to describe the anisotropic evolution of the subsequent yield surfaces under the loading path change and different yield definition is validated.

Experimental Material and Specimen
The material used in this study is the high-purity (99.89%) polycrystalline aluminum, and the metallographic photo shows an equiaxed grain structure and average grain size is proximately 80 µm in Figure 1. Its chemical composition and the basic mechanical properties from uniaxial tension test and pure torsion test are given in Table 1. These experiments are carried out by thin-walled tube specimens made of aluminum alloy materials subjected to combined axial-torsional loading, the geometry of the specimens is shown in Figure 2. The steel inserts are embedded in both ends of the specimens to ensure reliable clamping and resist the severe deformation.

Experimental Material and Specimen
The material used in this study is the high-purity (99.89%) polycrystalline aluminum, and the metallographic photo shows an equiaxed grain structure and average grain size is proximately 80 μm in Figure 1. Its chemical composition and the basic mechanical properties from uniaxial tension test and pure torsion test are given in Table 1. These experiments are carried out by thin-walled tube specimens made of aluminum alloy materials subjected to combined axial-torsional loading, the geometry of the specimens is shown in Figure 2. The steel inserts are embedded in both ends of the specimens to ensure reliable clamping and resist the severe deformation.

Experimental Procedure
The experiments are conducted by an MTS809 tension-torsional servo-controlled hydraulic testing machine (MTS Systems Corporation, Eden Prairie, MN, USA) at room temperature. An extensometer with a 25-mm gauge length is calibrated to measure the axial and torsion strain of the gauge section of specimens. The stress-strain curves of σ(ε) in uniaxial tension and τ(γ) in pure torsion are displayed in Figure 3. The stress-strain behaviour under uniaxial tensile loading exhibits the pronounced strain hardening effect in which a high strain-hardening rate is observable with the increase of axial strain. Unlike the axial component of flow stress, the stress-strain curve obtained from the pure torsion experiment shows the lower strain hardening behavior.

Experimental Material and Specimen
The material used in this study is the high-purity (99.89%) polycrystalline aluminum, and the metallographic photo shows an equiaxed grain structure and average grain size is proximately 80 μm in Figure 1. Its chemical composition and the basic mechanical properties from uniaxial tension test and pure torsion test are given in Table 1. These experiments are carried out by thin-walled tube specimens made of aluminum alloy materials subjected to combined axial-torsional loading, the geometry of the specimens is shown in Figure 2. The steel inserts are embedded in both ends of the specimens to ensure reliable clamping and resist the severe deformation.

Experimental Procedure
The experiments are conducted by an MTS809 tension-torsional servo-controlled hydraulic testing machine (MTS Systems Corporation, Eden Prairie, MN, USA) at room temperature. An extensometer with a 25-mm gauge length is calibrated to measure the axial and torsion strain of the gauge section of specimens. The stress-strain curves of σ(ε) in uniaxial tension and τ(γ) in pure torsion are displayed in Figure 3. The stress-strain behaviour under uniaxial tensile loading exhibits the pronounced strain hardening effect in which a high strain-hardening rate is observable with the increase of axial strain. Unlike the axial component of flow stress, the stress-strain curve obtained

Experimental Procedure
The experiments are conducted by an MTS809 tension-torsional servo-controlled hydraulic testing machine (MTS Systems Corporation, Eden Prairie, MN, USA) at room temperature. An extensometer with a 25-mm gauge length is calibrated to measure the axial and torsion strain of the gauge section of specimens. The stress-strain curves of σ(ε) in uniaxial tension and τ(γ) in pure torsion are displayed in Figure 3. The stress-strain behaviour under uniaxial tensile loading exhibits the pronounced strain hardening effect in which a high strain-hardening rate is observable with the increase of axial strain. Unlike the axial component of flow stress, the stress-strain curve obtained from the pure torsion experiment shows the lower strain hardening behavior.  The subsequent yield behaviors of polycrystal aluminum with the multiple-specimen method are studied at five different unloading cases (A→OA1, A→OA2, B→OB, C→OC, D→OD) under pre-cyclic loading paths by experiment and the crystal plasticity finite element scheme, and the unloading points OA1, OA2, OB, OC, or OD were defined as the center of the yield surface, respectively. The determined processes of the subsequent yield surfaces are mainly summarized as the following three steps. First (pre-cyclic loading process), each specimen is cyclically loaded to a steady stress state after 30 symmetrical tensile-compressive cycles at the 0.3% cyclic strain amplitude with the strain-controlled mode in Figure 4a. Then (unloading process), for each unloading case, the specimen is unloaded from different locations (A, B, C, or D) on a saturation hysteretic loop to the unloading end points (OA1, OA2, OB, OC or OD) under the stress-controlled mode respectively, as shown in Figure 4a (A→OA1, A→OA2, B→OB, C→OC, D→OD). It is worth noting specially that all unloading locations should be limited within the elastic range, or the lower yield strength caused by reverse loading will change the size of yield surface. Finally (reloading process), the specimen is reloaded along a specific reloading direction, which is represented by a radial ling in Figure 4b, to determine the yield points of subsequent yield surfaces with different target offset strain by combined tension-torsion proportional loading paths ranging from 0° to 180° at intervals of 15°. . It is worth noting specially that all unloading locations should be limited within the elastic range, or the lower yield strength caused by reverse loading will change the size of yield surface. Finally (reloading process), the specimen is reloaded along a specific reloading direction, which is represented by a radial ling in Figure 4b, to determine the yield points of subsequent yield surfaces with different target offset strain by combined tension-torsion proportional loading paths ranging from 0 • to 180 • at intervals of 15 • .
The definition of yield points adopted here is the unloading stiffness method, considering the unloading stiffness decreases with the increase of offset strain [5][6][7], The unloading stiffness E' corresponding to each offset strain is determined by the method of gradual reloading and unloading, as shown in Figure 4c. First, the specimen is reloaded from the unloading point to the stress state approach to the target offset strain ∆E p offset , then the specimen is unloaded back to the unloading point, and the unloading stiffness E is determined by fitting the linear portion of the unloading curve. By drawing a straight line Y = E (X − ∆E p offset ), which parallels the unloading stiffness E , the intersection of the straight line and the equivalent stress-strain curve is defined as the subsequent yield point represented by Y. By repeating this reloading-unloading process, the subsequent yield points corresponding to each given target strain 20 µε, 50 µε, 100 µε, 200 µε, 600 µε, and 1000 µε in each reloading direction can be determined completely.  The definition of yield points adopted here is the unloading stiffness method, considering the unloading stiffness decreases with the increase of offset strain [5][6][7], The unloading stiffness E' corresponding to each offset strain is determined by the method of gradual reloading and unloading, as shown in Figure 4c. First, the specimen is reloaded from the unloading point to the stress state approach to the target offset strain

Constitutive Framework
A crystalline plastic model is used which considers the typical 12 slip systems for a single crystal of face centered cubic (FCC) structure materials. The total deformation gradient within the crystal plasticity framework can be multiplicatively decomposed as: where F * is the elastic component arising from elastic stretching U e and rigid body rotation R L of the configuration according to the standard polar decomposition, which can be expressed as: Materials 2020, 13, 3069 6 of 20 F * = R L ·U e . F p is the plastic component describing the dislocation glide along the specific slip planes [28]. For a material point in a crystal grain, the Euler's velocity gradient tensor L can then be defined as the sum of the skew-symmetric gradient tensor L * and the plastic gradient tensor L p [29][30][31]: The Euler's velocity gradient tensor L can be further divided into a symmetric deformation velocity tensor, D = sym(L) = 1 2 (L + L T ), and the anti-symmetric rotation velocity tensor, transpose' matrix of the velocity gradient tensor. Further, we write: Subsequently, the deformation tensor and rotation tensor, respectively, can be decomposed into a lattice contribution elastic part and a plastic part [32]: where D * and D p are, respectively, the elastic part and the plastic part of stretching tensor. W p is the plastic spin tensor generated by crystal slip and W * is the lattice spin tensor, which can be defined as follows [33]: where R is the rigid body rotation of crystal lattice, . R T is elastic distortion. The plastic velocity gradient in the intermediate (relaxed) configuration, L p , is defined as the sum of the reference shear strain rate . γ (α) on all slip systems [34]: where the term (s (α) = m (α) n (α) ) represents the Schmidt tensor on the α-th slip system that is defined by the dyadic product of orthogonal unit vectors m (α) and n (α) in the reference configuration, with m (α) indicating the slip direction and n (α) being the normal to the slip plane.
. γ α is a reference shear strain rate on the α-th slip system, n is the total number of slip systems.
When the deformation gradient F * is applied, the crystalline lattice rotates and these vectors (m (α) and n (α) ) in the reference configuration transform to m (α) * and n (α) * respectively, in the current deformed configuration, such that: where it is assumed that m (α) and n (α) are orthogonal to each other.
The Schmid tensor on the α-th slip system can be written as the sum of asymmetric tensor p (α) and an anti-symmetric tensor ω (α) , such that: where the symmetric part and anti-symmetric part of the Schmid tensor on the α-th slip system can be defined as: The plastic deformation gradient can be expressed as: Similarly, the spin tensor, which represents the material rotation due to slip, can be expressed as: Further, the elastic stretching tensor and the lattice spin tensor can be written as:

Hardening Rules
In this work, considering the predominant cyclic strain preloading paths in the subsequent yield problem, while the effects of back-stress are essential to be considered for the case of cyclic loading. Accordingly, the conventional crystalline plastic theory has been developed to capture the cyclic plastic response of the alloy crystals [12]. According to the viscous regularization proposed by Hutchinson [35], an additional back-stress x (α) is here introduced into the viscoplastic constitutive model. The slip rate for a given slip system is given by: where . γ 0 is a reference shear strain rate on the α-th slip system, k is the power law exponent representing the rate sensitivity, τ (α) is the resolved shear stress on the α-th slip system of the single crystal, g (α) is the critical shear stress on the activated slip system α to govern the isotropic hardening of the crystal, and the back-stress x (α) is introduced here to characterize the nonlinear directional hardening of the crystal on the α-th slip system.
The evolution of the critical shear stress g (α) for a given activated slip system α is described as: where h αβ (γ) is the slip-plane hardening modulus which governs the interaction between various slip systems, γ (β) is the slip rate on the system β, n is the number of crystal slip systems, and γ is the accumulated shear strains on the all active slip systems. Moreover, Hutchinson suggested [36]: where q is a constant, whose value can be chosen in the range from 0 to 1, and Chang and Asaro suggested [37]: where h 0 is the initial hardening rate, τ 0 is the critical resolved shear stress and τ s is the saturation shear stress, which can be preliminarily determined according to elastic range of the experimental hysteresis loop.

Back Stress Evolution
According to Walker [38] and the Chaboche model [39], the evolution of back-stress .
where a is the initial hardening modulus of the slip system, c and d are the nonlinear hardening parameters, respectively, e 1 and e 2 are used to describe the law of the cyclic hardening saturation. The above material parameters which describe the cyclic plastic behavior are determined by the trial-and-error method combined with cyclic tests. The above constitutive relation has been implemented into the user material subroutine UMAT for the ABAQUS (Dassault Systems, Paris, France)/Standard module [10,32].

Finite Element Models at Different Scale
To treat the deformation behavior of micro-inhomogeneity at the grain scale, the finite element model (FEM) of polycrystal aluminum within crystal plasticity framework is built using a global finite element model (GFEM) as the same size of the thin-walled tube specimen used in the experiments. Two different representative models of the microstructure, which the GFEM consists of the same 76,000 C3D8 elements and different crystalline grains in the center gauge region, are used, as shown in Figure 5. The first one (GFEM) is shown in Figure 5a, where each grain as polycrystalline aggregate is composed of many elements with random shape, size, and random crystal orientations [40]. It is further demonstrated that this size is accurate enough of the statistical representation of the actual microstructure for obvious different number of grains in the center gauge region by sequential serial sectioning. Another one (GFEM) is built, as shown in Figure 5b, in which each cubic element standards for a grain [27,40]. Thus, the FEM model is divided into a great deal of grains resulted in a poor description of the grain shape. For a GFEM model, the macroscopic axial strain and torsion strain in each increment are calculated by outputting the axial displacement U 1 and torsion displacement U R1 at the reference point A and B simulating an axial-torsional extensometer of 25 mm fixed on the specimen's gauge section in the test. The model is constrained and loaded by coupling the nodes in fixed clamping end and load clamping end with their respective reference point, and the macroscopic loading conditions are introduced in the model by applying a axial load F and a torsion load T in the reference point of load clamping end in order to achieve the axial and torsional deformation history, as shown in Figure 5a.
In this research, the finite element model (FEM) for polycrystal aluminum is also built using a 3D cubic polycrystal RVE (cube of side 1 mm) with polycrystalline aggregate as a material point to simulate the gauge section of a thin-walled tube by the crystal plasticity model, as shown in Figure 6. The depicted RVE contains 27,000 elements and is divided into 125 equiaxed grains. The grains in the above FEM are given attributes with random morphology, size, and random crystal orientations, and the RVE is treated as a macroscopic isotropic body. The boundary conditions of loading and displacement are exerted to the RVE in three faces along the direction of x k (k = 1, 2, 3), respectively, shown in Figure 6b. In addition to using the equation constraints of ABAQUS software to constrain all nodes on the surface of the representative element body to meet the periodic boundary conditions [39,40], four nodes A, B, C, and D of the RVE are selected in the representative element body as reference points and the corresponding macroscopic constraints are applied to the element body in Figure 6b. In order to reflect the combined axial-torsional loading conditions as the experiments for the RVE model, macro axial force F and macro lateral torque T are applied to the reference point A on the loading surface relative to the fixed surface.
polycrystalline aggregate is composed of many elements with random shape, size, and random crystal orientations [40]. It is further demonstrated that this size is accurate enough of the statistical representation of the actual microstructure for obvious different number of grains in the center gauge region by sequential serial sectioning. Another one (GFEM) is built, as shown in Figure 5b, in which each cubic element standards for a grain [27,40]. Thus, the FEM model is divided into a great deal of grains resulted in a poor description of the grain shape. For a GFEM model, the macroscopic axial strain and torsion strain in each increment are calculated by outputting the axial displacement U1 and torsion displacement UR1 at the reference point A and B simulating an axial-torsional extensometer of 25 mm fixed on the specimen's gauge section in the test. The model is constrained and loaded by coupling the nodes in fixed clamping end and load clamping end with their respective reference point, and the macroscopic loading conditions are introduced in the model by applying a axial load F and a torsion load T in the reference point of load clamping end in order to achieve the axial and torsional deformation history, as shown in Figure 5a. In this research, the finite element model (FEM) for polycrystal aluminum is also built using a 3D cubic polycrystal RVE (cube of side 1 mm) with polycrystalline aggregate as a material point to simulate the gauge section of a thin-walled tube by the crystal plasticity model, as shown in Figure  6. The depicted RVE contains 27,000 elements and is divided into 125 equiaxed grains. The grains in the above FEM are given attributes with random morphology, size, and random crystal orientations, and the RVE is treated as a macroscopic isotropic body. The boundary conditions of loading and displacement are exerted to the RVE in three faces along the direction of 1,2,3 , respectively, shown in Figure 6b. In addition to using the equation constraints of ABAQUS software to constrain all nodes on the surface of the representative element body to meet the periodic boundary conditions [39,40], four nodes A, B, C, and D of the RVE are selected in the representative element body as

Parameters Calibration of Crystalline Plasticity Models
Material parameters for crystal plasticity are calibrated by a GFEM containing 36,000 elements and 36,000 crystalline grains in the center gauge region and the RVE consisting of 27,000 elements and 125 equiaxed grains by crystalline plasticity, respectively. Figure 7a shows that the cyclic stress-strain curve tends to saturate after 30 cycles during the tension-compression loading at 0.3% strain amplitude by experiment, in other words, the polycrystalline aluminum exhibits stable cyclic

Parameters Calibration of Crystalline Plasticity Models
Material parameters for crystal plasticity are calibrated by a GFEM containing 36,000 elements and 36,000 crystalline grains in the center gauge region and the RVE consisting of 27,000 elements and 125 equiaxed grains by crystalline plasticity, respectively. Figure 7a shows that the cyclic stress-strain curve tends to saturate after 30 cycles during the tension-compression loading at 0.3% strain amplitude by experiment, in other words, the polycrystalline aluminum exhibits stable cyclic hardening behavior. Therefore, the constitutive parameters of crystalline plasticity are calibrated by fitting the predicted stable hysteresis curve with the experimental data in Figure 7a. The constitutive parameters are further adjusted from the monotonic uniaxial curve in Figure 7b. hardening rate that describes the hardening rate of materials after they enter plasticity, a is the initial slipping hardening modulus and c is a constant leading to soften, e 1 and e 2 are the parameters to describe cyclic hardening rate and saturation rate of the stress-strain hysteresis curve, respectively. Further, the applied material parameters are adjusted repeatedly until the simulated stress-stra

The Number of Grain effects on the Macroscopic Behavior
Before proceeding further, the number of grains in the center gauge region of the GFEM should be enough to verify the convergence of the predicted results due to random grain orientation, size, and random morphology. However, the computational cost increases heavily with the number of grains, so an optimal number of grains should be discussed for the accurate predictions. For this reason, the center gauge region of the GFEM containing 64, 125, 216, 512, and 800 grains as shown in Figure 8a are studied, respectively. Each grain shown as a single color is a cube containing many elements. Figure 8b displays relatively finer grain size, which is slightly bigger than measured in the center gauge region of the GFEM containing 36,000 grains, each element represents one grain with a given orientation characterized by the color code.
The computationally obtained stress-strain cyclic curves for the GFEM with a different number of grains by crystal plasticity under cyclic tension-compression loading are given in Figure 9. It also predicts a slightly higher stress less than 2 MPa with 64 grains than that other grains. It is obvious that computational results of the stress-strain curves almost converge for the center gauge region of the GFEM containing more than 125 grains. From this sensitivity study, it is found that the GFEM with 125 grains are employed to show a converged macroscopic stress-strain response. This insensitivity of numerically computed mechanical behavior to the number of grains is resulted from the randomly arranged orientation and size of grains. The systematically predicted macroscopic stress-strain curves using the models with each grain containing many elements are more closely The model contains 13 parameters describing the mechanical behavior of an aluminum single crystal, such as the three elastic constants C 11 , C 12 , and C 44 ; the ten plastic constants: isotropic hardening constants τ 0 , τ s , h 0 ; kinematic hardening constants a, c, and d; cyclic hardening constants e 1 and e 2 ; the reference strain rate . γ 0 ; the rate sensitivity parameter k. The three independent elastic constants of single crystal aluminum that characterize the anisotropic elasticity tensor C 11 , C 12 and C 44 are obtained from the aluminum elastic properties as a function of elastic modulus E = 57 MPa, shear modulus G = 25 MPa, which are identified by matching the experimental uniaxial tensile and torsion stress-strain curves in the linear region for this alloy, and Poisson's ratio µ, which is prescribed as ν = 0.3. The reference strain rate . γ 0 is prescribed as 0.001/s −1 , which is assumed to the rate-independent material. The rate sensitivity parameter k is set to 200, since the pure aluminum is almost strain rate insensitive. The parameter d is set to 0, since creep is ignored during cyclic loading. Using these values, the remaining plastic parameters of crystalline plasticity model, namely, τ 0 , τ s , h 0 , a, c, e 1 , and e 2 are calibrated by fitting the stress-strain responses by the experiments under cyclic tension-compression loading, as shown in Figure 7a. The initial critical resolved value τ 0 and its saturated value τ s of shear stress are obtained from the elastic range from the peak points unloading to reverse yielding points in the initial and stable stress-strain hysteresis loop, respectively. For cyclic hardening is less pronounced, the material constant τ 0 is close to its saturated value τ s . The parameter h 0 is the initial hardening rate that describes the hardening rate of materials after they enter plasticity, a is the initial slipping hardening modulus and c is a constant leading to soften, e 1 and e 2 are the parameters to describe cyclic hardening rate and saturation rate of the stress-strain hysteresis curve, respectively. Further, the applied material parameters are adjusted repeatedly until the simulated stress-strain stabilized loop and tensile curve are well-matched with the experimental results in Figure 7. It is obvious that the computationally obtained cyclic stress-strain loop is in good agreement with the experimental results, verifying that the parameters used in crystal plasticity model are enough accurate. The finial elastic and plastic property parameters used in the crystal plasticity model are summarized and listed in Table 2.

The Number of Grain effects on the Macroscopic Behavior
Before proceeding further, the number of grains in the center gauge region of the GFEM should be enough to verify the convergence of the predicted results due to random grain orientation, size, and random morphology. However, the computational cost increases heavily with the number of grains, so an optimal number of grains should be discussed for the accurate predictions. For this reason, the center gauge region of the GFEM containing 64, 125, 216, 512, and 800 grains as shown in Figure 8a are studied, respectively. Each grain shown as a single color is a cube containing many elements. Figure 8b displays relatively finer grain size, which is slightly bigger than measured in the center gauge region of the GFEM containing 36,000 grains, each element represents one grain with a given orientation characterized by the color code.
Materials 2020, 13, x 12 of 21 approximated to that adopting the model with each element representing one grain. Therefore, the thin-walled finite element model in the center gauge region which contains a total of 36,000 elements and 36,000 grains can be adopted in the following computations.  The computationally obtained stress-strain cyclic curves for the GFEM with a different number of grains by crystal plasticity under cyclic tension-compression loading are given in Figure 9. It also predicts a slightly higher stress less than 2 MPa with 64 grains than that other grains. It is obvious that computational results of the stress-strain curves almost converge for the center gauge region of the GFEM containing more than 125 grains. From this sensitivity study, it is found that the GFEM with 125 grains are employed to show a converged macroscopic stress-strain response. This insensitivity of numerically computed mechanical behavior to the number of grains is resulted from the randomly arranged orientation and size of grains. The systematically predicted macroscopic stress-strain curves using the models with each grain containing many elements are more closely approximated to that adopting the model with each element representing one grain. Therefore, the thin-walled finite element model in the center gauge region which contains a total of 36,000 elements and 36,000 grains can be adopted in the following computations.

Results and Discussions
A family of the initial yield surfaces is determined by radial proportional stress-controlled paths with different ratio of shear stress  3 to tensile stress  ranging from 0° to 360° with 15° intervals with the different offset strain by using the GLBM, as shown in Figure 10. Firstly, we observe that the yield surfaces are close to isotropic expansion in the (  ,  3 ) plane with the increasing offset strains whose shape is smooth convex ellipses with two axes of symmetry: the tensile stress axis  and the shear stress axis  3 . The shape of the initial yield surfaces is independent of the yield definition choice, but the size of the initial yield surfaces is dependent on the various offset strain values. By comparison, for instance, the tensile and reverse compressive yield stress (38.8, −38.8 MPa) are higher than the shear yield stress and reverse shear (33.87, −33.87

Results and Discussions
A family of the initial yield surfaces is determined by radial proportional stress-controlled paths with different ratio of shear stress √ 3τ to tensile stress σ ranging from 0 • to 360 • with 15 • intervals with the different offset strain by using the GLBM, as shown in Figure 10. Firstly, we observe that the yield surfaces are close to isotropic expansion in the (σ, √ 3τ) plane with the increasing offset strains whose shape is smooth convex ellipses with two axes of symmetry: the tensile stress axis σ and the shear stress axis √ 3τ. The shape of the initial yield surfaces is independent of the yield definition choice, but the size of the initial yield surfaces is dependent on the various offset strain values. By comparison, for instance, the tensile and reverse compressive yield stress (38.8, −38.8 MPa) are higher than the shear yield stress and reverse shear (33.87, −33.87 MPa) for a specified offset strain of 1000 µε. The initial yield surface of the polycrystal aluminum may be described by the Tresca yield criterion. With reference to Figure 4, the different reloading paths in the (σ , τ 3 ) plane among 0° to 180° are adopted to probe the experimental subsequent yield surfaces at the beginning of the unloading points OA1 (and OA2, the reverse yielding point), OB, OC, and OD within the elastic domain as the assumed yield surface center corresponding to the four different positions A, B, C, and D, With reference to Figure 4, the different reloading paths in the (σ, √ 3τ) plane among 0 • to 180 • are adopted to probe the experimental subsequent yield surfaces at the beginning of the unloading points O A1 (and O A2 , the reverse yielding point), O B , O C , and O D within the elastic domain as the assumed yield surface center corresponding to the four different positions A, B, C, and D, respectively, after 30 symmetry tension-compression cycles at 0.3% strain amplitude under strain control. Furthermore, the experimental results are compared with the predicted ones by the crystal plasticity simulation using the GFEM containing 36,000 elements and 36,000 crystalline grains in the gauge section of a specimen and the RVE containing 27,000 elements and 125 equiaxed grains, and the differences of subsequent yield surfaces predicted between the two finite element models (FEM) with different scales and boundary conditions are discussed.
The experimental subsequent yield surfaces under the above five different unloading cases are plotted in Figure 11, the distinct features of the yield surface for anisotropic yielding can be observed. The subsequent yield surfaces expand rapidly towards the reverse loading direction and shrinks in the direction perpendicular to the pre-loading direction but has little influence in pre-deformation direction. This suggests a rapid hardening of slip under the reverse loading condition and a slow hardening in the preloading direction that dominates deformation under tension (or compression). This is the cause of the higher curvature of yield surface in the pre-loading direction than that in the other directions, and as a result, the subsequent yield surfaces exhibit the "sharp corner" in the pre-deformation direction and flat in its reverse direction. The anisotropic hardening characteristics of the subsequent yield surface strongly depend on the direction of the accumulated plastic strain, which is verified that the "sharp corner" of yield surfaces points to the pre-deformation direction corresponding to tension or compression, seen Figure 11, and their expansion and translation in the stress space are clearly observed. The yield surfaces are approximately symmetric about the stress axis σ located in the pre-loading direction. In Figure 11a,b, the effects of the unloading points O A1 and O A2 on the evolution of the yield surfaces are displayed. By comparison of the different unloading points, the results of the numerical and experimental study reveal that the yield surface distortion in the pre-loading direction becomes more evident owing to the unloading range increased (A-O A2 ) and the size of yield surfaces becomes shrinkage due to the small reverse plastic deformation at O A2 as the smaller offset strain of 20 µε and 50 µε, the yield surface in stress space is also translated observably.
As mentioned above, the size and shape of the subsequent yield surface depend on the yield definition and, in fact, similar experimental observations were reported in the literature [41,42]. The subsequent yield surfaces will exhibit pronounced kinematic and distorted hardening with a relatively small offset strain. The distorted evolution is characterized by the yield locus with a rounded triangle vertex in the loading direction and perfectly flat opposite to the loading direction. If a relatively large target deviation strain is chosen, for instance, in Figure 11, 200 µε, 600 µε, and 1000 µε, the subsequent yield surfaces exhibit approximately isotropic hardening.
As shown in Figure 12, the computationally obtained and experimentally measured subsequent yield surfaces are compared, and the trend of subsequent surface predicted is consistent with the experimental observations. It can be seen that the sizes of the yield surfaces have some differences between them, but the distinct features of the yield surfaces for anisotropic yielding can be accurately predicted. It demonstrates that the crystal plasticity model considering the microstructure heterogeneities and slipping mechanism within the crystal can accurately describe the anisotropic evolution of the subsequent yield surfaces caused by pre-loading deformation. As shown in Figure 12, the computationally obtained and experimentally measured subsequent yield surfaces are compared, and the trend of subsequent surface predicted is consistent with the experimental observations. It can be seen that the sizes of the yield surfaces have some differences between them, but the distinct features of the yield surfaces for anisotropic yielding can be accurately predicted. It demonstrates that the crystal plasticity model considering the microstructure To verify the predictive capabilities of the approach under various levels of the finite element models (FEM) from a micro material point (RVE) to macro specimen (GFEM) scale transition, the first scheme adopted here is the GFEM to simulate a same test specimen with the constrained and loaded conditions corresponding to the axial and torsional deformation history consistent with experiments. In the second scheme, the polycrystalline aggregate RVE model containing 125 grains is adopted to simulated a material point of specimen gauge section, and to reflect the continuity of deformation during loading, the RVE model is applied with the periodic boundary conditions, the macroscopic deformation is controlled by applying displacement and load on the master node of the model, as shown in Figure 6b. The subsequent yield surfaces at five unloading cases with different offset strain are presented with the GFEM and RVE models in Figure 12. As a compromise between accuracy and computational cost, the response of macroscopic deformation is believed to be better approximated when the GFEM is adopted in this study, whereas the RVE model provides a computationally cheap estimate for the polycrystalline response although it expects to be relatively inaccurate. The GFEM systematically predicts the yield surfaces which are slightly larger than the ones by RVE, as it corresponds to the same offset strain due to the development of lattice preferred orientations. As a consequence, there is a nearly indistinguishable difference in the predicted subsequent yield surfaces between the GFEM and the RVE model. The coincidence is attributed to the low sensitivity of finite element model sizes. The statistical analysis about the inhomogeneity of the Representative Volume Element (RVE) with micro-structure randomly generated is extremely important for the physical nature of macroscopic mechanical respond for polycrystalline material undergoing different deformation. In this work, subjected to the deformations of 20 uε and 1000 uε at the reloading direction of 180°, namely, the lower limiting value and upper limit value of yield definition in this work. The contour distribution of the accumulated equivalent plasticity strain and Mises stress in the RVE by The statistical analysis about the inhomogeneity of the Representative Volume Element (RVE) with micro-structure randomly generated is extremely important for the physical nature of macroscopic mechanical respond for polycrystalline material undergoing different deformation. In this work, subjected to the deformations of 20 µε and 1000 µε at the reloading direction of 180 • , namely, the lower limiting value and upper limit value of yield definition in this work. The contour distribution of the accumulated equivalent plasticity strain and Mises stress in the RVE by simulation using the crystal plasticity model are displayed in Figure 13, where heterogeneous distribution of the two variables can be clearly observed. To compare the inhomogeneity of equivalent plastic strain increment with different offset strain, a coefficient of variation (COV) was defined as the ratio of the standard deviation S of the accumulated effective plastic strain increment to its average absolute value [10]. Figure 14 illustrates the relation curves between the COV of the equivalent plastic strain increment and offset strain ranging from 20 με to 1000 με along various reloading paths from 0° to 180° with interval 30° at unloading point D-OD after 30 tension-compression cycles. The results show that heterogeneous distributions of strain and stress can be clearly observed in Figure 14 with increasing deformation as accumulated equivalent plasticity strain increments, respectively. For example, the COV of the equivalent plastic strain increment corresponding to the target offset strain of 20 με and 1000 με decreases from 1.453 to 0.26 at the reloading direction of 180° and the statistical-numerical distribution of the equivalent plastic strain increment are shown in Figure 15. These results demonstrate that with a very small offset strain, the representative volume element (RVE) presents very serious heterogeneity leading to anisotropy distortion of the yield surface due to a small amount of elements and grains generating new plastic strain in the RVE. However, if the larger offset strain is chosen, numerous elements and grains occur with new plastic deformation in the To compare the inhomogeneity of equivalent plastic strain increment with different offset strain, a coefficient of variation (COV) was defined as the ratio of the standard deviation S of the accumulated effective plastic strain increment to its average absolute value [10]. Figure 14 illustrates the relation curves between the COV of the equivalent plastic strain increment and offset strain ranging from 20 µε to 1000 µε along various reloading paths from 0 • to 180 • with interval 30 • at unloading point D-O D after 30 tension-compression cycles. The results show that heterogeneous distributions of strain and stress can be clearly observed in Figure 14 with increasing deformation as accumulated equivalent plasticity strain increments, respectively. For example, the COV of the equivalent plastic strain increment corresponding to the target offset strain of 20 µε and 1000 µε decreases from 1.453 to 0.26 at the reloading direction of 180 • and the statistical-numerical distribution of the equivalent plastic strain increment are shown in Figure 15. These results demonstrate that with a very small offset strain, the representative volume element (RVE) presents very serious heterogeneity leading to anisotropy distortion of the yield surface due to a small amount of elements and grains generating new plastic strain in the RVE. However, if the larger offset strain is chosen, numerous elements and grains occur with new plastic deformation in the RVE, which leads to the plastic deformation of material tending to be homogeneous. By comparison with the same offset strain, as shown in Figure 14, the higher heterogeneity in the preloading direction is revealed than that towards the reverse loading direction. Therefore, no violation of the convexity can be regarded as a requirement for choosing the proper specified offset strain.  Figure 14. Variation curves of COV (Here, COV was defined as the ratio of the standard deviation S of the accumulated effective plastic strain increment to its average absolute value) for the equivalent plastic strain increment vs offset strain along different reloading paths at unloading point D-OD after 30 tension-compression cycles.

Conclusions
In this work, the crystal plasticity model considering the back-stress associated with the finite element models (FEM) approaches including a micro material point (RVE) and macro specimen (GFEM) at different length scales of polycrystalline aggregates is adopted to validate the predictive capabilities for the subsequent yield surfaces of polycrystal aluminum, and the simulated results are assessed by their comparison with experimental observations. Conclusions are summarized as follows:

Conclusions
In this work, the crystal plasticity model considering the back-stress associated with the finite element models (FEM) approaches including a micro material point (RVE) and macro specimen (GFEM) at different length scales of polycrystalline aggregates is adopted to validate the predictive capabilities for the subsequent yield surfaces of polycrystal aluminum, and the simulated results are assessed by their comparison with experimental observations. Conclusions are summarized as follows: 1. The shape of the initial yield surfaces is independent of the choice of the yield definition, but the size of the initial yield surfaces is strongly associated with the various offset strain. The initial

Conclusions
In this work, the crystal plasticity model considering the back-stress associated with the finite element models (FEM) approaches including a micro material point (RVE) and macro specimen (GFEM) at different length scales of polycrystalline aggregates is adopted to validate the predictive capabilities for the subsequent yield surfaces of polycrystal aluminum, and the simulated results are assessed by their comparison with experimental observations. Conclusions are summarized as follows: 1. The shape of the initial yield surfaces is independent of the choice of the yield definition, but the size of the initial yield surfaces is strongly associated with the various offset strain. The initial yield surface of the polycrystal aluminum can be described by the Tresca yield criterion.
2. The directional hardening characteristics of the subsequent yield surfaces for the polycrystal aluminum are strongly associated with the direction of the accumulated plastic strain, which expands rapidly towards the reverse loading direction and shrinks in the direction orthogonal to the preloading direction but has little influence in the pre-deformation direction. Therefore, the subsequent yield surfaces exhibit the "sharp corner" in the pre-deformation direction and flat in the reverse direction.
3. Unlike the initial yield surfaces, the size and shape of the subsequent yield surfaces are obviously sensitive to the offset strain. The subsequent yield surfaces present pronounced kinematic hardening (Bauschinger effect) and distortion with a relatively small offset strain, which leads to a triangle vertex in the direction of loading but flat in the reverse loading direction. If a relatively large offset strain is chosen, the subsequent yield surfaces seem to the isotropic hardening.
4. The yield surface calculated by using the GFEM is believed to be better approximated whereas the RVE model requires a lower computational cost. There is a nearly indistinguishable difference regarding the simulation of subsequent yield surfaces between the two different scale models. We attribute this coincidence to the low sensitivity of finite element model scales.