Cold Rolling Texture Prediction Using Finite Element Simulation with Zooming Analysis

Cold rolling is widely employed in the manufacturing industry for the production of metal plates. In the cold rolling process, the thickness reduction of the metal plate under the recrystallization temperature generates severe anisotropy; this influences the subsequent forming processes. Therefore, the generation and prediction of metal plate anisotropy during cold rolling is a highly interesting research topic involving upstream studies of sheet metal forming. In this study, using the finite element method with zooming analysis, we established an efficient elastic–plastic analysis method to predict the metal plate texture after cold rolling. This method for cold rolling texture prediction was confirmed by comparing the experimental and simulation results of cold rolling for an S45C plate with a body-centered cubic lattice. Further, the numerical analysis method proposed in this study can contribute to the study of anisotropy as an alternative to experimental approaches.


Introduction
Cold rolling is a widely used process in the manufacturing industry for metal plate production [1]. In this process, the thickness reduction of the metal plate under the recrystallization temperature generates severe anisotropy [2], which influences the sheet metal formation. Wang et al. [3] studied the tension leveling process using different constitutive models and observed that the plastic anisotropy influences the residual curvature of the leveled metal strip. Wu et al. [4] studied cup drawing and hole expansion processes using the Hill48 yield function under the associated and non-associated flow rules. They concluded that the plastic anisotropy has a significant influence on the sheet metal forming results, such as in cup height and thickness reductions. Friedman and Pan [5] predicted the forming limit curves for sheet metal forming using different yield criteria and concluded that the plastic anisotropy can influence the predicted results. Because the anisotropy generated from the cold rolling process significantly affects the subsequent forming processes [6], many researchers have studied the texture evolution of cold-rolled plates using numerical and experimental methods. Bate and Fonseca [7] predicted the texture development of steels during cold rolling using crystal plasticity models. Yanagimoto [8] proposed a transition probability model for active slip systems for body-centered cubic (BCC) metals to predict the texture evolution of steels after 50% cold rolling reduction. Additionally, Morimoto et al. [9] proposed a deformation texture prediction model comprising the Orowan and Taylor models to predict the rolling texture of steels. The texture predictions of the abovementioned studies were all compared with experimental results and presented acceptable prediction accuracy.
However, there is still room to improve the accuracy and efficiency of texture prediction. For instance, Das [10] adopted the Bayesian neural network to efficiently calculate the texture evolution of steels during cold rolling. Fujita et al. [11] adopted a hybrid crystal plasticity approach combining the finite element (FE) method and fast Fourier transform method to study the through-process texture evolution during plate rolling. Taking into account the development of intragranular misorientations, Després et al. [12] adopted the grain fragmentation visco-plastic self-consistent model to calculate the cold rolling texture of steels. The calculation accuracy was improved compared with the results when using the visco-plastic self-consistent model. Our study aimed to establish an efficient elastic-plastic analysis method to predict the metal plate texture after cold rolling using the FE method with zooming analysis. The FE method can calculate the stress and strain distributions of cold-rolled plates accurately, while the zooming analysis can calculate the cold rolling texture efficiently. The zooming analysis can be widely applied in various types of FE analysis to improve the calculation efficiency [13][14][15][16].
In this study, the stress and strain distributions of the cold-rolled plate were calculated at the macro scale first, then at the micro scale to simulate the cold-rolled texture using the previously obtained strain distribution as an input. In this way, the texture calculation was only performed in the magnified area, thereby improving the calculation efficiency. The predicted rolling texture can reflect the anisotropy of a cold-rolled plate and clarify the generation of anisotropy caused by the cold rolling process. Moreover, if the predicted rolling texture is sufficiently accurate, this method can be used in the study of metal plate anisotropy as an alternative to experimental methods such as electron backscatter diffraction (EBSD) and X-ray diffraction.

FE Modeling of Cold Rolling Process
In this section, the FE modeling of the rolls and cold-rolled plate, modeling of the normal and tangent contact behaviors between the rolls and cold-rolled plate, and corotational formulation for the FE analysis are introduced in detail.

Modeling of Rolls and Cold-Rolled Plate
For the FE modeling of the cold rolling process, the rigid analysis body was used to model the rolls and the first-order brick element with eight nodes (details in Appendix A) was adopted to model the cold-rolled plate. Figure 1 shows a schematic of the 1 /4 FE model for the cold rolling process considering the symmetry of the rolls and cold-rolled plate. In the model, x is the rolling direction (RD), y is the transverse direction (TD), and z is the thickness direction. The 1 /2 surface of the upper roll is given by: where R r is the roll radius, y r is the half-width of the roll, and the line (x = x r , z = z r ) represents the upper roll axis.  [12] adopted the grain fragmentation visco-plastic self-consistent model to calculate the cold rolling texture of steels. The calculation accuracy was improved compared with the results when using the visco-plastic self-consistent model. Our study aimed to establish an efficient elastic-plastic analysis method to predict the metal plate texture after cold rolling using the FE method with zooming analysis. The FE method can calculate the stress and strain distributions of cold-rolled plates accurately, while the zooming analysis can calculate the cold rolling texture efficiently. The zooming analysis can be widely applied in various types of FE analysis to improve the calculation efficiency [13][14][15][16].
In this study, the stress and strain distributions of the cold-rolled plate were calculated at the macro scale first, then at the micro scale to simulate the cold-rolled texture using the previously obtained strain distribution as an input. In this way, the texture calculation was only performed in the magnified area, thereby improving the calculation efficiency. The predicted rolling texture can reflect the anisotropy of a coldrolled plate and clarify the generation of anisotropy caused by the cold rolling process. Moreover, if the predicted rolling texture is sufficiently accurate, this method can be used in the study of metal plate anisotropy as an alternative to experimental methods such as electron backscatter diffraction (EBSD) and X-ray diffraction.

FE Modeling of Cold Rolling Process
In this section, the FE modeling of the rolls and cold-rolled plate, modeling of the normal and tangent contact behaviors between the rolls and cold-rolled plate, and corotational formulation for the FE analysis are introduced in detail.

Modeling of Rolls and Cold-Rolled Plate
For the FE modeling of the cold rolling process, the rigid analysis body was used to model the rolls and the first-order brick element with eight nodes (details in Appendix A) was adopted to model the cold-rolled plate. Figure 1 shows a schematic of the ¼ FE model for the cold rolling process considering the symmetry of the rolls and cold-rolled plate. In the model, is the rolling direction (RD), is the transverse direction (TD), and is the thickness direction. The ½ surface of the upper roll is given by: where is the roll radius, is the half-width of the roll, and the line ( = , = ) represents the upper roll axis.   Figure 2 shows the schematics of the global node numbering and element numbering for the cold-rolled plate, where N x , N y , and N z are the numbers of elements in the x, y, and z directions, respectively. Nodes and elements are numbered in the following order of direction: y → z → x . The global node number for node (i nod , j nod , k nod ) and the element number for the element (i ele , j ele , k ele ) are given by: where i nod , j nod , k nod and i ele , j ele , k ele represent the respective number positions of the node and element in the x, y, and z directions. To ensure an accurate FE analysis in the bite area, the mesh size in the rolling direction is determined by the thickness reduction of the cold-rolled plate and the length of the contact region between the roll and the plate. Additionally, the mesh sizes in the transverse and thickness directions were set to values in the same order of magnitude as the mesh size in the rolling direction to ensure calculation convergence.  Figure 2 shows the schematics of the global node numbering and element numbering for the cold-rolled plate, where , , and are the numbers of elements in the , , and directions, respectively. Nodes and elements are numbered in the following order of direction: → → . The global node number for node ( , , ) and the element number for the element ( , , ) are given by: where , , and , , represent the respective number positions of the node and element in the , , and directions. To ensure an accurate FE analysis in the bite area, the mesh size in the rolling direction is determined by the thickness reduction of the cold-rolled plate and the length of the contact region between the roll and the plate. Additionally, the mesh sizes in the transverse and thickness directions were set to values in the same order of magnitude as the mesh size in the rolling direction to ensure calculation convergence.

Modeling of Contact
For the contact modeling, the surface of the roll was treated as the master surface and the surface of the cold-rolled plate was treated as the slave surface. A schematic of the contact is shown in Figure 3, where , , and are the tangent, transverse, and normal directions to the roll surface of the contact point in the local coordinate system, respectively. The normal and tangent contact behaviors are calculated in the local coordinate system and then converted to the global coordinate system when solving the equilibrium equations of the FE analysis, which will be mentioned in Section 4 (details in Appendix B).

Modeling of Contact
For the contact modeling, the surface of the roll was treated as the master surface and the surface of the cold-rolled plate was treated as the slave surface. A schematic of the contact is shown in Figure 3, where x , y , and z are the tangent, transverse, and normal directions to the roll surface of the contact point in the local coordinate system, respectively. The normal and tangent contact behaviors are calculated in the local coordinate system and then converted to the global coordinate system when solving the equilibrium equations of the FE analysis, which will be mentioned in Section 4 (details in Appendix B).  For normal contact behaviors, because the analytical rigid surface is adopted to model the roll surface as the master surface, the direct constraint elimination method can be applied by assuming a penetration of 0 to decide the location of the contact point of the cold-rolled plate in the direction. Subsequently, the rolling pressure can be calculated directly from the stress field of the contact element owing to the balance of internal and external forces in the normal direction [17].
For the tangent contact behaviors, we adopted the inverse tangent rule to calculate the friction between the roll and the cold-rolled plate in the and directions. The inverse tangent rule is given by: where is the friction between the roll and cold-rolled plate, with components and in the and directions, respectively. Further, is the relative displacement of the contact point to the roll surface in the contact plane, with and being the respective relative displacements in the and directions, respectively. The factor is the rolling pressure normal to the roll surface, is the friction coefficient, and is the parameter with an arbitrary value for numerical calculations. The recommended value range for the ∆ ⁄ ratio is 100-1000 [18].

Co-Rotational Formulation
In this study, the co-rotational formulation was applied to define the strain in the elastic-plastic FE analysis of cold rolling to ensure the objectivity of the strain tensor [19]. Figure 4 shows the three configurations in the co-rotational formulation, which are defined according to the right polar decomposition of the deformation gradient tensor , given by: where and are the rotation and right stretch tensor, respectively. As shown in Figure 4, the configuration before deformation is set as the reference configuration and is the material differential fiber before deformation. Meanwhile, the current configuration is that after deformation, becomes = considering the rotation and stretching in the current configuration. Thus, the deformation gradient tensor between the current configuration and reference configuration is defined as: For normal contact behaviors, because the analytical rigid surface is adopted to model the roll surface as the master surface, the direct constraint elimination method can be applied by assuming a penetration of 0 to decide the location of the contact point of the cold-rolled plate in the z direction. Subsequently, the rolling pressure f n can be calculated directly from the stress field of the contact element owing to the balance of internal and external forces in the normal direction [17].
For the tangent contact behaviors, we adopted the inverse tangent rule to calculate the friction between the roll and the cold-rolled plate in the x and y directions. The inverse tangent rule is given by: where f a is the friction between the roll and cold-rolled plate, with components f x and f y in the x and y directions, respectively. Further, u a is the relative displacement of the contact point to the roll surface in the contact plane, with u x and u y being the respective relative displacements in the x and y directions, respectively. The factor f n is the rolling pressure normal to the roll surface, µ is the friction coefficient, and A is the parameter with an arbitrary value for numerical calculations. The recommended value range for the ∆u a /A ratio is 100-1000 [18].

Co-Rotational Formulation
In this study, the co-rotational formulation was applied to define the strain in the elastic-plastic FE analysis of cold rolling to ensure the objectivity of the strain tensor [19]. Figure 4 shows the three configurations in the co-rotational formulation, which are defined according to the right polar decomposition of the deformation gradient tensor F, given by: where R and U are the rotation and right stretch tensor, respectively. As shown in Figure 4, the configuration before deformation is set as the reference configuration and dX is the material differential fiber before deformation. Meanwhile, the current configuration is that after deformation, dX becomes dx = FdX considering the rotation and stretching in the current configuration. Thus, the deformation gradient tensor between the current configuration and reference configuration is defined as:  The co-rotational configuration only considers the stretching, while the material differential fiber becomes = in this configuration. As expressed in the following equations, the strain increment tensor in the co-rotational configuration is defined by using the right stretch tensor , while the strain increment tensor in the current configuration is given by considering the rigid body rotation.
The strain increment tensor can then be calculated using the Green-Lagrange strain tensor , while in the current configuration can be calculated using the Green-Lagrange strain tensor and the rotation tensor , where and are given as: The elastic-plastic decomposition of the strain increment tensor in the co-rotational formulation is given by: where and are the elastic and plastic strain increment tensors observed in the co-rotational configuration, respectively. Meanwhile, the corresponding strain increment tensors observed in the current configuration are denoted without carets. In the next section, material constitutive modeling at both the macro and micro scales is introduced under the co-rotational and current configurations of the co-rotational formulation.

Macro-Scale and Micro-Scale Material Constitutive Modeling
To predict the cold rolling texture of the metal plate, a macro-scale phenomenological constitutive model and a micro-scale crystal plasticity constitutive model were incorporated into the current elastic-plastic analysis framework, as shown in Figure 5. The macro-scale model was applied to the elastic-plastic FE analysis of cold rolling to predict the stress and strain distributions of the cold-rolled plates, while the micro-scale model was coupled with the macro-scale model to predict the local cold rolling texture using zooming analysis. The details of the current analysis framework are described in this section. The co-rotational configuration only considers the stretching, while the material differential fiber dX becomes dx = UdX in this configuration. As expressed in the following equations, the strain increment tensor dε in the co-rotational configuration is defined by using the right stretch tensor U, while the strain increment tensor dε in the current configuration is given by considering the rigid body rotation.
The strain increment tensor dε can then be calculated using the Green-Lagrange strain tensor E, while dε in the current configuration can be calculated using the Green-Lagrange strain tensor E and the rotation tensor R, where E and R are given as: The elastic-plastic decomposition of the strain increment tensor in the co-rotational formulation is given by: dε = dε e + dε p co-rotational configuration dε = dε e + dε p current configuration (8) where dε e and dε p are the elastic and plastic strain increment tensors observed in the co-rotational configuration, respectively. Meanwhile, the corresponding strain increment tensors observed in the current configuration are denoted without carets. In the next section, material constitutive modeling at both the macro and micro scales is introduced under the co-rotational and current configurations of the co-rotational formulation.

Macro-Scale and Micro-Scale Material Constitutive Modeling
To predict the cold rolling texture of the metal plate, a macro-scale phenomenological constitutive model and a micro-scale crystal plasticity constitutive model were incorporated into the current elastic-plastic analysis framework, as shown in Figure 5. The macro-scale model was applied to the elastic-plastic FE analysis of cold rolling to predict the stress and strain distributions of the cold-rolled plates, while the micro-scale model was coupled with the macro-scale model to predict the local cold rolling texture using zooming analysis. The details of the current analysis framework are described in this section. Materials 2021, 14, x FOR PEER REVIEW 6 of 17 Figure 5. Schematic of the current analysis framework.

Macro-Scale Phenomenological Constitutive Model for Stress and Strain Calculation
We defined the macro-scale phenomenological constitutive model of a cold-rolled plate in the co-rotational configuration. For the elastic domain, this model is given by: where = ( ) (1 − 2 )(1 + ) ⁄ is the Lamé constant and = 2(1 + ) ⁄ is the shear modulus of elasticity. In both terms, is the Young's modulus and is Poisson's ratio.
For the plastic domain, the macro-scale phenomenological constitutive model comprises the von Mises yield criterion (Equation (10)), associated flow rule (Equation (11)), and Kim-Tuan isotropic hardening rule (Equation (12)) for an accurate description of the flow stress for large deformation [20]. The reason for adopting these three constitutive relations is because cold rolling is a relatively simple manufacturing process compared with some sheet metal forming processes; that is, for the work-hardening rule, there is no loading reversal during the cold rolling process. Therefore, in this study, it is not necessary to adopt some of the more complicated work-hardening rules (i.e., there is no need to consider the effect of material properties, such as the Bauschinger effect [21]). The isotropic hardening rule together with the adopted yield criterion and flow rule are sufficient to describe the macro-scale material responses during cold rolling.
where is the Cauchy stress tensor in the co-rotational configuration, is the plastic multiplier, and ̅ is the equivalent plastic strain. The parameters , , , , and ℎ of the Kim-Tuan isotropic hardening rule can be identified from the strain-stress curve obtained from the compression test in the thickness direction. It should be noted that this macro-scale constitutive model is defined in the co-rotational configuration of the formulation mentioned in Section 2.3; thus, to obtain the stress and strain distributions in

Macro-Scale Phenomenological Constitutive Model for Stress and Strain Calculation
We defined the macro-scale phenomenological constitutive model of a cold-rolled plate in the co-rotational configuration. For the elastic domain, this model is given by: where λ = (Eν)/[(1 − 2ν)(1 + ν)] is the Lamé constant and G = E/[2(1 + v)] is the shear modulus of elasticity. In both terms, E is the Young's modulus and ν is Poisson's ratio. For the plastic domain, the macro-scale phenomenological constitutive model comprises the von Mises yield criterion (Equation (10)), associated flow rule (Equation (11)), and Kim-Tuan isotropic hardening rule (Equation (12)) for an accurate description of the flow stress for large deformation [20]. The reason for adopting these three constitutive relations is because cold rolling is a relatively simple manufacturing process compared with some sheet metal forming processes; that is, for the work-hardening rule, there is no loading reversal during the cold rolling process. Therefore, in this study, it is not necessary to adopt some of the more complicated work-hardening rules (i.e., there is no need to consider the effect of material properties, such as the Bauschinger effect [21]). The isotropic hardening rule together with the adopted yield criterion and flow rule are sufficient to describe the macro-scale material responses during cold rolling.
whereσ is the Cauchy stress tensor in the co-rotational configuration, dλ is the plastic multiplier, and ε p is the equivalent plastic strain. The parameters σ y , K, t, a, and h of the Kim-Tuan isotropic hardening rule can be identified from the strain-stress curve obtained from the compression test in the thickness direction. It should be noted that this macro-scale constitutive model is defined in the co-rotational configuration of the formulation mentioned in Section 2.3; thus, to obtain the stress and strain distributions in the current configuration, the rigid body rotation must be considered. The macro-scale constitutive model is applied to the FE analysis of cold rolling to predict the stress and strain distributions for the cold-rolled plate (discussed in Section 4).

Micro-Scale Crystal Plasticity Constitutive Model for Lattice Orientation Calculation
The Taylor model [22] was adopted as the micro-scale crystal plasticity constitutive model to predict the cold rolling texture in the current configuration of the co-rotational formulation. This model is widely applied to calculate the texture evolution during metal forming processes for BCC [9], face-centered cubic [23], and hexagonal close-packed [24] metals owing to its simplicity and accuracy. To predict the texture evolution of the coldrolled plate using zooming analysis, the Taylor model, coupled with the macro-scale phenomenological constitutive model mentioned in Section 3.1, is given by: where τ (k) is the shear stress of the slip system k, dε p is the plastic strain increment tensor, and σ is the Cauchy stress tensor in the current configuration of the co-rotational formulation. Moreover, P (k) is the symmetric part of the Schmid tensor of the slip system k, dγ (k) is the slip increment of the slip system, e i is the orthonormal basis vector, a (k) is the unit normal vector of the slip plane, and b (k) is the unit slip direction vector. The yield condition for the slip system is given by: where τ s (k) is the critical resolved shear stress (CRSS) of the slip system. When τ (k) = τ s (k) , the slip of slip system begins. Considering Equations (13)-(15), the grain plastic work increment is given by: Based on the minimum slip principle [22], the grain plastic work increment dW p is minimized using the interior point legacy algorithm [25] under the constraints of Equations (14) and (15) to determine the slip for all slip systems. The lattice orientations can then be calculated using the slip values dγ (k) for all slip systems to predict the local cold rolling texture via zooming analysis (details in Section 4).

FE and Zooming Analysis Method for Cold Rolling Texture Prediction
In this study, the cold rolling process was analyzed as a static problem; thus, the static implicit method was applied to the FE analysis of cold rolling using the co-rotational formulation mentioned in Section 2.3. The equilibrium equation, based on the principle of virtual work using the co-rotational formulation, is given by: where Π is the first Piola-Kirchhoff stress tensor, V is the volume of the cold-rolled plate, t is the surface traction imposed on the cold-rolled plate, δu is the virtual displacement vector, and S t is the surface region where the surface traction is imposed. Here, X represents the coordinates of the reference configuration. Using the FE discretization for the cold-rolled plate mentioned in Section 2.1, the displacement vector u in an element is where N is the shape function of the element (details in Appendix A) and u (nod) is the nodal displacement vector of a node in the element. Thus, the FE form of the equilibrium equation for an element is given by: where V e and S e t represent the volume and surface subjected to the surface traction of the element, respectively. Dropping the nodal virtual displacement δu (nod) on both sides of Equation (20) gives: where T (nod) is the equivalent nodal force imposed on the node. In the FE analysis of cold rolling, T (nod) comprises the rolling pressure and friction mentioned in Section 2.2. Since Equation (21) is a nonlinear equation for the displacement vector u, the Newton-Raphson method was applied to solve for the u as an implicit method. By applying the Newton-Raphson method, we obtained the linearized form of Equation (21), which is called the element tangent stiffness equation, given by: where: The element tangent stiffness equation is used to assemble the corresponding global tangent stiffness equations [K]∆u = {T} − {Q} according to the global node and element numbering given by Equations (19) and (20). The equation is then used to solve ∆u for all nodes in the m th iteration of the Newton-Raphson method (Equation (24)).
The global tangent stiffness equations are stored using the skyline method and the modified Cholesky method is adopted to solve the equations because of its compatibility with the skyline method [26]. When solving Equation (21) using the Newton-Raphson method, after obtaining ∆u (m) in the m th iteration by solving Equation (24), the displacement u (m) in the m th iteration is given by: where β is a relaxation coefficient that promotes convergence of the calculation. For the convergence judgment of the Newton-Raphson method, the following residuals were used: (26) Figure 6 shows a flowchart of the FE analysis for cold rolling texture prediction using zooming analysis. For the calculation of the stress and strain distributions during the assembly of the global tangent stiffness equation, we adopt the return mapping algorithm as the stress integration scheme to calculate the stress and strain in the co-rotational configuration using the macro-scale constitutive model mentioned in Section 3.1. For each element, the strain is calculated using the selective reduced integration with eight Gaussian integration points and an integration point located at the element center to avoid volumetric locking [27]. For the calculation of the cold rolling texture, the local plastic strain increment obtained by the FE analysis, which is represented by the nodal plastic strain increment in the corresponding area, is selected to be substituted into Equation (14) for the zooming analysis. The texture calculation, as shown in Figure 6, can then be conducted with the micro-scale constitutive model mentioned in Section 3.2 using the nodal plastic strain increment as an input. The selection of the slip systems and the corresponding CRSS ratios for the texture calculation depends on the cold rolled materials, which is discussed in detail in Section 5.
Materials 2021, 14, x FOR PEER REVIEW 9 of 17 Figure 6 shows a flowchart of the FE analysis for cold rolling texture prediction using zooming analysis. For the calculation of the stress and strain distributions during the assembly of the global tangent stiffness equation, we adopt the return mapping algorithm as the stress integration scheme to calculate the stress and strain in the co-rotational configuration using the macro-scale constitutive model mentioned in Section 3.1. For each element, the strain is calculated using the selective reduced integration with eight Gaussian integration points and an integration point located at the element center to avoid volumetric locking [27]. For the calculation of the cold rolling texture, the local plastic strain increment obtained by the FE analysis, which is represented by the nodal plastic strain increment in the corresponding area, is selected to be substituted into Equation (14) for the zooming analysis. The texture calculation, as shown in Figure 6, can then be conducted with the micro-scale constitutive model mentioned in Section 3.2 using the nodal plastic strain increment as an input. The selection of the slip systems and the corresponding CRSS ratios for the texture calculation depends on the cold rolled materials, which is discussed in detail in Section 5.

Validation of Current Analysis Framework
To validate the proposed elastic-plastic analysis framework under the co-rotational formulation, we conducted a cold rolling experiment with a 28.33% thickness reduction for an S45C plate (Nippon Steel, Tokyo, Japan) with an 84 mm width and 6 mm thickness. The FE analysis and texture prediction for the cold rolling process were then conducted for comparison using the analysis method mentioned in Section 4.
The FE modeling parameters were matched with the experimental cold rolling conditions, while the detailed FE analysis conditions for the S45C plate are listed in Table  1. The ¼ plate was discretized using 3840 elements (i.e., 160, 8, and 3 grids along the rolling, width, and thickness directions were generated, respectively). It is worth mentioning that a case study was conducted to finally determine the appropriate mesh

Validation of Current Analysis Framework
To validate the proposed elastic-plastic analysis framework under the co-rotational formulation, we conducted a cold rolling experiment with a 28.33% thickness reduction for an S45C plate (Nippon Steel, Tokyo, Japan) with an 84 mm width and 6 mm thickness. The FE analysis and texture prediction for the cold rolling process were then conducted for comparison using the analysis method mentioned in Section 4.
The FE modeling parameters were matched with the experimental cold rolling conditions, while the detailed FE analysis conditions for the S45C plate are listed in Table 1. The 1 /4 plate was discretized using 3840 elements (i.e., 160, 8, and 3 grids along the rolling, width, and thickness directions were generated, respectively). It is worth mentioning that a case study was conducted to finally determine the appropriate mesh size under the premise of guaranteeing both the calculation efficiency and accuracy. Before cold rolling, the flow stress curve of the S45C plate was obtained first by conducting compression tests in the thickness direction using cylindrical S45C plate samples with a diameter of 4 mm and a height of 6 mm. Two compression tests were conducted and the obtained flow stress curves were identical (experimental results are shown in Figure 7). For the elastic deformation, the Young's modulus and Poisson's ratio were set to 210 GPa and 0.3, respectively. For the plastic deformation, the Kim-Tuan isotropic work-hardening rule (Equation (12)) was adopted to fit the experimental flow stress curve, and the results are shown in Figure 7. Moreover, the friction coefficient of the inverse tangent rule (Equation (3)) is assumed to be 0.3. samples with a diameter of 4 mm and a height of 6 mm. Two compression tests were conducted and the obtained flow stress curves were identical (experimental results are shown in Figure 7). For the elastic deformation, the Young's modulus and Poisson's ratio were set to 210 GPa and 0.3, respectively. For the plastic deformation, the Kim-Tuan isotropic work-hardening rule (Equation (12)) was adopted to fit the experimental flow stress curve, and the results are shown in Figure 7. Moreover, the friction coefficient of the inverse tangent rule (Equation (3)) is assumed to be 0.3.  The distributions of the rolling pressure and rolling friction in the contact area calculated by the FE analysis are shown in Figures 8 and 9, respectively. The neutral line of the friction in the direction was successfully predicted by the inverse tangent rule. Moreover, the position of the peak rolling pressure is near the neutral line, which is consistent with the results obtained by Richelsen and Tvergaard [28]. The experimental and simulation results of the width spread ratio, rolling force, and rolling torque are shown in Table 2, highlighting the accuracy of the proposed elastic-plastic FE analysis method. The distributions of the rolling pressure and rolling friction in the contact area calculated by the FE analysis are shown in Figures 8 and 9, respectively. The neutral line of the friction in the x direction was successfully predicted by the inverse tangent rule. Moreover, the position of the peak rolling pressure is near the neutral line, which is consistent with the results obtained by Richelsen and Tvergaard [28]. The experimental and simulation results of the width spread ratio, rolling force, and rolling torque are shown in Table 2, highlighting the accuracy of the proposed elastic-plastic FE analysis method.   For the texture prediction, because of the BCC lattice structure of the cold-rolled S45C plate, 48 slip systems (listed in Table 3) were considered for the Taylor model to calculate the cold rolling texture. Because the {110}<111> and {112}<111> slip systems contribute   [29].
A node at the central layer of the center part of the cold-rolled S45C plate was selected for the validation of the texture prediction. The texture of the corresponding area was obtained by EBSD measurements for comparison with the prediction results. At the central layer of the cold-rolled plate, the lattice rotation tensor dω (k) caused by the slip of slip system k is given by: where Q (k) is the unsymmetric part of the Schmid tensor of slip system k and is expressed as: The total lattice rotation tensor of the calculated node in the central layer is the sum of the lattice rotations of all slip systems. Considering all the lattices contained by the calculated node, the lattice orientations can be calculated using Equations (27) and (28) to predict the cold rolling texture. Figure 10a shows the EBSD results. The measurements obtained 4260 lattice orientations in a scan area of 2 mm 2 in the central layer of the center part of the cold-rolled plate. Comparing this with the texture prediction, we calculated the same number of random original lattice orientations using the Taylor model and the nodal plastic strain increment obtained using FE analysis. The predicted texture is shown in Figure 10b was also validated by Yanagimoto [8] using the transition probability model for active slip systems for BCC metals. It is worth mentioning that the pole figure with normalized density can be used for the prediction of directional R-values [30] and yield stresses [31], which are important manifestations of anisotropy; thus, the proposed analysis framework can contribute to the studies of the anisotropy generation and prediction related to the cold rolling process. It is also a potential alternative to experimental measurements of cold rolling texture.

Conclusions
We proposed an elastic-plastic analysis framework, which incorporated a macroscale FE analysis, using a phenomenological constitutive model for stress and strain calculations and a micro-scale zooming analysis using a crystal plasticity constitutive model for lattice orientation calculation. The proposed framework was applied to the elastic-plastic FE analysis of cold rolling using the static implicit method under a corotational formulation to predict the local cold rolling texture.
The cold rolling texture between the predicted results and EBSD measurements on the S45C plate was in good agreement, thereby validating the proposed elastic-plastic analysis framework. Therefore, the proposed method qualifies as an effective numerical method for future studies on anisotropy generation and prediction related to the cold rolling process. It is also a potential alternative to experimental measurements of cold rolling texture. Further improvement of the current analysis framework may include adding feedback from the micro scale to the macro scale (e.g., adopting an anisotropic constitutive model at the macro scale and calculating the anisotropic parameters using the predicted texture at the micro scale as the feedback) for multi-scale analysis, which will

Conclusions
We proposed an elastic-plastic analysis framework, which incorporated a macro-scale FE analysis, using a phenomenological constitutive model for stress and strain calculations and a micro-scale zooming analysis using a crystal plasticity constitutive model for lattice orientation calculation. The proposed framework was applied to the elastic-plastic FE analysis of cold rolling using the static implicit method under a co-rotational formulation to predict the local cold rolling texture.
The cold rolling texture between the predicted results and EBSD measurements on the S45C plate was in good agreement, thereby validating the proposed elastic-plastic analysis framework. Therefore, the proposed method qualifies as an effective numerical method for future studies on anisotropy generation and prediction related to the cold rolling process. It is also a potential alternative to experimental measurements of cold rolling texture. Further improvement of the current analysis framework may include adding feedback from the micro scale to the macro scale (e.g., adopting an anisotropic constitutive model at the macro scale and calculating the anisotropic parameters using the predicted texture at the micro scale as the feedback) for multi-scale analysis, which will be the focus of our future studies.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

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

Appendix A. Shape Function of the First-Order Brick Element
The first-order brick element with eight nodes (local node numbering is shown in Figure A1) was adopted to model the cold-rolled plate. The shape function of the first-order brick element in the natural coordinate system is given by: where g i , h i , and r i are the coordinates of node i in the natural coordinate system. Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

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

Appendix A. Shape Function of the First-Order Brick Element
The first-order brick element with eight nodes (local node numbering is shown in Figure A1) was adopted to model the cold-rolled plate. The shape function of the firstorder brick element in the natural coordinate system is given by: where , ℎ , and are the coordinates of node i in the natural coordinate system. Figure A1. Local node numbering of the first-order brick element in the natural coordinate system.

Appendix B. Local Coordinate System for Contact Calculation
As shown in Figure A2, when the angle between the axis of the local coordinate system and the axis of the global coordinate system is , the transformation from the local to the global coordinate system is given by: The contact behaviors are calculated in the local coordinate system and then converted to the global coordinate system using Equation (A2) when solving the equilibrium equations mentioned in Section 4. Figure A1. Local node numbering of the first-order brick element in the natural coordinate system.

Appendix B. Local Coordinate System for Contact Calculation
As shown in Figure A2, when the angle between the z axis of the local coordinate system and the z axis of the global coordinate system is θ, the transformation from the local to the global coordinate system is given by: The contact behaviors are calculated in the local coordinate system and then converted to the global coordinate system using Equation (A2) when solving the equilibrium equations mentioned in Section 4. Figure A2. Schematic of local coordinate system for contact calculation.