E ﬀ ects of Tension–Compression Asymmetry on Bending of Steels

: Stainless steels (SUS) and dual-phase (DP) steels have tension-compression asymmetry (TCA) in mechanical responses to full loading cycles. This phenomenon can signiﬁcantly inﬂuence sheet metal forming of such metals, however, it is di ﬃ cult to describe this behaviour analytically. In this research, a novel analytical method for asymmetric elastic-plastic pure bending using the Cazacu–Barlat 2004 asymmetric yield function is proposed. It only uses material parameters in tension along with an asymmetry coe ﬃ cient related to the yield function. Bending operations of SUS304 and DP980 are investigated as two case studies. In the pure bending for both SUS304 and DP980, moment–curvature diagrams are analytically obtained. Furthermore, linear and nonlinear springback behaviours of SUS304 are analytically investigated. Moreover, using the analytical model as a user-deﬁned material, a numerical model is developed for both steels under pure bending. In the V-bending case of SUS304 with and without TCA e ﬀ ects, the springback behaviours of the material are investigated numerically. In addition, considering friction e ﬀ ects, the analytical method is further modiﬁed for predicting springback behaviours in the V-bending of 16 types of SUS304 with various strengths are determined. All the analytical and numerical results have good agreement with those experimental results from literature for validation. both rolling directions, while the symmetric model shows a discrepancy between measured and calculated results.


Introduction
Over the last decade, due to the development of advanced materials and manufacturing technologies, sheet metal forming has become more challenging for industry sectors and researchers. This manufacturing process has a wide range of applications in automobile parts, electronic parts, food and drink cans, etc. Bending is a widely used forming process in sheet metal forming which is usually accompanied with an unwanted springback phenomenon. Springback is one of the most challenging issues in sheet metal forming. Accurate prediction of springback helps manufacturing industries to design appropriate tools to compensate for this effect and enhance the dimensional accuracy of the finished parts [1].
Anisotropy in plastic behaviour is a common phenomenon in metal forming. An initial anisotropy in yield points and the subsequent hardening behaviour can be caused by tension-compression asymmetry (TCA) [2,3]. Several metals have been reported to have asymmetrical yield surface due to the fact of their eccentricity in yielding or in their subsequent hardening. Hexagonal closed pack (HCP) metals, such as magnesium alloys, have strong asymmetry due to the fact of their limited slip systems and the strong basal texture [4]. Cubic materials, such as body-centred cubic (BCC) and face-centred cubic (FCC) metals, are usually considered as symmetrical materials. However, some cubic metals are reported to exhibit TCA too. Kuwabara et al. [5] designed a special testing apparatus and performed which show great agreement. By using the UMAT, the finite element model of the pure bending was devised, and the obtained numerical results were in good agreement with the analytical and experimental results. In the second case, V-bending, the effects of TCA on the springback behaviour of SUS304 plates were investigated numerically. The springback results from the simulations with and without TCA effects were compared with those experimental results from the literature. The results show that springback prediction improved significantly by considering the TCA of the material. Finally, the analytical method was modified for prediction of springback in V-bending. Based on this method, the springback values of 16 types of SUS304 with different strengths were calculated and compared with those from the experiments. The predicted springback behaviours were consistent with the experimental results. This paper is outlined in the following sections. Section 2 provides a detailed development of the proposed analytical method. Section 3 explains the numerical method and discusses the analytical and numerical results in detail, and, finally, Section 4 gives comprehensive conclusions based on the research findings obtained.

Analytical Method
The model materials used in this study were austenitic stainless-steel SUS304 and DP980 sheets with thicknesses of 0.3 mm and 1.2 mm, respectively, and widths of 1.5 mm and which were studied in the literature [6,7]. An in-plane uniaxial tension-compression test apparatus was used to prevent buckling of the specimens during the compression tests. In addition, a specially designed pure bending test apparatus used to measure the bending moment-curvature diagrams. Details of these devices can be found in the literature [6,7].
Some basic assumptions were considered in developing the analytical solution of pure bending as follows: • The plane section remains planar during pure bending; • The sheet thickness is assumed to be constant during the process; • The bending strain is proportional to its distance from the neutral surface.
The 1-, 2-and 3-directions were assumed to be the circumferential (θ), width (z), and radial (r) directions, respectively. The strain in the circumferential direction due to the pure bending was defined as: where ε b is the bending strain in the circumferential direction, y is the distance measured from the neutral plane and R is the bending radius. The neutral surface of pure bending is not always the same as the mid-surface. It can shift from the mid-surface for several reasons, some of which are applied tensile or compressive forces, TCA in yield points or the subsequent hardening. The TCA is only observed during the plastic bending and disappears for elastic bending. In this case, the neutral surface shifts to the outer or inner radius of the bending whenever the compressive yield point is lower or greater than the tensile yield point (RD and TD specimens in this study, respectively). Figure 1 shows the stress distribution through the thickness for an asymmetric material under pure bending. As it can be seen, the neutral surface with the length of L n and the radius of R n is not coinciding with the mid-surface with the length of L m and radius of R m . The parameter d shows the coordinate from the mid-surface to the neutral surface. In this paper, to consider the effects of neutral surface shift, a pseudo-neutral surface was defined as suggested by Lee et al. [16]. It is a surface where the strain equals the membrane strain and coincides with the real neutral surface in case of no applied tensile force. Hence, the general form of strain distribution of an asymmetric material for pure bending is given as:

Equivalent Stress and Strain
Usually, in the bending analysis of sheet metals, two kinds of loading conditions are assumed: uniaxial and plane strain-plane stress loading conditions. By applying these two loading conditions and considering the associated flow rule for plastic deformation, a relationship among stresses in three directions can be found, and the yield function for the specific loading condition can be written only as a function of bending stress, σ1 as described in Equation (3), where β is a function of loading condition. In fact, the yield function is assumed to be a homogeneous function of degree one with σ1 as the only independent variable. For asymmetric materials, β is a function of loading condition and the sign of the stress. Therefore, depending on the loading condition, it has different values in tension βt, and compression βc, as defined in Equation (4), By using the associated flow rule and applying equivalent plastic strain increment as the plastic multiplier, we have: From Equations (4) and (5), the bending plastic strain for an asymmetric material is given as: In this paper, to consider the effects of neutral surface shift, a pseudo-neutral surface was defined as suggested by Lee et al. [16]. It is a surface where the strain equals the membrane strain and coincides with the real neutral surface in case of no applied tensile force. Hence, the general form of strain distribution of an asymmetric material for pure bending is given as:

Equivalent Stress and Strain
Usually, in the bending analysis of sheet metals, two kinds of loading conditions are assumed: uniaxial and plane strain-plane stress loading conditions. By applying these two loading conditions and considering the associated flow rule for plastic deformation, a relationship among stresses in three directions can be found, and the yield function for the specific loading condition can be written only as a function of bending stress, σ 1 as described in Equation (3), where β is a function of loading condition. In fact, the yield function is assumed to be a homogeneous function of degree one with σ 1 as the only independent variable. For asymmetric materials, β is a function of loading condition and the sign of the stress. Therefore, depending on the loading condition, it has different values in tension β t , and compression β c, as defined in Equation (4), By using the associated flow rule and applying equivalent plastic strain increment as the plastic multiplier, we have: From Equations (4) and (5), the bending plastic strain for an asymmetric material is given as: Appl. Sci. 2020, 10, 3339 5 of 20

Plastic Deformation
Usually, the plastic flow of steels is concave down in both tension and compression loading. Generally, in order to capture the hardening behaviour of these materials, isotropic hardening rules, which can reproduce concave down curves, are applied. In this study, Voce and Swift hardening laws were considered and are expressed in Equations (7) and (8), respectively.
where σ t is the tensile yield point and B, C, K s and ε 0 are material parameters determined from the tensile SS curve. According to Equations (4), (6), and (7), and considering 1/β t = γ t and −1/β c = γ c , the bending stress during plastic deformation based on Voce hardening can be defined as: Similarly, the bending moment based on Swift hardening rule can be expressed as: In order to write the equilibrium condition for a sheet under bending (Figure 1), the force acting on the strip with a unit width and thickness of dy is considered. Therefore, from the equilibrium condition, the tension T and the moment of the force element M can be written as: where t is the thickness of the sheet during bending. Substituting Equation (9) into Equations (11) and (12) and considering the elastic and plastic deformations through the thickness the bending moment based on Voce hardening rule is defined as: where t 0 is the initial sheet thickness and y t and y c are the coordinates from the mid-surface to the interface of elastic and plastic zones in tension and compression regions, respectively, given as: Finally, by integrating Equation (13) the final expression for the bending moment under uniaxial loading is expressed as: and from the equilibrium of forces, we have: By integrating the above equation, the nonlinear relationship between neutral surface shift, d and the applied tensile force is defined as: In pure bending, T = 0 and the neutral surface shift, d, can be found by solving the abovementioned nonlinear equation.
Similarly, for the Swift hardening rule, the following expressions for bending moment and equilibrium of forces can be expressed as: Appl. Sci. 2020, 10, 3339 7 of 20

Yield Function
In this study, the Cazacu-Barlat 2004 yield function was applied [13]. It is an odd function of the principal values of the stress deviator S. Therefore, the yield criterion is sensitive to the sign of the stress σ 1 , and it can capture the tension-compression asymmetry of the material. The Cazacu-Barlat 2004 yield function can be expressed as: where τ y is the yield stress in pure shear, c is a material parameter and J 2 and J 3 are the second and third invariants of the stress deviator tensor, respectively. The material parameter c is expressed as: By assuming the uniaxial loading condition, the yield function is reduced to: The sign function of σ 1 is defined as: Therefore, we have: From the above equation, we have: where β is the function of stress sign, loading condition, and the material constant c. In the case of symmetry, the material constant c is degenerated to zero, and β has the same value of unity for both tension and compression loading. Therefore, for uniaxial condition, and considering 1/β t = γ t and −1/β c = γ c , we have: The values of c, γ t , and γ c , for DP980 and SUS304 in rolling-direction (RD) and transverse-direction TD, are calculated from Equations (22) and (28) and listed in Table 1. Based on the experimental stress-strain curves available in the literature [6,7], the tensile and compressive yield points, which were determined using the commonly used 0.2% offset strain method, are also listed in Table 1.

Elastic Unloading
In elastic springback, the unloading process, which causes an equal value of moment in the opposite direction, is assumed to be elastic. Therefore, as shown in Figure  Moreover, the magnitude of the springback is considered as the curvature change from the unloaded centreline curvature as shown in Figure 2. Hence, the equation for the elastic unloading of applied M is expressed as: where 1/Δr is the curvature change due to the unloading (springback) and Mu e is the elastic unloading moment.

Nonlinear Unloading
Generally, elastic unloading is applied to find the Δ(1/R). However, it was shown that the unloading is not entirely elastic, and there are small scale plasticity effects in springback [19,20]. In order to consider the effects of nonlinearity in unloading and springback, several methods were Hence, the equation for the elastic unloading of applied M is expressed as: where 1/∆r is the curvature change due to the unloading (springback) and M u e is the elastic unloading moment.

Nonlinear Unloading
Generally, elastic unloading is applied to find the ∆(1/R). However, it was shown that the unloading is not entirely elastic, and there are small scale plasticity effects in springback [19,20]. In order to consider the effects of nonlinearity in unloading and springback, several methods were suggested by various researchers. Sun and Wagoner [21] introduced a quasi-plastic-elastic (QPE) strain component in addition to plastic and elastic strain components. The QPE strain allows an additional strain recovery, which cannot be captured by elastic unloading. They used some aspects of two-yield surface approaches. Another recently proposed method is a piecewise-linear approach that adopts the concept of multiple yield surfaces [22]. Although these methods are successful in capturing the nonlinearity in unloading, they have the complexity of applying the multiple yield surface approaches which make them suitable for finite element analysis. The most common approach is to define a variable elastic modulus for the unloading [19,23,24]. In this method, the elastic modulus is the function of the accumulated plastic strain. This method is called the chord model, since the coefficients of the proposed function are determined by fitting the function to the measured chord moduli. Chord modulus is the slope of a chord or a straight line drawn between the unloading start point and the end point (σ = 0) of the unloading stress-strain curve. In this study, to describe the nonlinear elastic unloading analytically, the chord method was adopted.
According to Yoshida et al. [25], the variation of the effective elastic modulus with accumulated plastic strain can be expressed in terms of an exponential equation as described in Equation (30) where E 0 and E a are the initial and saturated elastic moduli, ε p is the equivalent plastic strain and ξ is a material parameter. As it can be seen in Figure 3, Equation (29) was fitted to the experimental data achieved from experimental stress-strain curves in the loading-unloading with different pre-strain values for SUS304TD specimens from Kuwabara et al. (2009a) [6]. The applied coefficients for the chord method are listed in Table 2, where T and C denote tension and compression, respectively.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 20 where E0 and Ea are the initial and saturated elastic moduli, ̅ is the equivalent plastic strain and ξ is a material parameter. As it can be seen in Figure 3, Equation (29) was fitted to the experimental data achieved from experimental stress-strain curves in the loading-unloading with different prestrain values for SUS304TD specimens from Kuwabara et al. (2009a) [6]. The applied coefficients for the chord method are listed in Table 2, where T and C denote tension and compression, respectively.  The equation for the unloading moment with the variable elastic modulus and the assumption of no neutral surface change during unloading, is given as  The equation for the unloading moment with the variable elastic modulus M VE U and the assumption of no neutral surface change during unloading, is given as where E a T , and E a C are saturated elastic moduli in tension and compression, ξ T and ξ T are the material parameters in tension and compression and ε P T and ε P c are equivalent plastic strains in tension and compression, respectively. Although in the above equation, the nonlinear relationship between elastic modulus and plastic strain is applied, the relationship between curvature change and unloading moment is linear. For simplicity, the neutral surface shift during the unloading, which is caused by the variation of elastic modulus, is neglected. By integrating the equation and substituting the equivalent plastic strain, the final equation for a nonlinear unloading moment with variable elastic modulus is derived as: For the equivalent plastic strains we have: where D, R, y t and y c are the neutral surface shift, bending radius and the coordinates from the origin to the interface of elastic and plastic zones in tension and compression regions, respectively. It is worth mentioning that the above parameters are constant during the unloading.

Pure Bending of DP980
Dual-phase 980 steel has significant applications in the advanced manufacturing of automobiles. They are frequently used in car body applications due to the fact of their high strength. In order to validate the analytical method and the derived equations, pure bending of DP980 sheets in rolling and transverse directions is studied in this section. The analytical results were compared with the experimental data from Maeda et al. [7].
As mentioned in Section 1, a few analytical studies are available in the literature, which are semi-analytical methods. In those studies, the uniaxial stress-strain curves of the materials were divided into several regions, and at each region, the curve was fitted with a specific equation [5][6][7]16,17]. The values of stress were found by using those piecewise fitted curves or directly from the stress-strain curves. Applying these semi-analytical approaches results in using too many parameters and limits the methods only to the uniaxial loading condition. However, continuum plasticity is based on the concepts of yield surface, flow rule and hardening evolution. Although these concepts can be applied in finite element simulations, to the best of the authors' knowledge, no analytical methods to date were developed for modelling the evolving plastic behaviours of asymmetric materials using continuum plasticity concepts.
One of the benefits of the analytical method presented in this study is that compressive hardening parameters are not necessary. Based on Equations (9) and (10), bending stress is the only function of the tensile hardening parameters and the compressive yield point. In order to fit the stress-strain responses of the material, the Swift hardening rule was applied, and Equation (10) fits to the experimental results in Figure 4. As it can be seen, the proposed model can capture the plastic flow in tension and compression successfully. The material parameters and hardening coefficients applied in the analytical study are listed in Table 3 for reference. The tensile Young's modulus and tensile and compressive yield points were measured from the experimental stress-strain curves provided in Figure 4, and the hardening parameters were determined by curve fitting. It is worth mentioning that the experimental stress-strain curves used in this study were considered as the reference monotonic tension-compression test data for the DP980. This implies that the material properties among the test specimens were constant and the material was selected from one batch.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 12 of 20 applied in the analytical study are listed in Table 3 for reference. The tensile Young's modulus and tensile and compressive yield points were measured from the experimental stress-strain curves provided in Figure 4, and the hardening parameters were determined by curve fitting. It is worth mentioning that the experimental stress-strain curves used in this study were considered as the reference monotonic tension-compression test data for the DP980. This implies that the material properties among the test specimens were constant and the material was selected from one batch.
(a) (b)  Pure bending moment-curvature of DP980 in RD and TD were determined based on Equation (19). The outputs based on symmetric and asymmetric models are compared with the experimental results in Figure 5. As it can be seen, the asymmetric model can successfully model the pure bending of DP980 in both rolling directions, while the symmetric model shows a discrepancy between measured and calculated results.   Pure bending moment-curvature of DP980 in RD and TD were determined based on Equation (19). The outputs based on symmetric and asymmetric models are compared with the experimental results in Figure 5. As it can be seen, the asymmetric model can successfully model the pure bending of DP980 in both rolling directions, while the symmetric model shows a discrepancy between measured and calculated results.
Pure bending moment-curvature of DP980 in RD and TD were determined based on Equation (19). The outputs based on symmetric and asymmetric models are compared with the experimental results in Figure 5. As it can be seen, the asymmetric model can successfully model the pure bending of DP980 in both rolling directions, while the symmetric model shows a discrepancy between measured and calculated results.

Pure Bending and Springback of Stainless Steel 304
The Voce hardening defined in Equation (7) is fitted to the experimental results in Figure 6. As it can be seen, both tensile and compressive fitted curves were in good agreement with the

Pure Bending and Springback of Stainless Steel 304
The Voce hardening defined in Equation (7) is fitted to the experimental results in Figure 6. As it can be seen, both tensile and compressive fitted curves were in good agreement with the experimental results in RD and TD [6], respectively. Moreover, the material parameters applied for fitting are listed in Table 4.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 20 experimental results in RD and TD [6], respectively. Moreover, the material parameters applied for fitting are listed in Table 4.  The calculated bending moment-curvatures based on Equation (16) are compared with the experimental results in Figure 7. As it is shown, for both RD and TD specimens, the curvatures considering the TCA were in great agreement with the experimental results. On the other hand, the calculated curves based on symmetric behaviour did not coincide with the experimental results. In this analytical method, the symmetric condition occurs when the asymmetry parameter of the yield function is zero (c = 0).  The calculated bending moment-curvatures based on Equation (16) are compared with the experimental results in Figure 7. As it is shown, for both RD and TD specimens, the curvatures considering the TCA were in great agreement with the experimental results. On the other hand, the calculated curves based on symmetric behaviour did not coincide with the experimental results. In this analytical method, the symmetric condition occurs when the asymmetry parameter of the yield function is zero (c = 0).
The calculated bending moment-curvatures based on Equation (16) are compared with the experimental results in Figure 7. As it is shown, for both RD and TD specimens, the curvatures considering the TCA were in great agreement with the experimental results. On the other hand, the calculated curves based on symmetric behaviour did not coincide with the experimental results. In this analytical method, the symmetric condition occurs when the asymmetry parameter of the yield function is zero (c = 0). The variation of neutral surface shift versus the curvature change is shown in Figure 8 for pure bending in RD and TD. As it can be seen, the natural surface shift is zero at the beginning due to the elastic bending. However, by increasing the bending curvature, elastoplastic bending occurs and the natural surface shift increases. Moreover, the neutral surface shifts with a high rate initially and slows down as the plastic deformation governs the bending behaviours. In addition, the neutral surface shifts to the inner and outer radius of bending under bending in TD and RD, respectively. The neutral The variation of neutral surface shift versus the curvature change is shown in Figure 8 for pure bending in RD and TD. As it can be seen, the natural surface shift is zero at the beginning due to the elastic bending. However, by increasing the bending curvature, elastoplastic bending occurs and the natural surface shift increases. Moreover, the neutral surface shifts with a high rate initially and slows down as the plastic deformation governs the bending behaviours. In addition, the neutral surface shifts to the inner and outer radius of bending under bending in TD and RD, respectively. The neutral surface shift in TD shows much lower values than those in RD, because the TCA in TD is much smaller than that in RD.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 14 of 20 surface shift in TD shows much lower values than those in RD, because the TCA in TD is much smaller than that in RD. The unloading paths based on constant and variable elastic moduli are compared with the experimental path for RD and TD specimens in Figure 9. As it can be seen, the unloading path based on Equation (32) with variable elastic modulus matched the experimental path greatly. Because the 1/Δr is equivalent to springback, it can be concluded that the predicted springback behaviours are improved significantly by applying the variable elastic modulus. The unloading paths based on constant and variable elastic moduli are compared with the experimental path for RD and TD specimens in Figure 9. As it can be seen, the unloading path based on Equation (32) with variable elastic modulus matched the experimental path greatly. Because the 1/∆r is equivalent to springback, it can be concluded that the predicted springback behaviours are improved significantly by applying the variable elastic modulus.
The unloading paths based on constant and variable elastic moduli are compared with the experimental path for RD and TD specimens in Figure 9. As it can be seen, the unloading path based on Equation (32) with variable elastic modulus matched the experimental path greatly. Because the 1/Δr is equivalent to springback, it can be concluded that the predicted springback behaviours are improved significantly by applying the variable elastic modulus.

Numerical Study of Pure Bending of DP980 and SUS304
The analytical solution based on Cazacu-Barlat 2004 and the Voce hardening rule was implemented in the UMAT subroutine in ABAQUS/Standard. The radial return algorithm [26], which applies an implicit integration technique, was applied and improved to solve the constitutive equations. In the first step, the pure bending tests were simulated, and the results were compared with the analytical results. For pure bending tests, a 2D strip with the same dimensions used in experimental tests was modelled. Due to the symmetric boundary condition about the plane along the centreline of the sheet, only half of the sheet sample was modelled along with the rigid body to model the clamp as a fixture to the sample, as shown in Figure 10. The clamp was modelled as a discrete rigid part and fused to the right edge of the strip using the tie constraint. The moment was applied at point A, which is the reference point of the clamp. Thirty elements through the thickness

Numerical Study of Pure Bending of DP980 and SUS304
The analytical solution based on Cazacu-Barlat 2004 and the Voce hardening rule was implemented in the UMAT subroutine in ABAQUS/Standard. The radial return algorithm [26], which applies an implicit integration technique, was applied and improved to solve the constitutive equations. In the first step, the pure bending tests were simulated, and the results were compared with the analytical results. For pure bending tests, a 2D strip with the same dimensions used in experimental tests was modelled. Due to the symmetric boundary condition about the plane along the centreline of the sheet, only half of the sheet sample was modelled along with the rigid body to model the clamp as a fixture to the sample, as shown in Figure 10. The clamp was modelled as a discrete rigid part and fused to the right edge of the strip using the tie constraint. The moment was applied at point A, which is the reference point of the clamp. Thirty elements through the thickness were considered, since it showed good convergence in the results. The 4 node bilinear plane stress quadrilateral (CPS4R) element was adopted to model the sheet. were considered, since it showed good convergence in the results. The 4 node bilinear plane stress quadrilateral (CPS4R) element was adopted to model the sheet. The bending moment-curvature curves of SUS304 and DP980 were obtained by using the proposed model by using UMAT subroutine and by using the default isotropic material in Abaqus (von-Mises), respectively. The results are compared with the analytical results in Figure 11. As it can be seen, the curves based on the default model of Abaqus did not coincide with the analytical curves. However, the analytical and simulation results using the UMAT subroutine were in good agreement for both RD and TD specimens. This shows that the model was implemented in UMAT successfully. The bending moment-curvature curves of SUS304 and DP980 were obtained by using the proposed model by using UMAT subroutine and by using the default isotropic material in Abaqus (von-Mises), respectively. The results are compared with the analytical results in Figure 11. As it can be seen, the curves based on the default model of Abaqus did not coincide with the analytical curves. However, the analytical and simulation results using the UMAT subroutine were in good agreement for both RD and TD specimens. This shows that the model was implemented in UMAT successfully.
The bending moment-curvature curves of SUS304 and DP980 were obtained by using the proposed model by using UMAT subroutine and by using the default isotropic material in Abaqus (von-Mises), respectively. The results are compared with the analytical results in Figure 11. As it can be seen, the curves based on the default model of Abaqus did not coincide with the analytical curves. However, the analytical and simulation results using the UMAT subroutine were in good agreement for both RD and TD specimens. This shows that the model was implemented in UMAT successfully.

Numerical Study of SUS304 under V-Bending
To further study the effect of TCA in the springback of the SUS304, the V-bending process was simulated, and the resultant springback was determined and compared with the experimental results. The V-bending simulations were performed on SUS304 plates with 0.3 mm thickness. The bending angle and bending radius were 90° 5.0 mm, respectively. In this simulation, a 2D sheet was modelled along with two discrete rigid parts as punch and die. Like the pure bending simulation, the CPS4R element type was chosen and thirty elements were used through the thickness to secure accuracy after mesh refinement. The V-bending simulation with three stages-before bending, bending and springback-is shown in Figure 12.

Numerical Study of SUS304 under V-Bending
To further study the effect of TCA in the springback of the SUS304, the V-bending process was simulated, and the resultant springback was determined and compared with the experimental results. The V-bending simulations were performed on SUS304 plates with 0.3 mm thickness. The bending angle and bending radius were 90 • 5.0 mm, respectively. In this simulation, a 2D sheet was modelled along with two discrete rigid parts as punch and die. Like the pure bending simulation, the CPS4R element type was chosen and thirty elements were used through the thickness to secure accuracy after mesh refinement. The V-bending simulation with three stages-before bending, bending and springback-is shown in Figure 12. In the V-bending process, the friction between sheet, punch and die is inevitable. It was reported by Matuszak [27], Trzepieciński et al. [28] and Ramezani et al. [29] that the average value of friction coefficient in steel sheet forming is in the range of 0.1 to 0.2, depending on the surface properties, geometrical features, velocity and lubricant condition, etc. Therefore, in this study, the average value of 0.15 was considered as the friction coefficient for the process. The predicted springback from the simulation is compared with the experimental results in Figure 13. As it can be seen, considering TCA in calculations improved the springback results significantly. The numerical springback values with TCA are closer to the experimental results, compared to those with tension. By applying TCA in calculations, the relative errors were reduced from −7.86% to −2.11% and 23.01% to 4.32% for TD and RD bending, respectively. In the V-bending process, the friction between sheet, punch and die is inevitable. It was reported by Matuszak [27], Trzepieciński et al. [28] and Ramezani et al. [29] that the average value of friction coefficient in steel sheet forming is in the range of 0.1 to 0.2, depending on the surface properties, geometrical features, velocity and lubricant condition, etc. Therefore, in this study, the average value of 0.15 was considered as the friction coefficient for the process. The predicted springback from the simulation is compared with the experimental results in Figure 13. As it can be seen, considering TCA in calculations improved the springback results significantly. The numerical springback values with TCA are closer to the experimental results, compared to those with tension. By applying TCA in calculations, the relative errors were reduced from −7.86 to −2.11% and 23.01 to 4.32% for TD and RD bending, respectively. of 0.15 was considered as the friction coefficient for the process. The predicted springback from the simulation is compared with the experimental results in Figure 13. As it can be seen, considering TCA in calculations improved the springback results significantly. The numerical springback values with TCA are closer to the experimental results, compared to those with tension. By applying TCA in calculations, the relative errors were reduced from −7.86% to −2.11% and 23.01% to 4.32% for TD and RD bending, respectively.

The Modified Analytical Method for Springback Prediction of SUS304 under V-Bending
In this section, a modified analytical method is developed for springback prediction of Vbending accounting for friction effects, noting that the friction has to be considered in the simulations, but it was always ignored in analytic modelling. Figure 14 shows the variation of the springback obtained from numerical, analytical and modified analytical models versus the friction coefficient for 3 bending radii of 3, 5 and 8 mm. As it can be seen, for R = 5 mm, the numerical springback values increase by increasing the friction coefficient from 0 to around 0.1. For greater values of the friction coefficient, the springback is approximately constant. It is interesting that as the friction coefficient increases, the numerical values

The Modified Analytical Method for Springback Prediction of SUS304 under V-Bending
In this section, a modified analytical method is developed for springback prediction of V-bending accounting for friction effects, noting that the friction has to be considered in the simulations, but it was always ignored in analytic modelling. Figure 14 shows the variation of the springback obtained from numerical, analytical and modified analytical models versus the friction coefficient for 3 bending radii of 3, 5 and 8 mm. As it can be seen, for R = 5 mm, the numerical springback values increase by increasing the friction coefficient from 0 to around 0.1. For greater values of the friction coefficient, the springback is approximately constant. It is interesting that as the friction coefficient increases, the numerical values get closer to the analytical values of springback. For example, for the bending radius of 5 mm and the friction coefficients of higher than 0.1, the numerical and analytical springback results are almost the same. This indicates that the springback behaviours in V-bending are sensitive to friction between the sheet and the die for lower friction coefficients. This indicates that the underlying parameter, which causes the difference between the numerical and analytical springback results, is the friction. Therefore, in this section, the analytical model is modified to take the frictional effects into account.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 17 of 20 get closer to the analytical values of springback. For example, for the bending radius of 5 mm and the friction coefficients of higher than 0.1, the numerical and analytical springback results are almost the same. This indicates that the springback behaviours in V-bending are sensitive to friction between the sheet and the die for lower friction coefficients. This indicates that the underlying parameter, which causes the difference between the numerical and analytical springback results, is the friction. Therefore, in this section, the analytical model is modified to take the frictional effects into account. Based on simulation results, it was observed that in the frictionless case, during the bending process, the sheet slipped towards the die and reached the bottom of the die before the bending process ended. This caused the sheet to bend extra at a lower degree (almost 45°) in the middle of the sheet and, consequently, forced it to unbend to the die radius at the end of the bending process. This behaviour causes non-uniform bending in the sheet which consequently reduces the springback. By increasing the friction coefficient, the sheet slip reduces and almost stops for higher friction. Based on simulation results, it was observed that in the frictionless case, during the bending process, the sheet slipped towards the die and reached the bottom of the die before the bending process ended. This caused the sheet to bend extra at a lower degree (almost 45 • ) in the middle of the sheet and, consequently, forced it to unbend to the die radius at the end of the bending process. This behaviour causes non-uniform bending in the sheet which consequently reduces the springback. By increasing the friction coefficient, the sheet slip reduces and almost stops for higher friction. Therefore, the radius reduction due to the slip was measured from the simulation results for different friction coefficients. As it can be seen in Figure 15, the radius reduction increases with increasing the coefficient friction and almost reaches zero for higher friction coefficients. This is in line with the variation of springback angle with the friction coefficient which indicates the relationship of radius reduction and springback. The variation of the radius reduction is different for TD and RD bending which shows that the radius reduction is related to material properties. In addition, the radius reduction is related to the bending radius, as it is different for different bending radii. To include the effects of friction in analytical solution, the springback caused by the radius reduction, which is called slip springback in this study, was deducted from the analytical springback results. Therefore, to find the analytical values of springback in the V-bending, the following relationship was defined: where ∆θ VB is the V-bending springback, ∆θ PB is the pure bending springback and ∆θ SS is the slip springback, respectively. Based on Equation (34), springback values of V-bending in RD and TD for the bending radius of 5 mm were calculated. According to Figure 15, for the bending radius of 5 mm and for friction coefficients higher than 0.1, the radius reduction is zero. Therefore, in this case, the slip springback was zero, and the pure bending and V-bending springbacks were the same. Figure 16 shows the experimental and calculated springback values of V-bending in RD and TD for 16 types of SUS304 with different strengths. The horizontal axis shows the material strength, σ2.9, which is the plastic flow stress when the strain reaches the value of 2.9% in the uniaxial tensile test. In order to estimate the properties of the material with different strengths, the material parameters listed in Table 4 were uniformly changed by keeping the asymmetry ratio constant until the desired strength was reached. While the analytical springback based on tension did not match with the experimental results, the analytical springback considering TCA showed great agreement with the experimental results in both rolling directions. The results based on tension showed a slight difference between springback in RD and TD, since the tensile flow stress was almost the same in both rolling directions. This clearly shows the importance of considering TCA in predicting springback of asymmetric metals. Based on Equation (34), springback values of V-bending in RD and TD for the bending radius of 5 mm were calculated. According to Figure 15, for the bending radius of 5 mm and for friction coefficients higher than 0.1, the radius reduction is zero. Therefore, in this case, the slip springback was zero, and the pure bending and V-bending springbacks were the same. Figure 16 shows the experimental and calculated springback values of V-bending in RD and TD for 16 types of SUS304 with different strengths. The horizontal axis shows the material strength, σ 2.9 , which is the plastic flow stress when the strain reaches the value of 2.9% in the uniaxial tensile test. In order to estimate the properties of the material with different strengths, the material parameters listed in Table 4 were uniformly changed by keeping the asymmetry ratio constant until the desired strength was reached. While the analytical springback based on tension did not match with the experimental results, the analytical springback considering TCA showed great agreement with the experimental results in both rolling directions. The results based on tension showed a slight difference between springback in RD and TD, since the tensile flow stress was almost the same in both rolling directions. This clearly shows the importance of considering TCA in predicting springback of asymmetric metals. uniformly changed by keeping the asymmetry ratio constant until the desired strength was reached. While the analytical springback based on tension did not match with the experimental results, the analytical springback considering TCA showed great agreement with the experimental results in both rolling directions. The results based on tension showed a slight difference between springback in RD and TD, since the tensile flow stress was almost the same in both rolling directions. This clearly shows the importance of considering TCA in predicting springback of asymmetric metals.

Conclusions
A novel analytic elastic-plastic method based on the Cazacu-Barlat 2004 yield function and isotropic hardening was developed to investigate the bending and springback behaviour of asymmetric materials. Compared to existing methods, the proposed analytical method can be an effective tool to simplify springback prediction of asymmetric sheet metals. This new method can be used for analytic study and numerical analysis by implementing it into the commonly-used commercial finite element analysis package Abaqus/Standard. As a case study, pure bending of stainless steel 304 and DP980 was investigated analytically and numerically. In addition, the V-bending of the stainless steel 304 plates was studied numerically, and the analytical method was modified. The following conclusions can be drawn based on the comparison of the analytical and numerical outputs of the experimental results: 1. The comparison of the analytical and experimental results shows that for pure bending, the new method can successfully capture TCA of the materials with less material parameters; 2. The linear and nonlinear springback of SUS304, determined based on the analytical method, showed great agreement with the experimental results; 3. Implementing the analytical method into the UMAT subroutine, the springback results from the V-bending simulation of SUS304 were a great match with the experimental results. By applying TCA effects in calculations, the relative error was reduced from −7.86 to −2.11% and 23.01 to 4.32% for TD and RD bending, respectively. 4. From simulations, it was observed that, in V-bending of SUS304, the springback increased by increasing the friction coefficient, depending on the bending radius and the material properties. In addition, it was observed that the numerical and analytical springback values were almost the same at higher friction coefficients, and the analytical method can be modified for lower friction coefficients.
5. The modified analytical method can successfully predict the springback behaviours of 16 types of SUS304 in the V-bending process, and the results were consistent the with experimental results.