Inﬂuence of Crystal Plasticity Parameters on the Strain Hardening Behavior of Polycrystals

: The effective mechanical properties of a polycrystal depend directly on the single-crystal properties of each grain and its crystallographic orientation with respect to the load axis. While the micromechanical approach has been used quite extensively to study the inﬂuence of grain shape and crystallographic texture on the resulting mechanical behavior of a polycrystal, the inﬂuence of the crystal plasticity parameters, which describe the constitutive behavior of the single crystal, requires to be investigated systemically because, typically, these parameters are ﬁtted to describe a given material behavior. In the current research, this gap is ﬁlled by systemically studying the effect of changes in crystal plasticity parameters on the effective mechanical properties of polycrystals. The numerical model employed here consists of a representative volume element of 100 grains, and the material properties are described by using a non-local crystal plasticity model. A proper homogenization technique was used to homogenize the micromechanical results to an effective macroscopic material response. The equivalent stress versus equivalent plastic strain curve was obtained numerically by introducing the Voce-type hardening law, mimicking the material behavior in uniaxial tensile tests. The four parameters of the Voce-type hardening law were ﬁtted to the macroscopic stress-strain curves, and the correlation between the crystal plasticity parameters and the Voce parameters has been studied, which is an efﬁcient way to study the inﬂuence of microscopic material descriptions on the macroscopic behavior of polycrystals.


Introduction
Micromechanical analysis is the study of heterogeneous materials at the level of the individual components making up these materials. The main goal of micromechanical modeling is to build a bridge between large-scale modeling and micromechanical phenomena [1]. Heterogeneous materials have, due to the presence of various phases, different physical and mechanical properties, leading to anisotropic material behavior. Therefore, micromechanical analyses are used to predict the macroscopic anisotropic behavior of these heterogeneous materials based on each phase's material geometry and properties. This procedure is known as homogenization [2][3][4][5][6], which aims at deriving a mean-field description of the macroscopic material behavior from microstructural information [7]. The homogenization technique was applied to analyze macroscopic responses of textured polycrystals [3], composites [4,6], and the calibration of the yield surface of different steels [5].
From a materials science point of view, establishing a link between microstructure features and mechanical properties is crucial to obtain an accurate mechanical response of the material on the macroscopic level. One key principle of micromechanical modeling is generating a representative volume element (RVE) that can capture important microstructure features such as grain size and grain shape distribution [6]. Among the various existing techniques, developing an RVE based on the weighted Voronoi tessellation concept is the most favorable technique due to the least input parameters [8][9][10][11]. Based on this concept, an artificial microstructure model is developed here to predict the strain hardening behavior of DP-Steel [8], a micromorphology model is expanded under heat treatment [9], and a two-phase composite is characterized [10].
Due to advances in computational power, the crystal plasticity finite element (CPFE) approach has successfully been applied to simulate the plastic deformation behavior of various crystalline materials in recent decades [12], from single crystal [13], polycrystal [14] to multicrystalline materials [15]. The most significant benefit of using the CPFE method is that it can consider the anisotropy in the individual grains' deformation behavior, which is a consequence of their crystallographic orientations with respect to the main loading axis. From a mechanical point of view, mesoscale approaches to the multiscale modeling of granular materials are fundamentally based on the development of crystal plasticity theories [16,17]. A universal characteristic of these theories is the explicit modeling of slip systems within the crystal lattice to constitute a model that naturally accounts for plastic slip and dislocation motion and their interaction. Even though the widespread implementation of a dislocation-based multiscale plasticity model is generally rather limited due to the tremendous need for computational resources, crystal plasticity theories provide a robust theoretical framework leading to the development of better phenomenological plasticity models [18,19].
Although the influence of deformation gradients has been ignored in some of the proposed crystal plasticity theories, size effects are crucial in some applications based on experimental results, such as polycrystalline nickel bending [20], single-crystal copper, and single-crystal aluminum microbending experiments [21,22], and polycrystalline copper twisting [23]. Advanced non-local constitutive models have been developed to consider the effect of the deformation gradient, which is responsible for different types of size effects, such as bending and twisting of a polycrystalline. It mostly occurs at grain boundaries between polycrystal grains whose misorientations cause a large deformation mismatch in their behavior. Hence, the produced internal stresses due to mechanical loads are characterized by non-local plasticity models. Most of these constituent models are derived from the principle of the geometrically necessary dislocation (GND) density tensor [24][25][26].
On the one hand, the mechanical response of polycrystalline materials using non-local crystal plasticity theory is very time-consuming and computationally ineffective. On the other hand, the large parameter space for the non-local crystal plasticity model makes it difficult to parametrize all the parameters [27]. Therefore, it will be beneficial to simulate the material response by using empirical hardening functions [28], which have the advantage of consisting of just a few parameters while precisely capturing the hardening behavior of the material during plastic deformation. For instance, the hardening behavior of singlecrystalline niobium was simulated by finding a proper model which could identify crystal plasticity parameters that matched well with experimental results [29]. Many classical work hardening laws were used to describe the mechanical behavior of different materials [30]. Based on experimental data, all parameters defined in the different work hardening laws have been identified precisely.
In the present study, micromechanical analyses are conducted by using the homogenization method and creating an RVE to investigate the effect of non-local crystal plasticity parameters on the deformation behavior of a polycrystalline solid. Then the relation of non-local crystal plasticity parameters with macroscopic hardening parameters has been obtained through a scale-bridging approach to represent the macroscopic response of the material.
In Section 2, a numerical model consisting of 100 grains is developed to simulate the polycrystals, and the implemented non-local crystal plasticity model describes their behavior. The homogenization technique is introduced to produce periodic boundary conditions and to derive the stresses and strains at the macroscopic level. In the last part of Section 2, a 4-parameter Voce-type hardening law is proposed to mimic the plastic behavior of the polycrystals by the non-local crystal plasticity model. In Section 3, the influence of the non-local crystal plasticity parameters is studied by changing them in a range, and the trends of the proposed hardening law parameters are obtained by changing the non-local crystal plasticity parameters.

RVE Generation
For this study, quasi-2D RVEs are created using the dynamic microstructure generator (DMG) method [31][32][33], as depicted in Figure 1. This in-house software couples a particle simulation method with a radical Voronoi tessellation algorithm. As input parameters, it requires the target grain size distribution determined via a log-normal distribution in the form of the average grain diameter and the standard deviation. The number and size of spheres are calculated to mimic the targeted grain size distribution with the prescribed distribution parameters. In the next step, the spheres are randomly placed into a large finite volume, which is compressed. During this process, the spheres are allowed to move freely under a repulsive potential to avoid overlapping. The updated sphere positions and the diameter of every sphere from all time steps are then forwarded to a radical Voronoi tessellation algorithm from the open-source software Voro++ [33] to create the RVE. The resulting grain size distribution of the RVE obtained from each time step is then compared with the target distribution, and the RVE with the minimum difference is selected. The geometry of the 2D RVE is then extruded for 1% of a side length and meshed with 8-node linear brick elements (C3D8) by using CUBIT [34]. Section 2, a 4-parameter Voce-type hardening law is proposed to mimic the plastic behav-ior of the polycrystals by the non-local crystal plasticity model. In Section 3, the influence of the non-local crystal plasticity parameters is studied by changing them in a range, and the trends of the proposed hardening law parameters are obtained by changing the nonlocal crystal plasticity parameters.

RVE Generation
For this study, quasi-2D RVEs are created using the dynamic microstructure generator (DMG) method [31][32][33], as depicted in Figure 1. This in-house software couples a particle simulation method with a radical Voronoi tessellation algorithm. As input parameters, it requires the target grain size distribution determined via a log-normal distribution in the form of the average grain diameter and the standard deviation. The number and size of spheres are calculated to mimic the targeted grain size distribution with the prescribed distribution parameters. In the next step, the spheres are randomly placed into a large finite volume, which is compressed. During this process, the spheres are allowed to move freely under a repulsive potential to avoid overlapping. The updated sphere positions and the diameter of every sphere from all time steps are then forwarded to a radical Voronoi tessellation algorithm from the open-source software Voro++ [33] to create the RVE. The resulting grain size distribution of the RVE obtained from each time step is then compared with the target distribution, and the RVE with the minimum difference is selected. The geometry of the 2D RVE is then extruded for 1% of a side length and meshed with 8-node linear brick elements (C3D8) by using CUBIT [34].
Based on the commercial finite element packaging software, ABAQUS [35], and implementing many subroutines, the material's mechanical response in terms of stress-strain curves is obtained.
From mesh sensitivity analyses performed on the generated RVE consisting of 100 grains with an average grain size of 50 µm and a standard deviation of 0.1 µm, the optimum mesh size of approximately 0.003 mm is adopted, which totally included 9796 elements. The 100-grain RVE is meshed with 8-node linear brick elements (C3D8) and has a thickness of one element with the size of 0.002 mm. Therefore, the simulations are assumed to be quasi-2D under plane stress conditions.  Based on the commercial finite element packaging software, ABAQUS [35], and implementing many subroutines, the material's mechanical response in terms of stressstrain curves is obtained.

Nonlocal Crystal Plasticity Model
From mesh sensitivity analyses performed on the generated RVE consisting of 100 grains with an average grain size of 50 µm and a standard deviation of 0.1 µm, the optimum mesh size of approximately 0.003 mm is adopted, which totally included 9796 elements. The 100-grain RVE is meshed with 8-node linear brick elements (C3D8) and has a thickness of one element with the size of 0.002 mm. Therefore, the simulations are assumed to be quasi-2D under plane stress conditions.

Nonlocal Crystal Plasticity Model
The material considered in this study is Ferrite, which is categorized as a bodycentered cubic (BCC) form of iron, and its behavior is characterized by a non-local crystal plasticity model developed by a user-defined material subroutine (UMAT) in ABAQUS. With the assumption of large displacements, the total deformation gradient tensor, F, is separated by multiplicative decomposition as where F e and F p are the elastic and the plastic part of the deformation gradient tensor. The rate of plastic deformation is related to the plastic part of the deformation gradient tensor as below where L p is the plastic part of the gradient velocity tensor and in the case of dislocation slip as the only deformation process, results in where . γ α is the slip rate and M α = d α n α defines the Schmid tensor in the intermediate configuration for the slip system α, which is defined by the slip direction d α and the slip plane normal n α . The symbol ⊗ denotes the dyadic product of two vectors resulting in a second-rank tensor. N counts the total number of slip systems. Based on large deformation theory and by introducing the stiffness tensor C, the elastic response can be formulated as where S is the second Piola-Kirchhoff stress tensor in the intermediate configuration. The Cauchy stress in the current configuration is defined as The plastic deformation mechanism here is governed by the slip mechanism where dislocations slip in well-designed slip systems. In the plastic regime, the flow rule and the strain hardening law are defined as below: .
where . γ 0 is the reference shear rate, p 1 the inverse value of the strain rate sensitivity, h 0 the initial hardening rate,τ sat the saturation slip resistance, and p 2 a fitting parameter. χ αβ is the cross-hardening matrix between crystallographic mobile dislocations and super GNDs.
The flow rule described in Equation (6) includes two back stresses,τ GNDi α and τ GNDk α which define the additional hardening caused by GNDs due to strain gradients [26]. This additional hardening can be separated into isotropic (τ GNDi α ) and kinematic (τ GNDk α ) hardening contributions.
In the case of treating F p as additional degree of freedom (DOF) to consider the nonlocal effect, [36,37], it is possible to calculate the dislocation density tensor in the reference configuration as following The net Burgers vector b can be determined with the help of the dislocation density tensor for an arbitrary unit area with a normal vector n, [38] as b = Gn.

of 18
The dislocation density tensor is projected to the global Cartesian coordinates of the system, and then geometrically necessary super dislocations are defined individually. Then, the far-field stress of the crystallographic GND population can be described with good accuracy [26]. In this way, the GND density tensor is separated into nine independent parts where d α and t α are permutations of the Cartesian unit vectors and b is the norm of the Burgers vector. ρ α is named as super GND density, of which the three first densities belong to the screw super GND densities, and the last six ones represent the edge super GND densities. Due to strong cross hardening produced by GNDs, the additional passing stress has to be considered [39] for mobile dislocations caused by super GNDŝ where c 1 is a geometrical factor [26] and µ is the shear modulus. With the assumption of small elastic strains, the resolved shear stress, τ α , and the back stress, τ GND α , within the intermediate configuration are written as where S GND is the internal stress in the intermediate configurations [26].
In the present study, the parameters c 1 , τ 0 , . γ 0 , p 2 ,τ sat and h 0 are changed in a specified range, and their influences are evaluated on the equivalent stress-equivalent plastic strain diagram. The equivalent stress is defined according to the von Mises yield criterion, which in a Cartesian coordinate system (along x, y, and z-direction) is defined as and the equivalent plastic strain is written in a manner consistent with the equivalent von Mises stress as follows

Homogenization Technique
In simulations, to preserve the periodicity during deformation, periodic boundary conditions are applied [8] such that from one side of the model, the DOFs of the selected nodes are coupled to those of their counterparts on the opposite side, and due to cross-contraction effects, lateral shrinkage is not constrained. The full description of the implemented formulations of the periodic boundary condition applied in this work can be found in [31].
Periodic boundary conditions in the Cartesian coordinate system have been used for displacement and DOF for the non-local effect, i.e., the plastic deformation gradient. In order to apply periodic boundary conditions, opposite surfaces must have exactly the same surface mesh, which requires a well-controlled meshing process.
The macroscopic loading condition of the RVE is quantified by transferring the mechanical boundary conditions of the model into macroscopic stress-and strain tensors. Therefore, to implement periodic boundary conditions, the deformation state is applied to the four reference nodes V 1 , V 2 , V 4 , and H 1 as visualized in Figure 2. We constrained displacements along x, y, and z directions at V 1 in order to hinder rigid body motion and along x and y directions at H 1 . Displacement along y and x directions, respectively, are fixed at nodes V 2 and V 4 , together with rotation around z-axis.
where ∆x, ∆y, and ∆z are the periodic box dimensions in the global Cartesian coordinate system. The macroscopic stress tensor is calculated by using the reaction force vectors at the four reference nodes and the current nodal position vectors x node of the reference nodes [5] as follows: where the symmetric tensor is defined as sym = 1/2[A + A T ] for tensor A and its transposed A T .

Empirical Hardening Relation
An empirical relation is proposed to mimic the work hardening behavior of the considered material in the present study,. From the many available empirical and mathematical descriptions of the strain hardening phenomenon, a Voce-type hardening relation is selected to project the strain hardening of the Ferrite as follows. Since the specified boundary nodes of the model track the kinetics of these reference nodes, the macroscopic strain tensor can easily be calculated from the nodal displacements of these nodes [31]. The mathematical description of the macroscopic strain tensor with the assumption of small macroscopic deformations and its corresponding macroscopic stress tensor has been formulated from the nodal displacements u node i and the reaction force vectors (F node ) of these four nodes, respectively [5]. Due to the periodic coupling, the reaction forces at these four nodes are the resultant of all forces on the considered surface representing directly the resultant of all forces needed to determine the stress tensor.
where ∆x, ∆y, and ∆z are the periodic box dimensions in the global Cartesian coordinate system. The macroscopic stress tensor is calculated by using the reaction force vectors at the four reference nodes and the current nodal position vectors x node of the reference nodes [5] as follows: where the symmetric tensor is defined as sym = 1/2[A + A T ] for tensor A and its transposed A T .

Empirical Hardening Relation
An empirical relation is proposed to mimic the work hardening behavior of the considered material in the present study,. From the many available empirical and mathematical descriptions of the strain hardening phenomenon, a Voce-type hardening relation is selected to project the strain hardening of the Ferrite as follows.
where σ ε p represents the equivalent stress and ε p is the equivalent plastic strain. a, b and c are hardening parameters, and σ y represents the yield limit. The most influence of parameter a is on the slope of the hardening during plastic deformation and an increase of a results in a higher level of equivalent stress. Parameter b changes the hardening slope with a higher rate compared to a. Parameter c affects both the yield limit and the hardening behavior such that higher values of c lead to a lower quantity for the yield limit and to a change of the hardening behavior from exponential to linear. The fitting procedure of the hardening relation defined in Equation (18) is described in a MATLAB environment [40]. The "fit" function has been chosen from MATLAB tools and uses the hardening function as a guess function; then, its parameters are identified by defining a lower and upper bound for parameters and making an initial guess.
The curve-fitting toolbox software uses the "nonlinear least-squares" formulation to fit a nonlinear model to data. A nonlinear model is defined as a nonlinear equation in the coefficients or a combination of linear and nonlinear in the coefficients.
In matrix form, nonlinear models are given by the formula where y is an n-by-1 vector of responses, f is a function of λ and x, λ is an m-by-1 vector of coefficients, x is the n-by-m design matrix for the model, and is an n-by-1 vector of errors.
Since estimating the coefficients using simple matrix techniques is difficult in nonlinear models, an iterative approach is required [40]. In the following fitted results, the error is less than 0.1%.

Results and Discussion
The considered non-local crystal plasticity parameters are changed in a specific and physically meaningful range to study their influence on the hardening behavior of the material and their relations with the hardening parameters defined in Equation (18). Based on the macroscopic response using the crystal plasticity model, the hardening parameters are obtained by the nonlinear least-squares method.
The default values of the parameters are listed in Table 1, and to study the effect of each parameter on the material response, the corresponding parameter is changed, and the others will keep unchanged and equal to the default value.  The results of changing the initial critical resolved shear stress, τ 0 , on the equivalent, stress-equivalent plastic strain curve are shown in Figure 3. Parameter τ 0 has been varied in the range of 50-70 MPa with steps of 2 MPa. It is evident that, with the increase in τ 0 , both the initial yield stress and the stresses during the hardening of the material will increase with approximately the same hardening slope. As a whole, changing the τ 0 parameter will affect the overall response of the material, improve the mechanical capability of the considered material, and shift the overall response to the higher stress level. Furthermore, based on the hardening function and using the fitting procedure, the material's behavior is identified as largely comparable with the equivalent stress-equivalent plastic strain graphs. Figure 4 visualizes the trend of changes in identified Voce-type hardening parameters with different τ 0 values. An increase of the τ 0 in the considered range leads to a rising trend in both a and σ y although the trend is different as a nonlinear and linear function, respectively (Figure 4a,b). A nonlinear trend is observed in b and c parameters by an increase in the τ 0 . Visibly, the wide range of change in parameters a and b by varying the τ 0 is highlighting the mechanical behavior observed in Figure 3 that the hardening behavior has been most affected. by an increase in the τ 0 . Visibly, the wide range of change in parameters a and b by varying the τ 0 is highlighting the mechanical behavior observed in Figure 3 that the hardening behavior has been most affected. by an increase in the τ 0 . Visibly, the wide range of change in parameters a and b by varying the τ 0 is highlighting the mechanical behavior observed in Figure 3 that the hardening behavior has been most affected.  In the following, the results of changing the saturation resolved shear stress, τ sat , on the equivalent, stress-equivalent plastic strains are shown. Here, parameter τ sat has been changed in the range of 300-600 MPa with steps of 50 MPa, of which the results are visualized in Figure 5. Visibly, with an increase in τ sat , stresses during material hardening will In the following, the results of changing the saturation resolved shear stress,τ sat , on the equivalent, stress-equivalent plastic strains are shown. Here, parameterτ sat has been changed in the range of 300-600 MPa with steps of 50 MPa, of which the results are visualized in Figure 5. Visibly, with an increase inτ sat , stresses during material hardening will increase, and the initial yield stress remains unchanged. In addition, the hardening slope of the equivalent stress-plastic strain curves goes up.
In the following, the results of changing the saturation resolved shear stress, τ sat , on the equivalent, stress-equivalent plastic strains are shown. Here, parameter τ sat has been changed in the range of 300-600 MPa with steps of 50 MPa, of which the results are visualized in Figure 5. Visibly, with an increase in τ sat , stresses during material hardening will increase, and the initial yield stress remains unchanged. In addition, the hardening slope of the equivalent stress-plastic strain curves goes up.
The equivalent stress-equivalent plastic strain curves represented in Figure 5 are fitted well, using the Voce-type hardening law with numerical results based on the crystal plasticity model. The trend of change in the Voce-type hardening law parameters with different saturation resolved shear stresses are visualized in Figure 6, where a and σ y decrease and b and c increase while the trend for the two latter parameters is quadratic in nature. The yield limit, σ y in the Voce-type hardening law is mostly unaffected by a change in the saturation resolved shear stress as also has been represented in Figure 5.  The equivalent stress-equivalent plastic strain curves represented in Figure 5 are fitted well, using the Voce-type hardening law with numerical results based on the crystal plasticity model. The trend of change in the Voce-type hardening law parameters with different saturation resolved shear stresses are visualized in Figure 6, where a and σ y decrease and b and c increase while the trend for the two latter parameters is quadratic in nature. The yield limit, σ y in the Voce-type hardening law is mostly unaffected by a change in the saturation resolved shear stress as also has been represented in Figure 5.
In the following, the results of changing the initial hardening rate, h 0 , from 300 to 600 MPa on equivalent stress-equivalent plastic strain curves are depicted in Figure 7. As the results highlight, with an increase in h 0 , stresses during the hardening of the material will increase approximately with the same slope, and the initial yield stress remains unchanged. Figure 7 also emphasizes that the proposed hardening function can describe the material response accurately.
Adopting different h 0 will cause various trends in the hardening parameters, see Figure 8. With a change in the h 0 in the defined range, a reduces at the initial level and remains constant for the larger h 0 as represented in Figure 6a. On the other hand, σ y illustrates a different behavior such that a very small increase occurs and then it decreases, mimicking a quadratic function. b and c hardening parameters both go up in quadratic form by increasing h 0 (see Figure 8b-d). Among the Voce-type hardening parameters, σ y and c are not changing significantly by various initial hardening rates. These two parameters are mostly affecting the yield limit, as it is verifiable in Figure 7. In the following, the results of changing the initial hardening rate, h 0 , from 300 to 600 MPa on equivalent stress-equivalent plastic strain curves are depicted in Figure 7. As the results highlight, with an increase in h 0 , stresses during the hardening of the material will increase approximately with the same slope, and the initial yield stress remains unchanged. Figure 7 also emphasizes that the proposed hardening function can describe the material response accurately.
Adopting different h 0 will cause various trends in the hardening parameters, see Figure 8. With a change in the h 0 in the defined range, a reduces at the initial level and remains constant for the larger h 0 as represented in Figure 6a. On the other hand, σ y illustrates a different behavior such that a very small increase occurs and then it decreases, mimicking a quadratic function. b and c hardening parameters both go up in quadratic form by increasing h 0 (see Figure 8b-d). Among the Voce-type hardening parameters, σ y and c are not changing significantly by various initial hardening rates. These two parameters are mostly affecting the yield limit, as it is verifiable in Figure 7.    Another crystal plasticity parameter and its influence that was studied is the exponent of strain hardening, p 2 . It is changed in the range of 8-12 with steps of 0.5, and the results are shown in Figure 9 together with the fitted behavior using hardening parameters. As p 2 increases, stress during hardening of the material decreases, but the initial yield stress remains unchanged, and the hardening slope decreases to some extent. Figure 10 shows the trend of changes in identified hardening parameters with different p 2 . With an increase in p 2 , a and σ y go up exponentially and like a quadratic function, respectively, but b and c diminish linearly. By adopting different values of p 2 , in the Voce-type hardening law, a, c and σ y will not vary drastically, as the yield limit was unaffected and the hardening rate did not change significantly compared to the change in other crystal plasticity parameters. 2 yield stress remains unchanged, and the hardening slope decreases to some extent. Figure  10 shows the trend of changes in identified hardening parameters with different p 2 . With an increase in p 2 , a and σ y go up exponentially and like a quadratic function, respectively, but b and c diminish linearly. By adopting different values of p 2 , in the Voce-type hardening law, a, c and σ y will not vary drastically, as the yield limit was unaffected and the hardening rate did not change significantly compared to the change in other crystal plasticity parameters. tively, but b and c diminish linearly. By adopting different values of p 2 , in the Voce-type hardening law, a, c and σ y will not vary drastically, as the yield limit was unaffected and the hardening rate did not change significantly compared to the change in other crystal plasticity parameters. The next crystal plasticity parameter in this study is the reference shear rate, γ̇0. It has been changed in the range of 0.001-0.005 with steps of 0.0005, of which the results are visualized in Figure 11. It is seen that with an increase in γ̇0, both, the initial yield limit and the stresses during the hardening of the material, decrease, but the reduction trend is The next crystal plasticity parameter in this study is the reference shear rate, . γ 0 . It has been changed in the range of 0.001-0.005 with steps of 0.0005, of which the results are visualized in Figure 11. It is seen that with an increase in . γ 0 , both, the initial yield limit and the stresses during the hardening of the material, decrease, but the reduction trend is not significant, and also the hardening slope remains unchanged. Therefore, it can be concluded that the overall behavior is not sensitive to a change in the reference shear rates. On the other hand, the hardening parameters could show results coming close to the numerical results depicted in Figure 11, and the trend of change in these parameters with different reference shear rates is visualized in Figure 12, where a and σ y decrease and b and c increase while the trend for all them is quadratic. By changing the reference shear rate, the macroscopic hardening parameters do not change in a wide range, as predicted from the stress-strain curve in Figure 11. The next crystal plasticity parameter in this study is the reference shear rate, γ̇0. It has been changed in the range of 0.001-0.005 with steps of 0.0005, of which the results are visualized in Figure 11. It is seen that with an increase in γ̇0, both, the initial yield limit and the stresses during the hardening of the material, decrease, but the reduction trend is not significant, and also the hardening slope remains unchanged. Therefore, it can be concluded that the overall behavior is not sensitive to a change in the reference shear rates. On the other hand, the hardening parameters could show results coming close to the numerical results depicted in Figure 11, and the trend of change in these parameters with different reference shear rates is visualized in Figure 12, where a and σ y decrease and b and c increase while the trend for all them is quadratic. By changing the reference shear rate, the macroscopic hardening parameters do not change in a wide range, as predicted from the stress-strain curve in Figure 11.  The influence of both, saturating critical resolved shear stress,τ sat , and geometrical factor, c 1 , was investigated on the equivalent stress versus equivalent plastic strain of the considered material. Therefore, their values are changed, and the material response is obtained, as illustrated in Figure 13. Asτ sat increases, stresses during the hardening of the material will increase, and the initial yield stress remains unchanged. In addition, the hardening slope of the equivalent stress-equivalent plastic strain curves increased with the value ofτ sat . On the other hand, adopting higher quantities for the geometrical factor results in an increase in the initial yield limit, the stresses during the hardening of the material, and the slope of the hardening behavior.
The trend of variation in the identified hardening parameters due to the simultaneous change of bothτ sat and c 1 parameters is visualized in Figure 14. It is seen that higher geometrical factors and lower values of saturating shear stress lead to an increase in the hardening parameter. Therefore, a is very sensitive to a higher contribution of c 1 and lower quantities ofτ sat . On the contrary, variation of the other hardening parameters, σ y , b and c, is mostly affected by lower c 1 and higherτ sat . The influence of both, saturating critical resolved shear stress, τ sat , and geometrical factor, c 1 , was investigated on the equivalent stress versus equivalent plastic strain of the considered material. Therefore, their values are changed, and the material response is obtained, as illustrated in Figure 13. As τ sat increases, stresses during the hardening of the material will increase, and the initial yield stress remains unchanged. In addition, the hardening slope of the equivalent stress-equivalent plastic strain curves increased with the value of τ sat . On the other hand, adopting higher quantities for the geometrical factor results in an increase in the initial yield limit, the stresses during the hardening of the material, and the slope of the hardening behavior.
The trend of variation in the identified hardening parameters due to the simultaneous change of both τ sat and c 1 parameters is visualized in Figure 14. It is seen that higher geometrical factors and lower values of saturating shear stress lead to an increase in the hardening parameter. Therefore, a is very sensitive to a higher contribution of c 1 and lower quantities of τ sat . On the contrary, variation of the other hardening parameters, σ y , b and c, is mostly affected by lower c 1 and higher τ sat .     Figure 15 visualizes a heatmap revealing the correlation between the crystal plasticity parameters and the Voce-type hardening parameters. The light yellow and dark blue colors represent the positive and negative linear relation between these parameters. The heatmap highlights that almost all crystal plasticity parameters have a positive/negative linear correlation with the yield limit, σ y . Among all Voce-type hardening parameters, a has the least linear correlation with the crystal plasticity parameters.
Crystals 2021, 11, x FOR PEER REVIEW 16 of 18 Figure 15 visualizes a heatmap revealing the correlation between the crystal plasticity parameters and the Voce-type hardening parameters. The light yellow and dark blue colors represent the positive and negative linear relation between these parameters. The heatmap highlights that almost all crystal plasticity parameters have a positive/negative linear correlation with the yield limit, σ y . Among all Voce-type hardening parameters, a has the least linear correlation with the crystal plasticity parameters.
The represented heatmap leads to a proper understanding of the correlation between the crystal plasticity parameters and the macroscopic material behavior, which improves the fitting of crystal plasticity parameters to experimental stress-strain curves. Figure 15. A heatmap representing the correlation between the crystal plasticity parameters and the Voce-type hardening parameters.

Conclusions
This study investigates the effect of different non-local crystal plasticity parameters on polycrystals' resulting macroscopic equivalent stress-equivalent plastic strain curves. The non-local crystal plasticity parameters are changed in a specific range that is physi- Figure 15. A heatmap representing the correlation between the crystal plasticity parameters and the Voce-type hardening parameters.
The represented heatmap leads to a proper understanding of the correlation between the crystal plasticity parameters and the macroscopic material behavior, which improves the fitting of crystal plasticity parameters to experimental stress-strain curves.

Conclusions
This study investigates the effect of different non-local crystal plasticity parameters on polycrystals' resulting macroscopic equivalent stress-equivalent plastic strain curves. The non-local crystal plasticity parameters are changed in a specific range that is physically meaningful. The macroscopic behavior of the material is fitted by an empirical strain hardening relation using the nonlinear least-squares method, and the trend of variations in the hardening parameters induced by a change in material properties has been obtained.
The crystal plasticity properties that varied in this study include the initial critical resolved shear stress, τ 0 , the saturation critical resolved shear stress,τ sat , the initial hardening rate, h 0 , the exponent of strain hardening, p 2 , the reference shear rate, . γ 0 , and the geometrical factor, c 1 . The influence of all these parameters on the effective stress-strain behavior of the polycrystal has been quantified by fitting the resulting polycrystal flow curves with a Voce-type hardening model.
Changing the τ 0 parameter led to an increase in both, the initial yield limit and the hardening rate of the polycrystal and, thus, generally to a shift in the material response towards a higher stress level. Varying the value of h 0 revealed no change in the initial yield stress but an increase in the hardening behavior with the same slope. By increasing p 2 , the initial yield stress did not change, but the hardening rate decreased significantly. Higher . γ 0 quantities reduced both, the initial yield stress and the hardening behavior of the material, although this reduction is very insignificant, and it can be concluded that the material behavior is unaffected by an increase in . γ 0 values in the defined range. It is noted here, however, that this parameter might influence the loading rate dependence of the strength. Results showed that, with an increase in c 1 , the initial yield limit and the hardening rate of the material increased significantly.
To represent the correlation between crystal plasticity and Voce-type hardening parameters, a heatmap was created that indicates that all crystal plasticity parameters strongly correlate with the yield limit and the hardening law. These partly nonlinear relationships between crystal plasticity parameters and the global mechanical properties of polycrystals are clarified in this work.