Effect of the yield surface evolution on the earing defect prediction

1. INTRODUCTION

The cup drawing process is employed to turn the sheet metals into a three-dimensional cup with a high production rate. However, this process involves complex loading mechanisms such as bending, unbending, compression, stretching, and sliding concurrently (Chung and Shah, 1992Chung, K., Shah, K. (1992). Finite element simulation of sheet metal forming for planar anisotropic metals. Int. J. Plast. 8 (4), 453-476. https://doi.org/10.1016/0749-6419(92)90059-L.; Brabie et al., 2013Brabie, G., Costache, E.M., Nanu, N., Chirita, B. (2013). Prediction and minimisation of sheet thickness variation during deep drawing of micro/milli parts. Int. J. Mech. Sci. 68, 277-290. https://doi.org/10.1016/j.ijmecsci.2013.01.028.). Further, the lubrication between the tools, blank holder force, clearance between the punch and die, and materials anisotropy substantially affect the product quality (Taherizadeh, 2010Taherizadeh, A. (2010). Numerical simulation of sheet metal forming using non-sssociated flow rule and mixed isotropic-nonlinear kinematic hardening model. Doctoral Thesis, University of Windsor, Ontario, Canada. ; Taherizadeh et al., 2015Taherizadeh, A., Green, D.E., Yoon, J.W. (2015). A non-associated plasticity model with anisotropic and nonlinear kinematic hardening for simulation of sheet metal forming. Int. J. Solids Struct. 69-70, 370-382. https://doi.org/10.1016/j.ijsolstr.2015.05.013.). One of the defects encountered in these processes is earing that results from the non-homogeneous material flow induced by the sheet’s anisotropy. These defects should be removed by the trimming operations after the forming. Therefore, predicting the earing defect may lead the automotive companies to optimize the blank shape, prevent the material loss, and increase the production rate. From this point of view, finite element (FE) analyses become vital for predicting the earing defect without conducting any physical tests. Constitutive equations embedded or implemented into the FE software are employed to define the behaviors of the material, such as anisotropy, hardening, etc. Since the earing phenomenon stems from the anisotropy-based uneven material flow, the plasticity model within the constitutive equations is substantial.

Plasticity modeling has three steps, namely yield function definition, flow rule definition, and hardening rule definition (Firat, 2008Firat, M. (2008). A numerical analysis of sheet metal formability for automotive stamping applications. Comput. Mater. Sci. 43 (4), 802-811. https://doi.org/10.1016/j.commatsci.2008.01.068.; Taherizadeh et al., 2010Taherizadeh, A., Green, D.E., Ghaei, A., Yoon, J.W. (2010). A non-associated constitutive model with mixed iso-kinematic hardening for finite element simulation of sheet metal forming. Int. J. Plast. 26 (2), 288-309. https://doi.org/10.1016/j.ijplas.2009.07.003.). All these three steps are associated with the yield surface concept. The yield function describes the boundaries of the yield surface, the flow rule defines the direction of the incremental plastic strain, and the hardening rule depicts the evolution of this surface. The effect of yield surface on the earing in the cup drawing process is extensively studied (Chung and Shah, 1992Chung, K., Shah, K. (1992). Finite element simulation of sheet metal forming for planar anisotropic metals. Int. J. Plast. 8 (4), 453-476. https://doi.org/10.1016/0749-6419(92)90059-L.; Yoon et al., 1998Yoon, J.W., Barlat, F., Chung, K., Pourboghrat, F., Yang, D.Y. (1998). Influence of initial back stress on the earing prediction of drawn cups for planar anisotropic aluminum sheets. J. Mater. Process. Technol. 80-81, 433-437. https://doi.org/10.1016/S0924-0136(98)00148-4.; Yoon et al., 2000Yoon, J.W., Barlat, F., Chung, K., Pourboghrat, F., Yang, D.Y. (2000). Earing predictions based on asymmetric nonquadratic yield function. Int. J. Plast. 16 (9), 1075-1104. https://doi.org/10.1016/S0749-6419(99)00086-8.; Yoon et al., 2006Yoon, J.W., Barlat, F., Dick, R.E., Karabin, M.E. (2006). Prediction of six or eight ears in a drawn cup based on a new anisotropic yield function. Int. J. Plast. 22 (1), 174-193. https://doi.org/10.1016/j.ijplas.2005.03.013.; Izadpanah et al., 2016Izadpanah, S., Ghaderi, S.H., Gerdooei, M. (2016). Material parameters identification procedure for BBC2003 yield criterion and earing prediction in deep drawing. Int. J. Mech. Sci. 115-116, 552-563. https://doi.org/10.1016/j.ijmecsci.2016.07.036.; Park and Chung, 2012Park, T., Chung, K. (2012). Non-associated flow rule with symmetric stiffness modulus for isotropic-kinematic hardening and its application for earing in circular cup drawing. Int. J. Solids. Struct. 49 (25), 3582-3593. https://doi.org/10.1016/j.ijsolstr.2012.02.015.). Singh et al. (2018)Singh, S., Basak, S., Prakash, L., Roy, G.G., Jha, M.N., Mascarenhas, M., Panda, S.K. (2018). Prediction of earing defect and deep drawing behavior of commercially pure titanium sheets using CPB06 anisotropy yield theory. J. Manuf. Process. 33, 256-267. https://doi.org/10.1016/j.jmapro.2018.05.003. conducted physical experiments and FE analyses on a cup drawing process of commercially pure titanium. Firstly, they determined the limiting drawing ratio and then predicted the ear formation and optimized the blank shape. The performance of a yield function is evaluated by the r-value and yield stress ratio predictions of different material orientations. The maximum cup height in the earing profile is correlated with the r-value distribution, while the earing magnitude depends on the yield stress ratio predictions of the yield function (Singh et al., 2018Singh, S., Basak, S., Prakash, L., Roy, G.G., Jha, M.N., Mascarenhas, M., Panda, S.K. (2018). Prediction of earing defect and deep drawing behavior of commercially pure titanium sheets using CPB06 anisotropy yield theory. J. Manuf. Process. 33, 256-267. https://doi.org/10.1016/j.jmapro.2018.05.003.).

Moreover, the yield surface constructed by the yield function should be convex for the uniqueness of the incremental plastic strain direction (Firat, 2008Firat, M., Kaftanoğlu, B., Eser, O. (2008). Sheet metal forming analyses with an emphasis on the springback deformation. J. Mater. Process. Technol. 196 (1-3), 135-148. https://doi.org/10.1016/j.jmatprotec.2007.05.029.; Mendiguren et al., 2015Mendiguren, J., Rolfe, B., Weiss, M. (2015). On the definition of an kinematic hardening effect graph for sheet metal forming process simulations. Int. J. Mech. Sci. 92, 109-120. https://doi.org/10.1016/j.ijmecsci.2014.12.005.). However, in some cases, a proper yield surface shape definition may not be sufficient for earing defect prediction. In addition, decreasing the solution time is another issue related to the complexity of the constitutive modeling. Vladimirov et al. (2010)Vladimirov, I.N., Schwarze, M., Reese, S. (2010). Earing prediction by a finite strain multiplicative formulation for anisotropic elastoplastic materials. GAMM Mitt. 33 (1), 116-129. https://doi.org/10.1002/gamm.201010009. investigated the effect of the Hill48 yield function on the earing prediction of two alloys, namely Al-5wt.%Mg and AA2008, with the multiplicative decomposition approach. Grillo et al. (2015)Grilo, T.J., Valente, R.A.F., Alves de Souza, R.J. (2015). Modelling non-quadratic anisotropic yield criteria and mixed isotropic-nonlinear kinematic hardening: analysis of forward and backward-Euler formulations. Int. J. Mater. Form. 8, 533-547. https://doi.org/10.1007/s12289-014-1176-9. proposed implicit and explicit stress update schemes in order to decrease the solution time for the earing prediction of AA2090-T3 aluminum alloy. Cruz-Gonzalez et al. (2022)Cruz-Gonzalez, C., Vargas-Arista, B., León-Mendez I., Guzmán-Flores I. (2022). On the application of the forming limit diagrams for quality control of blanks for wheelbarrow of ASTM A1008 carbon steel. Rev. Metal. 58 (2), e218. https://doi.org/10.3989/revmetalm.218. related the forming limit diagram to the cup drawing and pointed out the importance of forming limit curves for the forming processes having different loading mechanisms. Recently, Habraken et al. (2022)Habraken, A.M., Aksen, T.A., Alves, J.L., Amaral, R.L., Betaieb, E., Chandola, N., Corallo, L., et al. (2022). Analysis of ESAFORM 2021 cup drawing benchmark of an Al alloy, critical factors for accuracy and efficiency of FE simulations. Int. J. Mater. Form. 15, 61. https://doi.org/10.1007/s12289-022-01672-w. carried out an extensive study related to the effects of several parameters, including the yield functions, associated and non-associated flow rules, hardening models, friction coefficient, stress update schemes, etc., on the performance of earing prediction.

Although many studies were carried out for earing formation prediction in cup drawing processes, the studies were generally concentrated on the influence of yield surface shape. The evolution of the yield surface also affects the earing profile (Yoon et al., 2010Yoon, J.H., Cazacu, O., Yoon, J.W., Dick, R.E. (2010). Earing predictions for strongly textured aluminum sheets. Int. J. Mech. Sci. 52 (12), 1563-1578. https://doi.org/10.1016/j.ijmecsci.2010.07.005.; Vladimirov et al., 2010Vladimirov, I.N., Schwarze, M., Reese, S. (2010). Earing prediction by a finite strain multiplicative formulation for anisotropic elastoplastic materials. GAMM Mitt. 33 (1), 116-129. https://doi.org/10.1002/gamm.201010009.; Taherizadeh et al., 2015Taherizadeh, A., Green, D.E., Yoon, J.W. (2015). A non-associated plasticity model with anisotropic and nonlinear kinematic hardening for simulation of sheet metal forming. Int. J. Solids Struct. 69-70, 370-382. https://doi.org/10.1016/j.ijsolstr.2015.05.013.). When the sheet slides over the punch and die shoulder, the sheet metal zones corresponding to the shoulder are exposed to bending-unbending loads, also called cyclic or reverse bending (Firat, 2007Firat, M. (2007). U-channel forming analysis with an emphasis on springback deformation. Mater. Des. 28 (1), 147-154. https://doi.org/10.1016/j.matdes.2005.05.008.; Firat et al., 2008Firat, M., Kaftanoğlu, B., Eser, O. (2008). Sheet metal forming analyses with an emphasis on the springback deformation. J. Mater. Process. Technol. 196 (1-3), 135-148. https://doi.org/10.1016/j.jmatprotec.2007.05.029.). The reverse bending situation induces the Bauschinger effect, known as the decrease in yield stress, when the material is loaded in the reverse direction. The Bauschinger effect can only be introduced by the kinematic hardening rule, which defines a motion of the yield surface by back stress. Although many kinematic hardening rules were proposed in the past, they can be decomposed into the linear models (Prager, 1956Prager, W. (1956). A new method of analyzing stresses and strains in work-hardening plastic solids. J. Appl. Mech. 23 (4), 493-496. https://doi.org/10.1115/1.4011389.; Ziegler, 1959Ziegler, H.A. (1959). A modification of Prager’s hardening rule. Q. Appl. Math. 17 (1), 55-65. http://www.jstor.org/stable/43634629.), multilinear or multifield models (Besseling, 1958Besseling, J.F. (1958). A theory of plastic and creep deformations of an initially isotropic material Showing Anisotropic Strain-Hardening, Creep Recovery, and Secondary Creep. J. Appl. Mech. 25 (4), 529-536. https://doi.org/10.1115/1.4011867.; Mróz, 1967Mróz, Z. (1967). On the description of anisotropic workhardening. J. Mech. Phys. Solids. 15 (3), 163-175. https://doi.org/10.1016/0022-5096(67)90030-0.; Dafalias and Popov, 1976Dafalias, Y.F., Popov, E.F. (1976). Plastic internal variables formalism of cyclic plasticity. J. Appl. Mech. 43 (4), 645-651. https://doi.org/10.1115/1.3423948.) and nonlinear models (Armstrong and Frederic, 1966Armstrong, P.J., Frederic, C.O. (1966). A mathematical representation of the multiaxial Bauschinger effect. Central Electricity Generating Board, Report RD/B/N 731.; Ohno and Wang, 1993Ohno, N., Wang, J.-D. (1993). Kinematic hardening rules with critical state of dynamic recovery. Part 1: Formulations and basic features for ratcheting behavior. Int. J. Plast. 9 (3), 375-390. https://doi.org/10.1016/0749-6419(93)90042-O.). Among them, nonlinear models can successfully depict the Bauschinger effect and transient behavior. This study aimed to assess the earing defect prediction performance of the combined hardening model that involves isotropic and a nonlinear kinematic hardening model simultaneously on the cup drawing process of the AA6016-T4 aluminum alloy. A combined hardening model consisting of isotropic and Armstrong-Frederic kinematic hardening models was implemented to depict the hardening response of the material, while a sixth-order complete homogeneous polynomial yield function was implemented in order to ensure perfectly describing the shape of the yield surface. The plasticity model was implemented into the Marc FE software through a user-defined material subroutine called Hypela2. The analyses were conducted separately, assuming only the isotropic hardening and combined isotropic-kinematic hardening rules to show the pure yield surface evolution effect. The FE results of punch force - stroke and earing prediction were compared with the experimental results.

2. MECHANICAL TESTS

This section gives information about the experimental procedure of the benchmark study on AA6016-T4 aluminum alloy with a thickness of 0.98 mm (Habraken et al., 2022Habraken, A.M., Aksen, T.A., Alves, J.L., Amaral, R.L., Betaieb, E., Chandola, N., Corallo, L., et al. (2022). Analysis of ESAFORM 2021 cup drawing benchmark of an Al alloy, critical factors for accuracy and efficiency of FE simulations. Int. J. Mater. Form. 15, 61. https://doi.org/10.1007/s12289-022-01672-w.). The standard tensile tests were carried out separately using the Shimadzu tensile testing machine at the Tokyo University of Agriculture and Technology, the University of Aveiro, and the University of Liege. In this work, the uniaxial tensile tests conducted at the Tokyo University of Agriculture and Technology were regarded. These tests were conducted from the rolling direction to the transverse direction, with intervals of 15° degrees. Conventional uniaxial tensile test specimens having 50 mm gauge length and 12.5 mm gauge width were used. The tests were carried out at a constant strain rate of 10^-3 s^-1, and a mechanical extensometer was used for sensitive strain measurement. In addition, the biaxial yield stress was procured from the biaxial tensile test in compliance with the ISO 16842 (2014)ISO 16842 (2014). Metallic materials - Sheet and strip - Biaxial tensile testing method using a cruciform test piece. Norma ISO. standards conducted on cruciform specimens. The test standards and the geometry of the test specimens were explained elaboratively in ref (Habraken et al., 2022Habraken, A.M., Aksen, T.A., Alves, J.L., Amaral, R.L., Betaieb, E., Chandola, N., Corallo, L., et al. (2022). Analysis of ESAFORM 2021 cup drawing benchmark of an Al alloy, critical factors for accuracy and efficiency of FE simulations. Int. J. Mater. Form. 15, 61. https://doi.org/10.1007/s12289-022-01672-w.). Within the scope of this study, the uniaxial tensile test data was utilized in order to obtain the hardening parameters and the yield function calibration of the material. Table 1 shows the experimental r-value, yield stress ratio directionalities, and the balanced biaxial yield stress data.

Here, θ refers to the orientation in which the parameters were procured. The r-values indicate the resistance against the thinning, while the yield stresses are the stress values corresponding to the 0.2% plastic strain for the related orientation. Yield stresses in different orientations were evaluated through their normalized values. All the yield stress values were divided by the reference yield stress procured from the rolling direction (RD). Besides, the orientation of 90° corresponds to the transverse direction (TD), respectively. The combined hardening model parameters were obtained from the reversal shear test conducted by the University of Aveiro. These data were also provided by the ESAFORM 2021 benchmark organizing team (Habraken et al., 2022Habraken, A.M., Aksen, T.A., Alves, J.L., Amaral, R.L., Betaieb, E., Chandola, N., Corallo, L., et al. (2022). Analysis of ESAFORM 2021 cup drawing benchmark of an Al alloy, critical factors for accuracy and efficiency of FE simulations. Int. J. Mater. Form. 15, 61. https://doi.org/10.1007/s12289-022-01672-w.). The reversal shear test data performed in RD was considered as the reference data. Single element simulations were considered for the combined hardening parameter’s determination, and the inverse method was adopted.

As the main issue, a cup drawing process of the AA6016-T4 alloy was carried out at University of Porto using a hydraulic universal testing machine of 300 kN capacity. The schematic view of the ESAFORM 2021 benchmark cup drawing process (Fig. 1a) and the tools (Fig. 1b) are illustrated below.

The cup drawing process comprises a blank, die, punch, blank holder, and body stoppers. The blank is placed between the die and the blank holder initially. In the sequel, the blank holder force was applied to the sheet. With the motion of the punch, the blank was subjected to flow into the cavity of the die. The diameter and the thickness of the blank were measured through a micrometer. The diameter of the blank sheet was 107.54 mm, corresponding to a 1.79 drawing ratio, and the blank’s thickness was 0.98 mm. The blank holder force of 40 kN was applied to the sheet during the deformation. The rigid body stoppers were used to block the blank holder and maintain the distance between the die and blank holder when the blank’s periphery reached the cup wall region. The drawing was suspended when the punch stroke reached 54 mm vertical distance. The punch speed was constant and 0.5 mm/s during the drawing operation. The earing profile of the cup was measured relatively through a digital dial gauge micrometer with 0.001 resolution. When the cup height was measured, the cup was rotated simultaneously to measure the height in all orientations. The total cup height was also measured with a high-precision height gauge (Habraken et al., 2022Habraken, A.M., Aksen, T.A., Alves, J.L., Amaral, R.L., Betaieb, E., Chandola, N., Corallo, L., et al. (2022). Analysis of ESAFORM 2021 cup drawing benchmark of an Al alloy, critical factors for accuracy and efficiency of FE simulations. Int. J. Mater. Form. 15, 61. https://doi.org/10.1007/s12289-022-01672-w.).

3. PLASTICITY MODEL

This section gives information related to the implemented plasticity tools comprising the yield criterion, the flow rule, and the hardening rule. The co-rotational plasticity approach was implemented in the numerical simulations. According to this approach, the total strain increment can be disintegrated additively into the elastic and plastic components, Eq. (1).

Incremental elastic strain incloses the deviatoric and hydrostatic parts, Eq. (2), while an incremental plastic strain is solely deviatoric (.

where, e represents the deviatoric strain component, ε_m denotes the hydrostatic strain component, and the term I is the kronecker delta. Similarly, incremental stress involves the deviatoric shape-changing (S_ij) and hydrostatic volume changing (σ_m) parts.

A correlation expressed in the Eq. (4) can be established between the incremental deviatoric strain and deviatoric stress.

where, G is the shear modulus. In the three-dimensional stress space, the yield criterion represents the bound between the elastic and plastic regions, Eq. (5).

Gotoh (1977)Gotoh, M. (1977). A theory of plastic anisotropy based on a yield function of fourth order (plane stress state). In. J. Mech. Sci. 19 (9), 505-512. https://doi.org/10.1016/0020-7403(77)90043-1. proposed the first fourth-order polynomial yield function. This criterion has nine analytically adjustable parameters; however, the yield surfaces constructed by this criterion suffers from non-convexity. Later, Soare (2007)Soare, S.C. (2007). On the Use of Homogeneous Polynomials to Develop Anisotropic Yield Functions with Applications to Sheet Forming. PhD Thesis, University of Florida, Florida, USA. and Soare et al. (2008)Soare, S.C., Yoon, J.W., Cazacu, O. (2008). On the use of homogeneous polynomials to develop anisotropic yield functions with applications to sheet forming. Int. J. Plast. 24 (6), 915-944. https://doi.org/10.1016/j.ijplas.2007.07.016. introduced a new calibration process, including an optimization step for Gotoh’s criterion regarding the convexity of the yield surface. Moreover, the HomPol6 yield function was introduced by Soare (2007)Soare, S.C. (2007). On the Use of Homogeneous Polynomials to Develop Anisotropic Yield Functions with Applications to Sheet Forming. PhD Thesis, University of Florida, Florida, USA.. This yield function can accurately predict the r value and yield stress ratios in different material orientations owing to the higher flexibility provided by having 16 parameters. The basic polynomial form of the criterion makes the determination of the plastic strain increment simpler. Moreover, biaxial yield stress and r value can be captured by the HomPol6 yield function, and this criterion is given in Eq. (6).

c₁, c₂, c₆, and c₇ can be calculated by the analytical equations obtained by considering the yield conditions (c₁ and c₇) and experimental r data (c₂ and c₆) in rolling and transverse directions, separately.

Other parameters of HomPol6 were computed by employing an optimization procedure to minimize the errors between the experimental and predicted data sets of r value and yield stress ratio (Eq. (11)).

The stress acting in an arbitrary orientation having θ angle to the RD can be decomposed into the components as in the Eq. (12-14)Eq. 12, 13, 14.

The corresponding yield stress ratio and r value prediction in the arbitrary orientation is given in Eq. (15) and Eq. (16), respectively.

Correlatively, the predicted biaxial yield stress ratio and r value can be computed by using the Eq. (17) and Eq. (18).

Associated flow rule adopted in this work can be expressed as Eq. (19).

Here, dλ, f and dε^p are the proportionality factor, plastic potential and the plastic strain increment respectively (Firat, 2008Firat, M. (2008). A numerical analysis of sheet metal formability for automotive stamping applications. Comput. Mater. Sci. 43 (4), 802-811. https://doi.org/10.1016/j.commatsci.2008.01.068.). In this work, back stress curve was characterized by the Armstrong - Frederic type hardening model that is given in the Eq. (20), Armstrong and Frederic (1966)Armstrong, P.J., Frederic, C.O. (1966). A mathematical representation of the multiaxial Bauschinger effect. Central Electricity Generating Board, Report RD/B/N 731..

where, C and γ are the material coefficients. The first term at the right side of the expression is the linear part of the back stress in line with the Prager’s rule, while the second term represents the fading dynamic memory (Paul et al., 2010Paul, S.K., Sivaprasad, S., Dhar, S., Tarafder, M., Tarafder, S. (2010). Simulation of cyclic plastic deformation response in SA333 C-Mn steel by a kinematic hardening model. Comput. Mater. Sci. 48 (3), 662-671. https://doi.org/10.1016/j.commatsci.2010.02.037.; Zhang et al., 2011Zang, S.L., Guo, C., Thuillier, S., Lee, M.G. (2011). A model of one surface cyclic plasticity and its application to springback prediction. Int. J. Mech. Sci. 53 (6), 425-435. https://doi.org/10.1016/j.ijmecsci.2011.03.005.). dp is the equivalent plastic strain increment and expressed as Eq. (21).

4. MATERIAL CHARACTERIZATION

This section presents the calculated parameters of the yield function and the combined hardening model implemented in this study. Besides the constructed yield surface, the r value and yield stress ratio predictions of the HomPol6 criterion were presented. Moreover, the back stress and the isotropic hardening curves were demonstrated in this section.

The identification procedure of the HomPol6 parameters was explained in Section 3, Plasticity Model. Parameters c₁, c₂, c₆, and c₇ can be computed through explicit formulas, while the other parameters were obtained using the least square optimization method. The HomPol6 yield function has 16 coefficients, and the obtained parameters are listed in Table 2. Figure 2 illustrates the predicted yield surface.

According to the normality postulation, a yield surface should be convex (Firat, 2008Firat, M. (2008). A numerical analysis of sheet metal formability for automotive stamping applications. Comput. Mater. Sci. 43 (4), 802-811. https://doi.org/10.1016/j.commatsci.2008.01.068.). This condition yields that the direction of the incremental plastic strain should be outward from the surface for each loading path. Figure 4 shows that the HomPol6 yield function produces a convex elliptical yield surface and captured the balanced biaxial yield stress data. Since the yield stress in balanced biaxial loading is lower than the uniaxial tension in RD, HomPol6 produced an elliptical yield locus. The elliptical surfaces are convenient for fulfilling the convexity condition. Figure 3 shows the yield stress and r-value directionalities of HomPol6.

As it is seen in Fig. 3, the experimental r-values and yield stress ratios tend to decrease from RD to 45° while they increase from 45° to TD. These general trends are also consistent with the literature studies, (Kuwabara et al., 2017Kuwabara, T., Mori, T., Asano, M., Hakoyama, T., Barlat F. (2017). Material modeling of 6016-O and 6016-T4 aluminum alloy sheets and application to hole expansion forming simulation. Int. J. Plast. 93, 164-186. https://doi.org/10.1016/j.ijplas.2016.10.002.; Zhang et al., 2019Zhang, S., Lu, Y., Shen, Z., Zhou, C., Lou, Y. (2019). Prediction of ductile fracture for Al6016-T4 with a ductile fracture criterion: Experiment and simulation. Int. J. Damage Mech. 29 (8), 1-23. https://doi.org/10.1177/1056789519865771.). Only a slight decrease occurred in the yield stress ratio directionality between 45° and 60°. However, this slight decrease is not essential since the yield stress ratios are not changed significantly at 45° and 60° (0.908 for 45° and 0.898 for 60°). The sheets are manufactured by the rolling processes. Furthermore, AA6016-T4 includes a T4 temper involving the solution heat treatment and natural aging processes. Therefore, the r values and yield stress ratios are affected by metallurgical events. The HomPol6 yield function could successfully predict the variation of the plastic properties with respect to the different material orientations. In particular, the yield stress ratio directionality was completely captured. Although the r-values changed dramatically for different orientations, the HomPol6 criterion showed good agreement with the experimental outcomes. This enhanced prediction capability relies on a high number of parameters to describe the material’s anisotropic responses. In earing prediction, successful r-value and yield stress ratio predictions are recommended since the cup height location is related to the r-value prediction, whereas the magnitude of the ear is related to the yield stress ratio predictions.

The cup drawing simulations of AA6016-T4 aluminum alloy were conducted based on the isotropic hardening rule and combined hardening rule separately in order to show the pure influence of the combined hardening rule. Therefore, the parameters for the isotropic hardening rule and combined isotropic-kinematic hardening rule assumptions were obtained by using the curve fitting method. The RD was regarded as the reference direction to the hardening parameter’s determination, and The Levenberg - Marquardt algorithm, which is derived from steepest descent and Newton algorithm (Sapna et al., 2012Sapna, S., Tamilarasi, A., Kumar, M.P. (2012). Backpropagation learning algorithm based on Levenberg Marquardt algorithm. Comput. Sci. Inf. Technol. 393-398. https://doi.org/10.5121/csit.2012.2438.; Okkan et al., 2018Okkan, U., Serbeş, Z.A., Gedik, N. (2018). MATLAB ile Levenberg-Marquardt algoritması tabanlı YSA uygulaması: Aylık yağış-akış modellemesi. Dicle Univ. J. Eng. 9, 351-362. https://dergipark.org.tr/tr/pub/dumf/issue/36316/412696.), was adopted in MATLAB Curve Fitting Tool. The Swift law was assumed (Eq. (22)) for the hardening characterization.

where, K and p refer to the strength coefficient and hardening exponent, respectively, while ε₀ represents the initial strain value. The ε_p denotes the plastic strain calculated from the following expressions.

Here ε_tot and ε_e are the total strain and elastic strain, respectively, while E represents the Elasticity modulus. In Eq. (22), σ_true depends on the plastic strain value ε_p. The hardening parameters only for isotropic hardening rule assumption were given in Table 3. Figures 4 shows the comparison between the experimental and the generated flow curves.

The experimental hardening curve was procured from the uniaxial tensile test conducted in RD. In the benchmark study, the uniaxial tensile tests were conducted separately at Aveiro University, Tokyo University, and Liege University (Habraken et al., 2022Habraken, A.M., Aksen, T.A., Alves, J.L., Amaral, R.L., Betaieb, E., Chandola, N., Corallo, L., et al. (2022). Analysis of ESAFORM 2021 cup drawing benchmark of an Al alloy, critical factors for accuracy and efficiency of FE simulations. Int. J. Mater. Form. 15, 61. https://doi.org/10.1007/s12289-022-01672-w.). In the present study, the test results obtained from Tokyo University were utilized for the hardening curve determination. The material hardens when the yield strength of the material is exceeded. Correspondingly, the required force value producing the stress should increase in order to maintain the plastic deformation. The experimental and fitted (in compliance with the Swift law) curves increase for this reason and mostly overlap with each other.

The reversal shear test data provided by the ESAFORM 2021 benchmark organizers (Habraken et al., 2022Habraken, A.M., Aksen, T.A., Alves, J.L., Amaral, R.L., Betaieb, E., Chandola, N., Corallo, L., et al. (2022). Analysis of ESAFORM 2021 cup drawing benchmark of an Al alloy, critical factors for accuracy and efficiency of FE simulations. Int. J. Mater. Form. 15, 61. https://doi.org/10.1007/s12289-022-01672-w.) were employed in order to obtain the combined isotropic-kinematic hardening parameters. To this end, a representative unit volume element (RVE) was utilized, and an inverse method was adopted. Reversal shear test conditions were applied to the RVE (Fig. 5). By adjusting the kinematic and isotropic hardening parameters simultaneously, the reversal behavior of the material was tried to be estimated.

The total hardening curve equals to the sum of the reference yield stress, isotropic hardening, and kinematic hardening curves. The mathematical expression of the total hardening curve was given in Eq. (24). Besides, Table 4 summarizes the combined hardening parameters.

Here, the second and the third terms on the right side of the Eq. (24) represent isotropic hardening and kinematic hardening, respectively. Besides, σ_ref is the initial radius of the yield surface. Correspondingly, Q and b are the material coefficients associated with isotropic hardening that controls the proportional expansion of the yield surface, while C and γ are the material constants related to the backstress curve controlling the motion of the yield surface (Chaboche and Rousselier, 1983Chaboche, J.L., Rousselier, G. (1983). On the plastic and viscoplastic constitutive equations-Part I: rules developed with internal variable concept. J. Press. Vessel Technol. 105 (2), 153-158. https://doi.org/10.1115/1.3264257.; Zang et al., 2011Zang, S.L., Guo, C., Thuillier, S., Lee, M.G. (2011). A model of one surface cyclic plasticity and its application to springback prediction. Int. J. Mech. Sci. 53 (6), 425-435. https://doi.org/10.1016/j.ijmecsci.2011.03.005.). The third term is the derived form of the Eq. (20). Based on this expression, the backstress curve saturates at a constant stress level of C/γ. In addition, Fig. 6 demonstrates the estimated reversal shear response and predicted isotropic hardening and back stress curves.

As it can be seen in Fig. 8a, the combined hardening model accurately captures the transient behavior and the Bauschinger effect. Only the permanent softening was not predicted; that is not the primary issue in the cup drawing. The Bauschinger effect and transient behavior are essential. As for Fig. 8b, the backstress curve saturates at approximately 90 MPa, while the isotropic hardening curve has a consistent upward trend.

5. CUP DRAWING SIMULATIONS

The cup drawing process was analyzed in implicit Marc software coupled with the Hypela2 user-defined material subroutine. Blank was discretized with full integration constant dilatational hexahedral elements (Hex7 in Marc software), which are free from shear locking (Marc, 2018aMarc, MSC. (2018a). Volume A: Theory and User Information. MARC Analysis Research Corporation, USA.; Marc, 2018bMarc, MSC. (2018b). Volume B: Element library. MSC. Software Corporation, USA.). Hence, these elements ensure precise results in forming operations (Taşdemir, 2022Taşdemir, V. (2022). Finite element analysis of the springback behavior after V bending process of sheet materials obtained by differential speed rolling (DSR) method. Rev. Metal. 58 (2), e219. https://doi.org/10.3989/revmetalm.219.). Three thousand seventy-two elements were generated for the blank meshing. Die, punch and blank holder were assumed as rigid surfaces, and control nodes were assigned to the blank holder and punch in order to define the force and displacement boundary conditions. Rigid body stoppers were also modeled to maintain the distance between the die and the blank holder. Figure 7 illustrates the FE model of the cup drawing process and the generated mesh layout.

In the numerical approach, the blank was clamped between the blank holder and die, then the blank holder force of 40 kN was applied to the blank. Then the punch started moving in the vertical direction. When the punch stroke reached 54 mm, the simulation was suspended. The control nodes were attributed to each component. Thus, the components were fixed in the radial direction. The punch displacement and the blank holder force were specified through these nodes in the drawing (z) direction. The friction coefficient between the blank and the other tools was set to 0.1, and the Coulomb bilinear friction type was assumed. The segment-to-segment contact algorithm with the Lagrangian multiplier method was adopted to minimize the penetration.

6. RESULTS AND DISCUSSION

This section presents a summary of the results containing the prediction of ear formation and punch force-stroke responses. Since the ironing issue was reported by most participants of the ESAFORM 2021 benchmark study, the clearance between the punch and die was increased from 1.2 mm to 1.4 mm in the simulations in order to show the performance of the implemented plasticity model. Besides, the wear on the die radius region was also discussed and a higher clearance than the initially determined was indicated in the ESAFORM 2021 benchmark study (Habraken et al., 2022Habraken, A.M., Aksen, T.A., Alves, J.L., Amaral, R.L., Betaieb, E., Chandola, N., Corallo, L., et al. (2022). Analysis of ESAFORM 2021 cup drawing benchmark of an Al alloy, critical factors for accuracy and efficiency of FE simulations. Int. J. Mater. Form. 15, 61. https://doi.org/10.1007/s12289-022-01672-w.). In the first stage of the study, the cup drawing simulations were conducted using the HomPol6 yield criterion coupled with the isotropic hardening rule assumption. The convergence criterion was regarded as the absolute residual force of a constant value in order to restrict the sudden peaks in punch force response. Firstly, the results of isotropic hardening coupled with HomPol6 criterion are presented in Fig. 8.

The HomPol6 criterion coupled with isotropic hardening captured the experimental punch force-stroke behavior accurately (Fig. 8a). The sudden increase and decrease in punch force response stem from the contact definition of the tools. FE solver Marc is trying to satisfy the yield condition and restrain the penetration of the die tools simultaneously. The absolute residual force convergence method restricts the sudden changes in punch force. However, this model predicted lower ear magnitude when compared to the experimental earing profile as it is observed in Fig. 8b. It may result from the inadequate gap between the die and the punch at the cup wall region. The combined hardening rule is adopted in the second stage, coupled with HomPol6 criterion. Figure 9 demonstrates the comparison of the punch force-stroke response and the earing profile predictions.

It was seen from Fig. 9a that the combined hardening rule captured the punch force-stroke curve obtained from the physical experiment. There are only minor differences between the results obtained by the hardening rules, which shows that the hardening rule slightly affects the punch force response. Nonetheless, a distinct difference was observed in the earing prediction results as shown in Fig. 9b. The combined hardening model simultaneously increased the mean cup height level and earing magnitude. Thus, the combined hardening results turned out to be more consistent with the experimental results. Correspondingly, the combined hardening model showed better agreement when compared to the experimental results. In addition, at the vicinity of 270° degrees, two ears were observed in the physical test, which this model approximately captures. It is noteworthy that the kinematic part of the combined hardening model was represented by the Armstrong-Frederic model. This model involves a linear part and a memory term that endows a nonlinear feature to the backstress response. Including these two effects in the numerical simulations may have played an essential role in increasing the prediction performance of ear formation. Figure 10 demonstrates the numerically and experimentally drawn cups.

As it can be seen in Fig. 10b, there are eight prominent ears in the results of isotropic hardening. However, the ear formation predicted by the combined hardening rule was more consistent with the experimentally drawn cup (Fig. 10c). Figure 11 shows the mean cup height values of the experimental and numerical results.

Figure 11 shows that the mean cup height was noticeably increased when the combined hardening effect was included. Although the minimum cup heights remain stable, the maximum values of cup heights are dramatically increased. Figure 12a shows the cup height for distinct orientations, while Fig. 12b summarizes the mean cup height and the earing percentages. The earing percentage was measured by using Eq. (25).

According to Fig.12a, the increase in the cup height values is palpable. Besides, an apparent asymmetry in the experimental values is observed. For 0° and 270°, the cup heights predicted by the combined hardening rule were close, while primarily for 180°, the asymmetry in experimental height is especially salient. This asymmetry may arise from the misalignment of the blank; however, it increases the mean cup height, as seen in Fig. 12b. Considering the earing percentage, the combined hardening rule shows better agreement with the experiment. It should be noted that the sheet is exposed to bending-unbending or reverse situations along the punch and die shoulders. Therefore, a difference in yield stress is induced along these regions due to the reversal loading conditions. Since the yield stress is decreased, material flow is facilitated on these regions that affect the cup wall and the perimeter of the drawn cup. This effect is taken into consideration with the kinematic hardening, and the prediction performance is increased. In order to examine the reason for this enhanced prediction capability, the cup wall thickness values were investigated. Figure 13 compares the cup wall thickness values for 25 mm (H25) and 30 mm (H30) height from the cop bottom. The experimental outcomes were acquired from the study of Habraken et al. (2022)Habraken, A.M., Aksen, T.A., Alves, J.L., Amaral, R.L., Betaieb, E., Chandola, N., Corallo, L., et al. (2022). Analysis of ESAFORM 2021 cup drawing benchmark of an Al alloy, critical factors for accuracy and efficiency of FE simulations. Int. J. Mater. Form. 15, 61. https://doi.org/10.1007/s12289-022-01672-w..

The thickness values were notably decreased when the combined hardening rule was assumed, and the thickness values of the combined hardening assumption converged to the experimental profile. The improvement in the prediction capability for the combined hardening assumption may result from decreased thickness since the clearance value was increased to 1.4 mm. The cup height values are higher than 30 mm, which is the verge of the measured thickness profiles. Noting that the thickness predicted by the isotropic hardening assumption reached 1.4 mm approximately at 30 mm height, it may be concluded that the predicted thickness reaches the clearance threshold at the cup’s perimeter. Hence, the drawn cup based on the isotropic hardening rule was subjected to an ironing issue. On the other hand, the predicted thickness value was decreased with the combined hardening assumption that may prevent the drawn cup from being exposed to ironing issues.

7. CONCLUSIONS