Analysis of Side-Wall Wrinkling in Deep Drawing Processes

Wrinkling, developed due to compressive instability, is one of the main failure modes in deep drawing processes. Depending on the die geometry and blank holder force, it can be observed as flange wrinkling and/or side-wall wrinkling. In this study, side-wall wrinkling in cylindrical cup drawing process is investigated by using different constitutive models that are formed by using CPB06ex2, Hill’48, BBC2008-8p and BBC2008-16p yield criteria with isotropic hardening, and von Mises yield criterion with isotropic, kinematic and combined hardening. Numerical simulations are performed by implementing CPB06ex2, BBC2008-8p and BBC2008-16p models to a commercial finite element code through user subroutines. AA5042-H2 aluminium alloy is selected as a sheet material with 0.2083 mm thickness. In order to verify the developed models, the experimental results of the previous works for wrinkling profile and punch force evolution are employed. The results are compared and presented for the considered yield criteria.


Introduction
Deep drawing is one of the main sheet metal forming processes, which is employed especially in automotive and aerospace industries. In the deep drawing process, the sheet metal is clamped by the blank holder and then the punch draws the sheet into the die cavity in order to achieve the desired shape. During the process, the sheet may fail due to excessive tensile, shear or compressive stresses. Necking, tearing and earing are some examples of these failures of the sheet metal.
Wrinkling is one of the main failure modes arising in deep drawing processes. It is a kind of buckling occurred due to compressive instability. There are two types of wrinkling; flange wrinkling and side-wall wrinkling. Since side-walls are comparatively unsupported by forming tools, the required blank holding force to prevent side-wall wrinkling is much higher than that for flange wrinkling [1].
For deep drawing production, both types of wrinkling are undesired for stamping parts not only for visual appearance but also large wrinkles may damage dies and lead to formed sheet to be out-ofusage. For these reasons, prediction of the wrinkling is crucial to produce final parts with higher quality and less cost.
In literature, some analytical and numerical methods have been developed for prediction of the wrinkling formation. One of these is the energy method which has been widely used [1][2][3][4][5]. Using this method, Wang et al. [1] equated the external work done by membrane forces to the bending energy so that the critical buckling stress can be calculated. When the actual applied compressive stress becomes higher than the critical buckling stress, wrinkling is expected to occur. Besides, the bifurcation theory has also been employed to predict the occurrence of wrinkling [6,7]. Chu and Xu [7] examined flange wrinkling by using the general bifurcation theory to find out the critical condition for the onset of wrinkling. It was found that the bifurcation theory provides reliable results. However, instead of the complex plastic bifurcation theory, finite element simulation can also be preferred as a method due to its simplicity and reliability. For numerical simulations, several constitutive models have been proposed to estimate the behavior of the material during plastic deformation. The selection of appropriate constitutive model is crucial for modeling and predicting wrinkles at the correct time with the correct form.
In the present study, the capability of different constitutive models on simulation of side-wall wrinkling for cylindrical cup drawing process is investigated through finite element method [8]. The employed constitutive models are formed by applying CPB06ex2, Hill'48, BBC2008-8p and BBC2008-16p yield criteria with isotropic hardening, and von Mises yield criterion with kinematic, isotropic and combined hardening. Numerical simulations are performed by using the ABAQUS finite element code. CPB06ex2, BBC2008-8p and BBC2008-16p yield functions are implemented to the code through the developed user material subroutines whereas other constitutive models are already available in the program. Also, the predictions of the models are verified with the experimental data obtained from literature.

Yield Functions
CPB06ex2 Yield Criterion. According to the orthotropic yield criterion proposed by Plunkett et al. [9], the effective stress for materials with no tension-compression asymmetry can be calculated as where  indicates the degree of homogeneity and ( ) are the principal values of transformed tensors  and  , respectively [10]. The transformed tensors are determined by where s denotes deviatoric stress tensor and C and C  are the fourth-order tensors of material anisotropy coefficients. For plane stress condition, the components of the transformed tensor  are of the form [10]  ( In the above equations, The material anisotropy coefficients C and C  can be found by minimizing the error function. The error function is In the above equation, the values indicated by BBC2008 Yield Criterion. The BBC2008 plane-stress yield criterion is developed by Comsa and Banabic for highly orthotropic materials [11]. The yield function is of the form where k , r l 1 , r l 2 , r m 1 , r m 2 , r m 3 , r n 1 , r n 2 and r n 3 are defined as material parameters. For BCC alloys, k is equal to 3 while for FCC materials, k is equal to 4. Other quantities can be determined by the minimization of the following error function:

Achievements and Trends in Material Forming
where  y and  r are the normalized uniaxial yield stress and anisotropy coefficient along a direction that is inclined by the angle  from the rolling direction, respectively. Moreover, while b y denotes the bi-axial yield stress, b r is the anisotropy coefficient in case of bi-axial traction.

Implementations of Yield Functions
The BBC2008-8p, BBC2008-16p and CPB06ex2 yield functions are implemented to ABAQUS/Explicit finite element code employing the user material subroutine VUMAT. The details of the implementation are given in the following. The yield surface can be defined as where  and p  denote the effective yield stress and plastic strain, respectively. The elastic stressstrain relationship is of the form and e C is the fourth order elasticity tensor. From the normality rule, the plastic strain increment is calculated as where dλ is the plastic multiplier. In case of isotropic hardening, For the calculation of the stress state, the return mapping algorithm is preferred. In this elastoplastic algorithm, the trial elastic stress is calculated by using the following equation [12] e e old trial new If the trial stress is lower than the current yield stress, the deformation is elastic. Otherwise, the plastic deformation occurs and the new stress state can be determined with plastic corrector In many of the previous studies that were about the implementation of yield functions into a finite element code, the return mapping algorithm was used and it was stated that this algorithm provides reasonable results [10,13,14,15].

Finite Element Model
In this study, the simulation of side-wall wrinkling is performed by considering the Benchmark 4 that was proposed in Numisheet 2014 conference. For the verification, the experimental data of the Benchmark 4 reported by Neto et al. [16] is compared with the predictions of the developed constitutive models. Fig. 1 shows the sheet, punch, die and pressure pad, and the related dimensions are given in Table 1. The magnitude of the blank holding force was set equal to 8.9 kN. In addition, the Coulomb friction model was employed and the friction coefficient between the tools and sheet was taken from the specifications of the benchmark as 0.03. In numerical simulations, the quarter of the sheet was modelled by four node shell elements with reduced integration (ABAQUS S4R) while the punch, die and blank holder were assumed as rigid bodies. In the initial mesh, there were 1875 elements. The FE model of the cylindrical cup drawing process and meshing of the quarter sheet are shown in Fig. 2.

Results and Discussion
The simulation of side-wall wrinkling is investigated by employing different constitutive models for AA5042 aluminum alloy sheet. The mechanical properties and anisotropy coefficents of AA5042 are given in Table 2. C are given in Table 3. For plane stress, the coefficients 44 C , 44 C , 55 C and 55 C are equal to zero [18]. Also, for the BBC2008-8p and BBC2008-16p criteria, the material parameters are provided in Table 4 and Table 5, respectively.  Table 4. Material parameters of BBC2008-8p yield function [11] Material parameters of BBC2008-8p k s w  Table 5. Material parameters of BBC2008-16p yield function [11] Material parameters of BBC2008-16p k s w The comparison of the punch force variations with respect to punch displacements obtained by using the developed constitutive models and experimental data is shown in Fig. 3. The agreement between the numerical and experimental results is found to be satisfactory. However, the two versions of the BBC2008 model underestimate the punch force. Also, the most accurate result is provided by the CPB06ex2 based constitutive model.
In order to investigate the effect of different constitutive models on prediction of wrinkling, the side-wall wrinkling profiles were obtained at two different cross sections perpendicular to the drawn cup axis at 3.5 mm and 5.5 mm in the negative z direction. The measurements were taken from the bottom of the drawn cup, which had been formed by 20 mm punch displacement. In Fig. 4, the deformed cup and the used coordinate system are shown.
Comparing the predicted and experimental wrinkle numbers; 13 wrinkles were observed at both cross sections in the experiments, on the other hand, in numerical simulations, the CPB06ex2 and Hill'48 based constitutive models predicted 14 wrinkles while other constitutive models predicted 12 wrinkles. It is believed that experimental errors such as eccentricity may affect the number of wrinkles observed [16]. Consequently, it can be said that, the numerical results agree well with the experimental results.
In addition, the estimated side wall wrinkling profiles were compared with the results of experiments (Fig. 5-11).
In the plane mm z 5 .
3 − = , the numerical results show similar behavior with experiment, except the kinematic hardening model. At the regions in the transverse and rolling direction of the sheet, compared to von Mises yield criterion with isotropic and combined hardening, the BBC2008-8p and BBC2008-16p models estimate the wrinkling profiles more accurately. Moreover, the agreement between the predictions of the CPB06ex2 model and the experimental results is found considerably better.
In the plane mm 5.5 z − = , there are appreciable deviations between the results obtained numerically and experimentally. The estimated wrinkle amplitude is higher than the experimental one. The reason of that may be attributed to the use of the recommended constant friction coefficient value; the friction coefficient may vary as the punch moves. The contact pressure, the surface roughness of tools and sheet and the lubrication conditions can affect the friction coefficient value. Additionally, considering both cross sections, the von Mises with kinematic hardening model is unsuccesfull in predicting the side-wall wrinkling.

Conclusions
In this study, the capability of various constitutive models on the prediction of side-wall wrinkling is examined. Also, for punch force variation and wrinkling profile the experimental results of Numisheet 2014 Benchmark 4 are compared with the predictions of the developed constitutive models. Following conclusions have been obtained from this study: • The comparison between the experimental and numerical punch force variation shows that the two versions of the BBC2008 yield criterion underestimate the punch force variation. Moreover, the predictions of the von Mises with isotropic hardening, kinematic hardening and combined hardening, Hill'48 and CPB06ex2 models are similar to the experimental data but the closest result to the experiment is obtained by using the CPB06ex2 model.
• For side-wall wrinkling, at mm z 5 . 3 − = , the predictions of the constitutive models for wrinkling profile are in good agreement with experiments, except the von Mises with kinematic hardening model. The numerical results obtained by using the BBC2008-8p and BBC2008-16p yield criteria are very similar to the experimental results especially at the regions in the rolling and transverse direction. Also, the most accurate result is acquired employing the CPB06ex2 model comparing to the other models.
• According to the comparison of side-wall wrinkling profile in the plane mm 5.5 z − = , there are considerably deviation between the experimental amplitude of wrinkles and the numerical results. In addition of that, considering both cross sections, the kinematic hardening model does not give reasonable results for the prediction of side-wall wrinkling.
• The agreement between the predictions of the constitutive models for wrinkle number and the experimental results is found to be satisfactory.