Prediction of ductile damage and fracture in the single- and multi-stage incremental hole-flanging processes using a new damage accumulation law

Incremental hole-flanging (IHF) is a process in which a sheet with a pre-cut hole is flanged by the single point incremental forming (SPIF) process. Fracture prediction in IHF, such as SPIF, is associated with many challenges due to the deformation mechanisms. The purpose of this paper is to overcome the existing limitations and challenges and, thus, to predict accurately failure in single- and multi-stage IHF processes. To this end, the modified Mohr–Coulomb (MMC) criterion was implemented using an appropriate user subroutine in a finite element method (FEM) model. The AA6061-T6 aluminum alloy sheet, which has low formability and is fractured from its free edges in the IHF process, was examined as an example. Initially, a linear damage accumulation law, in which the prediction error is high due to the nonlinear stress and strain states in the IHF, was used to predict the fracture. Therefore, in the next step, a nonlinear damage accumulation function was utilized. While the nonlinear accumulation accurately predicts the single-stage IHF fracture, it is not able to predict the fracture well in the multi-stage IHF. It was observed that in multi-stage IHF, the rate of damage accumulation decreases with increasing the number of forming stages. Accordingly, a new nonlinear damage accumulation rule was developed. Experimental and numerical results indicated acceptable accuracy of the proposed nonlinear accumulation in the fracture prediction in the single- and multi-stage IHF processes.


Introduction
In conventional hole-flanging (HF) process, a sheet with a pre-cut hole that is fully clamped at the periphery area is formed to a smooth rounded flange by a punch and a die. In contrast, in incremental hole-flanging (IHF) process, the edge of the pre-cut hole is stretched and bent by the single point incremental forming (SPIF) process using a hemispherical tool [1]. In SPIF, a hemispherical tool moving in a predefined path by a CNC machine forms incrementally a sheet or tube clamped in non-deformation zones. So, in this process, no dedicated die is required [2]. Flanging has a wide range of applications in the household appliances and automobile industries. Providing material for thread cutting, hidden joints, rigidity, and strengthening of the edge of the sheet parts such as automobile front fenders and complex panels are some of these applications [3,4].
In IHF process, the gradual and localized forming of the sheet leads to obtaining flanges with deeper wall, compared to conventional HF. In addition, the lack of a dedicated die in IHF process significantly reduces the time and cost required to manufacture the process tool set, especially in multi-stage forming. As the formability in SPIF is limited by the fracture due to suppression of local necking [5], in IHF the necking is similarly suppressed and the part is failed by fracture [6].
Due to the local forming of the sheet in SPIF, and particularly in IHF process, severe thinning occurs in the parts. This limits the application of this process in industry. Multi-stage forming is an appropriate way to solve this limitation. In this technique, the forming area is expanded to a wider one, and severe deformation is generally prevented [7]. Moreover, by using a proper multi-stage strategy, the forming limits can be improved [8].
The first study on the multi-stage IHF was performed by Cui and Gao [1]. They presented three different strategies and showed that the strategy of reaching the final diameter by increasing the flange diameter in several stages improves the formability and enhances the forming limit. Centeno et al. [6] stated that the increase in sheet formability is achieved by multi-stage forming, when a certain range of process parameters is properly selected. An optimized two-stage strategy was proposed by Morales-Palma et al. [9] for the IHF process. Compared to other previously employed strategies, the one they presented improved the thickness distribution and production time. The effect of the pre-cut hole diameter on the four-stage IHF was analyzed by Hussain et al. [10]. It was shown that there is a threshold pre-cut hole diameter that has the highest formability in aluminum sheets. Besong et al. [11] reduced multi-stage IHF defects such as low geometrical accuracy and long process time and improved formability by a new method and tool.
The main difference between IHF and SPIF of parts without a pre-cut hole is the presence of a pre-cut hole in the sheet in IHF; as the diameter of the hole increases, the deformation mode changes. This discrepancy is because of a combination of in-plane stretching under plane strain conditions with bending around the perimeter of the flange [6]. In addition, it has been revealed that in IHF, different regions of the flange experience various strain and stress states [12,13]. This causes the deformation and failure mechanism in the IHF to be different from that in SPIF of a part without a pre-cut hole. Silva et al. [14] investigated the deformation mechanism of circular and square flanges and flanges with complex geometry. Variations of stress triaxiality and equivalent plastic strain in the single-stage IHF were studied by Martinez-Donaire et al. [15]. Concerning fracture in the IHF, Cristino et al. [16] investigated this phenomenon experimentally by tracing strain paths in different flange regions through a circular grid analysis. Using the obtained fracture toughness values by Ayada [17], and Rice and Tracy [18] criteria, the highest amount of damage was observed in the middle of the flange walls. Gandla et al. [19] assessed the fracture forming limit diagram (FFLD) of deep drawing quality steel sheets with six different ductile damage criteria in IHF. They found that the Ayada criterion provides the best prediction with the experimental FFLD.
As discussed in the literature on IHF, fracture usually occurs at the flange edges similar to that in the hole expansion process. Due to the complexities of edge fracture, it is an interesting topic to investigate. In recent years, some research has been done about edge fracture, especially in hole expansion of high strength metals such as dual-phase advanced high strength steels. Wang et al. [20] evaluated edge fracture in hole expansion process for DP780 steel using MMC ductile fracture criterion. They found that the shearing process to fabricate a pre-cut hole on the sheet has a significant effect on forming limit or hole expansion ratio (HER). Mu et al. [21] also studied the edge fracture by considering pre-exciting damage made by the cutting process. They used some ductile fracture models and showed that the Lou-Huh model has a better prediction for edge fracture in the hole expansion process. The edge fracture mechanisms of two dual-phase steel sheets were investigated by Teng and Chen [22]. They showed that edge cracking behaviors are essentially fractures with preexisting micro-cracks introduced by the shearing process. Mu et al. [23] improved the accuracy of edge fracture simulation by considering the effect of pre-damage field induced by the hole blanking on hole expansion process using a ductile fracture criterion.
Ductile fracture in metal forming processes is the result of severe plastic deformation due to the gradual deterioration of the material strength that causes loss of load-carrying capacity [24]. So far, various fracture models have been proposed to predict ductile damage and failure, which can be classified into the two general categories of coupled and uncoupled models. In the coupled models such as GTN [25,26] and Lemaitre [27], damage accumulation is involved in the constitutive equations, and the material yield surface is also considered. In contrast, in the uncoupled models, such as in most of the phenomenological models, the effect of damage on the yield surface is ignored, and the onset of fracture occurs when the damage indicator reaches a critical value [28,29].
From a microscopic point of view, the ductile fracture is the result of three stages: nucleation, growth, and accumulation of microvoids. Primary ductile fracture models based on microvoid growth were proposed by McClintock [30] and Rice and Tracey [18]. In these models, it was assumed that the main factor for the microvoid growth, and as a result the fracture, is stress triaxiality which is the ratio of the mean stress to equivalent stress. It was later shown that in addition to stress triaxiality, the Lode angle parameter that depends on the third deviatoric stress invariant must also be considered in the fracture mechanism in order to explain all the failure modes [31]. Hence, new criteria such as MMC [32], Mohr and Marcadet [33], Lou and Huh [34], and some practical models [35][36][37] have recently been proposed considering both the stress triaxiality and Lode angle.
Nonlinear stress and strain variation on the one hand and the fracture at the edges of the sheet perpendicular to the tool path in the IHF, especially in low formability sheets such as AA6061-T6 aluminum alloy, on the other hand, make the damage variation in the IHF as a complex phenomenon. As mentioned above, the deformation and fracture mechanism in IHF has been investigated by some researchers, but the interaction of damage and material deterioration leading to failure has not yet been accurately evaluated. In addition, although fracture prediction using various ductile fracture models has been done in IF processes, but it has not been studied in IHF, and the failure in IHF has been predicted only by using FFLD Although the edge fracture has been investigated in some forming processes such as the hole expansion operation, but damage and fracture mechanisms in IHF are different from those in the hole expansion process. Therefore, this paper aims to predict fracture in IHF, especially in the multi-stage forming, using ductile damage and fracture models. Due to the capability of the MMC ductile fracture criterion on damage and fracture prediction in various forming processes at a wide range of stress triaxiality, this criterion is used in this study. Initially, by using appropriate calibration tests, the MMC criterion was calibrated as a hybrid experimental-numerical method for AA6061-T6 aluminum alloy sheets. The calibrated criterion was then used to predict fracture in single-stage IHF using a linear damage accumulation law. The nonlinearity of stress and strain state in IHF prompted that linear accumulation could not accurately predict fracture in single-stage IHF. Therefore, an existing nonlinear accumulation law was used to predict the damage. The preliminary studies in this research showed that nonlinear accumulation could not predict failure in multi-stage IHF, especially in more than two stages. Therefore, there is a necessity for a new nonlinear accumulation law. Consequently, a new damage accumulation law is developed to improve fracture prediction in the multi-stage IHF. The most important goal of the current work is to improve fracture prediction in multi-stage IHF, which is done by developing a new nonlinear damage accumulation law for the first time.

MMC ductile fracture criterion
In this paper, the uncoupled phenomenological criterion MMC is used to model ductile fracture in AA6061-T6 aluminum alloy sheets. Most phenomenological models define damage based on stress and strain states. Stress state can be expressed based on the principal stresses or dimensionless stress parameters called stress triaxiality ( ) and normalized Lode angle parameter (NLAP) ( ) that are defined as follows [38]: (1) = m von in which m and ̃v on are mean stress and equivalent von Mises stress, and J 3 represents the third deviatoric stress invariant. Under plane stress condition, there is a relationship between two parameters η and , according to Eq. 3: The MMC ductile fracture criterion is a fracture function obtained by transferring the stress-based Mohr-Coulomb fracture criterion into the space of η , , and equivalent plastic strain ( p ) . This criterion can predict fracture in a wide range of , i.e., − 1 3 ≤ η ≤ 2 3 , and is defined as follows [32]: where A and n are strength coefficient and strain hardening exponent in the assuming work hardening law, respectively, and c 1 , c 2 , and c 3 are the three material parameters that must be calibrated through at least three fracture tests.

Damage accumulation law
The fracture function (as in Eq. 4) is defined under proportional loading conditions; thus, a damage accumulation law is required to consider damage variations under other conditions [39]. The damage accumulation law is usually a function of the variations of p with a weighting function that controls its increase rate. The weighting function depends on the stress state [40], and here it is assumed that this function is explained in terms of η and . Therefore, the damage accumulation law can be written as Eq. 5.
In the above equation, D is the damage indicator, a scalar variable, that describes the relative loss of the material ductility, and d p is the amount of increase in the equivalent plastic strain. At the beginning of the process, D is equal to zero, while at the instant of fracture, D = 1 and p = f , where f represents the fracture equivalent plastic strain. The weighting function f ( , ) is a constant under proportional loading, and at fracture, it can be written [40]: As a result, Eq. 5 can be reformulated as follows: Regarding Eq. 7, the value of D is a linear function of equivalent plastic strain; hence, this equation is called linear damage accumulation. Linear damage accumulation under nonproportional conditions cannot correctly explain the effects of p on fracture. Consequently, a nonlinear damage accumulation is required for nonproportional loading, as reported by Papasidero et al. [39]. So far, a number of nonlinear damage accumulation laws have been proposed, most of which are inspired by fatigue failure rules, such as the ones proposed by Cortese et al. [40] and Šebek et al. [41]. In this paper, the nonlinear damage accumulation presented by Xue [24], Eq. 8, is used.
where m is an exponent greater than zero. With m = 1 , Eq. 8 becomes the linear damage accumulation, Eq. 7. As m decreases relative to the unity, the rate of damage accumulation variation at the beginning of the process will be very high, and as the process progresses to fracture, it will continuously decrease. Conversely, if m is higher than unity, the reverse of the previous case will occur, and the damage accumulation will be faster in the final stages.
As mentioned, multi-stage forming improves formability over single-stage forming. The deformation in multistage forming is reduced to smaller parts according to the number of stages. This causes the stress and strain states to be different from those in single-stage forming. As reported by Mirnia et al. [42] and Talebi-Ghadikolaee et al. [43], in the multi-stage forming, the strain and stress states are highly nonlinear and will become more intense as the number of stages increases. This makes damage prediction in multi-stage forming much more complex than that in single-stage forming. Overcoming this complexity requires modifying the existing fracture models or damage accumulation laws. These modifications should be such that the effect of the number of forming stages is taken into account, so that the damage can be accurately predicted by changing the number of stages. Therefore, in this paper, according to the experimental observations and finite element examination of the IHF, Eq. 8 is developed in order to improve the damage prediction in multi-stage IHF in the form of Eq. 9 based on the number of stages (N) . In Sect. 5.4, it is described how the modifications are made.

Mechanical properties and fracture
behavior of AA6061-T6 sheet

Tension tests
In this research, AA6061-T6 aluminum alloy sheet with a thickness of 1 mm is used. In order to obtain the mechanical properties and anisotropy of the sheet, the uniaxial tension test according to ASTM E8-M standard was used in three directions of 0°, 45°, and 90° relative to the rolling direction. In addition to the uniaxial tension test (UT), three other tension tests, including notched tension (NT), plane strain tension (PS), and in-plane shear (SH) tests that involve different stress states, are considered to investigate the fracture behavior of AA6061-T6. Dimensions of tension test and prepared samples are shown in Figs. 1 and 2, respectively. All samples were prepared using a wire cut machine. Tension tests were performed on the SANTAM universal testing machine, and a 50-mm gauge length extensometer with the measurement range of 10 mm was employed to measure the exact displacement of the specimens. Each fracture test was repeated 3 times. The speed of the machine for performing the tension tests was set at 1 mm/min.

Finite element method (FEM) modeling of tension tests
Calibration of fracture criterion is performed by an inverse method based on the stress and strain history of tension tests. To do so, FEM simulations of tension tests were performed using the commercial FE code Abaqus/Explicit to obtain p f , as well as the variation of η and during the process until the onset of experimental fracture of the test. The expressed parameters were defined using a VUSDFLD user subroutine in the FEM software for each integration point. As shown in Fig. 3, all tension tests, except the SH test, were modeled fully without considering geometric symmetry. Only a half of the SH test was modeled with proper symmetry boundary conditions due to having two identical gauge sections. Tension tests were discretized using a 3D solid element C3D8R with a dimension of about 0.2 mm in the central section where the localized necking and failure is likely to happen and with a larger size in other sections. Through the thickness, 8 elements were considered.
The implemented mechanical properties and hardening law will be described in detail in Sect. 3.3. Hill's 48 yield criterion was assumed to express material anisotropy. As will be explained in Sect. 3.3, the simulation of tension tests was performed in two steps. In the first step, without considering the fracture criterion, the stress and strain states were investigated, and the MMC criterion was calibrated using the criterion. In the second step, the fracture was considered to validate the calibrated model. The element deletion technique was utilized to predict the onset and propagation of cracks. In this way, whenever the damage indicator D reaches unity for each element at all integration points, the respective element is removed and considered as the location of the fracture onset.

Plasticity model
As described in Sect. 3.1, the UT test was used in three directions of 0°, 45°, and 90° relative to the rolling direction to obtain the elastic-plastic properties and material  Table 1. The strain hardening behavior of the material was approximated by Eq. 10 through the linear combination of Swift hardening law and the behavior without hardening as described in [31].
According to Eq. 10, the material behavior up to the necking strain is expressed by Swift hardening law, and its constants k , 0 , and n called strength coefficient, pre-strain, and strain hardening exponent, respectively, are obtained by fitting the Swift equation with the experimental true stress-strain curve shown in beyond the necking strain. The parameter ̃U TS in Eq. 10 indicating no hardening behavior is named ultimate tensile stress at necking strain. The weighting factor Q , called the post-necking hardening parameter, is determined by the inverse method so that the force-displacement curve predicted numerically corresponds well to the experimental result. The parameter values of Eq. 10 are given in Table 1. The force-displacement curves obtained by experiment and predicted by the Swift model, as well as Eq. 10, are illustrated in Fig. 4b. According to Fig. 4b, compared to the Swift hardening, the calibrated hardening reduces the error between experimental and numerical force-displacement curves, especially near the fracture initiation region. Hill's 48 yield criterion is used to define the material anisotropy behavior according to Eq. 11.
where ij is the stress component and F, G, H, L , and N are the anisotropic parameters, calculated using the Lankford coefficients (Table 1) based on the relationships presented in [44]. These parameters are presented in Table 2.

MMC fracture criterion calibration
The MMC fracture criterion has three parameters that must be calibrated using appropriate tests that cover various stress states. To this end, four tests of UT, NT, PS, and SH were considered, as shown in Figs. 1 and 2. These tests were simulated using the material behavior described in the previous section, and their stress states were extracted. A comparison between the experimental and numerical simulation force-displacement curve of NT, PS, and SH tension tests is performed in Fig. 5. Figure 5a shows the experimental force-displacement curves obtained for the three PS tests. As can be seen from the figure, the difference between the results is not considerable which shows the repeatability of the PS experiments. For the other fracture tests, similar results were observed. According to the figure, the differences between the experimentally obtained and numerically predicted curves increase after the necking strain. It is due to the material softening because of plastic damage, which is not considered in uncoupled fracture criteria, including the MMC criterion. To improve the accuracy of predicting the force-displacement curve after the local necking, it is better to consider a coupled damage model to account for the softening effect of the material. In order to calibrate the fracture criterion using tension tests, the test loading paths are required. Using the results of FE simulation, the loading path, or in other words, the variation of η and with respect to p , for each test at the element with the highest p was obtained until the instant of the experimental fracture of the test. The element with the highest p was chosen because the failure in this element is probably higher than those of the other elements, so it designates the critical element. The onset of experimental fracture was determined when a sudden drop in the force-displacement curve was observed. The loading paths of the tension tests along with the critical element of each test are demonstrated in Fig. 6. As can be seen, the critical element for the UT, NT, and PS tests is in the center of the gauge within the thickness, while for the SH test, the critical element is on the outer surface away from the free edges. In addition, the loading paths are nonproportional for all the tests, or in other words, the values of η and are not constant until the fracture initiation. For this reason, the average value of each of these two parameters, given in Table 3, was calculated using Eqs. 12 and 13.

Fig. 3 FE model of tension tests
Regarding the stress state variation in tension tests, the loading history of the tests is considered for MMC criterion calibration, and the three parameters c 1 , c 2 , and c 3 can be obtained by minimizing the following function: where i fr is the fracture strain for the i th test and n t is the number of calibration tests, that is, 4 in this study. A proper code was scripted in the MATLAB program to minimize Eq. 14. The identified parameters of the MMC criterion are c 1 = 0.048 , c 2 = 227.67 , and c 3 = 0.932 . The 3D fracture locus of the calibrated MMC criterion is shown in Fig. 7.
Considering rupture, the tensile tests were again simulated by the calibrated criterion to assess its ability to fracture prediction. The fracture displacement values of the experimental and numerical results with the error percentage are presented in Table 3. The maximum error is observed in the SH test that is about 7%. It is worth mentioning that in all the simulations, the first deleted element indicating the onset of fracture was completely consistent with the critical element.

Experiment
AA6061-T6 square sheets with initial dimensions of 100 mm × 100 mm and pre-cut holes were used for experimental IHF tests. Figure 8 shows the die set comprising the main structure, backing plate, blank holder, and hemispherical forming tool mounted on a three-axis CNC milling machine. The backing plate has a hole with a 60 mm diameter; consequently, the diameter of the formed final flange will be the same. In IHF, the parts can be fabricated by different strategies. Strategies illustrate how a blank is deformed incrementally to convert to a final product. In other words, they are the paths that the forming tool passes to form parts. Strategies can be done in a single-stage or multiple stages. In multi-stage strategies, the parts are formed by several preforms in several stages. In this study, a multi-stage strategy with a gradual increase in wall angle was implemented, as shown in Fig. 9. In all the strategies, except the four-stage strategy, the wall angle of the first stage (ϕ) was 50°, and the increase rate of the angle in the next stages (α) was considered equal according to the number of stages until reaching the wall angle of 90°. In the four-stage strategy, the angle ϕ was considered 45° to obtain the identical value of 15% for α at the subsequent stages.
The process was first performed by placing the sheet on the backing plate, so that the pre-cut hole on the sheet was centered with the backing plate hole. The sheet was then clamped firmly to the backing plate by a blank holder to prevent the free edges of the sheet from moving. Finally, the forming tool moved on the sheet according to the predefined path in a three-axis CNC milling machine to form the sheet. The tool path was extracted using PowerMill software. The step-down and feed rate of the spiral tool path w adjusted to 0.5 mm and 1000 mm/min, respectively. Furthermore, the spindle, and consequently the tool, did not rotate. The pre-cut hole of the sheets was made by a wire cut machine with finishing low speed. Moreover, in order to eliminate or minimize the pre-damage resulting from the shearing process, the cut edges were thoroughly abraded with sandpaper. For reducing friction between the tool and sheet, hydraulic oil was used along the entire path of the tool.

FEM of IHF
The IHF process was modeled 3D, using commercial FE code Abaqus/Explicit, as shown in Fig. 10. Although solid elements were used in the tension test simulations to calibrate the fracture model (Sects. 2 and 3), in the IHF process modeling, the sheet was discretized with the S4R shell elements with five integration points in the thickness direction. As stated by Pack and Mohr [45] and Besson [46], fracture modeling is sensitive to the element type, i.e., solid or shell, and the element size.
On the other hand, using shell elements in processes with long computational time such as IF, is inevitable. Mirnia and Shamsari [47] showed that in the SPIF process, the difference between the modeling results of the shell elements with those of the solid elements is negligible. The mesh size of the sheets in the forming area was considered 0.5 mm, which was determined based on the convergence analysis of the mesh to p , so that the process simulation time was not too long, and that the effect of mesh size on fracture modeling was regarded. The forming tool and the backing plate were modeled as analytical  . 6 Variation of η and with respect to p at the critical element until fracture onset for tension tests and discrete rigid with R3D4 elements, respectively. Instead of the blank holder, the peripheral edges of the sheet were fully constrained, as were in the actual experiments. The contact between the sheet and the forming tool was defined according to Coulomb friction law with a friction coefficient of 0.1 based on the [47]. The mechanical properties and material behavior were the same as described in Sect. 3.3. Fracture modeling and field variables were similar to those for fracture tests. The tool path obtained by PowerMill software was converted into suitable boundary conditions for the forming tool using a code scripted in the MATLAB program so that the tool paths in the experiment and simulation are quite identical. An appropriate time scaling was used to reduce the processing time, and it was ensured that the ratio of kinetic energy to internal energy was less than 1%.

Validation of IHF simulation
In this section, the FE model of IHF was verified without considering fracture. Figure 11 shows a flange with a wall angle of 70° (for simplicity, 70° flange or similar expression is substituted in the remaining of the text), and a pre-cut hole diameter of 38 mm formed experimentally and numerically in the single-stage strategy. The thickness distribution along with the Fig. 7 The fracture locus for AA6061-T6 based on calibrated MMC criterion at the space of ( , , f ) Fig. 8 Experimental equipment of IHF process (a) forming tool, (b) blank holder, (c) backing plate, and (d) base structure Fig. 9 The schematic of multi-stage IHF final flange profile is shown in Fig. 11c, d, respectively. According to Fig. 11c, the highest thickness reduction occurred in the middle of the flange wall, indicating this area is more stretched than the others. The thickness near the flange edge is greater than that of the middle area of the wall, because the bending is the dominant mechanism in the edges, while in other areas stretching is the principal deformation mechanism. In addition, Fig. 11c illustrates that the FE simulation correctly predicts the trend of thickness variation as well as the location of the maximum thinning. The maximum difference between the thicknesses obtained experimentally and numerically in Fig. 11c is less than 3%.
The final profile of the 70° flange was measured using a coordinate measuring machine (CMM), and the results are shown in Fig. 11d. As it is seen, the final experimental and simulation profile of the 70° flange is well matched. The final flange wall angle, in other words, the spring back of the flange wall, obtained from experiment and simulation, is illustrated in Fig. 12. In the experiment, the spring back is 3.4° for the middle of the flange wall, and it is about 8.2° for the flange edges. In the simulated sample, the spring back of the middle of the flange wall and the flange edges are about 4.5° and 6.7°, respectively. As can be observed in Fig. 12, the flange edge has more spring back than the middle parts of the flange, which can be due to the non-constraining of one side of the flange edge. Figures 11 and  12 reveal that the results of FEM are considerably consistent with the experimental results and, thus, are reliable.

Damage variation and fracture prediction in single-stage IHF
Two flanges with pre-cut hole diameters of 38 mm and 43 mm were considered to investigate damage and fracture in the single-stage IHF process. Based on the experimental tests, the flange with a pre-cut hole diameter of 38 mm is formed successfully without fracture up to a wall angle of 70° and fails at greater angles. Therefore, fracture prediction was performed using the calibrated MMC criterion for the mentioned flange at two wall angles of 70° and 80°. In experimental tests, the 70° flange was fully formed (Fig. 11), but the 80°  Fig. 13. According to Fig. 13a, the 80° flange cracks at the edge, with the crack propagation path perpendicular to the tool path direction. First, the linear damage accumulation law of Eq. 7 was used to predict the damage. Applying this equation in FE simulation, the 70° and 80° flanges were fractured at tool heights of 8.33 mm and 8.19 mm, respectively. Figures 13b and 14a illustrate the predicted fracture location for 80° and 70° flanges, respectively. Thus, the MMC criterion with the linear damage accumulation law underestimates fracture intensively in both 70° and 80° flanges, where the error value in the 80° flange was approximately 42%. As a result, in the IHF, such as the IF of the parts without a pre-cut hole as reported by Mirnia and Shamsari [47], the linear damage accumulation cannot predict fracture due to strongly nonproportional variation of stress and strain in the IHF process. Hence, the nonlinear damage accumulation law of Eq. 8 was used to predict damage variation. In Eq. 8, there is another parameter (m) that needs to be identified. The IHF process itself was used to calibrate m. As stated in Sect. 2.2, by decreasing m more than unity, the damage accumulation rate is reduced at the end of the process, and the numerical fracture occurs later. On the other hand, m must be determined so that the 70° flange does not fracture in the simulation. In this regard, using iterative simulations, m has determined the maximum value that 70° flange would not fail in the simulation, as shown in Fig. 14b. Accordingly, the value of m was determined as 0.185.
Damage variations relative to p with linear and nonlinear accumulations for the critical element (element with maximum p ) are represented in Fig. 15. As shown in the figure, for linear accumulation, the gradient of the curve is almost constant, meaning that damage variation occurs at a constant rate. In contrast, for nonlinear accumulation, in which the damage variation rate is high at the beginning of the process, the gradient of the curve decreases in the descending order so that the damage accumulation rate decreases relative to linear accumulation. In addition, in both the cases, the damage value at the bottom surface of the sheet is higher than that at the top surface that is in contact with the tool. It indicates that the crack in the IHF starts from the bottom surface of the sheet and propagates to the top surface. Figure 16 shows the damage variations for all 5 integration points of the critical element. According to the figure, although the damage is 1 for the integration points 1 and 5, which are at the top and bottom of the sheet, respectively, for the internal integration points, including 2, 3, and 4, the damage is less than 1. As a result, the respective element has not been deleted in the simulation. It indicates that micro-cracks are formed at the top and bottom surfaces, but these microcracks are not propagated to the inner surfaces. Micro-cracks of the bottom and top surfaces are shown in Fig. 17.
For the 80° flange, the nonlinear accumulation was also used in the simulation, and the fracture was predicted at 15-mm tool height, which was about 5.5% more than that of the experiment. Figure 13c exhibits the cracked flange obtained from simulation. Moreover, as shown in Fig. 13, the crack location at the flange edges in the rolling direction of the sheet in the simulation is exactly the same as that in the experiment, indicating the accuracy of the calibrated fracture criterion for crack prediction in the IHF.
In order to further investigate the calibrated fracture criterion and nonlinear damage accumulation, flanges with a precut hole of 43 mm and a wall angle of 90° were produced, and the results are compared in Fig. 18. In the experiment, the flange is ruptured at the tool height of 10.42 mm, while the simulation predicted the fracture at the tool height of 9.86 mm with an error value of about 5.4%. In addition, as can be observed in Fig. 18a, b, cracking occurred both in the simulation and experiment at the edges of the flange located in the rolling direction of the sheet. Therefore, it can be concluded that the calibrated MMC fracture criterion with nonlinear damage accumulation is not only able to predict the failure instant in the single-stage IHF with acceptable accuracy but also can precisely detect the fracture position.

Damage variations and fracture prediction in the multi-stage IHF
As described in Sect. 4.1, a multi-stage strategy with a gradual increase in the flange wall angle with different numbers of stages was used. According to the experimental results, Fig. 12 The final 70° flange obtained from (a) experiment, (b) simulation the minimum pre-cut hole diameter of a flange with a wall angle of 90° and a final diameter of 60 mm that can be successfully formed is 44 mm, corresponding to the 5-stage strategy. Therefore, to investigate the failure of the parts formed by the multi-stage strategy, the pre-cut hole diameter of the sheets was considered to be 43 mm. Figure 19 shows  Table 4. According to the table, as the number of stages increases, except for the 4-stage strategy, the forming height increases slightly. The 4-stage strategy has the lowest fracture height among all the strategies. By comparing the fracture height of the multi-stage strategy presented at Table 4 with that of the single-stage strategy (Sect. 5.2), it can be seen that the fracture height of the single-stage strategy is higher than that of the 2-, 3-, and 4-stage strategies, which indicates that the multi-stage strategies employed in this research are not appropriate. Thus, for the material used in this research, i.e., AA6061-T6, the multi-stage IHF with a gradual increase in the flange wall angle except the 5-stage strategy does not improve formability, which is in contradiction with what was stated in the literature about the advantage of multi-stage strategy over single-stage strategy. To improve formability, another strategy with a different tool path should be used, for example, the strategies introduced by Cui and Gao [1], which are not discussed in this article. Fracture prediction for the multi-stage strategies was performed using the nonlinear accumulation of Eq. 8, and the formed flanges and the corresponding failure heights are represented in Fig. 20 and Table 4, respectively. According to the table, while in experiments, the failure for all parts formed by the multi-stage strategies occurred in the last stage with a 90° wall angle; in the FE simulation, the parts fractured in the second stage. It indicates that although the calibrated nonlinear accumulation accurately predicted  fracture in the single-stage strategy, it is not able to predict fracture in the multi-stage strategy, and even the higher the number of stages, the greater the error value. To examine this issue, the strain path and state, as well as the damage variations for the elements shown in Fig. 21, are compared at different stage numbers without considering the failure and element deletion, which are shown in Figs. 22 to 24. Figure 22 demonstrates the variation of p for the respective elements. It can be seen from the figure that for all elements, p is higher in the first stage for the single-stage strategy due to more deformation compared to that in other strategies, while for all multi-stage strategies, it is the same, because the amount of deformation in the first stage for all the strategies is identical due to the forming of the 50° flange. Similarly, in the second stage, the 2-stage strategy has the highest p , which results from more deformation in this stage that forms the 90° flange. By increasing the number of stages, the amount of p also increases, and in each stage, the strategy that forms the 90° flange has more p than other strategies. Finally, the maximum total p is related to the 5-stage strategy that has the highest number of stages. Despite all strategies form flanges with the same pre-cut diameter and final wall angle, their p is not equal. As the number of stages increases, the value of p for each stage is added to the previous stages, resulting in an increase in the total p . Increasing p leads to raising the damage according to Eqs. 4 and 8, as shown in Fig. 23.
According to Fig. 22e, f, element 8199, located near the flange edge, has higher p , compared to the other elements, and the lowest value of p corresponds to element 8187, which is located near the radius of the flange (Fig. 21a, b). By evaluating Fig. 21e, f, it can be seen that with increasing the number of stages, the increase rate of p decreases. Hence, it can be concluded that the damage accumulation rate decreases by increasing the number of stages. Figure 23 illustrates the strain path for the bottom surface of the respective elements for all the strategies. In Fig. 23, α is the ratio of the major to minor in-plane strains. According to the figure, the strain state for all elements is intensively nonlinear and constantly changing. Element 8187 is under a pure shear state in the singlestage strategy, but it experiences a uniaxial pressure state for multi-stage strategies. The strain state of element 8193 varies between biaxial tension and plane strain states, and the more the number of stages, the greater the number of changes. For element 8199, the final strain state in all strategies is the uniaxial tensile state, and the strain path changes from a plane strain state in the single-stage strategy to a pure shear state in the 5-stage strategy by increasing the number of stages. On the one hand, in all elements, as the number of stages increases, the strain path shifts to the left side, and the minor strain tends to the negative strain. On the other hand, increasing the number of stages increases the number of changes in strain paths. These changes reduce the rate of damage accumulation by increasing the number of stages. It can be seen from Fig. 23 that the deformation mechanism in IHF is different from that in SPIF of parts without a pre-cut hole. While the strain state in SPIF of parts without a pre-cut hole is usually in the form of a plane strain state [42], in IHF different regions of the flange experience various strain states.
Damage variations are investigated using the accumulation law of Eq. 8 for the respective elements in all strategies in Fig. 24. In this investigation, element deletion and crack initiation are not considered in order to fully assess damage variations. As shown in Fig. 24   Results of Figs. 22 and 24 revealed that in multi-stage IHF, by changing the number of stages, the rate of damage accumulation also changes, and the more the number of stages, the lower the rate of damage accumulation. Therefore, to accurately predict damage and failure in multistage IHF, a damage accumulation law is required that does not have a constant accumulation rate with different numbers of stages, while in Fig. 24, it is observed that the rate of damage accumulation according to Eq. 8 is constant for different stage numbers. Consequently, the damage accumulation of Eq. 8 cannot accurately predict damage in multi-stage IHF. Figure 25 shows the maximum value of p and damage based on the nonlinear accumulation of Eq. 8 for different numbers of stages. The figure clearly shows the nonlinear variations of p and the linear variations of the damage with increasing the number of stages.

Fracture prediction by the proposed damage accumulation law
As discussed in the previous section, the nonlinear damage accumulation of Eq. 8 could not accurately predict failure in the multi-stage IHF. Hence, a new accumulation law is proposed in the form of Eq. 9 to improve the fracture The derivation procedure of Eq. 9 and Eq. 10 is as follows: first, for considering the number of stages (N) on damage accumulation and reducing damage accumulation rate by increasing N , the exponent m in Eq. 8 is divided by N . By simulation, it was revealed that N needs an exponent. The exponent of 1 4 showed the best result for every number of stages. For the single-stage strategy (N = 1) , Eq. 10 will be equal to Eq. 8, indicating that this equation can be applied to all the strategies. In addition, it can be noted that in Eq. 10, there is only one parameter (m) that needs to be calibrated, similar to Eq. 8, which is obtained through the single-stage strategy. According to Eq. 10, by increasing N , the exponent of the ratio of equivalent plastic strain to fracture strain ( p ∕ f ) decreases, which causes the damage accumulation rate to be higher at the beginning of the process and lower at the end of the process.
The results of the prediction of Eq. 10 along with the prediction error are given in Table 4. Regarding Table 4, the error value of the proposed accumulation law for all the multi-stage strategies is less than 10%, which significantly reduces the error value compared to Eq. 8. The average error obtained by Eq. 1 is about 9%. Figure 26 shows the crack location predicted by Eq. 10. This indicates that the proposed damage accumulation law is capable to predict fracture in multi-stage IHF very well.
The damage variation obtained from Eq. 10 as a function of p for multi-stage strategies for the first removed element is shown in Fig. 27. According to the figure, at the beginning of the process, the 5-stage strategy has the highest damage accumulation rate, while the lowest  rate is related to the 2-stage strategy. But, as the process progresses, the damage accumulation rate decreases for all the strategies, such that at the end of the process, as opposed to the start, the 5-stage strategy will have the lowest damage accumulation rate. Therefore, based on the proposed damage accumulation, by increasing the number of stages, the rate of damage accumulation decreases, and fracture is delayed. Figure 28 compares variations of the damage accumulation obtained from Eqs. 8 to 10 for multi-stage strategies. According to Fig. 28, the maximum and minimum differences between the two accumulation laws are related to the 5-stage strategy and the 2-stage strategy, respectively. So, adding the numbers of stages increases the differences between the two accumulation laws. This increase due to increasing the number of stages makes the proposed damage accumulation provided a good fracture prediction for different numbers of stages.

Conclusions
This study was conducted to numerically predict the damage and fracture in single-and multi-stage IHF processes. AA6061-T6 sheet was implemented in this study, and failure was predicted using the MMC ductile fracture criterion. Due to the nonlinear nature of stress and strain states in the IHF, a nonlinear damage accumulation law was employed to evaluate the damage in the single-stage IHF. The damage mechanism in multi-stage IHF was analyzed in detail, based on which a new nonlinear damage accumulation law was developed to precisely predict fracture in the IHF of AA6061-T6 sheets. The proposed damage accumulation law can be very effective for reducing time and cost in multi-stage flanging, especially for trial manufacturing. Based on the performed investigations, the following results are extracted: -Linear damage accumulation due to nonproportional and nonlinear strain paths in IHF does not accurately predict damage and fracture in this process, and its prediction error is very high. Therefore, a nonlinear damage accumulation was employed. Based on the simulation of the 70° flange, an additional nonlinear parameter of damage accumulation was obtained. The error of fracture prediction in single-stage IHF with nonlinear damage accumulation is less than 6%. -Although the nonlinear accumulation has a correct prediction of the instant and location of the fracture in the single-stage IHF, it cannot accurately predict the dam- age in the multi-stage IHF. In addition, as the number of stages in the multi-stage IHF increases, the prediction error value increases too. -By examining the strain path and state in single-and multi-stage IHF, it is found that with increasing the number of stages, not only the variation of the strain path increases but also the strain states and the deformation conditions change, compared to the single-stage IHF. It is observed that as the number of stages increases, the minor strain tends to the negative side, and the strain variation curve shifts to the left hand. -Variations of p and the major and minor in-plane strains in the multi-stage IHF indicate that the rate of damage accumulation changes and decreases with increasing the number of stages. Therefore, to correctly predict fracture in the multi-stage IHF, a damage accumulation law is required that has different trends and rates for various stage numbers. Hence, a new damage accumulation law was developed in which the effect of the number of stages was considered.
The proposed damage accumulation significantly improves the damage prediction in the multi-stage IHF compared to the existing nonlinear damage accumulation laws.
For future work, the proposed damage accumulation law can be evaluated for other strategies and geometries in IHF and for other incremental forming processes.
Author contribution All authors contributed equally to this work.
Funding This work was supported by the Babol Noshirvani University of Technology (grant number BNUT/955130011), Seyyed Emad Seyyedi.

Declarations
Ethics approval Not applicable.