A Review on Sheet Metal Forming Evaluation of Advanced High Strength Steels (AHSS)

1. Introduction

The 2023 report from the European Environmental Agency on climate change underscores that the transportation sector constituted a substantial contributor to greenhouse gas emissions, comprising nearly 25% of the total emissions across European Union member states in the year 2021 [1]. European Union government policies are committed to achieve a 23% reduction in greenhouse gas emissions from the transport sector in 2030, compared to 2005. In vehicles, a 10% weight reduction corresponds to fuel consumption saving of 6.5% [2]. The ultimate goal of mitigating CO2 gas emissions through the reduction of fossil fuel emissions in vehicles has prompted automobile manufacturers to prioritize lightweight construction as a main concern [3,4,5]. The CO2 emissions of newly registered cars are substantially decreasing, with noteworthy year-on-year declines of 12% for both 2020 and 2021. This trend is partially attributed to the reductions of Body in White (BiW) components overall weight, which are the main contributors of the overall weight of vehicles [6]. Along with environmental concerns, passenger safety may never be compromised. In order to comprise lightweight construction, safety, and structural performance, vehicle’s material list have been continuously shifting through the employment of a higher variety of materials, e.g., through the usage of aluminum alloys and Advanced High Strength Steels (AHSS).

The use of AHSS have become increasingly popular in the automotive industry, due to their appropriate mechanical properties, namely: high structural strength, high stiffness, high work hardening capacity, low cost, high performance and mass optimization capacity [7], while maintaining good levels of ductility [8,9]. Figure 1 presents the relationship between Ultimate Tensile Strength (UTS) and Total Elongation of AHSS compared to conventional steels. The enhanced high strength to density ratio of these materials renders them a competitive contender of aluminum and magnesium alloys, as thinner sheets can be used. Based on the material requirements of different automotive components conventional High Strength Steels are typically applied to panels (both outer and inner) and structural parts, while AHSSs are engineered for anti-intrusion barriers such as bumpers and impact beams [4], as illustrated in Figure 2.

Sheet metal forming operations are widely employed in the automotive 2nd tier [12,13]. Table 1 summarized AHSS types, characteristics and common applications in the field of automotive industry. The manufacturing of Body in White components heavily relies on various sheet metal forming processes such as stamping and deep drawing. Although the major advantages of AHSS employment, the forming process of such materials is increasingly complex. This can be attributed to their multiphase microstructural composition, high yield strength as well as the thinner sheet thickness considered. Thus, the usage of AHSS intensified and introduces forming challenges almost never encountered before. These forming challenges are depicted in Figure 3. The significant yield strength leads to higher forming forces, energy consumption and higher stresses [14,15]. Thus, requiring the selection of appropriate press machines. Moreover, AHSS exhibit greater tendency to wrinkle due to the lack of adequate hold-down and a reduction of sheet thickness is frequently observed. Springback and dimensional control inconsistencies are accentuated, requiring the need of tooling shape modification and/or adoption of progressive die compensation techniques [16,17,18,19]. On the other hand, the multi-phase microstructure of AHSS leads to batch-to-batch inconsistencies. Simultaneously, die designs are getting increasingly complex. Sharp edges and quick curvature changes are frequently introduced to meet aesthetic concerns [20].

In accordance with the aforementioned forming challenges associated with manufacturing AHSS, the evaluation of the sheet metal forming process solemnly through the analysis of the traditional forming limit diagrams (FLD) can lead to inadequate and misleading interpretations. The material susceptibility to different fracture typologies increased with the employment of AHSS. Thereby justifying the need for more accurate fracture modeling of sheet metal forming operations [22]. Once the factors affecting failure are understood, engineers can minimize the risk of failure. This, in turn can improve product quality and efficiency, reduce scrap rates and the number of experimental tryouts [23].

This review paper introduces an assessment of the applicability of the Forming Limit Diagram (FLD) in the post-processing of Advanced High-Strength Steel (AHSS) forming processes. Following this introduction, the limitations of the FLD are discussed. The subsequent discourse emphasizes the critical necessity for advanced and precise fracture/damage models within the realm of AHSS forming processes, thereby accentuating their significance in the contemporary analysis of post-processing results for AHSS.

Table 1. Summary of different types of 1st generation of AHSS: microstructure, mechanical properties and typical applications. Adapted from [4,7,14,24,25,26,27].

AHSS	Microstructure	Mechanical Properties							Characteristics	Typical Applications
		${\mathit{R}}_{\mathit{P}0.2}$ (MPa)		${\mathit{R}}_{\mathit{m}}$ (MPa)		A80 %		$\mathit{t}$ (mm)
Dual Phase (DP)	Ferrite + Martensite	Docol DP500							Conventional: Lower Yield Strength Low Yield-to-Tensile Strength Ratios High Initial Work Hardening No Yield Point Elongation Significant Bake Hardening Good Uniform Elongation High Energy Absorption Good Cold Formability	Body Panels; Front and Rear Longitudinal Rails; Supporting Structure
		290-380		490-600		24		0.5-2.1
		Docol DP600
		330-430		590-700		20		0.5-2.1
		Docol DP800
		440-550		780-900		14		0.5-2.1
		Docol DP1000
		590-900		980-1130		8		0.5-2.1
									Higher Yield Strength Grades:High Yield-to-Tensile Strength RatiosSome Yield Point ElongationLower Bake HardeningReasonable Uniform ElongationGood Energy Absorption	Passenger Safety Cage Components (limited by axial or transverseBending): Rockers, Pillars, Pillar Reinforcements,Roof Rails and Cross Members
Transformation Induced Plasticity (TRIP)	Ferrite Matrix + Retained Austenite Hard phase of Bainite and/or Martensite	TRIP700							Excellent Combination of Strength and Ductility High Total Elongation Considerable Strain Hardening Capacity High Energy Absorption Excellent Formability High Fatigue Endurance Resistance	Cross Members; Longitudinal Beams; B-Pillar Reinforcements; Sills; Bumper Reinforcements
		400-520	690-800		24		1
		TRIP800
		450-570	780-910		21		1
Martensitic Steel (MS)	Martensitic Matrix + Small Amounts of Ferrite and/or Bainite	Docol 900M							High Yield Strength High Tensile Strength Good Hardness High Toughness	EV Battery Protection; Sill (Rocker) Reinforcements; Roof Reinforcements; Door Beams; Bumper Beams; Seat Structures; Rocker Panel Inners
		700-1000	900-1100		3		0.5-2.1
		Docol 1100M
		860-1100	1100-1300		3		0.5-2.0
		Docol 1500M
		1220-1520	1500-1750		3		0.5-2.1
		Docol 1700M
		1350-1700	1700-2000		3		1-2.1
Complex Phase Steel (CP)	Ferrite/Bainite Matrix + Small Amounts of Martensite; Retained Austenite and Pearlite	Docol 600CP							High Yield Strength Very High Ultimate Tensile Strength Excellent Uniform Ductility High Energy Absorption Excellent True Fracture Strain High Fatigue Strength Good Impact Strength High Residual Deformation Capacity Good Bendability	Seat Flange; Door Bar; Tunnel Stiffener; Rear suspension Bracket; Fender Beam; Seat Tracks; Seat Recliners
		350-500	600-740		16		0.7-2.5
		Docol 800CP
		780-950	980-1140		6		0.5-2-1
		Docol 1200CP
		900-1100	1180-1350		5		0.5-2.1

Table 1. Summary of different types of 1st generation of AHSS: microstructure, mechanical properties and typical applications. Adapted from [4,7,14,24,25,26,27].

AHSS	Microstructure	Mechanical Properties							Characteristics	Typical Applications
		${\mathit{R}}_{\mathit{P}0.2}$ (MPa)		${\mathit{R}}_{\mathit{m}}$ (MPa)		A80 %		$\mathit{t}$ (mm)
Dual Phase (DP)	Ferrite + Martensite	Docol DP500							Conventional: Lower Yield Strength Low Yield-to-Tensile Strength Ratios High Initial Work Hardening No Yield Point Elongation Significant Bake Hardening Good Uniform Elongation High Energy Absorption Good Cold Formability	Body Panels; Front and Rear Longitudinal Rails; Supporting Structure
		290-380		490-600		24		0.5-2.1
		Docol DP600
		330-430		590-700		20		0.5-2.1
		Docol DP800
		440-550		780-900		14		0.5-2.1
		Docol DP1000
		590-900		980-1130		8		0.5-2.1
									Higher Yield Strength Grades:High Yield-to-Tensile Strength RatiosSome Yield Point ElongationLower Bake HardeningReasonable Uniform ElongationGood Energy Absorption	Passenger Safety Cage Components (limited by axial or transverseBending): Rockers, Pillars, Pillar Reinforcements,Roof Rails and Cross Members
Transformation Induced Plasticity (TRIP)	Ferrite Matrix + Retained Austenite Hard phase of Bainite and/or Martensite	TRIP700							Excellent Combination of Strength and Ductility High Total Elongation Considerable Strain Hardening Capacity High Energy Absorption Excellent Formability High Fatigue Endurance Resistance	Cross Members; Longitudinal Beams; B-Pillar Reinforcements; Sills; Bumper Reinforcements
		400-520	690-800		24		1
		TRIP800
		450-570	780-910		21		1
Martensitic Steel (MS)	Martensitic Matrix + Small Amounts of Ferrite and/or Bainite	Docol 900M							High Yield Strength High Tensile Strength Good Hardness High Toughness	EV Battery Protection; Sill (Rocker) Reinforcements; Roof Reinforcements; Door Beams; Bumper Beams; Seat Structures; Rocker Panel Inners
		700-1000	900-1100		3		0.5-2.1
		Docol 1100M
		860-1100	1100-1300		3		0.5-2.0
		Docol 1500M
		1220-1520	1500-1750		3		0.5-2.1
		Docol 1700M
		1350-1700	1700-2000		3		1-2.1
Complex Phase Steel (CP)	Ferrite/Bainite Matrix + Small Amounts of Martensite; Retained Austenite and Pearlite	Docol 600CP							High Yield Strength Very High Ultimate Tensile Strength Excellent Uniform Ductility High Energy Absorption Excellent True Fracture Strain High Fatigue Strength Good Impact Strength High Residual Deformation Capacity Good Bendability	Seat Flange; Door Bar; Tunnel Stiffener; Rear suspension Bracket; Fender Beam; Seat Tracks; Seat Recliners
		350-500	600-740		16		0.7-2.5
		Docol 800CP
		780-950	980-1140		6		0.5-2-1
		Docol 1200CP
		900-1100	1180-1350		5		0.5-2.1

2. Failure, damage and fracture in sheet metal forming

Failure and fracture are two distinct terms that are often used interchangeably, although having different meanings. In the realm of sheet metal forming community, failure refers to the point at which the material can no longer withstand the applied stress and loses its ability to maintain its shape, due to the loss of load-carrying capacity [28]. Thickness instability is typically referred to as failure. Fracture, on the other hand, relates to the point at which material separation occurs, due to excessive stress. Damage is defined as an intermediate state, a product property that accumulates between processes. It represents the undesired evolution of one or more structural characteristics, e.g., micro voids, that hinder an engineering capability [29].

The ductile fracture phenomenon can be divided into two modes: (i) Ductile Normal Fracture; (ii) Ductile Shear Fracture [30]. The former is caused by void nucleation, growth and coalescence leading to damage growth, eventually fracture where a crack appears. The latter is a result of shear band localization which may occur through the thickness or in the plane. According to Mattiasson et al. [30], there are types of fracture scenarios, as illustrated in Figure 4. The first one addresses the commonly observed outcome in sheet metal forming operations. Constitutive failure models are adequate to provide proper results in this scenario. However, for the other two scenarios, a fracture constitutive model is required to accurately capture fracture occurrence.

Formability is defined as the ability of sheet metal to be formed into the desired shape without necking or cracking [31]. Thus, formability evaluation should comprise failure mechanisms on account of localized necking and fracture mechanisms not preceded by necking. The FLD lies under the assumption that localized necking precedes damage growth and failure. In this context, necking instability defines the process limits [22]. AHSS require the need of modelling fracture occurrence without necking. Such eventualities can arise under specific stress states such as in-plane shear deformation and nearly pure bending loading [20].

3. Traditional forming limit diagram (FLD)

Formability evaluation is usually performed through the FLD. The simplicity of representing the principal strains, ${\epsilon }_1$ and ${\epsilon }_2$, in a graphical way is of great interest for industry. The critical feature of the FLD is the FLC, which describes the strain limits at the onset of necking: the strains located above the FLC fail in a critical state, where necking instability has occurred; the strains behind the FLC are in a safe region, characterized by remaining formability. Therefore, the failure locus for necking instability under various loading conditions, ranging from uniaxial tension to biaxial tension, is depicted by the FLC.

The FLC is affected by different factors. Even for the same material grade, variations in the FLC can be noticed. The reason for that lies in the dependency of factors such as forming temperature [32], material thickness [33], forming speed [33], strain rate [32] and friction conditions [34]. Liu et al. [35] investigated the strain rate effect through a comparison between AHSS and three different kinds of conventional steels. They concluded that the AHSS are much more dependent on strain rate variations rather than the conventional steels.

3.1. FLC experimental determination

The FLC is determined in accordance with the ISO 12004-2_2008 standard [36] using either Nakazima [37] or Marciniak [38] test procedures.

The Nakazima testing method involves the use of specimens with varying cutout widths which are clamped using a circular blank-holder. Subsequently, the specimens are drawn until fracture through a drawing die by the rectilinear movement of a hemispherical punch. Draw beads in tools are used to restrain the blank between the die and blank-holder and also to avoid stress concentration at the transition from strip section to the blank [39]. Nakazima tests should be executed in frictionless conditions between the contact zone of the punch and specimen pair [40]. As an out-of-plane stretching method, Nakazima test is dependent on the tool geometries, through-thickness pressure and lubrification conditions [22].

On the other hand, the Marciniak in-plane test involves the deformation of specimens under balanced biaxial conditions. By modifying the washer and specimen geometries, the FLC can be generated [22]. An adapted design of the punch geometry is adopted in order to minimize the friction in the specimen central region.

Irrespective of the performed method, a wide range of deformation modes can be achieved by modifying the specimen geometry. Figure 5 portrays the influence of specimen geometry on deformation modes. To cover a broad range of deformation modes, at least five different specimen widths must be tested [36]. The FLC shape is strongly dependent on strain path, hence on specimen geometry. The more the specimen tested, the more accurate the FLC formability prediction.

The experimental methods to determine the DLC are usually time-consuming, low cost-effective and limited by experimental conditions, even though they are more trustable than other approaches [43]. Many researchers select the Nakazima test procedure due to the simplicity of tools required [44].

3.2. Limit strains determination methods

The accurate assessment of the limit strains consists of a major step for the experimental determination of the FLC. The onset of localized necking is an instantaneous moment, difficult to detect and, even more, to measure. Frequently, the limit strain methods calculations are divided into: (i) position-dependent methods; (ii) time-dependent methods; (iii) time-position-dependent methods.

The position-dependent methods are characterized by the individual, well established single and unique evaluation of the onset of necking. ISO 12004-2:2008 [36] recommends a position-dependent approach. As a standard, it is intended to provide practices and unify test condition procedures in order to gather consensus between different laboratories. According to this method, the strain variations are analyzed prior to the occurrence of a crack. The procedure for limit strains determination in accordance with ISO standard are provided by Martínez-Donaire et al. [45]. An alternative position-dependent method was proposed by Zhang et al. [46] to provide more stable and precise FLC determination.

The time-dependent methods rely on the analyses of the strain, strain rate or strain acceleration to determine the onset of through-thickness necking. In the work developed by Merklein et al. [47] the necking initiation occurs in the highest value of a linear regression coefficient curve, selected based on the rate of change of strain rate. Iquilio et al. [48] developed a more realistic method to precisely calculate the limit strains of a stainless steel with an heterogeneous strain distribution. Another time-dependent methods were proposed by Huang and Shi [49], Volk and Hora [50] and Hotz et al. [51]. Overall, these techniques rely on mathematical manipulation of strain rate rather than physical behavior of materials, which can be unreliable due to noise in data signals and may not be applied to all sheet metals [52].

Time-position-dependent methods are hybrid approaches which consider either the temporal evolution of strain, position or displacement of points on a line drawn perpendicular to the fracture direction [52]. Martínez-Donaire et al. [53] inferred the initiation of necking when the slope in the first special derivative of the vertical displacement remained constant within the necking region. Furthermore, Min et al. [52] established a time—position-dependent method to identify the onset of localized necking under in-plane deformation based on the measured surface geometry of the test specimen.

Although the availability of a wide number of limit strain methods there is no consensus about the method to be chosen, since the overall methods possesses their specific applications. Therefore, a continuous search for more accurate, innovative methods to incorporate specific loading conditions, typology of fracture is ongoing to handle with ISO standard limitations, namely: (i) unfit for significant strain gradients across the sheet thickness, e.g., when using small punch radii or corner radii in forming dies; (ii) unsuitable in stretch-bending operations; (iii) inapt for multiple necking implying variable strain distribution; (iv) inappropriate for sudden fracture leading to limit strain imprecisions particularly under plane strain or biaxial testing conditions.

3.3. Grid marking and deformation determination methods

Materials employed in sheet metal forming operations have a reduced sheet thickness, particularly the AHSS where the sheet thickness commonly used is less than 1.5 mm. As a result, surface strain sensitivity significantly increases and requires more advanced strain measurements methods for an accurate formability evaluation. The FLC determination is influenced by the strain measurement technique employed, implying concerns regarding grid marking type and method, which also plays a key role in addressing major and minor strain calculations. Grid marking is performed prior to the formability test experiment and followed by strain measurement. Figure 6 displays the main tasks involved in the FLC determination.

Grid marking is defined as the process of printing line patterns on the surface of the area of interest on the sheet sample [54]. There are mainly four grid methods for creating grid patterns. Circular or square grid types are the most commonly used for assessing formability. The choice between the available methods must take into considerations parameters such as durability, quality, resolution, contrast, pattern accuracy and cost. Nevertheless, one must guarantee the independence of process and forming conditions, such as lubrification and friction, from grid marking [54]. Ozturk et al. [54] introduced the main differences between the overall grid marking methods and a qualitative comparison of the grid marking methods was appointed. Laser etching process require the study of power speed, focus and width parameters [55]. Guk et al. [56] utilized a diode laser marking system to generate grid patterns on multiphase steels: the scanning speed and laser power influence on the mark depth was investigated. Furthermore, Yildiz and Yilmaz [57] after applying four different grid marking methods, namely: serigraphy, electro-chemical etching, photo-chemical etching and laser concluded that laser is best suited due to ease of application and stability.

The circular or square patterns are printed before deformation and transform into ellipses or quadrilaterals after forming operation, respectively. The major and minor strains are determined by comparing previous and subsequent patterns dimensions [58]. ISO 12004-2:2008 [36] stipulates that the dimensions of the grids are related to the thickness of the material, e.g., the circle diameter on circular patterns has a significant impact on strain calculation [59]. Speckle patterns are required to be used simultaneously with Digital Image Correlation (DIC) systems for strain measurement. Li et al. [60] sprayed speckle pattern on the top surface of the sample and one must guarantee the adhesion of the paint to the surface before and after deformation.

The most practical method of measuring the strains in a formed part is to utilize measuring grid pattern on the part surface. This can be done either manually or by use of an automatic method [55]. Manual techniques are characterized by several drawbacks such as high measurement time requirements, less accuracy and low resolution. Even though, Bandyopadhyayet et al. [61] adopted a manual strain measurement method to build the stress based forming limit diagram of a tailor welded blank material. Additionally, the deformed grids of a Inconel 718 was assessed by means of a stereo microscope [62].

The alternative automatic methods are the Single Point Strain Measurement as well as the Full Field Strain Measurement [43]. Both techniques require less time, and have better resolution, hence contributing with more reliable results. The main components in the automatic strain measurement system are the image acquisition hardware and image processing software [63]. ASTM E2218-02 standard [64] provides guidelines for strain measurement.

Single Point Strain Measurement is characterized by the measurement of one grid element separately. Grid Pattern Analyser (GPA) and Forming Measurement Tool Innovations (FMTI) are two available commercial software for automatic single-point strain analysis built on image processing. The image acquisition demands machine vision cameras to capture the deformed grids after forming and specialized lenses. Afterwards, the software fits the ellipse to the deformed grids and calculates major and minor principal strains [59].

Full Field Strain Measurement methods such as Digital Image Correlation (DIC) possess the advantage of measuring larger areas considerably quicker. This method, applicable either in two or three dimensions, require a significant amount of material images prior and after forming operation with a geometrical relationship between them to derive gradients and correlate them in order to calculate the major and minor strains [65]. As previously mentioned, DIC system must be used in association with image processing software. The commonly used are namely: German Argus or ARAMIS system of GOM company, the American Grid Analyser Model 100U of the FMTI system Inc., the VIC-2D(3D) from the Correlated Solutions Inc., the Dante Dynamics as well as the ICASOF [59,63]. In spite of the high measurement accuracy provided by these systems, they are characterized by an expensive cost [66]. Open source DIC alternatives are a cost effective solution such as: Ncorr [67], ALDIC [68], py2DIC [69], µDIC [70] and Multi-DIC [71]. Commonly these alternatives require a huge memory of images taken from a high-speed image acquisition system [59]. The DIC Full Field Strain Measurement Systems have been broadly used to build the FLC [72,73]. Min et al. [74] and Huang and Shi [49] compared standard, spatial and temporal methods to determine the onset of localized necking by the use of DIC data. The introduction of these systems able to record and memorize the whole forming process has shifted research attention towards time-dependent methods as alternatives to position-dependent methods.

Additionally, DIC may be used to calibrate ductile damage models, evaluation of anisotropic materials as well as to obtain strain distribution either in tensile tests or in a forming operation. In the work conducted by Gkolfinopoulos and Chijiwa [75], DIC analysis alongside GOM commercial software and a full frame CMOS camera enabled to investigate strain concentration in the tensile specimen necking region before fracture. The damage model parameters were numerically validated with the aid of DIC analyses.

4. Limitations of FLD

Forming limit diagrams play a crucial role in assessing the formability of metal sheets by defining the maximum strain to failure in terms of major and minor strains. The primary objective is to identify the onset of necking in limited regions of the examined sheet metals [76]. Industrial sheet metal forming manufacturers are interested in the phenomenological outcome observed, as such the FLD is well established as a failure criterion. However, FLD has some limitations with respect to forming behavior modelling accuracy. Those limitations are depicted in Figure 7.

On one hand, the FLD as a necking based failure criterion is intended to detect through-thickness necks and splits which lead to a global loss of carrying capacity of the material. Moreover, the range of study of the traditional FLD is limited to the deformation paths illustrated previously in Figure 5. Therefore, the points below the left diagonal are not covered. On the other hand, local instabilities as well as direct fracture on material surface are not covered by the FLD analysis. The forming industry trend towards high strength/low density alongside with sharp geometrical features accentuates forming challenges. Thus, the need for more advanced models able to model forming behavior, failure and fracture is a current demand. The AHSS tend to exacerbate these forming challenges in comparison to the conventional mild steels, as such direct fracture on material surface is more likely to occur.

To overcome FLD limitations, the special characteristics of the AHSS must be taken into consideration. On accordance, more advance fracture and damage models have been developed to meet the special characteristics of these kind of materials. Usually, these models require experimental characterization at different stress states [77]. Different specimen geometries and loading conditions are adopted. Section 6 presents the state of the art of the testing methods employed for a fully material characterization. Section 7 addresses the current trends in failure and damage models.

4.1. Non-linear strain paths

The main disadvantage of the traditional FLD is that the different experimental points in the standard FLC are all determined under the assumption of simple proportional loading states. Therefore, it fails to provide proper results for non-linear strain paths histories in sheet metal forming processes [78]. In multi-stage forming processes, or even at deep drawing operations under complex loading conditions, the probability of loading paths deviation from linearity is more commonly observed [79]. Furthermore, the traditional LD can also fail to provide adequate precision of the forming process because the material exhibits a different strain hardening behavior when a strain path change occurs compared to monotonic loading [80]. Figure 8a) illustrates some possible non-linear strain paths observed in sheet metal forming operations. Da Rocha and Jalinier [81] stated that a biaxial solicitation preceded by a uniaxial loading leads to an increased formability, as opposed to a biaxial solicitation followed by a uniaxial loading, as illustrated by Figure 8b). Thus, the disregard for this shift in loading may result in both conservative and optimistic outcomes with respect to the traditional FLD. Figure 8c) presents a deformation path exhibiting high non-linearity leading to failure before the forming limit is reached under the standard FLC.

Manopulo et al. [83] used a theoretical approach based on an extended Modified Maximum Force (eMMF) criterion to account non-proportional loading. Stress-based models are effective in modelling fracture occurrences in multi-forming operations that are characterized by non-linear strain paths [82]. However, these models may be limited to deformation histories within the Nakazima range. Non-linearities that arise from predeformation in shear or compression can be effectively modeled by damage accumulation models [20]. Furthermore, the Generalized Forming Limit Curve (GFLC) offers advantages such as its independence from the material model and its ability to represent multilinear strain paths. The GLFC was developed by Volk et al. [84] and it is currently implemented in AutoForm R10. Section 6.2 presents experimental techniques specially designed to address non-linear strain paths generation.

4.2. 3D stress states

Typically, the process of thin sheet forming operates under the plane stress mode [85]. Under this context, as deformation leads to thinning, the traditional FLD predicts accurately the forming behavior. However, when strain states occur below and within the left diagonal (${\epsilon }_1$=${\epsilon }_2$), the presence of a compressive hydrostatic stress component in the thickness direction mitigates the risk of necking. Thus, even under severe thinning, necking will not occur, and the material will deform until the fracture limit is reached. For an accurate prediction of the forming behavior in ironing and coining operations, in-plane shear deformation modes as well as surface defects, the consideration of 3D stress states under the left diagonal of the FLD diagram is required.

The strain increments generated in this region can be analyzed by means of the stress triaxiality. The stress triaxiality is a reliable indicator to determine whether a material is subjected to volumetric tension or a compressive stress state. Hence, failure models considering stress triaxiality have garnered attention as they can identify strain increments arising from a compressive stress state, thereby not contributing to the risk of necking.

4.3. Shear fracture

Within the realm of AHSS, there is a heightened frequency of shear fracture occurrence on tight radii [86,87,88]. Shear deformation inherently maintains the thickness of the sheet [20] and the FLD fails to encompass the second quadrant beyond pure shear conditions [89]. The fracture limit in the low triaxiality regime exhibits a greater sensitivity to the Lode angle [20]. In fact, the pure shear fracture limit can be significantly lower than the ductile fracture limit [90]. Li et al. [89] have effectively showcased the limitations of the FLD in accurately predicting shear fracture, particularly due to the relatively lower pure shear fracture limit when compared to the uniaxial tension limit, as depicted in Figure 9.

The determination of fracture limits in pure in-plane shear remains an area of ongoing investigation, as there is currently no established standardized experimental method for characterizing shear fracture modes. The main difficulty associated with the achievement of pure shear fracture limits regards with the crack opening modes which generally occurs near the transition region between uniaxial tension and in-plane shear conditions [91]. Furthermore, the double-notched shear (DNS) test is commonly used to determine fracture strains induced by crack opening through in-plane shear. It should be noted that this test does not strictly fulfill the conditions of pure shear due to the presence of bending effects and the development of localized necking at the end of the testing process [92]. Roth and Mohr [93] proposed 600 distinct specimen geometries and demonstrated the value of using three instead of one type of shear specimen. Further insights on specimens implying shear fracture are given in section 6.1.1.

4.4. Edge cracks

The FLD is also not reliable in accurately predicting edge cracks on shear cut surfaces during component forming. The preparation of blank edges presents a wide range of possibilities, making it challenging to describe the reduced formability on the edge [20]. The stress state experienced at the blank edge is typically uniaxial, implying that it should undergo necking before fracturing under tension. However, due to the shearing process, the material properties and surface characteristics of blank edges are significantly altered. Consequently, the occurrence of fracture at the edge can vary depending on microscopic features or defects before the material reaches its actual necking limit [94]. Pre-forming is an important forming process in the prevention of edge cracking in flanging and edge stretching processes. Various process parameters in shear cutting have been identified to impact the residual formability of metallic materials. Nasheralahkami et al. [95] investigated the influence of die clearance on sheared edge quality, while Matsuno et al. [96]examined the impact of cutting edge geometry on residual forming capacity. Shih et al. [97] explored the effect of shear rake angle. Standardized experimental approaches are still necessary for a successful inclusion of this effect into the modelling [20].

Frómeta et al. [98] demonstrated the reliability of the Essential Work of Fracture (EWF) methodology in assessing the fracture toughness of CP and DP steel grades to properly address issues associated with edge cracking. More recently, Feistle et al. [99] stated that it is imperative to use the GFLC to correctly take into account the introduced pre-forming in the formation of the Edge-Crack-Sensitivity-Factor and to be able to use it in the finite element simulation. The pre-deformation should not be considered as pre-damage, which would increase the edge crack sensitivity when using multi-phase steels [99].

4.5. Bending influence

The traditional FLD determination is based on in plane deformation without taking a bending component into account [100,101]. Hence, in situations where the strains on the outer layer surpass those of the membrane layers, the FLD fails to provide proper forming behavior results [20,102]. Although known and not being a concern for many years, the advances towards materials with higher strength/lower ductility as well as lower post-necking plasticity such as the AHSS accentuates the impact of bending on forming behavior.

Under superimposed bending, the FLD usually gives conservative results. This behavior arises from the fact that bending introduces a partially compressive stress state throughout the material thickness, which stabilizes the sheet and delays or prevents the onset of necking. When a material is subjected to bending in addition to other loading conditions, the conventional FLD often provides inaccurate results. In these contexts, the materials exhibit a higher deformation capacity than indicated by the FLD, as major strain values located above the safety limit, are in reality, below the failure limit [20,103]. Research conducted by Kitting et al. [104] involved stretch-bending tests conducted on DP800 and CP800 steels under plane-strain and biaxial stretching conditions. The results revealed a significant underestimation of the forming potential of these materials by the traditional FLD in combined stretching and bending operations. However, it should be noted that under severe stretch bending conditions, materials may experience surface failure, even if the strains remain below the FLC [20].

Kitting et al. [103] proposed a phenomenological approach that offers the ability to predict stretch-bending formability by using the conventional FLC prediction method. Neuhauser et al. [101] conducted angular stretch bending tests on DP600 steel and presented a methodology for quantitatively assess the formability of steel sheets subjected to stretching with superimposed bending. They also demonstrated how to integrate these findings into the existing FLD framework. Despite the well-known relationship between instability and bending radius, this area of research has not yet reached full industrial maturity. Its development is crucial for achieving a comprehensive analysis of formability [20].

4.6. Ironing and coining

Through thickness compression in negative hydrostatic stress states is encountered during ironing and coining operations. As a result, the FLD is unable to correctly predict the component criticality, indicating premature failure [94]. Borrego et al. [105] evaluated the hole flanging operation through the traditional FLD analysis. Their findings revealed that strain path curves for successful hole flanged parts near the formability limit exceeded the limits defined by the FLC, yet no failure was observed. Alternatively, the Limiting Forming Ratio (LFR) is a more suitable measure of formability in conventional hole flanging operations. This parameter directly quantifies the maximum stretching of the material at the hole tip [105].

5. Constitutive Material Models

In sheet metal forming applications, an accurate representation of the material’s plastic behavior is crucial for obtaining reliable results in finite element simulations [106]. The accuracy of finite element simulation strongly depends on the extent to which a material constitutive model can characterize the real material properties [107]. The plastic behavior of a material in a general stress state is described by three elements [108]:

• Yield Criterion;
• Associated Flow Rule;
• Hardening Rule.

Figure 10 presents an overview of the required material model parameterization concerning sheet metal forming simulation. Currently, constitutive models with anisotropic yield criteria based on associated flow rule and isotropic hardening rule are extensively used in this field. However, the increased adoption of lightweight materials like AHSS and the growing demand for accurate simulation of forming behavior necessitate the development of more accurate constitutive material models. The selection of appropriate yield criteria and hardening rules as well as the fitting of these models through suitable material testing procedures is of significant importance in sheet forming simulations. Even though, the adoption of most advanced and comprehensive material model formulations is not always favored for predicting forming behavior. Factors such as cost, ease of implementation and accuracy must be considered. Cost concerns mainly pertain to the mechanical testing required for model calibration which explains the widespread utilization of Hill’s 1948 model. Furthermore, the chosen model should gain industrial acceptance, with emphasis on computational efficiency, ease of implementation in numerical simulation codes and user-friendliness. Ultimately, the accuracy of the model’s predictions concerning yield locus, uniaxial yield stress and the coefficient of plastic anisotropy in uniaxial conditions must be guaranteed. Banabic [108] offers a comprehensive and historically detailed analysis of various material models, while a recent review paper written by Hou et al. [43] outlines the latest trends in material models addressing sheet metal forming applications.

5.1. Yield criterion

A yield criterion serves as a defining boundary that marks the transition from elasticity to plastic flow in response to various stress combinations [109]. Thus, accurately characterizing the yield surface is of significance, as sheet metal forming processes typically operate within the plastic region. Although high-strength sheets generally exhibit average r-values of 0.8-1.0, certain steels exhibit significant planar anisotropy r-value [110]. Anisotropic yield criteria are favored in material models for sheet metal forming operations during the production stage. The yield locus depicted in Figure 11a) represents different forming states. The precision of the yield locus is a key indicator of the accuracy of failure/fracture prediction, making it vital for the sheet metal forming community to establish highly precise yield criteria. This precision is essential for effectively anticipating failures/fractures and preventing costly production mishaps [111].

Isotropic yield functions such as the Tresca [114], the Huber-Mises-Hencky [115,116,117], the Drucker [118] and the Hershey [119] yield criteria are not commonly employed in the forming industry. The Hill 1948 [120] anisotropic quadratic function remains widely used due to its ease of parameter determination from uniaxial tension experiments in three directions. However, it exhibits limited capability in accurately predicting anisotropy in sheet metals [110]. Several studies have concluded that only non-quadratic functions can accurately capture the influence of anisotropy on the yield surface [108]. Therefore, the non-quadratic Hill 1979 [121] yield function was proposed, even though it is applicable only when the principal stress directions align with the orthotropic axes [108]. In 1990, Hill 1990 [122] provided an excellent description of the anisotropy coefficient and uniaxial yield stress variation within the sheet plane. Afterwards, Hill 1993 [123] offered a more user friendly expression, albeit with limited applicability. The emergence of lightweight materials, such as AHSS has shifted the anisotropic yield criteria towards more advanced formulations in the 21st century. Noteworthy examples include Barlat [124,125,126,127], Banabic-Balan-Comsa (BBC) [128,129], Cazacu-Barlat [130] and Vegter [131,132] which offer improved descriptions of the yield surface and closely track the planar variation of uniaxial yield stress and the coefficient of plastic anisotropy [108]. Figure 11b) depicts the yield locus of various yield criteria for a DP980 steel. The Barlat yld2000 model was used by Gutierrez et al. [133] for the analytical determination of forming limit curves in two third generation steels. While these yield criteria successfully accounted steels anisotropy, they underestimated the yield stress under shear conditions. Yoshida et al. [110] developed a more advanced yield function that provides a more accurate description of anisotropy influence.

Table 2 summarizes the most commonly used yield criteria for sheet metal forming simulations and the corresponding mechanical property parameters. All models require measurement of the material’s plastic anisotropy, serving as an indicator of sheet metal formability. Hill 1948 only requires standardized uniaxial tension along the rolling, diagonal and transverse direction. Additionally, Barlat yld2000 and BBC2005 require biaxial tension. In its initial formulation the Vegter material model also necessitated plane strain and shear tests in order to calibrate 17 parameters. However, Vegter 2017 [131] only requires standardized tensile tests.

The research conducted by Gösling [134] demonstrated that Barlat yld2000 exhibits superior accuracy in predicting erring compared to the Hill 1948 material model, which is influenced by the yield strength in the diagonal direction. A comparative analyses of three different yield criterion, namely Hill 1948, Barlat yld2000 and Barlat yld2000-Var [135] was performed by Cheng et al. [136]. Among these, the third mentioned model provided more precise results addressing the hardening behavior, the forming limit curve, especially under biaxial tension stress state, and the stretch formability.

Regarding the future trends in yield criteria for sheet metal forming applications of AHSS, several directions can be identified. Firstly, there is a need for the development of new yield criteria aimed at improving flexibility and reducing the number of parameters and the required experimental tests [107]. Secondly, standardized shear and plane strain experimental tests need to be developed. Lastly, there is a growing interest in the development of crystal plasticity material models. According to Hou et al. [43], crystal plasticity material models are advanced multi-scale modelling approaches that capture the deformation behavior of polycrystalline materials by considering slip, twinning and phase transformation processes.

5.2. Hardening rule

The hardening laws play a crucial role in describing how the yield surface evolves under plastic deformation and various loading conditions. These laws can be categorized in to four main families [137]: (i) isotropic; (ii) kinematic, (iii) anisotropic; for example distortional hardening and (iv) combined. Isotropic hardening implies that the yield surface expands without changing its shape as plastic strain increases. Kinematic hardening maintains the shape of the yield surface while translating it in the direction of strain. Distortional hardening refers to the distortion of the yield surface shape during plastic deformation [138]. Combined hardening models, such as kinematic + isotropic or kinematic + distortional hardening have also been proposed. Figure 12 illustrates the representation of the yield locus for these four different families of hardening rules.

Typically, the isotropic hardening is selected to reproduce the hardening behavior of metals in sheet metal forming processes, due to its simplicity. Commonly used isotropic hardening models include Ludwik [141], Swift law [142] and the Voce law [143]. However, experimental observations indicate that the hardening behavior differs under various loading conditions and directions, even under proportional loading conditions [43,141,144,145]. Therefore, the anisotropic hardening behavior can not be described by the anisotropic yield functions under isotropic hardening [145]. Moreover, cyclic loading and unloading occur frequently during sheet flow through draw beads and die radii. Under such loading and continuous reverse loading conditions, the stress-strain relationship is significantly influenced by the Bauschinger effect [107], which refers to a material property of stress/strain characteristics under non-proportional loadings [43]. Isotropic hardening laws fail to accurately capture the Bauschinger effect. The aforementioned limitations result in inaccuracies representing the forming behavior, particularly in cases where springback, tearing and earing occur. Kinematic hardening models have been proposed to address the Bauschinger effect, including Prager’s model [146], Armstrong and Frederick model [147], Chaboche’s model [148] and the model of Yoshida and Uemori [149]. Rosenschon and Merklein [150] have stated that a kinematic hardening law leads to a significantly better approximation of the cyclic hardening behavior rather than a pure isotropic formulation. Distortional hardening models such as the Homogeneous Yield Function-based Anisotropic Hardening model (HAH) offer an advantage in that the expression and parameter calibration of isotropic and anisotropic hardening terms are independent from each other [43]. Zhu et al. [139] demonstrated that the HAH-2d model reasonably captures the Bauschinger effect, permanent softening effect in reverse loading and latent hardening effect in cross loading., while the predictive accuracy for cross-loading softening requires further improvement. Additionally, combined kinematic + distortional hardening models have been proposed to study the material response under non-proportional loadings [151,152,153].

Butuc et al. [154] performed a comparative analysis of different hardening models in predicting the FLC of a DC06 steel. The right-hand side was well predicted by all evaluated models. The left-hand side was overestimated by the Swift law and combined Swift/Armstrong-Frederik law. Dizaji et al. [155] adopted a similar approach but focused on the influence of different hardening models on ductile fracture criteria. Their results showed that the combined Chaboche-Zeigler hardening model [156] and the kinematic Zeigler-Prager equation [146] provided more accurate predictions for the fracture location in square cup drawing using the Brozzo et al. [157] and Cockcroft and Latham [158] ductile fracture criterion, respectively. Nevertheless, the isotropic hardening model delivered best results than the kinematic hardening model for the Ayada et al. [159] ductile fracture criterion.

5.3. Flow rule

Yield criteria can be formulated with associated flow rule (AFR) and non-associated flow rule (NAFR). AFR models such as Barlat yld2000 capture yield stress and strain rate ratio using a linear combination of two functions, suitable for high-strength steels. However, third generation AHSS such as Q&P steels, exhibit a strength differential effect that requires more advanced AFR yield functions leading to higher complex formulations and computational costs [43]. NAFR offers a simpler approach for modelling anisotropic yield and plastic flow in metallic materials. It uses separate yield and plastic potential functions with simpler forms and analytical parameter calibration. Stoughton’s non-associative flow model [160], based on Hill 1948 function, accurately predicts yield stress in strongly anisotropic materials. NAFR’s usage and applications are increasing as more research investigates its effectiveness [161,162,163,164,165].

6. Material characterization methods

In the realm of sheet metal forming processes, the material undergoes various types of loading. Accurate experimentation plays a significant role in accurately calibrating models which predict crack initiation on macroscopic features. The fracture models rely on the correct assessment of the fracture strain. These fracture models require experimental testing under different stress states to fully characterize fracture occurrence. Different sample geometries, test configurations and/or loading paths generate different levels of stress states which significant influence the fracture typology. In order to represent any stress state, two dimensionless parameters, namely the Lode Angle and the Stress Triaxiality need to be represented. The loading path leading to fracture, represented in terms of the equivalent plastic strain evolution as a function of the Stress Triaxiality, $\eta$ and Lode Angle parameter, $\theta$ significantly influences the precise assessment of fracture initiation location. Therefore, a wide range of stress states need to be covered by employing diverse experimental tests. Lian et al. [166] proposed a hybrid ductile damage modelling with a hybrid ductile experimental-numerical parameters calibration procedure specifically tailored to account for the damage characteristics of modern high-strength steels. Figure 13 presents the methodology followed by the authors in order to calibrate material and fracture parameters based on experimental tests. Although the knowledge of the strain to fracture is the main ingredient of damage models [167], they used a phenomenological model to characterize the onset of damage: the crack initiation strain, rather than the fracture strain.

Currently, most material characterization methods are based under the assumption of proportional loading histories. Figure 14 presents some conventional experimental techniques employed for the assessment of relevant data such as fracture strain, $\eta$ and $\theta$. However, the achievement of proportional loading conditions at the specific material point where fracture initiates proves to be challenging. Thes sample geometry and test configuration often induce a certain level of deviation from linearity during ductile fracture. Furthermore, non-linear deformation paths commonly arise in multi-forming operations. As such, method able to characterize the stress state under non-proportional loading are needed. The current trends on behalf of material characterization methods rely on industrial standardization of fracture limits determination and the enhancement of material assessment techniques for both proportional and non-proportional loading histories. Table 3 provides a review on the material characterization methods employed to calibrate fracture/damage models parameters.

6.1. Proportional loading

6.1.1. Pure shear ($\eta =0;\overline{\theta }=0$)

The accurate assessment of the shear fracture limit is currently one of the most questioned subjects in sheet metal fracture modelling. Over the years, a large amount of sample geometries and methods have been proposed to address this issue. One hand, the ability of metallic materials to undergo shear deformation facilitates the attainment of elevated levels of plastic strains. On the other hand, makes it more difficult to achieve pure shear deformation paths at the crack initiation location.

Miyauchi [179] introduced the initial planar shear test documented in literature. The specimen comprises two symmetrically sheared regions and it is subjected to an axial load [180]. Rauch and G’Sell [181] adopted a modified version of this test, which was further refined by Genevois [182]. In Genevois’s test, the rotation of the specimen under the influence of shearing forces is impeded by a highly rigid frame. The specimens employed in this test possess a simple rectangular shape which deforms into a parallelogram. Bouvier et al. [180] conducted an optimization of the mechanical test, focusing on the impact of sample geometry and setup configuration on stress and strain uniformity. Subsequently, the ASTM-B831 [183] standard test method was established for shear testing, proposing a novel sample geometry and loading conditions. Alternative shear specimen designs have also been proposed such as the symmetric double-shear specimen with two identical gage sections [184]. However, the application of these specimens is limited due to the interference of edge effects when calculating shear stress [43]. It is important to note that the aforementioned approaches were originally devised for the determination of plastic hardening and are not generally suitable for characterizing formability limits caused by fracture [185]. Furthermore, they are affected by fracture occurrence at the free edges.

The primary prerequisite for determining the fracture limit strain in shear tests is the ability to achieve pure shear deformation paths, thereby ensuring a consistent zero triaxiality throughout the entirety of the experiment. To achieve this goal, it is essential to assure homogeneous stress and strain distribution while also preventing crack initiation at free edges, which can prove to be challenging, as most of the materials exhibit higher strain to fracture under uniaxial tension rather than pure shear. The importance of it relies on the crack initiation likelihood occurrence on higher stress triaxialities. Iosipescu [186] was among the first to endeavour the development of a pure shear test with constant triaxiality evolution. The subsequent ASTM D5379 [187] standard was formulated based on Iosipescu’s work and it is predominantly employed for composite materials. However, both approaches fail to adequately address the issue of edge cracking, similar to the specimens proposed by Bao and Wierzbicki [188] and Shouler and Allwood [189]. Dunand and Mohr [190] introduced a novel specimen incorporating localized sheet thickness reduction. Nevertheless, the initiation of fracture in this specimen could be reliably guaranteed to commence precisely at the centre of the sheared zone, thus compromising the achievement of triaxiality constancy. Figure 15 illustrates the four aforementioned proposed specimens.

Significant developments have been made in addressing the need for more accurate pure shear strain deformation paths through adapted in-plane torsion tests. These tests offer the capacity to achieve maximum strains up to 1.0 [191]. In the context of in-plane torsion testing, a circular sheet sample is clamped at both the outer rim and in the centre. By inducing planar rotation of the outer fixture, a simple shear is applied to the circular region situated between the clamps [192]. Yin et al. [185] determined the fracture strains for DP1000 using an in-plane torsion test, with a grooved specimen free of slits. This test configuration facilitated fracture propagation under constant zero triaxiality and effectively mitigated incipient cracking arising from edge effects. Grolleau et al. [193] proposed an alternative clamping system featuring a highly rigid device for in-plane torsion testing using grooved specimens. The system allowed for full optical access to the specimen’s gage section, enabling DIC measurements. It provided an almost perfectly proportional loading path for simple shear loading making it suitable for cyclic loading as well as testing ultra-high strength steels. Moreover, Roth and Mohr [170] recently reported a geometry shape optimization procedure which effectively addressed the strain to fracture under pure shear. The smiley shear specimen, depicted in Figure 16a), is distinguished by a series of notches that define the contours of the shear gage sections. As observed in Figure 16d) a more linear strain path evolution is achieved.

6.1.2. Uniaxial tension ($\eta =1/3;\overline{\theta }=1$)

ASTME-8 outlines the procedure for conducting uniaxial tensile tests on metals at room temperature, enabling the determination of mechanical properties [194]. The amount of plastic deformation achievable through the uniaxial tensile test is typically limited to around 20-30%, which is noticeably lower when compared to other stress states such as shear. Uniaxial tensile tests are usually performed in three directions: rolling, transverse and diagonal, to investigate material anisotropy. The mechanical properties obtained from uniaxial tensile testing include Young Modulus, Yield Strength, Tensile Strength, Uniform and Total Elongation and r-values. When measuring the fracture strain under uniaxial tension using the conventional flat dog-bone specimen, the stress state undergoes changes during the loading process. To measure the limit strains under uniaxial tension, a tensile specimen with a central hole can be used instead, which offer the advantage of maintaining a more constant linear strain path until fracture occurs near the hole boundary [195]. However, careful selection of the hole manufacturing process as well as hole diameter dimension is relevant, as inadequate choices can lead to significant errors in the estimated strains to fracture for uniaxial tension [170]. Conversely, the work conducted by Behrens et al. [77] reveals that the central hole specimen yields non-linear strain path histories.

6.1.3. Notch tension

In order to encompass a wider spectrum of stress triaxiality values along the fracture locus, the utilization of notch tension is frequently employed. By adjusting the notch angle and size on notched flat specimens, a diverse range of stress triaxiality conditions can be attained [196]. The presence of notches ensures the development of a localized neck at the centre of the specimen, perpendicular to the loading axis [195]. For instance, notch tension tests have been proposed to calibrate the parameters of the extended Modified Mohr Coulomb model (eMMC) [171,173].

6.1.4. Plane strain tension ($\eta =1/\sqrt{3};\overline{\theta }=0$)

Fracture occurrence during plane strain tension is primarily attributed to the load carrying capacity in the out-of-plane shear direction. Various approaches have been employed to determine the strain limits under plane strain tension. Jia and Bai [173] utilized uniform plane strain specimens as well as butterfly-shaped specimens. Santos et al. [197], on the other hand, employed double-notched specimens to assess the deformation mode under plane strain tension. However, it should be noted that with these test configurations, the stress state leading to fracture may not remain constant. The point at which plane strain tension limit is typically reached is on the material surface during bending operations. To address this, Roth and Mohr [170] introduced a modified plane strain bending configuration able to deliver proportional loading histories (Figure 16d). Their V-bending experiment involves placing a rectangular sheet material coupon (Figure 16b) on top of two parallel rollers. The sheet specimen is then subjected to loading by the two parallel rollers through a thin, knife-like tool which remains stationary, resulting in the formation of a shar V-bend.

6.1.5. Equi-biaxial tension ($\eta =\frac{2}{3};\overline{\theta }=-1$)

Fracture testing of sheet materials under equi-biaxial tension is well standardized [170]. Usually either bulge test or punch test is used. The hydraulic bulge test involves the deformation of a flat specimen through the application of hydraulic pressure. The outer edge of the sample is firmly clamped with a circular die, while the internal area is subjected to increasing fluid pressure during the test, causing deformation until fracture occurs [168]. The ISO 16808 standard [198] specifies a method for determining the biaxial stress-strain curve using an optical measuring system. The bulge test allows for a much larger achievable strain prior to necking, owing to the nearly equi-biaxial stress state, compared to the uniaxial tensile test. The punch test, on the other hand, utilizes a miniature disc sample subjected to compressive force applied through a hemispherical punch indenter on the upper surface of the disc [199]. The small punch deformation test procedure for metallic materials is covered by the standard ASTM E3205-20 [200]. The configuration used to perform both tests is schematically presented in Figure 14. The punch test offers the advantage of avoiding the evacuation of excess fluid after fracture, although the presence of tool friction can influence the experimental results. Roth and Mohr [170] proposed an adaptive punch testing procedure, wherein the distinctive feature is that the punch remains stationary throughout the experiment while the die and clamping plate move downwards. The specimen geometry is presented in Figure 16c). This approach reduces the required focal depth of the DIC camera system, enabling shorter object distances and ultimately enhancing the spatial resolution of the acquired surface strain fields. Furthermore, the location of fracture initiation remains fixed throughout the entire experiment.

6.1.6. Butterfly test

Mohr and Henn [201] devised a butterfly-shaped specimen that facilitated the examination of fracture behaviour under low and intermediate stress triaxialities, ranging from 0.0 and 0.6. The generation of different stress triaxialities can be accomplished without altering the design or dimensions of the specimen, but rather by adjusting its orientation. The experimental procedure is conducted using a biaxial testing apparatus comprising a fixed portion, a vertically movable portion and the corresponding specimen. By adjusting the angle between the specimen axis and the vertical axis, different stress triaxialities can be achieved. The gage section of the specimen exhibits a uniformly reduced sheet thickness compared to its boundaries and it is designed to induce crack formation at the centre of the specimen. Bai [202] introduced modifications to the design of the gage section, resulting in a non-uniform sheet thickness within the gage section.

Dunand and Mohr [190] further improved the initial specimen design attempting to increase homogeneity of the stress state and strain fields at the specimen centre. Additionally, the configuration of the testing setup enabled the application of combined normal and shear loads, as illustrated in Figure 17a). The experimental investigation focused on TRIP780 steel, encompassing a range of loads from pure shear to transverse plane strain tension. The findings revealed that under shear loading, both stress triaxiality and lode angle parameter exhibited relative constancy. However, a linear increase in stress triaxiality after the onset of through-thickness necking was observed for tension loading. Moreover, it was noticed the highly sensitivity to imperfections in the specimen geometry.

Despite the improvements made to the initial butterfly specimen design proposed by Mohr and Henn [201], achieving uniform stress and strain distribution throughout the entire experiment, especially under shear-dominated loading conditions, remained challenging. Therefore, Peshekhodov [203] further enhanced the butterfly design with the aim of ensuring strain localization at the specimen’s gage section regardless of the applied load direction. The new specimen geometry was divided was divided into three areas, as illustrated in Figure 17b). The comparative FEM analysis between this new geometry and the one proposed by Bai [202] demonstrated an increased likelihood of fracture initiation at the desired central area under arbitrary stress states. Furthermore, a more precise definition of the stress state could be achieved by specifying the load application angle. The implemented apparatus, based on a single material characterization technique, facilitated the generation of a wide range of stress triaxialities, from pure shear to plane strain tension. This approach accurately modelled the fracture behaviour of a DP600 steel using the CrachFEM and MMC fracture models.

Behrens et al. [77] adopted a similar experimental and numerical approach to obtain fracture data for a CP800 steel and calibrate the MMC fracture model. Similarly, the butterfly specimen developed by Peshekhodov [203] demonstrated fracture initiation at the specimen centre under various stress triaxialities, ranging from shear to uniaxial tension. Even though, the butterfly specimen did not deliver entirely proportional strain paths, improvements compared to other experimental tests are undeniable. In order to solve the problem of non-linear strain path evolution, Stockburger et al. [204] proposed an innovative adaptive experimental configuration, as shown in Figure 17c). The key feature if this setup was the continuous real-time adjustment of the loading angle to maintain constant stress states. To achieve this, a Phyton script was incorporated into the numerical simulation, compensating for deviations in stress triaxiality exceeding 0.03. A comparative analysis between the calibration of the MMC fracture model using the new adaptive test methodology and the fixed loading angle approach revealed significant enhancements in achieving linear strain paths. Furthermore, as the strain paths became more consistent, the equivalent plastic strain at fracture increased.

6.2. Non-proportional loading

In the context of sheet metal forming operations, the presence of non-linear strain path histories at fracture initiation sites has been a frequent phenomenon. Hence, attention within the sheet metal forming research community has also been directed towards enhancing experimental methodologies able of generating controlled non-linear strain paths, further improving the accuracy of material fracture and failure modelling. Aisvaran [205] identified the main test procedures to generate non-linear strain paths: (i) Nakazima test with modified punch geometry; (ii) two step drawing; (iii) in-plane biaxial test with cruciform specimen; (iv) Bulge test with stepped dies; (v) cruciform specimen in conventional Nakazima test setup. Table 4 summarizes the main approaches for generating NLSPs tested on AHSS.

The Nakazima test with modified punch geometry offers the advantage of achieving non-linear strain paths in a single step procedure. Saxena et al. [78] designed five novel punch geometries to replace the conventional hemispherical punch in order to change the strain path during deformation in the necking initiation site. The Nakazima test, employing a cruciform specimen serves as an alternative technique to generate NLSPs. An innovative procedure proposed by Jocham et al. [210] incorporates a modified blank-holder with adjustable draw bead height, simplifying the need for a single specimen geometry to produce the FLC [210], as opposed to the two step drawing approach requiring varying specimen geometry [206]. He et al. [211] introduced an innovative test procedure designed to facilitate continuous non-linear biaxial tensile deformations of sheet metals via bulging with stepped dies. This bulging process lead into a distinctive change in curvature radius relative to the bulging height and as the difference between successive elliptical sections increases, the presence of NLSPs becomes more pronounced [111]. The main distinctive feature of the in-plane biaxial test with cruciform specimen is that it provides a frictionless way of investigating the strain path change in a cruciform specimen. The design of this specimen continues to undergo refinement within the scientific community, yet standardization has not been established.

7. Damage and fracture models

The incorporation of damage or fracture models holds substantial significance within the realm of sheet metal forming simulation, enabling the anticipation of fracture occurrences in industrial contexts as well as facilitating a comprehensive understanding of damage behaviour. These models can be systematically classified into two principal categories, namely coupled and uncoupled, each grounded either in micromechanical or phenomenological foundations. Figure 18 provides an overview of the classification of ductile fracture models and showcases various models within their respective family backgrounds. The fundamental divergence between coupled and uncoupled models resides in the former’s integration of damage mechanisms alongside elastic and/or plastic responses, as opposed to the latter, which distinctly segregates damage evolution from material behaviour. Furthermore, the domain of coupled models can be further bifurcated into micro-mechanical based models and continuum damage models (CDM). The former considers the deterioration of material during the loading process, attributed to damage propagation, while the later finds its underpinning in a thermodynamic framework and is formulated upon a phenomenological basis.

Damage and fracture model or fracture criteria operate by employing an accumulating damage parameter, which serves as a fundamental basis for assessing fracture. By using equation (1), a general fracture model can be expressed.

$$D{\int }_0^{{\epsilon }_f}f\left(\sigma ,\epsilon ,\dot{\epsilon },T,\dots \right)\overline{{d\epsilon }^p}\le 1$$

(1)

where ${\epsilon }_f$ is the equivalent plastic strain at fracture, $\overline{{d\epsilon }^p}$ is the equivalent plastic strain increment, and $f$ is a scaling function which depends on internal variables such as stress state $\sigma$, strain rate $\dot{\epsilon }$, and temperature $T$. $D$ is seen as a damage indicator [28].

In order to accurately depict the normalized stress state, it is necessary to consider both stress triaxiality $\eta$, as well as lode angle parameter $\overline{\theta }$. It is well established that the strain to localization decreases when the stress triaxiality increases. Conversely, a lower stress triaxiality value results in an increased sensitivity to the lode angle parameter [212]. The stress triaxiality is very important to characterize the stress state of metallic materials and consists in the relative degree of hydrostatic stress in a given stress state which is defined according to equation (2).

$$\eta =\frac{{\sigma }_m}{{\sigma }_{eq}}=\frac{\frac{1}{3}({\sigma }_1+{\sigma }_2+{\sigma }_3)}{\sqrt{\frac{{({\sigma }_1-{\sigma }_2)}^2+{({\sigma }_2-{\sigma }_3)}^2+{({\sigma }_1-{\sigma }_3)}^2}{2}}}$$

(2)

where ${\sigma }_m$ corresponds to the hydrostatic stress, while ${\sigma }_{eq}$ refers to the Von Mises equivalent stress. ${\sigma }_1$, ${\sigma }_2$ and ${\sigma }_3$ are the principal normal stresses. $\eta$ controls physically the size of a given void, i.e., controls the micro-void growth during ductile fracture. For isotropic materials under plane stress, the triaxiality alone is enough to define any possible stress state relevant for fracture characterization [213]. On the other hand, $\overline{\theta }$ controls the shape of a given void [29] and is required to model three dimensional stress states. To determine the lode angle parameter, it is worth mentioning the lode angle, $\theta$ which is defined through the normalized third stress invariant, $\xi$:

$$\theta =\frac{1}{3}\arccos \left(\xi \right)=\frac{1}{3}\arccos {\left(\frac{r}{q}\right)}^3=\frac{1}{3}\arccos {\left(\frac{{(\frac{27}{2}.J_3)}^{1/3}}{\sqrt{3J_2}}\right)}^3$$

(3)

where $r$ and $q$ are two stress state parameters and $J_2$ as well as $J_3$ are the second and third deviatoric stress invariants, respectively. Therefore, $\overline{\theta }$ is defined on account of the given transformation:

$$\overline{\theta }=1-\frac{6\theta }{\pi },with\overline{\theta }ϵ[-1;1]$$

(4)

This section presents a brief historical review on ductile fracture and damage models. For a more comprehensive overview, the reader should refer to the work conducted by Cao [212] as well as Tekkaya et al. [29].

7.1. Comparison between different fracture and damage models

There are a large number of ductile fracture models. Therefore, it can be difficult to choose the most appropriate one to employ in a particular sheet metal forming process. Table 5 presents a brief and simplified qualitative comparison between the three main families of ductile fracture models.

For the fracture or damage model to garner acceptance within real industrial forming applications, a wise equilibrium must be attained between ease of implementation and forming behavior accuracy. In this realm, the investment time in both model formulation as well as its implementation in finite element software must be as low as possible. On the other hand, the accuracy of the model should be as high as feasible to correctly predict the component forming behavior. Generally, these parameters work in an inverse manner.

In these contexts, phenomenological uncoupled models are commonly preferred, e.g., where the components exhibit relatively straightforward geometries, characterized predominantly by linear strain paths. These models suit to optimize process parameters aiming of minimizing damage occurrence. The major advantage of such models consists on their ease of implementation and use. The simpler calibration procedure as well as the comparatively less number of parameters contribute to their quicker implementation. Additionally, since damage softening effects are not considered, mesh dependency is avoided, facilitating the simulation process: less CPU time and more convergence stability. However, it is imperative to acknowledge that the principal limitation of this category of models pertains to their applicability in scenarios involving complex loading paths outside the identification zone, as well as those relating to large plastic deformation. Thus, special care should be taken for their application outside the identification domain, due to their phenomenological foundation [212].

In many industrial applications, the loading path frequently exhibits intricate characteristics. Therefore, coupled models may be seen as an appropriate choice, since they account for softening effect due to damage accumulation, as opposed to uncoupled models. Indeed, the coupled models predictive damage accuracy is significantly higher. These models can capture continuous degradation due to damage. Although important when the material experiences very large plastic deformation in multi-step forming processes, the complexity of the simulation process is unavoidably raised, i.e., softening effect introduces mesh dependency. As a result, finer mesh induces faster damage accumulation. In order to mitigate mesh dependency, two validated approaches can be pursued: (i) integration of non-local methodologies, which may be grounded in either integral formulations or implicit/explicit formulations, requiring careful calibration; (ii) delineation of damage model parameters tailored to specific element dimensions, ensuring that during computations, the element sizes remain beneath this predetermined threshold [29]. This mesh dependency results in an increased CPU time in comparison to uncoupled models. However, to handle the complex strain path the material may undergoes, it is required the use of phenomenological models [29]. No purely micro-mechanical model has been shown to be able of capturing both damage mechanisms at high and low stress triaxialities in real multi-forming operations, due to their construction and current computational efficiency [212]. Thus, the micro-mechanical based models are discarded for industrial forming processes implementation. With a significantly high physical background, their complexity is high. The CPU time as well as the calibration procedure dissuades its implementation. Conversely, the CDM models require comparatively less effort in the calibration process as a smaller number of damage parameters are used. Additionally, these models couple elasticity with damage evolution.

Therefore, each family of models has their own advantages and drawbacks. The effort one is willing to take in model accuracy, neglecting simplicity, or vice versa, is the main factor in model selection. In this context, the advanced uncoupled phenomenological models are very popular in the sheet forming industry. Even though, their comparatively diminished damage prediction accuracy, they have been used and demonstrated satisfactory results even for applications involving complex stress states. Moreover, these ductile fracture models facilitate the incorporation of strain rate dependencies through accurately constructed strain-rate dependent terms. Among all phenomenological fracture models, it is worth mentioning the modified Mohr-Coulomb (MMC) [214], the all-strain based MMC (eMMC) [173] and the Hosford-Coulomb model (HC) [175]. All of them showed good performance modelling sheet metal forming fracture of AHSS. Table 6 summarizes the implementation of ductile fracture models in some representative studies reported in literature.

Alternatively, the adoption of ‘modular’ models has been progressively gaining attention. Even though coupled models generally exhibit better predictive accuracy, they encounter limitations in addressing conditions characterized by low stress triaxialities and notably shear-dominated loading conditions. Therefore, the conceptual framework of the ‘modular’ models is to upgrade these models by adding terms which increase the predictive accuracy for specific situation (e.g., the Lode dependent term; a strain threshold for damage accumulation which depends on the stress triaxiality). When thoughtfully formulated, these diverse models may converge towards consistent predictions concerning damage location [212].

7.2. Phenomenological fracture models

In these models, fracture is assumed to occur when a critical stress or strain rate is reached. Damage variable is just a ‘warning’ for the approach of fracture [212]. Although most of the uncoupled models encountered in the literature have a phenomenological background, the first developed fracture models were rooted in physical principles, incorporating considerations of void growth. Examples of such models are the McClintock [220] and Rice-Tracey [221] models.

In the realm of phenomenological fracture criteria, the concept of damage is approached as an integration of a stress-based function along the loading trajectory. Initially, early research effort solely focused on considering either stress triaxiality, e.g., the Oyane model [222], or principal stress, e.g., Cockcroft-Latham model [158]. However, it was observed that when subjected to complex loading paths, the predictive accuracy of these models in simulating forming behavior is diminished. The primarily stems from the fact that the critical damage threshold is contingent upon the specific loading trajectory. Recognizing this limitation, Bao and Wierzbicki [223] highlighted the necessity of incorporating the lode angle parameter dependence into these criteria to establish a unique definition of the normalized stress state. Thus, the recently proposed models, named advanced uncoupled phenomenological models comprehensively account for both stress triaxiality and lode angle influence. These criteria usually define a fracture locus where the strain to fracture is a function of stress triaxiality and lode angle [29]. Nevertheless, it is prudent to acknowledge that while such models enhance predictive accuracy within complex stress states, further refinements are necessary to address non-linear strain paths. In such cases, the equivalent strain to fracture no longer represents the actual strain at fracture, but rather serves as a weighting function which accounts for the stress state [212].

7.2.1. Modified Mohr-Coulomb model (MMC)

The Modified Mohr-Coulomb model [171] was proposed by Bai and Wierzbicki. The model was formulated by the transformation of the stress-based Mohr-Coulomb (MC) criterion to the corresponding strain-based model, where the spherical coordinate system is derived in terms of $\overline{{\epsilon }_f}$, $\eta$ and $\overline{\theta }$. The damage indicator can be formulated as in equation (5).

$$D\left({\overline{\epsilon }}_f\right)={\int }_0^{{\overline{\epsilon }}_p}\frac{d{\overline{\epsilon }}_p}{f\left(\eta \overline{\theta }\right)}$$

(5)

The fracture locus is calculated according to:

$${\overline{{\epsilon }_f}=\left\{\frac{A}{c_2}\left[c_{\theta }^S+\frac{\sqrt{3}}{2-\sqrt{3}}\left(c_{\theta }^{ax}-c_{\theta }^S\right)\left(\sec \left(\frac{\theta \pi }{6}\right)-1\right)\right][\sqrt{\frac{1+c_1^2}{3}}.\cos \left(\frac{\theta \pi }{6}\right)+c_1(\eta +\frac{1}{3}sen\left(\frac{\theta \pi }{6}\right)\left)\right]\right\}}^{\frac{-1}{n}}$$

(6)

There are a total of six parameters ($A$, $n$, $c_1$, $c_2$, $c_{\theta }^s$, $c_{\theta }^c$)that need to be determined. Parameter $A$ and $n$ represent the Swift hardening law parameters, while the remaining parameters are material constants. These constants are determined through a minimum of three fracture experimental tests conducted at different stress triaxiality ranges.

The engineering applicability of the MMC fracture model for a TRIP690 steel was demonstrated by Li et al. [89]. They accurately predicted the fracture location and magnitude of punch travel for two deep-drawing operation setups: square punch and circular punch. Furthermore, the results indicated that the MMC model outperformed the FLC criterion in dealing with NLSPs and shear-induced fracture. Nevertheless, Behrens et al. [216] highlighted the significance of considering no damage accumulation effects for stress triaxialities below -0.33 in the MMC formulation. Their study detected premature shear fracture initiation of a DP600 steel workpiece in a deep drawing operation without the aforementioned condition. Luo and Wierzbicki [215] investigated the accuracy of the MMC model in predicting fracture behavior in stretch forming simulator tests. The comparative evaluation between different element types revealed that both 3D and shell elements are sufficiently adequate to fit the experimental fracture data. The fracture location/mode, fracture wall stress and fracture wrap angle are particularly well predicted by the 3D element type. Moreover, the MMC model can detect variations in experimental friction conditions. More recently, Fu et al. [224] also proved the effectiveness of the MMC model for predicting ductile fracture both in stretch-bending and Nakazima tests for DP980 and DP1180 steels.

7.2.2. Strain based Modified Mohr-Coulomb (eMMC)

The eMMC model is derived by converting the stress state variables: $\eta$ and $\overline{\theta }$ from the MMC model to the strain state. By considering the stress ratio, $\beta$ and the strain ratio, $\alpha$ of in-plane principal strain increments, the model can be expressed in a space based solely on strain. In this space, the fracture strain is affected by either the strain ratio or a new plastic strain based parameter, $\varphi$ [173].

$$\alpha =\frac{{d\epsilon }_2}{{d\epsilon }_1},\beta =\frac{{\sigma }_2}{{\sigma }_1}$$

(7)

The eMMC formulation is based on the premise of monotonic loading and the plastic flow rule of Mises-Levy. This formulation enables the establishment of a correlation between $\alpha$ and $\beta$:

$$\alpha =\frac{2\beta -1}{2-\beta },\beta =\frac{2\alpha +1}{2+\alpha }$$

(8)

One significant benefit of transforming the MMC into the eMMC lies in the decoupling of the fracture model from plasticity, allowing for the study of fracture independent of plasticity [173]. This is achieved by the introduction of the $\varphi$ parameter, denoted as:

$$\varphi =\frac{180°}{\pi }atan2(-\dot{{\epsilon }_2^p,}\dot{{\epsilon }_1^p)+90°}$$

(9)

The $\varphi$ parameter represents the effect of in-plane anisotropy under all possible loading conditions. Furthermore, it denotes the angle to the positive minor strain direction in the forming fracture limit diagram (FFLD) under proportional loading.

A calibration study conducted by Jia and Bai [173] demonstrated good accuracy in predicting the eMMC fracture locus for a TRIP780 steel. Their study combined experimental and numerical methods and produced precise estimates for material strength, plastic flow behaviour, fracture initiation location and fracture propagation modes.

7.2.3. Hosford-Coulomb (HC)

Mohr and Marcadet [175] proposed the Hosford-Coulomb phenomenological model, which is motivated by micromechanics and assumes that fracture onset is closely linked to the formation of a primary of secondary band of localization. Despite being a phenomenological model, the underlying physical significance of the stress-based localization criteria is preserved. The HC model converts the stress space of the localization criterion into a mixed strain-stress space, incorporating the equivalent plastic strain, stress triaxiality and lode angle parameter. The HC model was specifically developed by Mohr and Marcadet to cater to the distinctive characteristics of AHSS. The accuracy of the model was confirmed by these authors under proportional loading and at low strain rates for three AHSS materials (DP590, DP780, TRIP780).

According to the HC model, fracture occurs when the following conditions are met:

$${\int }_0^{{\overline{\epsilon }}_f}\frac{d{\overline{\epsilon }}_{pl}}{{\overline{\epsilon }}_f^{pr}[\eta ,\overline{\theta }]}=1$$

(10)

In the mixed strain-space, the criterion is defined as in equation (11)

$${\overline{\epsilon }}_f^{pr}=k^{-1}\left[{\sigma }_f\left[\eta \overline{\theta }\right]\right]$$

(11)

where in the Modified Haigh-Westergaard stress space, $\overline{\sigma }$ takes the form of:

$$\overline{\sigma }={\overline{\sigma }}_f[\eta ,\overline{\theta }]=\frac{b}{{\left\{\frac{1}{2}({\left(f_1-f_2\right)}^a+{\left(f_2-f_3\right)}^b+{\left(f_1-f_3\right)}^a)\right\}}^{\frac{1}{a}}+c(2\eta +f_1+f_3)}$$

(12)

The model parameters are referred as $a$, $b$ and $c$.

Compared to the MC criterion, the HC model exhibits greater suitability for handling biaxial loading [175]. Roth and Mohr [225] proposed an empirical extension to the HC model which considers the influence of strain rate on fracture initiation in two AHSS materials (DP590 and TRIP780). This extension introduces an additional scalar parameter to the model formulation to account for strain rate effects. A similar approach was adopted by Erice et al. [176] who verified that the HC model provides a satisfactory approximation of the impact the stress triaxiality and lode angle parameter on three other AHSS (DP980, CP980 and CP1180).

The calibration procedure of the aforementioned studies was addressed through numerical simulations using solid elements. However, when it comes to calibrating the HC using shell elements, it becomes necessary to establish the Domain of Shell-to-Solid Equivalence (DSSE) in order to achieve convergence in the initiation of ductile fracture relative to the mesh size [178]. To address this issue, Pack and Mohr [217] introduced the DSSE-HC model, which combines the DSSE concept with the Hosford-Coulomb model. The DSSSE-HC model can capture fracture initiation without prior necking, as well as through-thickness necking failure. Notably, the DSSE approach offers the advantage of distinguishing between membrane and bending loading conditions, setting it apart from the FLD/FLC criterion. Fracture is assumed to occur when the loading path in a shell element either departs from the DSSE boundary or satisfies the HC fracture criterion. The DSSE-HC model was successfully calibrated for a DP780 steel and applied to accurately predict the displacement to fracture in both stretch bending and ‘V-bending’ operations. However, as highlighted by Pack et al. [178], further enhancements are necessary to mitigate the pronounced mesh size dependence of fracture predictions when employing a combined ductile fracture model and a necking based forming limit criterion.

Furthermore, the efficacy of the HC model in investigating the impact of interrupted loading conditions, as often encountered in servo press forming operations, has been demonstrated by Prasad et al. [177].

7.3. Continuum damage models (CDM)

CDM models are classified as phenomenological due to their foundational development hinging upon a rigorous thermodynamic framework, thereby ensuring a non-negative dissipation [212]. The degradation of material stiffness is addressed through the introduction of an internal scalar damage variable, which interacts with the stress tensor to account for isotropic damage phenomena. To impart the macroscopic response with damage effects, the concept of effective stress is introduced based on the principle of strain equivalence. This notion of effective stress allows a unified integration of damage evolution with the material behavior law, enabling the consideration of damage effects on both elasticity and plasticity [29].

The first contributor towards the formulation of models grounded in thermodynamical assumptions was Kachanov [226]. Subsequent improvements in this theoretical domain were introduced by Chaboche [227] and Lemaitre [228]. Indeed, among CDM models, the Lemaitre damage model is widely used, wherein stress triaxiality and equivalent stress are the pivotal factors governing damage propagation. Notably, Lemaitre’s model underwent enhancements, aiming to tailor their applicability to metal forming applications. These adaptations encompassed the augmentation of the damage dissipation potential, with Lemaitre and Desmorat [229] introducing a strain threshold beyond which damage remain constant [29]. To further enhance model accuracy, particularly at low stress triaxialities, Cao et al. [230] devised the Lemaitre Enhanced Model (LEL) by including a parabolic lode dependent term. The main modifications focused on [212]: (i) accounting for the influence of the Lode Angle parameter on damage evolution; (ii) improving the influence of stress triaxiality on damage threshold; (iii) enhancing the weakening function at low stress triaxialities.

Another notable CDM model is the Generalized Incremental Stress State Dependent Damage Model (GISSMO) [213]. In this model, beside the damage variable, D an instability function, F dependent of the critical strain, ${\epsilon }_{crit}\left(\eta \right)$ is introduced. Once F reaches 1, the coupling between damage and the stress tensor is activated [29].

7.3.1. GISSMO model

Andrade et al. [213] formulated the GISSMO model. The model aims to predict fracture of ductile materials by depicting a fracture curve dependent on the stress triaxiality ratio. The incorporation of a damage accumulative parameter, which is dependent on the plastic strain increment allows the accurate modelling of non-proportional loading histories. In this context, the damage variable with respect to time is given by equation (13)

$$\dot{D}=\frac{n}{{\epsilon }_f\left(\eta \right)}.D^{(1-\frac{1}{n})}.{\dot{\epsilon }}^p$$

(13)

This expression serves as an evolutionary equation which effectively addresses non-proportional loading histories. Additionally, the authors have introduced an instability function denoted as F able to understand when a strain localization arises, and therefore to accelerate material degradation when those situations occur. Using this approach, the eventual non-linearity affecting the onset fracture behavior is taken into consideration. To deal with mesh dependency, a triaxiality-dependent regularization strategy is adopted adjusting the fracture curve to the corresponding element size.

One of the main advantages of the GISSMO model lies in its ease of numerical implementation in comparison to other CDM models. Its feature of allowing coupling with any plastic material model ensures high flexibility. The resulting curves from the GISSMO model are presented in Figure 19. The calibration procedure involves the determination of several parameters, including the fracture curve as a function of stress triaxiality, ${\epsilon }_f\left(\eta \right)$, the damage exponent, $n$, the curve for critical strain, ${\epsilon }_{crit.}\left(\eta \right)$ and the fading exponent, $m$. The GISSMO model is implemented in the commercial finite element code LS-DYNA.

7.4. Micro-mechanics based models

The micro-mechanics based models consider the underlying physical meaning attained to the damage. Hence, the macroscopic feature is a result of micro-mechanical void, void growth and coalescence of neighboring voids. The material deterioration is linked with the fracture occurrence. The micro-mechanics based models can be subdivided into: (i) Gurson-like approach, (ii) non-linear homogenization based approach.

The effect of void growth on material degradation was introduced by Gurson [231]. The widely recognized Gurson-Tvergaard-Needleman (GTN) model [232] represents a further advancement of the Gurson framework. The GTN is a porous plasticity theory, accounting for the evolution of the void volume fraction: interaction, nucleation and final coalescence of voids, where the void growth depends on the trace of the plastic strain rate [29]. Multiple extensions have been proposed to enhance the GTN model’s accuracy, exemplified by the development of the GLD model which accounts for void shape effects. Stress triaxiality and lode angle dependent terms can be integrated into these models. As the main drawback of the micro-mechanics based models, particularly the GTN model is the low accuracy predicting shear-dominated loading cases, the incorporation of these stress state parameters increases the model’s accuracy. The stress triaxiality allows the void nucleation process to take place earlier at high stress triaxialities [212]. Xue [233] and Nahshon and Hutchinson [234] added a lode-dependent term in order to more accurately model the void growth [29]. However, these modifications are purely phenomenological.

The non-linear homogenization based approach models can describe the change of void shape and orientation, which have been showed important at low stress triaxialities [212].

7.4.1. Gurson-Tvergaard-Needleman (GTN) model

The GTN model postulates that damage growth is an outcome of material degradation occurring as the loading process takes place. Void evolution plays a crucial role in shaping the yielding behaviour, and is appropriately incorporated into the GTN model through a modification of the Von Mises yield potential [235]:

$$\varphi ={\left(\frac{{\sigma }_v}{{\sigma }_Y}\right)}^2+2q_1{f}^{*}\cosh \left(\frac{3}{2}{q}_2\frac{{\sigma }_H}{{\sigma }_Y}\right)-\left(1+q_3{f}^{{*}^2}\right)=0$$

(14)

where ${\sigma }_v$ represents the Von Mises equivalent stress, ${\sigma }_Y$ denotes the matrix material yield stress, and ${\sigma }_H$ refers to the hydrostatic component. $q_1$, $q_2$ and $q_3$ are material constants. $f^{\mathrm{*}}$ represents the effective volume void fraction, capturing the effects of void coalescence.

The evolution of voids is mathematically defined as follows:

$${\dot{f}=\dot{f}}_{growth}+{\dot{f}}_{nucleation}$$

(15)

where ${\dot{f}}_{growth}$ and ${\dot{f}}_{nucleation}$ are determined through equations (16) and (17).

$${\dot{f}}_{growth}=(1-f){\dot{\epsilon }}_{kk}^{pl}$$

(16)

$${\dot{f}}_{nucleation}=\frac{f_N}{S_N\sqrt{2\pi }}{\dot{\overline{\epsilon }}}^{pl}{e}^{{(-\frac{1}{2}\frac{{\overline{\epsilon }}^{pl}{\epsilon }_N}{S_N})}^2}$$

(17)

The volume dilatation rate is denoted by ${\dot{\epsilon }}_{kk}^{pl}$, whereas ${\overline{\epsilon }}^{pl}$ represents the equivalent plastic strain, and ${\dot{\overline{\epsilon }}}^{pl}$ is the rate of equivalent plastic strain. The standard deviation is represented by $S_N$.

The effective volume void fraction, $f^{\mathrm{*}}$ accounts for the micro-mechanics at work, including the presence of micro-voids the growth of voids, and the coalescence of neighboring voids leading up to the fracture at the macroscopic level [111]. $f^{\mathrm{*}}$ is determined as per the equation (18):

$$f^{\mathrm{*}}\left(f\right)=\left\{\begin{array}{c}f\\ f_C+K(f-f_C)\end{array}\right\}\begin{array}{c}f\le f_C\\ f>f_C\end{array}$$

(18)

where $K=\frac{f_{u-f_C}^{\mathrm{*}}}{f_f-f_C}$.

In this context, the void volume fraction, critical void volume fraction and void volume fraction at macroscopic failure are denoted by $f$, $f_c$ and $f_f$, respectively. The modified void volume fraction at macroscopic failure is represented by $f_u^{*}=f^{*}\left(f_f\right)$. The calibration process for the GTN model requires determination of a total of eight parameters, where a few parameters can be found in existing literature ($q_1$, $q_2$ and $S_N$) [218]. Moreover, $f_f$ is frequently assumed as 10% [236]. In the study conducted by Cui et al.[237], $f_f$ was computed through observations of the fracture morphology using SEM technology. Employing response surface methodology, five other parameters ($f_0$, $f_N$, $f_c$, ${\epsilon }_N$, $S_N$) were determined. In addition, Achineethongkham and Uthaisangsuk [238] identified damage model parameters through the use of metallographic analyses, local strain measurement and representative volume element (RVE) simulations.

The analysis of local crack initiation in Nakazima forming process on DP600 and TRIP600 steel was carried out by Uthaisangsuk et al. [235] through the application of the GTN damage model. Their findings indicated that this microstructure-based damage model provided satisfactory predictions of the fracture behavior of both multiphase steels in Nakazima stretching tests. Nevertheless, the influence of stress triaxiality on the fracture behavior has to be taken into account due to the void formation mechanism, which might be described by modifying model parameters [235]. Doig and Roll [218] stated that the GTN model is appropriate to treat real industrial cases. In their study, the GTN model accurately predicted the localization of cracks in a cross die part. However, similar to the FLD criterion, the coupled model overestimated the drawing depth. As opposed to the work conducted by Cui et al. [237], whom observed that the fracture position and limit strains in various loading conditions of a bulge test were almost identical between experimental and numerical simulation results of a DP590 steel. Furthermore, the GTN model was employed for evaluating the forming behavior in two-stage forming processes, as concluded by Achineethongkham et al. [238]. The authors concluded that the GTN model provided reliable predictions for the fracture states of samples subjected to NLSPs.

8. Conclusions

The shift towards the employment of AHSS in BiW components has been leading to forming challenges almost never encountered when forming conventional mild steels. Shear cracks, edge cracks and surface cracks are more frequently noticed within sheet metal forming industrial operations. The conventional FLD is unable to accurately model and account for these fracture occurrences.

The state of the art presented in this paper addresses the current trends to improve the forming behavior and fracture modelling accuracy of sheet metal forming processes. The FLD limitations and main drawbacks dealing with fracture occurrences are analyzed as well as the academic directions followed to overcome those limitations. The GFLC is an upgrade to the conventional FLD, currently implemented in AutoForm R10, aiming to improve non-linear strain paths modeling accuracy. However, further enhancements are required to deal with fracture modelling. Material modelling is significant for the forming behavior accuracy. Thus, a comprehensive review on material modelling is given. The directions followed by academic community with respect to material characterization methods and techniques is presented. Further improvements particularly for pure shear stress state are required in order to develop a standardized specimen. Furthermore, damage and fracture models are subjected to analysis.

The formability evaluation should be attained by the analysis of failure on account of through-thickness necks and splits as well as direct fracture. The consideration of a central framework is necessary to account for all relevant phenomenological outcomes. In this context, the advanced uncoupled phenomenological models are a promising approach. In order to garner industrial acceptance, the simplicity of use and implementation are favored. GISSMO model as well as ‘modular’ models have also been shown to provide accurate results. Even though, the analysis of relevant case studies, covering a wide range of stress states and complex strain paths is required to further enhance fracture modelling accuracy of such models when applied on AHSS in real industrial applications.