Flange Wrinkling in Deep-Drawing: Experiments, Simulations and a Reduced-Order Model †

: Flange wrinkling is often seen in deep-drawing process when the applied blankholding force is too small. This paper investigates the plastic wrinkling of ﬂange under a constant blankholding force. A series of deep-drawing experiments of AA1100-O blanks are conducted with different blankholding forces. The critical cup height and wrinkling wave numbers for each case is established. A reduced-order model of ﬂange wrinkling is developed using the energy method, which is imple-mented to predict the ﬂange wrinkling of AA1100-O sheet by incrementally updating the ﬂange geometry and material hardening parameters during the drawing process. A deep-drawing ﬁnite element model is developed in ABAQUS/standard to simulate the ﬂange wrinkling of AA1100-O blanks under constant blankholding force. The predicted cup height and wave numbers from the ﬁnite element model and reduced-order model are compared with the experimental results, which demonstrates the accuracy of the reduced-order model, and its potential application in fast prediction of wrinkling in deep-drawing process.


Introduction
Reliable prediction of failure such as wrinkling, earing and tearing is necessary for the virtual design of sheet metal forming processes such as deep-drawing and stamping. To handle the complex material and geometric features of a real component, analytical methods have given rise to numerical methodologies, such as finite element (FE) analysis, which can provide helpful information in a less expensive way than trial-and-error with physical prototypes. However, simulations of the manufacturing processes can be quite time-consuming, particularly when a quick evaluation is needed. Even worse, such a modeling approach is oftentimes incompatible with the needs of real-time process control, including the recent efforts on smart and autonomous manufacturing, e.g., Industry 4.0 [1]. Feedback process control requires sensing, a fast and efficient evaluation, decision-making and actuation, oftentimes while the part is being processed. In the case of drawing and stamping, this processing cycle is typically a fraction of a minute. This in turn leads to the development and implementation of reduced-order models with high efficiency and reasonable accuracy. On the other hand, higher fidelity numerical models can be used for process controller design and selection [2].
Wrinkling and tearing are two of the main failure modes commonly encountered in deep-drawing. As the material of the flange is drawn towards the die cavity, compressive hoop stresses develop. These will lead to wrinkling of the thin annular sheet (the flange of the blank). To prevent this, the flange is "held-down" by the blankholder. When the applied blankholding force (BHF) is insufficiently small, the flange may still wrinkle at Many of the analytical plastic wrinkling models were developed based on a fixed flange ratio (i.e., the ratio of inner to outer diameter of the flange) and a simple stress state in the flange, estimated from perfectly plastic Tresca model. It is important to note that the flange ratio changes along the drawing process, and accordingly so does the stress state. Therefore, estimation of the stress state based on a more realistic hardening model and updating the stress state along the drawing process can be crucial to the wrinkling prediction. Furthermore, most of the analytical models were developed based on springloaded blankholders. For typical pneumatic or hydraulic-type blankholders and cushions, a constant BHF is often exerted to the flange, which must be appropriately taken into account in the modeling. The fact that the restraint remains constant mathematically complicates the stability problem to be solved; recently, an analytical treatment of the elastic wrinkling problem has been proposed by the authors, considering the post-wrinkling equilibrium of the sheet [14,15]. This paper investigates the plastic wrinkling of the flange under a constant BHF during the drawing process. Section 2 presents results of the deep-drawing experiments conducted under a range of BHFs. Section 3 develops a 2-D analytical, or reduced-order wrinkling model based on energy method. The model is then used to predict the wrinkling of Plastic wrinkling of the flange in relation to cup drawing process has been extensively studied analytically and numerically. Geckeler [3] developed a simple 1-D wrinkling model regarding the critical wrinkling stress and wave number for no blankholder case, which was extended by Senior [4] to a more general case accounting for the effect of a blankholder, using the energy method. Yu and Johnson [5] extended the 1-D model to 2-D case for both elastic and plastic wrinkling of the flange based again on an energy criterion. Cao and Boyce [6] proposed an energy method to determine both elastic and plastic wrinkling of the flange by incorporating numerical and analytical work. An alternative method for plastic wrinkling analysis is using the bifurcation functional [7,8]. This method has been widely adopted in the analysis of structural instabilities [9][10][11]. Based on this idea, Triantafyllidis and Needleman [12] studied the plastic wrinkling of a thin sheet with normal anisotropy. Following the same idea, Chu and Xu [13] investigated the flange wrinkling of a thin sheet with more complicated anisotropy.
Many of the analytical plastic wrinkling models were developed based on a fixed flange ratio (i.e., the ratio of inner to outer diameter of the flange) and a simple stress state in the flange, estimated from perfectly plastic Tresca model. It is important to note that the flange ratio changes along the drawing process, and accordingly so does the stress state. Therefore, estimation of the stress state based on a more realistic hardening model and updating the stress state along the drawing process can be crucial to the wrinkling prediction. Furthermore, most of the analytical models were developed based on springloaded blankholders. For typical pneumatic or hydraulic-type blankholders and cushions, a constant BHF is often exerted to the flange, which must be appropriately taken into account in the modeling. The fact that the restraint remains constant mathematically complicates the stability problem to be solved; recently, an analytical treatment of the elastic wrinkling problem has been proposed by the authors, considering the post-wrinkling equilibrium of the sheet [14,15]. This paper investigates the plastic wrinkling of the flange under a constant BHF during the drawing process. Section 2 presents results of the deep-drawing experiments conducted under a range of BHFs. Section 3 develops a 2-D analytical, or reduced-order wrinkling model based on energy method. The model is then used to predict the wrinkling of AA1100-O sheet in our deep-drawing tooling. To verify the reduced-order model, Section 4 presents the FE simulations of the flange wrinkling under different BHFs. Predictions from the reduced-order model and the simulation results are compared with the experiments. This article is a revised and expanded version of a paper entitled: "Industry 4.0 in stamping: A wrinkling indicator for reduced-order modeling of deep-drawing processes", which was presented at the 30th International Conference on Flexible Automation and Intelligent Manufacturing (FAIM 2020-2021), Athens, Greece [16]. In comparison to our previous work that focused on wrinkling theory and simulations, this paper includes a series of wrinkling experiments of AA1100-O blanks, which serves to validate the proposed theory. The modelling work is also revisited and updated to reproduce the measured punch force-displacement responses and the wrinkling and ironing phenomena.

Material and Tooling
The material of this study is commercially-pure aluminum, AA1100-O, of 0.51 mm thickness. The sheets are received in the -H24 temper, cold-rolled. They are then heattreated at 343 • C for 5400 s (1.5 h), followed by air-cooling. The mechanical and forming properties of this material have been extensively characterized in earlier works [17,18]. Here, circular AA1100-O blanks of 35 diameter and 0.51 mm thick are produced via waterjetting.
To probe the behavior of this sheet against wrinkling in forming, a set of deep-drawing experiments are conducted with AA1100-O blanks under a wide range of blank holding forces (BHF). The circular cup drawing tooling used in this work is modular. The interchangeable parts consist of a die insert, a flat-ended punch and a blankholder inset, all made of A2 tool steel, heat-treated to 44 (punch) and 58 HRC (die), and with a surface finish Ra better than 0.8 µm. The tooling is integrated with an Instron 8872 servohydraulic testing machine of 25 kN capacity, as shown in Figure 2.
AA1100-O sheet in our deep-drawing tooling. To verify the reduced-order model, Section 4 presents the FE simulations of the flange wrinkling under different BHFs. Predictions from the reduced-order model and the simulation results are compared with the experiments.
This article is a revised and expanded version of a paper entitled: "Industry 4.0 in stamping: A wrinkling indicator for reduced-order modeling of deep-drawing processes", which was presented at the 30th International Conference on Flexible Automation and Intelligent Manufacturing (FAIM 2020-2021), Athens, Greece [16]. In comparison to our previous work that focused on wrinkling theory and simulations, this paper includes a series of wrinkling experiments of AA1100-O blanks, which serves to validate the proposed theory. The modelling work is also revisited and updated to reproduce the measured punch force-displacement responses and the wrinkling and ironing phenomena.

Material and Tooling
The material of this study is commercially-pure aluminum, AA1100-O, of 0.51 mm thickness. The sheets are received in the -H24 temper, cold-rolled. They are then heattreated at 343 °C for 5400 s (1.5 h), followed by air-cooling. The mechanical and forming properties of this material have been extensively characterized in earlier works [17,18]. Here, circular AA1100-O blanks of 35 diameter and 0.51 mm thick are produced via waterjetting.
To probe the behavior of this sheet against wrinkling in forming, a set of deep-drawing experiments are conducted with AA1100-O blanks under a wide range of blank holding forces (BHF). The circular cup drawing tooling used in this work is modular. The interchangeable parts consist of a die insert, a flat-ended punch and a blankholder inset, all made of A2 tool steel, heat-treated to 44 (punch) and 58 HRC (die), and with a surface finish Ra better than 0.8 μm. The tooling is integrated with an Instron 8872 servohydraulic testing machine of 25 kN capacity, as shown in Figure 2. A solid model of the tooling is shown in Figure 3a, and the dimensions of the tooling are shown in Figure 3b. Visible in Figure 3a are the three pneumatic actuators that are used for applying the constant BHF; however, in this study, due to the very low forces involved and to eliminate any influence of friction inside the actuators, the latter were A solid model of the tooling is shown in Figure 3a, and the dimensions of the tooling are shown in Figure 3b. Visible in Figure 3a are the three pneumatic actuators that are used for applying the constant BHF; however, in this study, due to the very low forces involved and to eliminate any influence of friction inside the actuators, the latter were replaced by a series of steel weights, added evenly on top of the BH (self-weight of the BH is 42 N). In this way, precise and unambiguous small forces could be applied to the blank during drawing. The punch-die clearance is about 1.27 times the sheet thickness [19].
replaced by a series of steel weights, added evenly on top of the BH (self-weight of the BH is 42 N). In this way, precise and unambiguous small forces could be applied to the blank during drawing. The punch-die clearance is about 1.27 times the sheet thickness [19].  For each experiment, the tooling and the blank are lubricated with Multidraw PL 61 SE by Zeller+Gmelin. The blanks are centered in the tooling with a centering ring. Proper alignment is ensured before each experiment. The drawing experiments are conducted under displacement control, with the punch velocity set to 0.1 mm/s. The punch force and displacement are automatically recorded during the experiments. Figure 4 shows the drawn cups, which correspond to a series of BHF ranging from zero to 176 N. Note that since the self-weight of the BH is 42 N, the only one smaller BHF that can be achieved with our tooling is zero. The flange of that cup wrinkles almost immediately after forming commences, as will be detailed later. As the BHF increases, the flange wrinkling is delayed, with more wrinkling waves appearing, and of smaller wrinkling amplitudes. When the BHF exceeds a certain value, the wrinkling can be fully suppressed (e.g., the last 2 cups in Figure 4). The average cup heights at the onset of wrinkling and the observed wrinkling waves corresponding to each BHF case are presented in Table  1.  For each experiment, the tooling and the blank are lubricated with Multidraw PL 61 SE by Zeller+Gmelin. The blanks are centered in the tooling with a centering ring. Proper alignment is ensured before each experiment. The drawing experiments are conducted under displacement control, with the punch velocity set to 0.1 mm/s. The punch force and displacement are automatically recorded during the experiments. Figure 4 shows the drawn cups, which correspond to a series of BHF ranging from zero to 176 N. Note that since the self-weight of the BH is 42 N, the only one smaller BHF that can be achieved with our tooling is zero. The flange of that cup wrinkles almost immediately after forming commences, as will be detailed later. As the BHF increases, the flange wrinkling is delayed, with more wrinkling waves appearing, and of smaller wrinkling amplitudes. When the BHF exceeds a certain value, the wrinkling can be fully suppressed (e.g., the last 2 cups in Figure 4). The average cup heights at the onset of wrinkling and the observed wrinkling waves corresponding to each BHF case are presented in Table 1.

Punch Force-Displacement
The punch force-displacement curves of all the experiments are shown in Figure 5. They all show a consistent ascending slope at the beginning of the process. Depending on the BHF, the responses show some difference. For the fully drawn cups, e.g., the 176 N  The punch force-displacement curves of all the experiments are shown in Figure 5. They all show a consistent ascending slope at the beginning of the process. Depending on the BHF, the responses show some difference. For the fully drawn cups, e.g., the 176 N case, the punch force develops a maximum and then drops to nearly zero, due to the competition between the diminishing flange (geometric softening) and the material work-hardening [20]. For the wrinkling cases, e.g., the 91 N case, the punch force develops a lower maximum force (compared to the fully drawn ones). As the punch force drops to some level, it starts to increase drastically again, due to ironing between the punch and the die of the wrinkles that have already formed in the flange. For the zero N case, the wrinkling occurs as the punch force is increasing. To avoid damage to the tooling (e.g., galling), the wrinkling experiments were terminated before the wrinkles were fully drawn through the die, see However, the sudden increase of the punch force is an indication of the onset of wrinkle ironing, which corresponds to the moment when the radial wave front reaches the inner die radius. If the wrinkling is assumed to occur only when the blank is in contact with the BH and the radial wave front is at the outer radius of the die (see Figure 2, otherwise a fully drawn cup without wrinkling may be impossible for this tooling), the draw depths required to move the radial wave front from the outer die radius to inner die radius over the die fillet can be estimated as h = 3.57 mm (this value depends on the dimensions of the tooling, see Appendix A for detailed derivation). The draw depth of the onset of flange wrinkling can be estimated by subtracting 3.57 mm from that at the onset of ironing of the wrinkles. This estimation method allows the instantaneous identification of the onset of wrinkling during the experiments, see Appendix A for more discussion.

Reduced-Order Modeling of Flange Wrinkling
In this section, an analytical model of the occurrence of wrinkles in the flange of a drawn cup is proposed. In comparison to a full-scale numerical simulation, the order of this model is much lower, and hence it can be profitably used for process design and control.
Consider an annular flange with an inner radius of a and an initial outer radius of R , as shown in the schematic of the circular cup drawing tooling used in this analysis in Figure 6. A constant BHF is applied to the flange during the drawing. Estimation of the pre-wrinkling stress and strain states in the flange is crucial for the wrinkling analysis. This section presents the approximation of the stress and strain states in the flange. As the

Identifying the Onset of Wrinkling from Force-Displacement Curves
Establishing the onset of wrinkling from experiments is quite challenging, because it is difficult to capture the deformation of the flange, which is covered by the BH during the drawing process. Even if the motion of the BH was available, identifying the onset of wrinkling would be difficult, because the motion of the BH does not show a sudden increase but a gradual change. Therefore, it is necessary to establish an alternative way of identifying the onset of flange wrinkling. We chose to do this from the punch forcedisplacement curve, Figure 5. For the low BHF cases, such as the zero N or 42 N, there is an obvious deviation of the responses between the wrinkled and the fully drawn case, which is due to the onset of flange wrinkling. The instant that deviation occurs can be approximately considered as the onset of flange wrinkling. For the medium or relatively high BHF cases that still show wrinkling, e.g., 91 N or 111 N, the onset of wrinkling may occur even after the force maximum. Therefore, the deviation of the F-d response from the fully drawn case, and consequently the different force maximum, may be due to the difference between BHF's; the onset of wrinkling is thus more difficult to identify. However, the sudden increase of the punch force is an indication of the onset of wrinkle ironing, which corresponds to the moment when the radial wave front reaches the inner die radius. If the wrinkling is assumed to occur only when the blank is in contact with the BH and the radial wave front is at the outer radius of the die (see Figure 2, otherwise a fully drawn cup without wrinkling may be impossible for this tooling), the draw depths required to move the radial wave front from the outer die radius to inner die radius over the die fillet can be estimated as h = 3.57 mm (this value depends on the dimensions of the tooling, see Appendix A for detailed derivation). The draw depth of the onset of flange wrinkling can be estimated by subtracting 3.57 mm from that at the onset of ironing of the wrinkles. This estimation method allows the instantaneous identification of the onset of wrinkling during the experiments, see Appendix A for more discussion.

Reduced-Order Modeling of Flange Wrinkling
In this section, an analytical model of the occurrence of wrinkles in the flange of a drawn cup is proposed. In comparison to a full-scale numerical simulation, the order of this model is much lower, and hence it can be profitably used for process design and control.
Consider an annular flange with an inner radius of a o and an initial outer radius of R b , as shown in the schematic of the circular cup drawing tooling used in this analysis in Figure 6. A constant BHF is applied to the flange during the drawing. Estimation of the pre-wrinkling stress and strain states in the flange is crucial for the wrinkling analysis. This section presents the approximation of the stress and strain states in the flange. As the flange is drawn-in, the inner radius remains the same while the outer radius decreases monotonically and eventually equals the inner radius a o . It is convenient to define the following geometric parameters (see Figures 3b and 6): • The initial inner-outer radius ratio of the flange The current inner-outer radius ratio of the flange ρ = a o /b. Note that ρ increases from ρ o to 1 for the whole flange drawing process.

•
The radial coordinate r normalized by the current outer radius s = r/b (note that ρ ≤ s ≤ 1).

Strain State
For circular cup drawing, if the thickness strain is assumed to be much smaller than the in-plane strains, the only nonzero displacement component in the flange is the radial displacement u , which furthermore depends on the radial position only. The radial and hoop strains can thus be expressed as: Neglecting the elastic deformation and thinning, and considering plastic incompressibility, we can obtain: The radial displacement can be then found as:

Strain State
For circular cup drawing, if the thickness strain is assumed to be much smaller than the in-plane strains, the only nonzero displacement component in the flange is the radial displacement u r , which furthermore depends on the radial position only. The radial and hoop strains can thus be expressed as: Neglecting the elastic deformation and thinning, and considering plastic incompressibility, we can obtain: du r dr The radial displacement can be then found as: where C is a constant to be determined. Substituting (3) into (1), the strain components can be obtained as: Using the boundary condition at the outer periphery of the flange, we have: where b and R b are the current and initial (i.e., blank) outer radius, respectively. Then C can be determined and Equation (4) can be rewritten as: The normalized strain distribution are plotted in Figure 7a

Stress State
The stress state in the flange is obtained from the radial equilibrium condition. Again, assuming the thickness change in the flange can be neglected during the drawing process, we have: For Tresca yielding, from plastic work compatibility, the equivalent plastic strain can be found as: where ε e,b is the equivalent plastic strain at the outer periphery of the flange, which undergoes uniaxial compression in the hoop direction.

Stress State
The stress state in the flange is obtained from the radial equilibrium condition. Again, assuming the thickness change in the flange can be neglected during the drawing process, we have: ∂σ r ∂r For Tresca yielding, Equation (8) becomes: or using the normalized radial coordinate, s, where σ e is the equivalent stress. Neglecting the frictional forces from the die and blankholder and using the boundary condition at the outer periphery: The radial stress can be integrated as: The hoop stress can then be obtained as: It can be seen that the stresses depend on the equivalent stress-strain relationship of the material. For example, for the perfectly plastic model, the equivalent stress is a constant σ o , and the stress components can then be explicitly expressed as: For the power-law hardening model, i.e., σ e = K(ε e ) N , the stress state can be obtained as: where Y b is the uniaxial stress at the outer periphery, which is the same as the equivalent stress at that location.
It can be seen from Equations (14) and (15) that the stresses (normalized by the equivalent stress at the outer periphery) don't depend on the deformation of the flange for Tresca yielding and perfectly plastic or power-law hardening materials. For the equivalent stress-strain relationships that have more complicated expressions, numerical integration method can be used for Equation (12). The stress distributions in general depend on the deformation of the flange.
For our AA1100-O sheet, the equivalent stress-strain relationship can be fitted with a Voce hardening model: where the fitting parameters are shown in Table 2 [17]. Since it is difficult to obtain an explicit expression from Equation (12) with the Voce model, numerical integration method is adopted to obtain the stress distributions in the flange. Figure 8a, b compare the radial and hoop stresses (normalized by the uniaxial/equivalent stress at the outer periphery) vs. the normalized radial position for the perfectly plastic, power-law hardening and Voce hardening cases. Since the stress distribution from the Voce hardening model depends on the deformation of the flange, here the stress distribution when the equivalent strain at the outer periphery reaches 0.04 is shown as an example, which corresponds to a draw-in of 0.7 mm for the 35 mm blank used here. The radial stresses do not show significant difference for the three models. The hoop stress from the Voce model is the largest in amplitude.

Plastic Instability Using the Energy Method
The well-established energy method [21] is used here to analyze the wrinkling of the flange in the deep-drawing of AA1100-O sheet under a constant BHF. Assume a wrinkling mode in the following form:

Plastic Instability Using the Energy Method
The well-established energy method [21] is used here to analyze the wrinkling of the flange in the deep-drawing of AA1100-O sheet under a constant BHF. Assume a wrinkling mode in the following form: where A represents the normalized wrinkling amplitude and n the wrinkling wave number. For a given angle, the wrinkling amplitude is assumed to increase linearly from zero at the die entrance radius, a o (see Figure 3b), to maximum at the outer periphery, R b . Note the maximum amplitude of the wrinkling mode is where b is the current outer radius of the flange (see Figure 6). In this way, the wrinkling criterion is constantly being updated during drawing. The bending strain energy due to this wrinkling mode can be expressed as [22]: where D = E r t 3 9 and E r = is the Reduced Modulus [23], which is a function of the Young's Modulus and Tangent Modulus. Note that the Tangent Modulus varies radially through the flange because of the equivalent strain variation as shown in Equation (7), so that D also varies radially.
Substituting (17) into (18), the bending energy can be expressed as: where F = 1 ρ E r (s)(M 1 (s, n, ρ) + M 2 (s, n, ρ))ds and: If the drawing is performed under constant BHF Q (i.e., a constant force is applied to the flange during the process), the potential energy increase of the blankholder due to wrinkling is: The work carried out by the pre-wrinkling membrane stresses is: where σ r and σ θ are evaluated from Equations (12) and (13) based on the Voce model as shown in Equation (16). Numerical integration method is adopted to obtain the expressions for the stresses in the flange. The work carried out by the membrane stresses can be rewritten as: where H = −4 1 ρ σ r (s)3s + σ θ (s)n 2 s 1 − ρ s 2 ds.
Flange wrinkling is assumed to occur when the following condition is satisfied: Substituting Equations (19), (21) and (23) into (24), the wrinkling criterion can be expressed as: where E r,b and Y b are the Reduced Modulus and equivalent stress at the outer periphery, respectively, and ϕ = Qa o E r,b t 3 is the normalized BHF.

Model Implementation
Note that in Equation (25) the flange ratio ρ = a o /b monotonically increases during the drawing process. As a result, parameters such as the Reduced Modulus E r,b and equivalent stress Y b at the outer periphery also change accordingly. However, the maximum wrinkling amplitude γ max is unknown a priori, which may depend on the post-wrinkling analysis [14,15]. Here for simplicity γ max is assigned to be t (i.e., the sheet thickness), which is supposed to be a reasonable approximation.
To determine the onset of flange wrinkling from Equation (25), ρ is increased incrementally. For a specific ρ, the equivalent strain at the outer periphery is updated as ln (ρ/ρ o ); E rb and Y b are then updated using the material hardening model as shown in (16). Once ρ is known, the right-hand side (RHS) of Equation (25) is a function of the wave number n only. Maximizing the RHS with respect to an integer n gives the most possible wave number for a given ρ. With the RHS being a maximum for the wave number n, Equation (25) is checked again. If the left-hand side (LHS) of Equation (25) is smaller than or equal to the right-hand side, i.e., LHS ≤ RHS, wrinkling is assumed to occur, otherwise, ρ is increased and the updating process is repeated until the flange drawing process is finished, i.e., the flange is completely consumed. A flow chart showing this process is presented below in Figure 9. ln (ρ/ρ ); E and Y are then updated using the material hardening model as shown in (16). Once ρ is known, the right-hand side (RHS) of Equation (25) is a function of the wave number n only. Maximizing the RHS with respect to an integer n gives the most possible wave number for a given ρ. With the RHS being a maximum for the wave number n, Equation (25) is checked again. If the left-hand side (LHS) of Equation (25) is smaller than or equal to the right-hand side, i.e., LHS ≤ RHS, wrinkling is assumed to occur, otherwise, ρ is increased and the updating process is repeated until the flange drawing process is finished, i.e., the flange is completely consumed. A flow chart showing this process is presented below in Figure 9. The step-by-step updated results are shown in Figure 10. The LHS (dash line) and RHS (solid lines) are plotted against the draw-in (normalized by initial flange width, 4.85 mm). Results for three different BHFs of 40, 80, and 125 N are presented. The points where the solid and dashed line intersect indicate the onset of wrinkling. It can be seen that as the BHF is increased, the onset of wrinkling is delayed. When the BHF is greater than 125 N, intersection between the solid and dash line becomes impossible and thus no flange wrinkling is predicted.

Finite Element Model
A FE model is developed in ABAQUS/Standard (i.e., implicit) to simulate the flange wrinkling in deep-drawing. A schematic of the model is shown in Figure 11a, where the punch, die and BH are modelled as rigid bodies. The blank is modelled with shell elements (S4R) and one quarter of the mesh is shown in Figure 11b. A coarse mesh (680 elements) is adopted in the central part of the blank since the deformation is small and relatively uniform. A refined mesh of 19 (radial) × 216 (circumferential) elements is used in the annular region, which includes the wall and flange of the drawn cup (the projected initial positions of the punch and die are identified by the red and blue lines). Note that in the simulation a full blank is used because the wrinkles may not have quarter-symmetry. Contact between the blank, die and blankholder is modelled as exponentially soft contact. The soft contact is ensured to be sufficiently "hard" so that a further increase of the contact stiffness does not influence the simulation results significantly (zero pressure at clearance of 0.02 mm and 30 MPa at clearance of 0 mm). The BHF is kept constant during the drawing simulations. The material behaviour is modelled with von Mises yield function and isotropic hardening (Voce hardening model in Table 2).
(a) (b) Figure 10. Evolution of the LHS and RHS of Equation (25) with the normalized draw-in during the drawing process, for different BHFs.

Finite Element Model
A FE model is developed in ABAQUS/Standard (i.e., implicit) to simulate the flange wrinkling in deep-drawing. A schematic of the model is shown in Figure 11a, where the punch, die and BH are modelled as rigid bodies. The blank is modelled with shell elements (S4R) and one quarter of the mesh is shown in Figure 11b. A coarse mesh (680 elements) is adopted in the central part of the blank since the deformation is small and relatively uniform. A refined mesh of 19 (radial) × 216 (circumferential) elements is used in the annular region, which includes the wall and flange of the drawn cup (the projected initial positions of the punch and die are identified by the red and blue lines). Note that in the simulation a full blank is used because the wrinkles may not have quarter-symmetry. Contact between the blank, die and blankholder is modelled as exponentially soft contact. The soft contact is ensured to be sufficiently "hard" so that a further increase of the contact stiffness does not influence the simulation results significantly (zero pressure at clearance of 0.02 mm and 30 MPa at clearance of 0 mm). The BHF is kept constant during the drawing simulations. The material behaviour is modelled with von Mises yield function and isotropic hardening (Voce hardening model in Table 2).

Finite Element Model
A FE model is developed in ABAQUS/Standard (i.e., implicit) to simulate the flange wrinkling in deep-drawing. A schematic of the model is shown in Figure 11a, where the punch, die and BH are modelled as rigid bodies. The blank is modelled with shell elements (S4R) and one quarter of the mesh is shown in Figure 11b. A coarse mesh (680 elements) is adopted in the central part of the blank since the deformation is small and relatively uniform. A refined mesh of 19 (radial) × 216 (circumferential) elements is used in the annular region, which includes the wall and flange of the drawn cup (the projected initial positions of the punch and die are identified by the red and blue lines). Note that in the simulation a full blank is used because the wrinkles may not have quarter-symmetry. Contact between the blank, die and blankholder is modelled as exponentially soft contact. The soft contact is ensured to be sufficiently "hard" so that a further increase of the contact stiffness does not influence the simulation results significantly (zero pressure at clearance of 0.02 mm and 30 MPa at clearance of 0 mm). The BHF is kept constant during the drawing simulations. The material behaviour is modelled with von Mises yield function and isotropic hardening (Voce hardening model in Table 2). Friction plays a crucial role in the simulation of deep-drawing process, as it influences the global response and local thickness significantly [24]. Without a direct measurement of the friction coefficient in our deep-drawing process, we chose to identify the friction coefficient with an inverse method, i.e., by matching the measured and simulated punch force-displacement for the fully drawn cup. The friction coefficient that produced the best match is µ = 0.17. This friction coefficient is then used for simulations of all BHF cases.
A geometric imperfection is introduced to the blank, to trigger wrinkling of the flange. Specifically, an out-of-plane coordinate offset of the nodes in the undeformed flange region is introduced. The offset value is assumed to increase linearly, from zero at the die opening radius to maximum at the periphery of the blank along radial direction. Since the wrinkling wave number is not known a priori, the offset amplitude at the periphery of the blank is assigned a random value within ±1.5% of the initial thickness of the blank. The random distribution of the imperfection amplitude is similar to the one used in [18], see Figure 20 there, where the out-of-plane amplitude was amplified 30 times to better visualize the imperfection. Figure 12 compares the simulated and measured punch force-displacement responses for two different BHF cases. These results are typical of the remainder of the cases studied. The ascending part of the responses are all matched, which is mostly due to the appropriate friction coefficient adopted [18,[25][26][27][28]. The post-maximum force responses typically involve the onset of wrinkling and post-wrinkling behaviour. As seen in the figures, the postmaximum force responses and particularly the onset of ironing are well matched for these two cases. Similar agreement is observed for the rest of the cases, with the exception of 0 N and 42 N cases, for which the simulations overpredict the maximum forces by about 13%, and also overpredict the draw depth at the onset of wrinkle ironing by about 15%. ment of the friction coefficient in our deep-drawing process, we chose to identify the friction coefficient with an inverse method, i.e., by matching the measured and simulated punch force-displacement for the fully drawn cup. The friction coefficient that produced the best match is μ = 0.17. This friction coefficient is then used for simulations of all BHF cases.

Simulation Results
A geometric imperfection is introduced to the blank, to trigger wrinkling of the flange. Specifically, an out-of-plane coordinate offset of the nodes in the undeformed flange region is introduced. The offset value is assumed to increase linearly, from zero at the die opening radius to maximum at the periphery of the blank along radial direction. Since the wrinkling wave number is not known a priori, the offset amplitude at the periphery of the blank is assigned a random value within ±1.5% of the initial thickness of the blank. The random distribution of the imperfection amplitude is similar to the one used in [18], see Figure 20 there, where the out-of-plane amplitude was amplified 30 times to better visualize the imperfection. Figure 12 compares the simulated and measured punch force-displacement responses for two different BHF cases. These results are typical of the remainder of the cases studied. The ascending part of the responses are all matched, which is mostly due to the appropriate friction coefficient adopted [18,[25][26][27][28]. The post-maximum force responses typically involve the onset of wrinkling and post-wrinkling behaviour. As seen in the figures, the post-maximum force responses and particularly the onset of ironing are well matched for these two cases. Similar agreement is observed for the rest of the cases, with the exception of 0 N and 42 N cases, for which the simulations overpredict the maximum forces by about 13%, and also overpredict the draw depth at the onset of wrinkle ironing by about 15%. Figure 13 compares the cups from simulations and experiments for the 65 N and 111 N cases. They simulations reproduce the experimental observation that the number of waves increases with the BHF, e.g., see Figure 4. This agreement is accomplished despite the complication of the simulations of post-wrinkling, e.g., contact and shell element type issues. The contours shown in Figure 13 indicate the sheet thickness. While the bottom of the cup and the lower parts of the wall have thinned-down (initial thickness is 0.51 mm), the upper part of the cup, including the wrinkles has become thicker, by over 30% in some locations.  Figure 13 compares the cups from simulations and experiments for the 65 N and 111 N cases. They simulations reproduce the experimental observation that the number of waves increases with the BHF, e.g., see Figure 4. This agreement is accomplished despite the complication of the simulations of post-wrinkling, e.g., contact and shell element type issues. The contours shown in Figure 13 indicate the sheet thickness. While the bottom of the cup and the lower parts of the wall have thinned-down (initial thickness is 0.51 mm), the upper part of the cup, including the wrinkles has become thicker, by over 30% in some locations.

Determination of the Onset of Wrinkling in FE Analysis
Determining the onset of wrinkling in the simulations can be ambiguous due to the existence of imperfections. In other words, it is not easy to pinpoint the onset of wrinkling by monitoring the out-of-plane displacement of selected nodes in the flange throughout the simulation; such curves show a gradual increase. Alternatively, here for each BHF case, the zero-imperfection case is also simulated, which does not show wrinkling, as expected. The vertical displacements of the BH for the perfect and imperfect cases are extracted and compared, as demonstrated in Figure 14 for the 91 N case. The onset of wrinkling is identified as the "bifurcation" point for the two cases. This process is conducted for all BHF cases. Flange wrinkling is observed in the simulations for BHF smaller than 190 N. Then, the critical draw depth and corresponding wave numbers at the onset of wrinkling are extracted from the simulations. Figure 15 compares the simulated critical cup heights with the reduced-order predictions and experiments for different BHFs. The cup heights are Flange wrinkling is observed in the simulations for BHF smaller than 190 N. Then, the critical draw depth and corresponding wave numbers at the onset of wrinkling are extracted from the simulations. Figure 15 compares the simulated critical cup heights with the reduced-order predictions and experiments for different BHFs. The cup heights are normalized by the one of the fully drawn cup, which in our case is 11.2 mm using volume constancy and assuming no change in thickness during forming. They all show that as the BHF increases, the critical cup height increases. When the BHF is greater than 156 N, no wrinkling is detected from the experiments. The reduced-order model provides a critical BHF of about 125 N, beyond which no wrinkling is predicted. Below that force, both experiments and the reduced-order model match very well. The FE results show mixed agreement. While at low BHF the FE results match the experiments and the reduced-order model very well, at larger BHFs it predicts a slower growth of the critical draw depth with the BHF. The difference can be partially attributed to the different methods used to identify the onset of wrinkling in FE analysis and experiments, including factors such as the limitation of shell element FE model in capturing the onset of ironing. The FE model predicts a critical BHF of about 190 N. Interestingly, the empirical recommendation for BHF in our case (i.e., 1.5% the yield stress x the original flange area) is 207 N [20]. In that sense, both modelling approaches provide limits that are consistent with the recommended one; furthermore, the experiments do not refute that value, either. The FE results are sensitive to the imperfection amplitude used. The one selected for this work (within ±1.5% of the initial thickness of the blank) provides the best match with both the punch force-displacement and the wrinkling predictions.   Figure 16 compares the simulated wave numbers with the reduced-order predictions and experiments for different BHFs. They agree quite well when the BHF is below about 125 N, all showing increase of the wave number with the BHF. For BHF larger than 121 N, wrinkling did not occur all around the periphery of the cups, i.e., local wrinkling occurred on a certain portion of periphery, as can be seen in Figure 4. The total number of the wrinkles is estimated for the whole periphery based on the local wave number density. Therefore, there is some uncertainty of the estimated wave numbers for the relatively larger BHFs. Furthermore, for the FE results, as the wave number increases, the prediction may be limited by the discretization (i.e., the element size) adopted. A more refined mesh is required to capture small wavelengths.

Conclusions
Flange wrinkling is a commonly encountered failure mode in deep-drawing processes, which can be prevented by applying a sufficiently large BHF. This paper investigates the flange wrinkling through combined experimental, analytical, and numerical ef-

Conclusions
Flange wrinkling is a commonly encountered failure mode in deep-drawing processes, which can be prevented by applying a sufficiently large BHF. This paper investigates the flange wrinkling through combined experimental, analytical, and numerical efforts. • Both the critical cup height (i.e., when wrinkling commences) and wave number of the wrinkled cups increase with the BHF. • For relatively low BHFs, wrinkling occurs before the punch force reaches a maximum. The onset of wrinkling is considered as the bifurcation of the punch force-displacement response from that of the fully drawn cups. • For relatively high BHF, wrinkling occurs after the maximum of the punch force. A method is proposed to identify the onset of wrinkling from the punch forcedisplacement curve, which involves identifying the onset of ironing from the post-force maximum regime first and then subtracting the draw depth due to the draw-in of wrinkling wave front from the onset of wrinkling to ironing.

•
The reduced-order model involves incrementally updating the geometry of the flange, the stress/strain state in the flange and material hardening parameters during the drawing process.

•
The friction coefficient used in the FE model is identified from an inverse method for the fully drawn cups. The FE model also successfully reproduces the measured punch-force displacement response for the pre-wrinkling and post-wrinkling regimes, and the post-wrinkling configurations of the wrinkled cups.

•
The wave numbers from the experiments, reduced-order model and FE model are overall in good agreement. • For our AA1100-O blanks, the critical BHF that suppresses wrinkling is about 156 N.
The reduced-order model and FE model predict minimum BHF of 125 N and 190 N to prevent flange wrinkling, respectively. The commonly adopted empirical equation recommends a BHF of 207 N, which is found to be conservative.
As a closing remark, we want to state that one of the major objectives of this work is to establish an FEA framework that predicts the appearance of wrinkling; when exactly wrinkling occurs (i.e., at what draw depth) is important, but secondary. As we show in Figures 12 and 15, the imperfection used strikes a very good balance between reproducing the experimental punch force-displacement curves, and the range of BHF for which wrinkles appear. Furthermore, as shown in Figure 13, the imperfection used also predicts the wrinkled configuration very well.
Author Contributions: Conceptualization, K.C. and Y.P.K.; methodology, K.C. and Y.P.K.; software, K.C.; validation, K.C., A.J.C. and Y.P.K.; formal analysis, K.C.; resources, K.C.; data curation, A.J.C.; writing-original draft preparation, K.C.; writing-review and editing, K.C., A.J.C. and Y.P.K.; project administration, Y.P.K.; funding acquisition, Y.P.K. All authors have read and agreed to the published version of the manuscript. Assuming the cylindrical shell fully wraps the punch, then the critical draw depth can then be found as: Substituting our tooling geometry (R = 2 mm, R = 10.65 mm, R = 10 mm) leads to ∆h = 3.57 mm. If the difference of the punch and die radii is negligible, Equation (A3) can be simplified as: which turns out to be ∆h = 3.36 mm for our tooling. Note that Equation (A3) gives a slightly larger number than Equation (A4) and thus an earlier occurrence of the wrinkling, which serves as a more conservative estimate of the onset of wrinkling.
The key assumption of this approach lies in associating the appearance of the second rise in the punch force after the force maximum, with the wrinkling wave front arriving at the inner radius of the die, and thus the wrinkles beginning to be ironed. However, FEA shows that as the wave front travels over the die entrance radius, the wrinkles tend to be reduced or even disappear, due to bending. This somewhat limits the accuracy of this approach. If the difference of the punch and die radii is negligible, Equation (A3) can be simplified as: which turns out to be ∆h w = 3.36 mm for our tooling. Note that Equation (A3) gives a slightly larger number than Equation (A4) and thus an earlier occurrence of the wrinkling, which serves as a more conservative estimate of the onset of wrinkling.
The key assumption of this approach lies in associating the appearance of the second rise in the punch force after the force maximum, with the wrinkling wave front arriving at the inner radius of the die, and thus the wrinkles beginning to be ironed. However, FEA shows that as the wave front travels over the die entrance radius, the wrinkles tend to be reduced or even disappear, due to bending. This somewhat limits the accuracy of this approach.