An analytical investigation on the wrinkling of aluminium alloys during stamping using macro-scale structural tooling surfaces

Structural surface texturing is believed to be a promising approach to modify tribological and thermal performances of tooling for sheet stamping processes. However, a fundamental study on the surface texturing design and resulting material deformation is currently lacking. In this paper, an advanced analytical buckling model specifically for the utilisation of textured tools at macro-scale, comprising dislocation-driven material model, isotropic yield criteria, bifurcation theory and Donnell-Mushtari-Vlasov (DMV) shell structure theory was established. The developed analytical buckling model was validated by cylindrical deep drawing experiments. Further finite element (FE) simulations with the implementation of material model via userdefined subroutine were also used to validate the bucking model for large surface texture designs. Effects of theoretical assumptions, such as yield criterion, boundary condition and test-piece geometry, on the accuracy of model prediction for wrinkling were investigated. It was found that the von Mises yield criterion and hinged boundary condition exhibited more accurate predictions. In addition, the DMV shell theory made this model most representative for large structural texturing designs. Furthermore, the implementation of induced shear strain component has an important effect on precisely predicting the wrinkling occurrence. The advanced analytical models developed in this study combine various classical mechanics, structure stability and material modelling together, which provides a useful tool for tooling engineers to analyse structural designs.


Introduction
Tooling surface performances, such as triblogical and thermal, are believed to play an imporant role in controlling material deformation and post-strengthening for hot stamping sheet metals [1]. In recent years, hot stamping of high strength aluminium alloys has become a popular and promising lightweight strategy for automotive OEMs (Original Equipment Manufacturers) due to the enviromental and fuel economy concerns [2,3]. Considering the highly reactive and formability behaviour of high strength aluminium alloys at elevated temperatures, tooling requirments have become more stringent [4,5]. Tooling technologies with low interfacial friction coefficient, good wear resistance and thermal properteies are urgently required. Structual surface texturing has been extensively investigated and applied in particular engineering areas due to its abilities of modifying tribological and thermal performances [6]. Surface texture can be roughly divided into micro-and macro-scales according to the texture dimension. Micro-scale surface texture, normally at µm or nm scale, is believed to be able to improve tribological properties of tool surfaces with the combination of lubricants [6]. Costa and Hutchings [7] investigated micro-scale patterned tool surface tribological performance using a strip drawing process, a significant reduction of friction coefficient was obtained as these surface textures were believed to exhibit a reservior function for lubricant to improve the lubrication effect. In addition, wear resistance of tool surfaces can also be increased as test-piece material debris during interface sliding are stored in the texture hollows. The great loss of hard particles at the tool/blank interface results in a significant improvement in wear resistance of tooling. Compared with the extensive research on micro-texturing of tool surfaces, limited research has been performed for macro-scale surface texturing. Franzen et al. [8] investigated the tribological performance of macroscopic textured tool surfaces using strip drawing process. The textured macro-scale tool surfaces can increase the interficial friction coefficient compared to flat tool surfaces, which might be a potentional approach to replace conventional draw beads for sheet stamping.
To this extent, most research focused on fundamental triblogy studies of textured tool surfaces. However, effects of surface texturing on draw-ability and thermal properties of hot stamping aluminium alloys have not been thoroughly studied. Zheng et al. [9] investigated the hot draw-abilty of aluminium alloys with a one-dimensional macro-scale textured tool surface both experimentally and numerically. A significant improvement of draw-ability was achieved due to a more uniform temperature field within the blank. The reasons are that, the utilisation of tool surface texture can reduce contact area between blank and tool material 3 significantly. Moreover, the interfacial heat transfer coefficient at the texture hollow was also decreased due to the vacant space [10]. However, the texture effect was only examined for blanks under a one-dimensional stress state. For hot stamping complex-shaped sheet components with a two-dimensional contary sign stress state, the compressive stress may cause material located above the texture hollow to buckle due to the absence of blankholding force constraint [11], as shown in Fig. 1 [12]. The modelling of flange wrinkling phenomena in deep drawing processes has been reviewed in detail by Kadkhodayan et al. [13]. To predict the induced buckling phenomena experienced when using macro-scale textured tools, Zheng et al. [12] estabilshed analytical buckling models based on Senior's beam [14] and Yu's plate structure buckling assumptions [15] for cylindrical deep drawing tests. However, the assumption of reduced Young's modulus which is used to reflect inelastic buckling is believed to be inaccurate [13]. In addition, the beam and plate assumptions may not be suitable for tool texture designs with larger dimensions.
To address the above disadvantages, the work in this paper aims to establish an advanced buckling model based on classical elastic-plastic mechanics of solids and futher extended to large surface textured tool designs. The use of constitutive relationship of an elastic-plastic solid can improve the calculation accuracy to cover the deformation behaviour at the plastic range. In this study, a more advanced analytical buckling model specifically for the metal forming processes using macro-textured tool design was developed based on elastic-plastic solid mechanics, bifurcation theory [16] and DMV shell structure theory [17]. The accuracy of the developed buckling model was validated by experimentation and FE simulations implemented with a dislocation-driven material model. Effects of related analytical assumptions and solutions, such as yield criteria and boundary condition, are discussed in detail. In addition, application ranges of the developed buckling model were also investigated.
The model can provide a prediction tool for stamping tool designers. The combination of a series of classic and robust theories provides a feasible approach to modelling actual forming processes, and can be a foundation for establishing buckling models of metal forming processes in hot forming conditions in the future.

Macro-textured tool design concept and modelling
Forming dies of sheet stamping processes are normally composed of the punch, die and blankholder. Fig. 2a shows the schematic diagram of a tool set with textured tool surfaces.
According to the contact condition between textured tools and sheet blank, textures can be divided into a groove (texture hollow) and surface (texture top). Correspondingly, the blank material in contact with the surface was defined as Zone S, while material located on the texture hollow (groove) was defined as Zone G. For sheet drawing processes with a twodimensional stress state generated within the flange material, flange material normally experiences a tensile stress in the material drawn-in direction and a compressive stress in the circumferential direction. Hence, Zone G material located between grooves of tools may buckle during draw-in as a result of the absence of blankholding force constraint, as shown in According to the structure stability theory, buckling can be avoided when the structure dimension is controlled within a certain magnitude. Hence, dimensions of texture features need to be defined. A cylindrical deep drawing design was used as the fundamental case study for the establishment of the buckling model. In this study, surface textures were only modelled on the blankholder, while the surfaces of punch and die remained flat. Fig. 3 shows the plane view of a quarter of the macro-textured blankholder. The narrow sectors indicate surface and wide sectors indicate groove. Since the length scale (radial direction) of the groove is parallel to material flow direction which has no significant effect on flange wrinkling, the length scale effect was ignored in this paper. The sizes of surface and groove of textures was defined by the subtended radial angles. Arc angles, and , of texture features were used respectively as shown in Fig. 3. The surface arc angle was fixed at 2.5°.
Then, a normalised ratio reflecting the texture feature variation, is defined as in Eq. (1). (1)

Fig. 3.
A plane view drawing of macro-textured blankholder used for cylindrical deep drawing. 6

Constitutive equations of elastic-plastic solids
To model the buckling phenomenon of Zone G flange material which is plastically deformed at room temperature, constitutive relationships of elastic-plastic solids needs to be obtained.
The incremental total strain can be expressed as the sum of elastic and plastic strain in a tensor expression as in Eq. (2) [18]. ( where represents the total strain increment, and represent elastic strain increment and plastic strain increment respectively.
The elastic strain increment is given in the generalised Hooke's law, described by Eq.
The plastic stress increment is given by Hill's theory [19] using the plastic potential .
In terms of the plastic loading, the plastic stress increment is given in Eq. (4).
where is a scalar factor associated with deformation history, material stress state, and work hardening given as Eq. (5) [20]. (5) where in this equation, is tangent Young's modulus. The magnitude of equals to the instantaneous slope of the effective stress-strain curve given in Eq. (6).
Where is the effective stress and is the effective strain.
Substituting Eqs. (3)-(6) into Eq. (2), the constitutive relationship between the stress and strain increments in the plastic loading is fully defined and expressed in Eq. (7).

(7)
With regard to sheet stamping processes, the plane stress state normally applies, . By applying the principal stress state, an explicit expression of Eq. (7) can be obtained for a particular plastic potential and is given in Eq. (8).
In this paper, two classical yield criteria applicable for isotropic materials, the von Mises and Tresca yield criteria, are considered due to the isotropic nature of the material resulting from the annealing heat treatment performed for the test-piece material. In addition, the selection of isotropic yield criteria can simplify the calculation significantly.

The von Mises yield criterion
The classical von Mises yield criterion yield function is given in Eq. (9).
where for the principal stress states.
For the plane stress state cases, , Eq. (9) can be further modified to Eq. (10).
where , are the deviatoric stress and .

The Tresca yield criterion
Another classical isotropic yield criterion, for the plane stress state , is that of Tresca yield criterion expressed by Eq. (16). It should be noted that the subscripts of the stress terms are arbitrarily defined instead of corresponding to the principal stress sequence. The sequence of principal stress is determined according to the specific metal forming process.

Material modelling in cold forming condition
To obtain instant magnitudes of the above stiffness moduli and material flow stress at different stages of draw-in, a material model of a selected material is required. In this paper, a physical-based dislocation-driven material model is applied [21]. The list of equations is given below from Eqs. (18)- (23).
Eq. (18) is the traditional power-law viscoplastic flow formulation considering damage effect on viscoplastic flow. represents the plastic strain rate, is the initial yield stress, is the hardening stress. and is the factor reflecting damage. Eq. (19) represents the evolution of material hardening , which is governed by the evolution of normalised dislocation density . is equal to, , is the initial material dislocation density, and is the maximum dislocation density. Eq. (21) represents damage evolution for the uniaxial formulation, which is a function of plastic strain rate . Eq. (22) is Hooke's law for a simple uniaxial state. Details of the illustation of each equation is given by Lin et al. [22]. The material constants of this dislocation-drvien material model can be deterimined using uniaxial tensile experiments of a particular alloy. Then, such a material model can be used for both developing the buckling model and numerical simulations via user-defined subroutine. 10

Bifurcation function
To overcome the inaccuracy of reduced modulus used in beam and plate buckling models, the bifurcation function proposed by Hutchinson et al. [16] was utilised. In addition, Zone G material can be considered as a quasi-shallow shell, as defined in DMV shell theory. This shell is thought to be accurate when the texture hollow, groove, is sufficiently large, while the assumption may not be accurate for small texture designs which needs to be further validated using experiments. Hence, the accurate application range of the DMV shell theory assumption needs to be validated experimentally. Before experiencing deformation, the assumed shell is believed to remain as a relatively flat geometry, the displacement components are varying functions of the coordinates. The bifurcation function described in Eq. (23) is utilised here to evaluate buckling occurrence for the Zone G material [16]. In this paper, as the flange material experiences plastic deformation when drawn in, the wrinkling of the Zone G material is assumed to be plastic buckling.
where is the incremental stiffness matrix related to stress increments as mentioned above, is the incremental bending strain, is the incremental stretching strain, is the inplane principal stresses, is the displacement increment of buckling with a direction normal to the shell middle surface and denotes the area of the shell middle surface over which wrinkles occur. The bifurcation function can be divided into three terms for simplicity, where represents the bending strain energy term, represents the membrane stress term and represents the potential energy of the edge stress. In terms of this function, when for all admissible fields, bifurcation will not occur and buckling cannot happen. While for , bifurcation becomes possible indicating experiencing buckling for some non-zero fields [16].

Incremental strain and stress
According to the basic relationships of the DMV shell theory, incremental bending and stretching strains due to buckling from the uniform membrane state can be calculated using Eqs. (24) and (25) [13]: where is the mathematical expression of deflection and , represents the in-plane displacement fields.
Considering the geometry feature of cylindrical cup deep drawing process, a cylindrical coordinate system is established. Then the displacement fields become functions of instant radius and polar angle , which are , and for this cylindrical coordinate system. Substituting the displacement fields into Eqs. (24) and (25), the detailed expressions of incremental strains are given in Eqs. (26) and (27).
The incremental stress resultants and bending moments at buckling can be calculated using the above incremental strains. These resultants are given by Eqs. (28) and (29).
In terms of the in-plane stress components, the tensile and compressive stresses of the flange in cylindrical deep drawing process are normally expressed in Eqs. (30) and (31) respectively when neglecting the friction term due to the non-blankholding force constraint of Zone G [23].
where is the material flow stress, is outside flange radius and is the instant flange radius for a selected infinitesimal small unit. 12

Boundary conditions
To obtain analytical expressions for the above incremental strains, mathematical functions using different boundary condition of the displacement fields are required. Selecting a vertical cross-section view of one unit (two surfaces and one groove) for tool textures, Fig. 4 shows schematic illustrations of (a) the hinged and (b) the built-in boundary conditions [24].
Each Zone G material clamped by tools located above groove and between two adjacent surfaces can be assumed as a shell sector with constraints at two ends. Zheng et al. [12] has established the mathematical expressions of deflections for Zone G material based on the above boundary conditions utilising a geometry assumption of two-dimensional plate, a detailed discussion and demonstration of the accuracy of the expressions were also described.
The built-in boundary condition represents a fully clamped tool constraint at an extremely high blankholding force, however, this kind of blankholding force may increase the radial stress severely, and the material cannot be drawn in. Regarding the hinged boundary condition, sheet metal can be drawn into the die with no displacement in the vertical direction, and , which means that wrinkling will not occur for Zone As can be seen from the expressions for and , there exists a component which represents the induced shear stress component resulting from deformation. In deep drawing processes, and can be regarded as the principal stresses corresponding to and mentioned above. Chu and Xu [25] has concluded that this shear stress component is vitally important for predicting wrinkling occurrence. In this paper, it is assumed that to conduct the calculation [13].

Calculation procedure
To capture the instant forming stage when flange wrinkling occurs, a relationship between the actual deep drawing process and developed buckling model needs to established. In this study, a simple method [26] using the relationship between the hoop strain of flange material (actual forming process) and the corresponding stiffness matrix (buckling model) was applied, as shown in Fig. 5a. In this method, the outside edge of Zone G material can be assumed as a uniaxially compressed slab ( ) according to Eq. (30). Hence, the deformation type of this slab is assumed to be identical to a uniaxial compression test. In addition, the arc shape of material is assumed to remain unchanged during draw-in. When the flange material is drawn from the initial stage to a certain stage, the instant diameter reduction, , defined in Eq. (38) corresponds to an effective strain of uniaxial stress-strain curve, as described in Eq.
(39). Finally, the above relationship can be used to capture the onset of flange wrinkling during material draw-in. (1) Input a radius reduction , calculating corresponding "effective strain" using Eqs.
(2) Obtain the instant tangent modulus using dislocation-driven material model in Eq.
(40), the plane stress state and the corresponding displacement fields. 15 (3) Obtain the instant stiffness matrix and incremental strains and according to step (2) results.
(4) Substitute stiffness matrix and incremental strains in the bifurcation function. If , wrinkling occurs and output the critical radius value; If , return to step (2).
(40) Tresca yield criterion can also be established. Table 1 lists the process parameters needed for the model calculation.

Quasi-static uniaxial tensile test
Quasi-static uniaxial tensile tests were performed to obtain mechanical properties of testpiece material, AA6082. The as-received material was in T6 condition with a thickness of 1.5 mm, and the chemical composition is given in

Deep drawing test
Cylindrical deep drawing test results as in the author's previous publication [12] were used here to validate the developed buckling model, as shown in Fig. 6a. Fig. 6b shows the necessary tool dimensions for the calculation of the analytical buckling model. The test-piece material used throughout the deep drawing tests was the same as the uniaxial tensile tests.
Texture features were only machined on the blankholders. Different texture ratios were used, , for the deep drawing tests. The blankholding forces used were 10 kN and 40 kN, and the forming speed was 75 mm/s. All the tests were performed in cold stamping condition with the use of Omega-35 lubricant.

FE simulation programme
The minimal texture ratio reflecting the largest groove in experimentations was 1/6 of which FE simulations were performed. Additional simulations were computed for existing and theoretical texture ratios, utilising Pamstamp. Fig. 7 shows the FE model of cylindrical deep drawing, as well as the loading and boundary conditions.
The dimensions of tools were consistent with experimental set-ups. The mesh of blank is extremely important for guaranteeing the accuracy of simulation [28]. In this study, the mesh of test-piece material was selected as shell element, Belytschko-Tsay, with a fixed size at 4 mm. For different draw ratios, the number of test-piece elements was 3456, 4012 and 4320 respectively which were sufficiently large to guarantee simulation accuracy. The mesh size 18 was fine enough for utilising subroutine in Pamstamp software as numerically determined.
Five integration points at the mid-plane of the element was used. The tool mesh was default solid elements and treated as rigid bodies. Table 3 lists the computation conditions for FE simulations. All the simulations were performed in cold stamping condition.

Material modelling
The developed dislocation-driven material model of AA6082 in Section 2.3 was validated by the uniaxial tensile test result in Section 3.1. The test-piece material was considered as isotropic after annealing treatment, and only mechanical data (stress-strain) in the rolling direction was used. Fig. 8a shows the comparisons between experimental results and material model predictions. The symbols represent experimental data of the uniaxial tensile test, and the solid lines represent fittings of material model using three different strain rates. As can be seen in this figure, the level of flow stress can be regarded to be independent of strain rate, 19 which is consistent with elastic-plastic mechanism of aluminium alloy in cold stamping condition. A good agreement was obtained indicating that the developed material is capable of reflecting material deformation in deep drawing process. Table 4 list the material constants of the dislocation-driven material model. The material model was implemented into the FE simulations via user-defined subroutine. Typical computed contours were shown in Fig. 8b.

FE model verification
where is the measured or predicted thickness at a particular position, is the initial testpiece thickness.   Fig. 10 gives experimental and numerical definitions of the onset of Zone G wrinkling. The onset of wrinkling is defined when the height between the peak and lowest points of wrinkles reaches 1 mm. However, the occurrence of wrinkling is difficult to be captured precisely considering the stroke fluctuations of the hydraulic press ram and uncertainties in the measurement. However, such a limitation can be overcome using a numerical approach, as shown in Fig. 10b. In this figure, indicates the instantaneous flange radius when the testpiece experienced flange wrinkling.

Boundary condition assumption
The mathematical expressions of wrinkles which determine the incremental strain tensors, are different using different boundary condition assumptions, in Eqs. (32) and (33). Fig. 11 shows force is unrealistic in actual forming processes, as it might cause fracture due to induced high radial tensile stresses.
23 Fig. 11. Effect of boundary condition assumption on the buckling model prediction accuracy using BHF 10 kN and 40 kN at a forming speed of 75 mm/s.

Yield function effect
The stiffness matrix used in the bifurcation function depends on the chosen yield criterion. . For the Tresca yield function in Fig. 12b, there is some difference between the model calculations and experimentations for a texture design: , while closer agreement can be found in Fig. 12a using the von Mises yield criterion. In general, the analytical buckling model using Tresca yield criterion has a higher prediction than using the von Mises yield criterion. While when the texture ratio becomes greater, both yield criteria have good capabilities to predict the onset of wrinkling within reasonable errors. Tresca for texture ratios: and at a blankholding force 10 kN and forming speed of 75 mm/s. 25

Application range of buckling model
As discussed above, the analytical buckling model was established on the assumption of DMV shell theory of Zone G material, which is believed to be more accurate for large texture features. The geometry of Zone G material was determined by texture ratio and test-piece draw ratio. Fig. 13a shows the wrinkling occurrence of Zone G material for different draw ratios using different tool texture ratio. . However, the previous models were believed to be inaccurate in reflecting the geometry effects of large texture features. Another important factor affecting wrinkling is the test-piece draw ratio. Fig. 13b shows the wrinkling occurrence of Zone G material using different test-piece draw ratio. The solid lines represent results from the analytical buckling model predicted results using the hinged boundary condition and the von Mises yield criterion, while the square symbols represent experimental results using a tool texture ratio of at a blankholding force of 10 kN. The forming speed used was 75 mm/s. The circular symbols represent further validated numerical results using a series of texture ratios to minimize the experimental tests. With increasing draw ratio, the critical reduction of test-piece diameter decreases, which means that, at a given texture ratio, the deep drawing process with a large draw ratio is easier for wrinkling to occur. The reason is similar to that of texture ratio effect, that is, the wider unsupported flange arising with larger diameter blanks is less resistant to buckling. Most of the numerical results are higher 26 than analytical results for the three large tool textures. One possible reason for this may be due to the wrinkling determination in the simulation. Wrinkling occurrence is defined when the largest height difference reaches 1 mm. However, the wrinkling occurrence detected in the analytical buckling model may be earlier than that defined in the numerical simulation.
Based on above discussions, the application range of developed buckling model is for smaller texture ratios at least lower than 1/6 and larger draw ratios.
27 Fig. 13. Application ranges of buckling models: (a) Different texture ratio and (b) Different draw ratio.

Conclusions
In this paper, an advanced buckling model that can be used to predict flange wrinkling behaviour in sheet stamping processes with the utilisation of textured tool surfaces was established. A series of classic mechanics theories, such as bifurcation and DMV theories, are integrated in this buckling model with advanced dislocation-driven based material model.
Experimental and numerical approaches were used to validate the accuracy of the proposed buckling model by evaluating the wrinkling occurrence using a cylindrical deep drawing tests utilising textured tool surfaces on the blankholders. For the buckling model accuracy, the effects of different boundary condition assumption due to the tool surface constraint and yield function selection on the wrinkling occurrence were investigated. It can be concluded that the hinged boundary condition reflects the actual tool surface constraint more accurately.
Moreover, using the von Mises yield criterion in the analytical buckling model exhibits better agreement with experimental results than using the Tresca yield criterion. It is also possible that the developed buckling model can be improved by introducing anisotropy yield criteria in the future. 28 Corresponding experiments were performed using different tool designs. The buckling model based on DMV shell theory can accurately predict the smaller texture ratio designs as well as greater draw ratios. While the DMV shell theory was not suitable when the tool texture ratio was increased to a lager magnitude.