Effects of yield point and plastic anisotropy on results of elastic–plastic finite element analysis of tension leveling

Tension leveling is an important industrial process to eliminate the flatness defects and residual stresses of metal strips to provide high-quality sheet metals for subsequent sheet metal forming. The finite element (FE) method can be applied to elucidate the effects of process parameters on the quality of sheets after tension leveling for various materials. In our previous investigation, an accurate FE model has been established for the elastic–plastic FE analysis of tension leveling. In this study, we further studied the effects of the yield point and plastic anisotropy on tension leveling using the FE model established in our previous investigation. Aiming at improving the accuracy of simulation, a modified constitutive model was developed to describe the anisotropic hardening of materials under cyclic loading. The modified constitutive model was implemented into Abaqus/Standard as a user-defined material subroutine to simulate the development of the anisotropy in materials during tension leveling. The modified model was also applied to the FE analysis of sheet metal forming processes to demonstrate its simulation capability and accuracy.


Introduction
Tension leveling is widely applied in industry to eliminate the flatness defects and residual stresses of metal strips [1,2]. The finite element (FE) method can be applied to elucidate the effects of process parameters on the quality of sheets after tension leveling for various materials [3]. The tension leveling experiment involving different process parameters for different materials can be conducted efficiently using the FE method, and the quality of leveled sheets under different leveling conditions can be checked by the numerical results. Subsequently, the effects of process parameters on the quality of different leveled sheets can be clarified by analyzing the numerical results with different process parameters. Recently, with aims to clarify the factors affecting tension leveling and achieve simulations of tension leveling with high accuracy, many researchers conducted related studies. Zhao et al. [4] conducted the FE analysis of tension leveling to determine the relationship between the intermesh of the leveler rolls and the working curvatures of leveled strips under a specific leveling arrangement. Zhang et al. [5] conducted tension leveling tests to investigate the effects of tension leveling elongation and intermesh of leveler rolls on tension leveling. They found that both factors can influence the microstructure and mechanical properties of leveled steel strips, indicating the possibility of improving steel strip properties by adjusting the elongation and intermesh during tension leveling. Garcia et al. [6] applied different work hardening models to the FE analysis of tension leveling and concluded that the introduced hardening rule is an important input in the FE analysis of tension leveling since it influences simulation results. Although previous investigations have clarified some factors affecting tension leveling such as the arrangement and intermesh of leveler rolls and the work hardening properties of materials, further investigations must be conducted (1) to reveal much more factors that affect tension leveling to provide a deeper insight into this process, (2) to improve the accuracy of tension leveling simulation, and (3) to realize consistency analysis of tension leveling and subsequent sheet metal forming considering the tension leveling history.
Regarding much more factors that affect tension leveling, the yield point phenomenon and plastic anisotropy are worth studying because they are remarkable characteristics of some metal materials, but their effects on tension leveling have not been well studied. The yield point phenomenon is a common property of metal materials such as mild steels [7], aluminum alloys [8], and titanium alloys [9]. Metal materials with yield point phenomenon exhibit a characteristic stress drop in the stress-strain curve at the end of the elastic region [10]. Regarding the true stress-strain behavior, it is characterized by a yield drop and the subsequent strain hardening [7], as shown in the enlarged area of Fig. 1. The effects of the yield point on processes involving cyclic loading have been investigated by many researchers. Žerovnik et al. [11] conducted uniaxial cyclic loading tests to study the effects of the yield point on the stress-strain responses and strain fields of cyclically loaded specimens. They also proposed a constitutive model for the yield point phenomenon in cyclic plasticity and applied it to the simulation of a cyclically loaded console beam to reveal the effects of the yield point phenomenon on the accumulated plastic strain after several loading cycles [12]. Kim et al. [13] conducted roller leveling to study the reduction of yield point during leveling. Plastic anisotropy is a significant factor affecting the forming of sheet metals [14]. Recently, several studies have been conducted to determine the effects of plastic anisotropy on sheet metal forming. Wu et al. [15] combined the yld91 model with the Hill48 model under the non-associated flow rule with the anisotropic hardening to elucidate the influence of the evolution of plastic anisotropy during the hole expansion of JSH590R sheets. Kim et al. [16] developed several constitutive models based on the non-associated flow rule using the Hill48, Yld2000-2d, and Yoshida polynomial functions to simulate hole expansion tests. They found that the anisotropy of the plastic potential function plays an important role in the thinning behavior prediction. Hu et al. [17] coupled the poly4 and Hosford yield functions under the associated flow rule to evaluate the anisotropic yielding behavior of aluminum alloy sheets during cup drawing to achieve an accurate FE simulation. Materials are also subjected to cyclic loading and demonstrate plastic anisotropy during tension leveling, but their effects on this process have not been well studied.
In our previous investigation [18], an FE model considering the anisotropy and cyclic plasticity of materials using the Hill48 yield criterion with the non-associated flow rule and the mixed hardening rule has been established for the elastic-plastic FE analysis of tension leveling. In this study, we further studied the effects of the yield point and plastic anisotropy on tension leveling using the FE model established in our previous investigation. On the basis of the effects of the plastic anisotropy on the results of the FE analysis of tension leveling, the constitutive model in our previous investigation [18] is modified to simulate the development of the anisotropy in materials during tension leveling accurately.

FE implementation of tension leveling
The implementation procedure of this study consists of (1) the FE modeling of tension leveling and (2) the numerical case studies using the FE model of tension leveling to study the effects of yield point and plastic anisotropy. The arrangement of leveler rolls used for tension leveling [19] in this study is shown in Fig. 2. The FE model of tension leveling is based on this arrangement in Abaqus/Standard, where leveler rolls are modeled as analytical rigid bodies and leveled strips are discretized as second-order brick elements. The dimensions of a leveled strip were 1300 mm length (x) × 200 mm width (y) × 0.5 mm thickness (z), and the corresponding mesh size was set to 5 mm (x) × 10 mm (y) with four elements in the thickness direction. For the material properties of the leveled strips, the user-defined material (UMAT) subroutine of Abaqus/Standard was used to implement the material properties (e.g., the Bauschinger effect, plastic anisotropy, and yield point) investigated in this study with the fully implicit backward Euler return mapping algorithm [18]. The static implicit step was applied to simulate this process. For all simulated cases, the Young's modulus and Poisson's ratio of the leveled strips were set to 210 GPa and 0.3, respectively. The friction coefficient between the leveler roll and the leveled strip was set to 0.15 according to the recommended value for cold rolling production line. For each simulated case, the tension was set to 40 MPa and the intermesh of the leveler rolls was set to a constant value. The FE analysis of springback was conducted by removing the leveler rolls after tension leveling to predict the final shape of leveled strips (residual curvatures) after elastic recovery. The details of the numerical case studies to investigate the effects of yield point and plastic anisotropy on tension leveling will be discussed in Sects. 3 and 4, respectively.

Effects of yield point on finite element analysis results of tension leveling
As mentioned in Sect. 1, the yield point phenomenon is a common property of metal materials, and the effects of the yield point on processes involving cyclic loading have been investigated by many researchers [11][12][13]. Materials are also subjected to cyclic loading during tension leveling and the effects of the yield point on this process have not been investigated; thus, we conducted FE analysis using the constitutive model applied in our previous investigation [18] (the Hill48 yield criterion with the non-associated flow rule and the mixed hardening rule, briefly summarized in Appendix 1) to study the effects of the yield point phenomenon on tension leveling. In this section, two types of virtual material with and without a yield point are adopted to study the effects of the yield point phenomenon on tension leveling. The virtual materials are both considered as isotropic using the constitutive model in our previous investigation [18] the stress-strain responses during cyclic loading are shown in Fig. 1. The only difference between these two materials is in the presence or absence of a yield point. The corresponding work hardening parameters of the constitutive model in our previous investigation [18] are shown in Table 1. The FE model mentioned in Sect. 2 was adopted to conduct the FE analysis of tension leveling. The intermeshes were set to  Table 2 (all results are the averages of the measured residual curvatures and total elongations at the center and edge of the leveled strip). Considering that the stress-strain responses are the same for these two materials except the yield point, the difference between these two materials indicates that the yield point phenomenon has an effect on tension leveling. To evaluate the effect in detail, the total elongation of each layer, which is calculated as the average of the center and edge parts of the leveled strip when the leveled strip passes through the leveling unit with an intermesh of 2.5 mm, was measured numerically for the two materials. The elongation history results of the two materials are compared in Fig. 3. It is shown that not only the yield point leads to a smaller elongation at the exit of the leveling unit than in the simulated case without the yield point (as shown in Table 2), but also the elongation histories in the leveling unit for materials with and without a yield point are slightly different (especially for the central layer). Taken together, it can be concluded from the FE analysis results that the yield point has influence on the elongation of the leveled strip during tension leveling.  Table 1 Parameters of work hardening rule with and without yield point yield point  200  -40  3000  50000  1300  900  without yield point  160  10  15  50000 1300 900

Effects of plastic anisotropy on finite element analysis results of tension leveling
As mentioned in Sect. 1, plastic anisotropy is a significant factor affecting the forming of sheet metals, and several recent studies have been conducted to determine the effects of plastic anisotropy on sheet metal forming [15][16][17]. To clarify the effects of plastic anisotropy on tension leveling, we conducted FE analysis using the Hill48 yield criterion with the associated flow rule and the mixed hardening rule, and the tension leveling results of virtual materials with different anisotropies were evaluated in our previous study [20]. The corresponding constitutive model and parameters are shown in Table 3 (the reader is referred to our previous study and Appendix 1 for more details about the corresponding constitutive model, note that the Hill48 yield criterion with the non-associated flow rule becomes the Hill48 yield criterion with the associated flow rule when the anisotropic parameters are the same for both the yield function and the plastic potential function). The FE analysis of tension leveling using the FE model mentioned in Sect. 2 was conducted for those materials with different anisotropic parameters to further study the effects of plastic anisotropy on the FE analysis results of tension leveling. Note that in each simulated case, only one of the anisotropic parameters was adjusted compared to the isotropic situation ( F = G = H = 1.0, L = M = N = 3.0 ) to study the effects of each single anisotropic parameter on tension leveling. In all simulated cases, the intermesh was set to 10 mm.
For the residual curvatures in the x direction after tension leveling (regarded as curl, obtained by conducting springback simulation after tension leveling, shown in Fig. 4), when the anisotropic parameters F, G, and H change, the values of curl also change. For L, M, and N, the effects are not significant since those anisotropic parameters are related to shear stresses. It is concluded that plastic anisotropy affects the strain and stress distributions during the FE analysis of tension leveling, especially for the central layer [20]. Figure 5 shows the plastic elongation in the x direction of the central layer after the tension leveling of the strips with different anisotropic parameters, indicating that the different anisotropic parameters F, G, and H lead to different elongations of the leveled strips, resulting in different curl values. The effects of plastic anisotropy on the residual curvatures in the y direction (regarded as gutter, obtained by conducting springback simulation after tension leveling, shown in Fig. 6) after tension leveling are also evaluated in this investigation. Similar to the results of curl values, when F, G, and H change, the values of the gutter also change, but the trend of the values of the gutter is opposite to that of the curl values when F, G, and H increase. For L, M, and N, no significant effects are observed. Plastic anisotropy can both affect the curl and gutter of the leveled strips after tension leveling; thus, it is an important factor affecting tension leveling.

Application of modified anisotropic constitutive model to finite element analysis of tension leveling
In our previous investigation [18], the Hill48 with the nonassociated flow rule and the mixed hardening rule was adopted as the constitutive model to conduct the FE analysis of tension leveling (details in Appendix 1). However, the Hill48 yield criterion applied in that constitutive model is with constant anisotropic parameters identified from the initial yielding surface (Eqs. 6-7). Thus, the constitutive model in our previous investigation [18] can only capture the anisotropy of initial yield stresses of materials. For steel sheets, besides the anisotropy of initial yield stresses, the evolution of the subsequent hardening also shows anisotropic behavior [21]. It is concluded in Sect. 4.1 that the anisotropic parameters affect the FE analysis results of tension leveling since the FE results vary with the anisotropic parameters. This finding implies the necessity of applying an anisotropic yield function with variable anisotropic parameters in a constitutive model to capture the anisotropic evolution of the work hardening behavior of materials during the FE analysis of tension leveling. Thus, the constitutive model applied in our previous investigation [18] was modified considering the anisotropic hardening under cyclic loading in the FE analysis of tension leveling. The details of the modified constitutive model are shown as follows. Work hardening rule: (1) here f ( ) is the yield function. , and denote the stress tensor considering the back stresses, the Cauchy stress tensor and the back stress tensor, respectively. The subscript m represents six components (xx, yy, zz, xy, xz, yz) of the stress tensor. The subscript n represents four directions of the metal sheet, where 0 represents the rolling direction, 45 represents the diagonal direction of the rolling and transverse directions, 90 represents the transverse direction, and ND represents the thickness direction. g( ) is the plastic potential function. n ( p ) is the flow stress in the n direction considering the back stresses. p is the equivalent plastic strain. C 1,m , C 2,m and m can be calculated using Eqs. (7-9). 0 n , Q n , b n , C 1,n , C 2,n and n are the parameters of the work hardening rule and can be identified from the cyclic stress-strain curve in the n direction. Compared with the parameters of the back stresses in the mixed hardening rule of the constitutive model before modification ( C 1 , C 2 and 1 in Eq. A5), which are only identified from the cyclic stress-strain curve in the rolling direction, the parameters of the back stresses in the mixed hardening rule after modification ( C 1,m , C 2,m and m in Eqs. [4][5] are independent in different directions, and are identified from the cyclic stress-strain curves in four directions; thus, the anisotropic hardening under cyclic loading can be captured.
Yield criterion: the anisotropic parameters of the yield function and can be calculated from flow stresses considering back stresses in the rolling, diagonal, transverse and thickness directions (Eqs. [11][12][13][14]. Compared with the anisotropic parameters in the yield function of the constitutive model before modification (Eq. A6), which are all identified from the directional initial yield stresses, all the anisotropic parameters of the yield function in the constitutive model after modification are identified from the directional flow stresses and become evolutionary during the working hardening. In this case, both the initial yield stress anisotropy and the subsequent anisotropic hardening can be captured. Flow rule: here p is the plastic strain tensor. d is the plastic multiplier.
potential function and can be calculated from R-values in the rolling, diagonal and transverse directions (Eqs. [17][18][19][20]. Note that it is possible to modify this plastic potential function with variable anisotropic parameters to capture the changes in R-values during forming when necessary. In this study, the changes in R-values are not considered. For the work hardening rule of the modified constitutive model (Eqs. 1-9), the mixed hardening rule consisting of the isotropic and kinematic hardening rules was modified considering the cyclic anisotropic hardening behaviors in the rolling, diagonal, transverse, and thickness directions. For the yield criterion, the Hill48 yield function was modified with the variable anisotropic parameters shown in Eqs. (10)- (14) to capture the anisotropic hardening behaviors in the rolling, diagonal, transverse, and thickness directions. For the flow rule shown in Eqs. (15)- (20), the classical Hill48 plastic potential function with constant anisotropic parameters was adopted. For the constitutive parameter identification, flow stresses in the rolling, diagonal, transverse, and thickness directions are required to determine the anisotropic parameters of the yield criterion and the parameters of the work hardening rule. The R-values in the rolling, diagonal and transverse directions are required to calculate the anisotropic parameters of plastic potential functions. The modified constitutive model considering the anisotropic hardening was implemented in Abaqus/Standard via user subroutine UMAT with the fully implicit backward Euler return mapping algorithm, and the successful application of the algorithm for anisotropic hardening was also reported by Wu et al. [22]. The FE analysis of tension leveling was conducted using the FE model mentioned in Sect. 2 with the modified constitutive model. The material of the leveled strip is SPCC and the identified constitutive parameters for SPCC are shown in Table 4. The parameters of flow stresses and the plastic potential function are identified from experimental data reported by Yoshida [19] and Huh [21], and the flow stress in the thickness direction is calculated as the average of the flow stresses in the rolling, diagonal and transverse directions. The tensile and cyclic loading results in the rolling, diagonal, transverse and thickness directions obtained from the modified constitutive model for SPCC are shown in Fig. 7, indicating the capability of the modified model to capture the anisotropic hardening behaviors of SPCC under monotonic and cyclic loading conditions. After conducting the FE analysis of tension leveling using the modified constitutive model, the residual curvatures in the x direction were measured and compared with the simulation results obtained using the constitutive model in our previous study [18] (details in Appendix 1, referred to as previous constitutive model) and the experimental results [19]. The comparisons are shown in Fig. 8, the mean absolute error for all measured residual curvatures decreases from 0.1041 to 0.0489 m −1 after applying the modified constitutive model, indicating that the modified constitutive model can provide more accurate FE analysis results of tension leveling than the previous constitutive model because of the consideration of anisotropic hardening. The modified constitutive model can also be applied to simulate sheet metal forming processes involving plastic anisotropy under cyclic loading (an example of applying the modified constitutive model to the accurate simulation of a two-step drawing process is shown in Appendix 2). The FE analysis of sheet metal forming processes considering the forming history (e.g., the forming strain history and residual stresses from the previous forming process) has a higher accuracy than that without considering the forming history [23]. Thus, the modified constitutive model can also contribute to the improvement of the accuracy of the FE analysis of sheet metal forming processes via a continuous simulation from leveling to forming.

Conclusions
The FE analysis of tension leveling was conducted using an established FE model in Abaqus/ Standard with a static implicit analysis step, and the user subroutine UMAT was adopted to implement constitutive models to the FE analysis of tension leveling with the fully implicit backward Euler return mapping algorithm [18]. The effects of the yield point on the FE analysis results of tension leveling were confirmed by FE analysis using constitutive models with and without a yield point. The effects of plastic anisotropy on the curl and gutter of the leveled strips after tension leveling were also confirmed by the FE analysis of tension leveling with different anisotropic parameters of the Hill48 yield function. Since the anisotropic parameters can affect the results of the FE analysis of both curl and gutter after tension leveling,  the constitutive model in our previous investigation [18] was modified using variable anisotropic parameters of the Hill48 yield function and different parameters for different directions of the mixed hardening rule to capture the anisotropic hardening behavior of materials under monotonic and cyclic loading conditions. The modified constitutive model was applied to the FE analysis of tension leveling, and more accurate simulation results were obtained because of the consideration of anisotropic hardening. The modified model can also contribute to the improvement of the accuracy of the FE analysis of sheet metal forming via a continuous simulation from leveling to forming. The constitutive modeling method and its FE implementation demonstrated in this study can be adopted as an effective tool to investigate the effects of material properties on metal forming processes.

Appendix 1: Hill48 yield criterion with non-associated flow rule and mixed hardening rule
The Hill48 yield criterion (Eq. A1) with the non-associated flow rule (Eqs. A2-A3) and the mixed hardening rule (Eqs. A4-A5) of the constitutive model applied in our previous investigation are briefly summarized as follows: represents the stress tensor considering the back stresses, where = -, and and denote the Cauchy stress tensor and back stress tensor, respectively. p , d , p and Y are the plastic strain tensor, plastic multiplier, equivalent plastic strain and yield stress in the rolling direction of the sheet metal, respectively. 0 , Q, b, C 1 , 1 and C 2 are the parameters of the mixed hardening rule, and can be identified from the cyclic stress-strain curve in the rolling direction of the sheet metal. The anisotropic parameters of the yield and plastic potential functions are identified as:

Appendix 2: Application of modified constitutive model to FE analysis of two-step drawing process
The schematic of the two-step drawing process for a DP600 steel sheet is shown in Fig. 9, where the drawing direction is the z direction and the lower surface of the undeformed sheet is the zero position of the z-axis. This process is part of Benchmark 2 of NUMISHEET 2014 to predict the sprung shape of a rectangular drawing cup after a two-step drawing process [24]. In the first drawing, the tooling radius R of the upper die and punch is 12 mm and the punch stroke is Fig. 9 Schematic of two-step drawing process for DP600 steel sheet Fig. 10 Dimensions of tools [24] 32.4 mm. In the second drawing, the tooling radius R of the upper die and punch is 8 mm. The deformed DP600 steel sheet in the first drawing process is drawn again, and the punch stroke of the second drawing is 47.5 mm. A blank holder force of 595 kN is provided by the lower die in the two-step drawing process. Note that all dimensions of the lower die, upper die, and punch are the same for the first and second drawing processes except for the tooling radius of the upper die and punch (tool dimensions shown in Fig. 10).
The dimensions of the DP600 steel sheet are 300 mm (x) × 250 mm (y) × 1.0 mm (z) with the rolling direction aligned along the 300 mm length [24]. After the two-step drawing process, the rectangular drawing cup is trimmed along the centerline of the y direction, and the sprung shape of the trimming edge in the x − z plane is measured. The FE model of the two-step drawing process was established on the basis of the tool and sheet dimensions, where the tools were modeled as rigid bodies and the DP600 steel sheet was discretized as the first-order brick element. The mesh size of the DP600 steel sheet was 2.5 mm (x) × 2.5 mm (y) with three elements in the thickness direction. The friction coefficient between the tool and the sheet was set to 0.1. The modified constitutive model mentioned in Sect. 4.2 was applied to describe the material behaviors via UMAT during the FE analysis, and the constitutive parameters of DP600 identified from material data reported by Stoughton [25] are shown in Table 5 (note that only the cyclic stress-strain curve in the rolling direction is available [25]; thus, C 1 , C 2 and in the diagonal, transverse and thickness directions are assumed to be equal to the corresponding values in the rolling direction). The static implicit step was adopted to simulate the two-step drawing process, the FE analysis of springback was conducted by removing the punch and dies after the twostep drawing to predict the sprung shapes of the rectangular cup after elastic recovery, and the comparison of the sprung shapes of the simulation and the experiment [26] after the two-step drawing is shown in Fig. 11, showing the capability of the modified constitutive model to simulate sheet metal forming accurately.