Crystal Plasticity Finite Element Simulation Study on Thread Rolling Pure Aluminum

In order to study the effects of grain orientation on thread rolling deformation, the mechanical performance was analyzed from a microscopic perspective. The Voronoi method was employed to establish a polycrystalline model and a user-defined material subroutine (UMAT) was adopted. Meanwhile the crystal plasticity finite element method (CPFEM) was used to simulate the deformation behavior of the grain on thread rolling process. The simulation results show that the deformation degree of the grain is different due to the different initial grain orientation and the representative element strain curve is more accurate after homogenization. It is also indicated that the initially activated slip systems cause the movement of the other slip systems. By analyzing the rotation of the grain orientation and the evolution of the texture through the pole figure, it can be seen the CPFEM results show high consistency with the test results of the electron back scatter diffraction (EBSD).


Introduction
With the increase in the miniaturization of metal materials, the higher requirements are put forward for high-quality products. When the size is reduced to the micro-scale, the size effects [1,2] are significant. In the plastic deformation process of polycrystalline materials, the grains of random orientation will appear to be preferentially oriented, and the material will show anisotropy in the macro-scale [3]. The grains concentrated in the same orientation positions will form texture [4,5]. The appearance of texture will have an important impact on the mechanical properties of the material. It is of great significance to study the grain orientation and dislocation slip in the micro-scale in the forming process of metal materials.
The material is regarded as isotropic in the general finite element simulation, and the macroscopic deformation behavior of the material is studied. A new simulation method needs to be introduced when it comes to micro-scale. The finite element simulation shifts from macro-scale to micro-scale, from the general finite element simulation to the crystal plasticity finite element simulation [6,7]. CPFEM is based on the continuum theory to homogenize discrete dislocations and analyze the effects of dislocation slip, which reflects the microcosmic characteristics of materials [8][9][10]. The biggest advantage is that it can describe the orientation change of each grain during the deformation process and the evolution of different slip systems. Through drawing pole figures that can describe the grain orientations, the deformation behaviors of crystalline materials are explained clearly [11][12][13][14][15]. The crystal materials are divided into the single crystal and the polycrystalline. The deformation ways include slip and twinning. In polycrystalline materials, the deformation method is much more complicated than the single crystal due to the difference in grain size, shape, orientation, and the existence of grain boundaries. In general, the dislocation slip is major, while the twinning is difficult. The critical shear stress required for the twinning is much greater than the slip. The twinning begins only when the slip is difficult and the critical shear stress of the twinning is reached [16][17][18]. In the present paper, the plastic deformation of material was simulated and discussed based on the dislocation slip mechanism without considering the twinning deformation. The mechanisms involved in the macro and micro-deformation of materials have been studied by many scholars using CPFEM. Chen et al. [19] used CPFEM to study the copper rolling, and the results showed that the motion of the slip system is different. Pi et al. [20,21] used CPFEM to study textures appearing during the deformation of pure aluminum. It was concluded that the textures of simulation were consistent with the experimental results. Hama and Takuda [22,23] analyzed unloading behavior of magnesium alloy during the tensile deformation process. They found that only the base slip systems were activated during unloading. Alankar et al. [24] and Ha et al. [25] simulated the deformation of pure aluminum by using a dislocationdensity-based crystal plasticity model. This model showed a reasonable attempt at incorporating microstructural considerations into the model along with the kinematics of crystal plasticity to form a more physically realistic prediction. Li et al. [26] studied the rotation of grains during uniaxial deformation. The results showed that the initial orientation determined the degree of grain rotation.
Based on the above literature reviews, it is clear that the molding method is relatively simple and has limitations. In the present paper, the biggest problem in the calculation process is the convergence of the analysis. The thread rolling involved the contact problem requires reasonable parameter settings to make the analysis normal. The polycrystalline model is established based on CPFEM and the Voronoi polyhedron theory, which is used to study the effects of the grain orientations during the inhomogeneous thread rolling deformation and evolution of slip systems as well as rotation of the grain under such cyclical and complicated stress conditions.

Kinematics model
The theory of crystal plasticity follows the pioneering work of Talor [27]. The plastic constitutive equations are initially proposed by Hill and Rice [28]. The deformation is split into elastic deformation and plastic deformation. The elastic deformation is caused by lattice distortion and rigid rotation, and the plastic deformation is a series of shear behaviors caused by the dislocation slip. The total deformation gradient can be decomposed into the plastic component of sliding along the slip system and the elastic component of lattice distortion and rigid rotation. The total deformation gradient is given as follows: Where and * are the plastic and the elastic deformation gradient, respectively. The elastic constitutive relation can be expressed as follows: Where is the second Piola-Kirchhoff stress; C e is the fourth-order anisotropic elasticity matrix; E is the Lagrangian elastic strain tensor; σ is the Cauchy stress tensor. ̇ is the rate of change of ; ( ) −1 is the inverse matrix of ; ̇( ) is the shear strain rate of the slip system α. The relationship between ̇a nd ̇( ) is given as follows: Where ( ) is the slip direction vector of the slip plane in the slip system α; ( ) is the normal direction vector of the slip plane in the slip system α. The summation range includes all the slip systems activated.

Constitutive model
Assuming that the crystal slip follows the Schmid's law, the relationship between and the resolved shear stress is given as follows: Where m * (α) and S * (α) are the unit normal direction vector of slip plane and the unit slip direction vector in the slip system α, respectively. ρ 0 and ρ represent the mass density of the reference state and the current state, respectively.
Crystal plasticity theory is completed by relating the shear strain rate ̇( ) to the stress . The generalized Schmid's law states that slip occurs if the resolved shear stress = : reaches a critical value. is the Schmid tensor. In the rate-dependent crystal plasticity model, the exponential equation describes the relationship between the shear strain ratė( ) and the resolved shear stress τ α of the slip system α. The formula is given as follows: Where is the critical resolved shear stress of the slip system α; ̇0 ( ) is the reference shear strain rate; is the strain rate sensitive coefficient. = 0 and = ∞ correspond to viscoelasticity and rateindependent materials, respectively. The critical shear rate ̇ is given as follows: Where ℎ is the slip hardening modulus; ̇( ) is the shear strain rate of the slip system . represents the total number of slip systems; The above formula sums all activated slip systems. When α = , ℎ = ℎ , ℎ is called self-hardening modulus; When α ≠ , ℎ is called latent hardening modulus. The formula for calculating the self-hardening modulus is given by Peirce [29,30] et al. and Asaro et al. [31] as follows: The formula for calculating the latent hardening modulus is as follows:γ Where ℎ 0 represents the initial hardening modulus; represents the cumulative shear strain rates of all slip systems; 0 represents the initial critical shear stress; is the saturation value of the initial critical shear stress; is the latent hardening parameter. The above constitutive model is compiled by FORTRAN language into the UMAT for the CPFEM simulation, and the combination of finite element method and crystal plasticity constitutive theory is realized.

Finite element model
In this paper, the model is established by the NEPER [32] software. The NEPER needs less code than Python programming. By running the corresponding module, and importing the GEO file generated into the ABAQUS, the polycrystalline models can be established successfully. The models are shown in Figure 1. After generating a polycrystalline model with 50 grains, each grain is given a random orientation. On the actual thread rolling experiment, the length of the aluminum bar changed from 10mm to 21mm, and the diameter changed from 3mm to 2mm. The thread rolling model is shown in Figure 2. The deformation of the polycrystalline aluminum with 50 grains is simulated. The size of the upper and lower plates is 200mm×20mm; the thickness is 5mm; the diameter of the aluminum bar to be formed is 3mm; the length is 10mm; the total number of nodes is 8236; the total number of elements is 5302; the element type is C3D10. The lower plate is fixed on the Y and Z axis; the upper plate is fixed on the Y axis and given a displacement in the negative direction of the Z axis; the upper and lower plates move in opposite directions on the X axis. During the processing, through the relative movement of the upper plate and the lower plate in the horizontal (X) and vertical (Z) directions, the metal bar is uniformly thinned and formed. In the present paper, the pure aluminum is taken as the research material. The elastic coefficients are C 11 = 108,000 MPa, C 12 = 62,000 MPa, and C 44 = 28,300 MPa. The parameters of the constitutive model are shown in the following Table 1.  In this paper, the simulation of thread rolling forming is obviously different from the uniaxial tension and compression process. During the thread rolling forming process, the stress should have a cyclical phenomenon duo to the rotation of the polycrystalline model. Figure 3 shows the stress-strain curve of the element 2032 (E2032) in the selected grain 20(G20). It is clear that the stress-strain curve is basically the same as the ordinary uniaxial loading in the deformation trend, but it is far more complicated. The stress peak value of this curve appears when E2032 is in contact with the plate. At the next moment, this element rotates to the next position and the stress value decreases. The bar becomes thinner and longer during the reciprocating motion. The length changed from 10mm to 12mm, and the diameter changed from 3mm to 2.5mm in the early stage of deformation. The bar no longer rotates due to the large deformation in the later stage of the simulation deformation.

Stress-Strain Curve Analysis
When the polycrystalline material is deformed, it can be seen that the difference of the deformation among grains is obvious. It can more truly reflect the mechanical properties of the polycrystalline material by extracting the mechanical response of the representative element. The volume average processing method is used in the polycrystalline region, and the mechanical response based on the homogenization theory is more accurate. After homogenization, the strain of the polycrystalline representative element with 50 grains is obtained in Figure 4. The strain curve obtained by calculating the extracted data is shown in Figure 4, which reflects the strain situation of the single grain and the homogenized representative element. The three curves in Figure 4 represent the change of the strain in G20, G35, and all homogenized grains. G20 and G35 are adjacent. Time is started counting from 1s. The equivalent strain value was zero before 0.04s, and the plastic deformation has not yet proceeded. Although there is a large difference in the strain among grains, the strain of the representative element is stable due to the coordinated deformation among grains. It is representative and accurate to better reflect the deformation behavior of the entire polycrystalline model.

Evolution of the slip systems
According to the crystal plasticity theory, whether the slip system is activated is controlled by the shear strain rate of the slip system. The pure aluminum is a face-centered cubic (FCC) structure with 12 slip systems shown in Table 2. In order to explain the dislocation slip mechanism of the thread rolling much better, the evolution of the slip system is specifically analyzed. The activation of the slip system in the grain inside and the grain boundary is different. The slip system is initially activated at the grain boundary, which causes the movement of the slip system at the grain inside, and finally forms a deformation band throughout the entire grain. The evolution of the slip system in different grains is shown in Figure 5.  Figure 5 shows the variation of the selected five slip systems in the different grains. Region A is located in the middle of the polycrystalline model, and region B is located in the end of the polycrystalline model. G20 and G35 belong to region A; G22 and G45 belong to region B. It can be seen that the movement of each slip system is completely different. The contact relationship was established before 1s. The slip systems did not start to move, and the shear strain rate was zero. The upper plate begins to press down, and the polycrystalline model performs the elasticity deformation at the initial stage of deformation. As the deformation continues, the slip system b1 of G20 begins to move along the [101] direction. The initially activated slip system causes the rotation of the slip plane and the slip direction, coupled with the interaction among adjacent grains, which results in the activation of slip systems a1, c3, and d2. The slip systems a1 and c3 move along the positive directions [01 ̅ 1] and [101 ̅ ] respectively, and the slip system d2 moves along the opposite direction [1 ̅ 01 ̅ ]. The shear stain rate of the initially activated slip system b1 is reduced under the influence of the later activated slip systems, and then moves along the opposite direction. With the deformation continuing, the shear strain rates of the activated slip systems are in the increase. It can be seen that the slip system c2 has not been activated during the entire deformation process. The slip system a1 of G35 is initially activated, which moves along the positive direction [01 ̅ 1]. As the deformation continues, the slip system c2 moves along the positive direction [110]. Then the slip is restricted and begins to move along the opposite direction [1 ̅ 1 ̅ 0]. The shear strain rates of b1 and d2 has obvious fluctuation. Finally, the slip system b1 moves along the positive direction [101], while the slip system d2 moves along the opposite direction [1 ̅ 01 ̅ ]. The increase rate of the slip system b1 is larger and the movement on this slip system is easier than the other slip systems, while the slip system c3 has not been activated from beginning to end. The movement of the slip systems is approximately the same in the grains of region B and region A. It can be seen that all the five slip systems are activated in G22 and G45 of region B. This is because the grains at the end of the polycrystalline model need less resistance due to coordinated deformation among grains and the slip system is easier to activate. The above shows that the movement of the slip systems in different grains is different due to the influence of grain location and adjacent grains. The easy slip plane in some grains may be the hard slip plane in the other grains. It is observed that some slip systems have not been activated in the deformation process. It can be concluded that the shear strain rates of the slip systems

Pole figures analysis
The lattice rotation can be well expressed by the pole figure. The pole figure of polycrystalline model with different numbers of grains is drawn before the deformation. The X, Y, and Z axes represent the transverse, rolling and normal directions respectively. The grain orientation distribution is uniform. The initial (100) pole figures are shown in Figure 6.  Figure 6 shows the pole figures with 20, 50, and 100 grains before the deformation. It can be seen that the distribution of poles is more uniform and the orientation is more random as the number of crystal grains increases. When the number of crystal grains is larger than 50, the grain orientation distribution is roughly uniform. The random orientation is imported into the polycrystalline model, and the orientation information is re-extracted after deformation. The pole figure after deformation is drawn for further analysis. The (100) pole figures of G20 and G35 before and after the thread rolling deformation of the orientation are shown in Figure 7. It is observed that the initial grain orientation is different in the pole figures of G20 and G35. The orientations of G20 rotates by discrete 10°-15° along the rolling direction (RD), and the rotation angle of the G35 along the RD is 15°-25°. The orientations of G20 diffuses 25°-35° along the normal direction (ND) and the orientations of G35 diffuses 15°-25° along the transverse direction (TD) and the ND. The dispersion of the G20 orientation pole figures is higher than G35. It can be seen that the rotation of each grain is inconsistent due to the difference in the initial orientation and position of the grains. It is shown that each grain is a single orientation and the poles are centrally distributed before the deformation. The grains rotate in the deformation process. The grains are not severely deformed due to coordinated deformation among grains, and the rotation appears near the initial orientation. It is concluded that the degree of internal deformation of the same grain is different, and the orientation is no longer single. In addition, the degree of the pole's dispersion is increased, and the rotation of different grains is different. The data obtained by the thread rolling simulation is converted into Euler angle, and the pole figures are drawn by using MATLAB. Figure 8 shows the changes of the pole figures before and after the thread rolling deformation. The maximum density value is 3.8 after deformation in Figure 8d, which is larger than 3.3 before deformation in Figure 8b. It can be seen that there is no obvious texture after deformation, and the grain orientation distribution is basically random. On the one hand, it is because the number of the grains is small, and the degree of deformation is relatively small, so the degree of grain orientation redistribution is not large. On the other hand, it is caused by this thread rolling deformation method. It is different from uniaxial loading. As the polycrystalline model continues to rotate, the stress state of the grains is periodic. The different grains will always reach the same stress state at different times, and there will be no obvious preferred orientation among the grains.  Taking into account the difficulty of calculation, we can increase the number of grains and adopt a large deformation rolling simulation in order to clearly analyze the evolution of texture and verify the correctness of the simulation results. The EBSD test is executed after rolling experiment. The {001} <100> and {112} <111> rolling textures are included in Figure 9a measured in the EBSD test. Figure  9b shows the results obtained by calculation of the polycrystalline model with 512 grains. The rolling experiment can be regarded as plane-strain compression in the simulation, and the finite element model is shown in Figure 9c. The initial orientation is random, and the texture appears after rolling deformation. It can be seen that the distribution of the grain orientation and the symmetry of the texture exist a high consistency. It is concluded that the experimental results are basically consistent with CPFEM simulation results.

Conclusions
This study aimed at thoroughly understanding the deformation mechanisms during the thread rolling forming. The deformation behavior is simulated by CPFEM. The study aimed at learning the inhomogeneous deformation caused by the grain orientation. The strain curve of polycrystalline representative element was obtained by the average volume method. The deformation of each grain was studied by analyzing the evolution of slip systems and the rotation of grain orientation. The pole figures were drawn to analyze the evolution of texture. The following conclusions are drawn: 1. CPFEM is used in thread rolling forming to study the orientation behavior of grain. The volume average processing method is used in the polycrystalline region. The strain curve of representative element is more accurate after homogenization.
2. The activation and movement state of each slip system are determined by the various initial grain orientations and grain locations. The grains at the end of the polycrystalline model need less resistance due to coordinated deformation and the slip system is easier to activate.
3. The pole figure shows that the grains are rotated. The pole figure obtained by the experiment is basically consistent with the pole figure by the simulation, which verifies the accuracy of CPFEM.