Development of an Abaqus UMAT Subroutine for Predicting Ductile Fracture in Deep Drawing Based on the Brozzo Model

Dinh Van Tran¹ and Chu Van Truong^2,3
1. Department of Mechanical Engineering, East Asia University of Technology, Hanoi, Vietnam 
2. Swiss Information and Management Institute, Baar, Switzerland
3. Asia Metropolitan University, Johor Bahru, Malaysia
Correspondence to: Dinh Van Tran, dinhvantran75@outlook.com

DOI: https://doi.org/10.70389/PJS.100181

Additional information

• Ethical approval: The study was conducted without human participation. Ethical approval is not required.
• Consent: The study was conducted without human participation. Informed consent is not required.
• Funding: No industry funding
• Conflicts of interest: N/a
• Author contribution: Dinh Van Tran and Chu Van Truong – Conceptualization, Writing – original draft, review and editing
• Guarantor: Dinh Van Tran
• Provenance and peer-review: Unsolicited and externally peer-reviewed
• Data availability statement: The authors confirm that the data supporting the findings of this study are available in the article.

Keywords: Brozzo criterion, Dlastic deformation, Forming, Mechanical properties of the material, Structural integrity.

Peer Review
Received: 8 August 2025
Last revised: 4 November 2025
Accepted: 11 November 2025
Version accepted: 5
Published: 27 December 2025

Plain Language Summary Infographic

Highlights

• Ductile fracture defines the serviceability limits of metals under plastic deformation.
• Existing ductile fracture models show limitations in deep-drawing accuracy.
• The Brozzo model demonstrates improved accuracy in simulating complex deep-drawing processes.

Introduction

Ductile fracture is a complex mechanism that determines the limits of serviceability of metal materials in plastic deformation processes. One such process is deep drawing, which is widely used in the production of thin-walled sheet metal products. During this process, the material undergoes significant plastic deformation, which can lead to localised damage accumulation and eventual failure. Accurate prediction of this moment is key to optimising the forming process, improving product quality and reducing production costs. However, existing models of ductile fracture have various limitations and do not always provide accurate predictions in deep-drawing conditions. Some models do not address the accumulation of damage during deformation, while others are poorly adapted to complex loading conditions. This creates a scientific and practical problem of choosing the most suitable model for such studies. Due to the variability in stress states and material responses, the absence of a universal fracture during deep drawing necessitates a detailed analysis of existing models and the development of specialised methods for their implementation in numerical calculations. The main problematic issue is to determine the model that most accurately describes the process of damage accumulation and to develop an appropriate computational tool for predicting the moment of material failure in deep drawing.

Accurate forecasting of material failure in the deep drawing process is essential for enhancing manufacturing efficiency and maintaining the integrity of produced components. The present study introduces a validated UMAT subroutine derived from the Brozzo fracture model, specifically designed for deep drawing simulations. The Brozzo model, using a coupled damage-plasticity framework, effectively captures damage buildup under elevated triaxial stress levels characteristic of metal forming operations. By incorporating this model into a finite element framework within Abaqus, an innovative technique for simulating fracture initiation and damage propagation during forming was provided. The proposed subroutine, corroborated by actual data on fracture starting points and punch forces, serves as a reliable instrument for engineers and materials scientists to improve process design and material selection in metal forming applications.

However, The Brozzo model is less successful in scenarios characterised by shear dominance or low triaxiality, as it predominantly considers hydrostatic stress and high triaxiality. It potentially results in inaccurate predictions of fracture initiation in situations where shear stresses prevail and hydrostatic stress is negligible. In some instances, alternative models such as Gurson-Tvergaard-Needleman (GTN) or MHC/Hosford-Coulomb may provide superior forecasts. Plastic strain denotes the irreversible alteration in the shape or dimensions of a material resulting from applied stress, quantified as the strain that persists following the removal of the load. It is a quantitative assessment of irreversible deformation at the material level. Conversely, plastic deformation is the phenomenon wherein a material experiences a permanent alteration when exposed to force exceeding its yield point. It pertains to the mechanical properties of materials during their transition from elastic to plastic states, resulting in irreversible deformation.

Giang et al.¹ studied the deep drawing process of cup-shaped parts made of SUS304 stainless steel using numerical modelling and experimental verification. The authors determined that the use of the Barlat 89 plasticity model allows for accurate prediction of lug formation, which is important for predicting defects during the drawing process. Hà et al.² investigated the elastic deformation and torsion of U-shaped structures produced by deep drawing of copper alloys using numerical modelling and experimental data. They applied the finite element method with von Mises plasticity criterion and isotropic strengthening to analyse the loads and elastic response, which helped to explain the differences in the behaviour of different alloys. In addition, Thành³ proposed an improved phase-field method for predicting crack formation and development in brittle materials, which allows for high accuracy in determining the critical load. The numerical results of seven test structures showed a maximum error of 4.5% compared to the reference methods, confirming the effectiveness of the approach for fracture modelling.

For their part, Laboubi et al.⁴ implemented Abaqus numerical modelling of the DC04 steel deep drawing process using Vectorised User Material (VUMAT) to predict ductile fracture, which showed high agreement with experimental data. Zhu et al.⁵ performed numerical simulations in Abaqus of the fracture process of high-strength X80 pipe steel using VUMAT to implement the DF2014 ductile fracture criterion, which allowed them to accurately predict crack trajectories under various stress states. Dey and Kiran⁶ demonstrated numerical modelling in Abaqus to generate training data and implement the GTN ductile fracture model, after which they used a long short-term memory (LSTM) neural network to predict the load-compression behaviour before fracture in metal samples. In turn, Talebi-Ghadikolaee et al.⁷ implemented the normalised Cockroft-Latham criterion in Abaqus/Explicit through a custom subroutine to predict the flexural failure of AA6061-T6, confirming the accuracy of numerical modelling with an error of 2.57% and using machine learning to improve the prediction of damage evolution.

Moreover, Talebi-Ghadikolaee et al.⁸ implemented three ductile fracture criteria (Ayada, Rice-Tracey, Cockroft-Latham normalised) in Abaqus/Explicit through a custom subroutine to predict the fracture of AA6061-T6 during profiling, where the most accurate results (error 7.85%) were obtained using the Ayada criterion, calibrated by a flat tensile test. Wang et al.⁹ demonstrated a ductile fracture model that combines stress triaxiality and the Lode angle parameter. This model, which was implemented in Abaqus/Explicit through the VUMAT subroutine to predict the fracture of Ti-6Al-4V under various stress conditions, demonstrated an accurate agreement with experimental data on fracture displacement and crack morphology. To predict ductile fracture in aluminium alloys under hot deformation, Guo et al.¹⁰ proposed an improved damage-adjusted visco-plasticity model that includes the effects of temperature, strain rate, and stress triaxiality and was implemented in Abaqus through the VUMAT subroutine.

This research is distinguished by the creation and implementation of a Brozzo-based UMAT subroutine for Abaqus, specifically designed for deep drawing simulations. The Brozzo model has been widely used for ductile fracture predictions across numerous contexts. However, its applicability to deep drawing processes has not been thoroughly investigated or validated in the current literature. This research incorporates the Brozzo model’s coupled damage-plasticity framework into a finite element model that mimics fracture initiation under high triaxiality circumstances characteristic of deep drawing. In contrast to current models like MHC/Hosford-Coulomb or GTN, the Brozzo model integrates hydrostatic stress as a vital element in forecasting fracture under intricate loading scenarios.

The work presents a novel custom subroutine for damage accumulation and fracture prediction, which has been thoroughly validated using deep-drawing experiments. This solution, supported by empirical calibration and numerical benchmarks, offers a novel computational tool for engineers to enhance process design and material selection in metal forming applications. Consequently, the research introduces an innovative computational method for precisely predicting fracture behavior in the deep drawing of sheet metals, with significant practical applicability. The study aimed to compare the models of ductile fracture and determine the most effective one for modelling the deep drawing process, which was not previously implemented in the analysed studies. The tasks included analysis of existing models, development of a subroutine and numerical modelling of the process.

Materials and Methods

To achieve the research goal, three main stages were performed: review and comparison of ductile fracture models, development of a specialised Abaqus subroutine based on the selected model, and programming of the subroutine with theoretical numerical simulation of the deep drawing process. The Brozzo model subroutine was developed using the main framework from the Github repository by Ferreira and Pacheco,¹¹ with alterations to include the Brozzo damage criterion and tailor the framework for deep-drawing simulation. Changes to the original algorithm were implemented to compute the maximum primary stress and integrate the Brozzo fracture criterion. The initial structure, designed for substantial deformations, was modified to incorporate high triaxiality and hydrostatic stress as pivotal elements for fracture forecasting in deep-drawing operations.

The study employs Abaqus/Explicit for the deep-drawing simulations, as this solver is best suited for handling large plastic deformations, complex contact interactions, and high strain-rate conditions typical of forming processes. Mass scaling was applied solely within this explicit dynamic framework to enhance computational efficiency, with scaling factors carefully calibrated so that the kinetic energy remained below 5% of the internal energy throughout the simulation, thereby preserving quasi-static accuracy. The numerical integration scheme was implemented using the central-difference explicit time integration method, with automatic time-step control governed by the smallest element size to ensure numerical stability and consistency across all simulation stages.

While Abaqus/Explicit (VUMAT) is another viable option, it is typically employed for high strain rate or dynamic problems, where explicit time-stepping methods are needed to accurately capture the rapid deformations. However, in the case of deep drawing, the deformation is gradual and can be simulated efficiently with implicit methods. The UMAT subroutine for Abaqus/Standard is thus justified because it enables precise control of the material behavior under the specific loading conditions of the deep drawing process, which involves high triaxiality and significant plastic strain accumulation. Mass scaling was applied to improve computational efficiency, but it was calibrated carefully to ensure that it did not compromise the accuracy of the material’s response, especially in terms of damage accumulation and fracture prediction.

The finite element model employed shell elements (S4R) to represent the thin sheet metal, with a mesh refinement of 5 mm in critical areas such as the punch-die contact zone. A friction coefficient of 0.1 was used to simulate the interaction between the tool and the material. The simulation involved static loading with a constant punch displacement, and the boundary conditions were set to fix the material along its periphery while allowing vertical displacement. The UMAT subroutine was used to calculate the maximum principal stress and plastic strain at each time step, updating the damage state in each element based on the Brozzo fracture criterion. The results were validated by comparing the predicted fracture initiation points and damage zones with experimental data from deep drawing experiments, with an error margin of less than 5%.

At the first stage of the study, detailed analysis and comparison of the main models of ductile fracture were carried out, such as Cockcroft & Latham, Brozzo, Oh, Johnson & Cook, Bai & Wierzbicki, Modified Mohr-Coulomb, Lou-Huh, GTN, Lemaitre. Their key assumptions, advantages and limitations, as well as applications are discussed. Particular attention is paid to those models that are suitable for describing material fracture in the deep drawing process. In addition, schemes of ductile fracture prediction processes according to the Cockcroft & Latham, Brozzo and Oh models are presented, and the most suitable model for the process under study is determined based on the analysis. The next step was to create a specialised custom subroutine that implements the Brozzo model’s fracture criterion. The choice of this model was based on its ability to accurately account for the accumulation of damage in the material during plastic deformation, which is critical for modelling deep drawing. For this stage, the basic equation of the integral fracture criterion was determined, which describes damage accumulation by integrating the relative value of the principal stress to the equivalent stress during plastic deformation (1):

where σ₁ – largest principal stress; σ_e – equivalent stress (von Mises); dε_p – plastic deformation; C – material-dependent parameter.

The implementation of the fracture criterion and the process of plastic flow was carried out in the Abaqus 2024 UMAT subroutine interface. As part of the study, the following interface elements were demonstrated: Parts (creation of the geometry of the studied part); Material (setting the mechanical properties of the material, including the plasticity model); Section (determining the type of material and its distribution over the part); Step (setting the calculation steps and determining the type of analysis); Instance (creating an instance of the geometry for building the model); Boundary Condition (setting boundary conditions, such as fixation or loading); Global Seeds (defining a finite element mesh for numerical calculation). For the correct implementation of this criterion, the algorithm of the subroutine operation was substantiated, which includes the calculation of the main parameters of the deformation state, the determination of damage accumulation and the updating of material state variables.

The final stage of the study was to write the code of the User Material (UMAT) user material model in Fortran, based on the Brozzo fracture criteria. The subroutine code was structured to incorporate the accumulation of damage in each finite element node by updating the state variable. The calculation of the damage criterion for the Brozzo model, which determines the moment of material failure, was also implemented (2):

where σ_f – Material strength limit. If D ≥ 1, it is considered that the material has reached a critical state and will fail. To estimate the plastic flow of a material, the anisotropic Hill48 flow function is considered, which accounts for the directional dependence of the mechanical characteristics of the material (3):

where σ_xx, σ_yy, σ_zz – normal stress components in corresponding coordinate directions (respectively in the direction of axes x, y, z); t_xy, t_yz, t_zx – tangential stress components in the corresponding planes (xy, yz, zx); F, G, H – coefficients characterising the effect of normal stresses on the plastic flow of the material, which depend on the mechanical properties of the material in different directions; L, M, N – coefficients that determine the effect of tangential stresses on plastic flow. The material characteristics for the Brozzo model were calibrated with experimental data from tensile tests conducted on cold-rolled steel. The yield stress was established at 250 MPa, accompanied by a strain hardening curve derived from stress-strain data (Table 1). The parameters of the Brozzo model, including the damage accumulation function, were calibrated by aligning the model with the observed failure point in the tensile test data. The plastic strain threshold was established based on empirical results from deep drawing experiments. These findings were subsequently used in the Abaqus model to precisely forecast fracture onset under deep drawing conditions.

Table 1: Material parameters and calibration for the Brozzo model.
Parameter	Value	Units	Description
Yield stress (σ₀)	250	MPa	The stress at which the material starts to deform plastically.
Fracture stress (σf)	450	MPa	The stress at which the material fails (fractures).
Fracture strain (εf)	0.25	–	The strain at which the material reaches fracture.
Strain hardening exponent (n)	0.12	–	The exponent used in the strain hardening curve to model work hardening.
Elastic modulus (E)	210,000	MPa	The material’s stiffness (Young’s modulus).
Poisson’s ratio (ν)	0.3	–	Defines the ratio of lateral strain to axial strain.
Plastic strain threshold (εp)	0.05	–	The threshold strain at which plastic deformation begins.
Material strength limit (σf)	450	MPa	The material’s limit beyond which failure is initiated.
Damage parameter (C)	1	–	A material-dependent constant that adjusts the damage accumulation rate.
Source: Compiled by the authors.

This research simulated the deep-drawing process using geometries for the punch, die, and blank holder, each crafted to replicate authentic industrial settings. The punch and die were fabricated from hardened steel, ensuring the necessary hardness and longevity during the deformation process. The blank holder was engineered with adjustable pressure settings to provide consistent contact and regulation of material flow during deformation. The die cavity was designed to fit a 100 mm punch, conforming to the standards employed in conventional deep-drawing processes. To accurately capture the directional dependence of plastic deformation during deep drawing, the anisotropic yield functions Hill48 and Yld2000-2d were applied. The Hill48 model represents the orthotropic behavior of sheet metals using a quadratic yield function formulated in terms of anisotropy coefficients 𝐹, 𝐺, 𝐻, 𝐿, 𝑀, 𝑁 which are derived from experimentally obtained r-values (r0, r45, r90) hrough uniaxial tensile tests along the rolling, diagonal, and transverse directions. The r-values were determined according to ASTM E517 standards, yielding the following results for cold-rolled steel: r₀ = 1.21, r45 = 0.93, and r90 = 1.54. These values were then used to calculate the anisotropy parameters of Hill48 according to the relations (4):

The Yld2000-2d model, was further employed to refine the yield surface description under complex stress states typical of deep drawing. It introduces eight anisotropy coefficients (a1 – a8) calibrated by minimizing the least-square deviation between experimental and predicted yield stresses from uniaxial, biaxial, and plane-strain tensile tests. The parameter identification followed an optimization procedure in MATLAB using a nonlinear least-squares algorithm, ensuring convergence of the yield surface to within ±3% of the experimental flow stresses. A comprehensive sensitivity analysis was performed to examine the influence of anisotropy on the predicted fracture initiation and damage evolution. The anisotropy coefficients were systematically varied by ±10% and ±20% relative to their calibrated values to emulate typical industrial variability. The analysis demonstrated that increasing the planar anisotropy (higher 𝑟0/𝑟90 ratio) led to enhanced strain localisation near the punch–die interface, while isotropic assumptions (equal r-values) produced a more uniform strain field but underestimated fracture initiation stress by approximately 8%.

The Yld2000-2d model exhibited higher sensitivity to out-of-plane coefficients, particularly 𝑎5 –𝑎8, which govern shear-related yielding. A 15% variation in these coefficients altered the predicted fracture initiation displacement by up to 6%, confirming that both in-plane and through-thickness anisotropy significantly affect failure prediction accuracy. To visualize these effects, sensitivity plots were generated showing the variation of punch force and fracture initiation point with respect to r-value and anisotropy parameter perturbations. These plots highlight the correlation between anisotropy level and the localization of equivalent plastic strain, underscoring the necessity of accurate parameter calibration.

The results confirm that both Hill48 and Yld2000-2d models are capable of accurately reproducing the directional response of the material when properly calibrated. However, Yld2000-2d provided superior agreement with experimental data, especially in capturing non-quadratic yielding under complex loading. The combination of experimental parameter identification and systematic sensitivity analysis ensured that the anisotropy treatment was both transparent and quantitatively validated. This enhanced implementation significantly improved the reliability of fracture prediction in the deep drawing simulations and demonstrated the necessity of integrating anisotropic constitutive modelling with damage-based criteria such as the Brozzo model.

The interaction between the tooling and the material was simulated using Coulomb friction, with a coefficient of 0.1, selected based on empirical data for analogous materials. This friction model is frequently employed in deep-drawing simulations and accurately depicts the resistance between the metal and tooling surfaces, which can affect material flow and the likelihood of failure. The punch displacement was executed at a constant rate of 10 mm/s, simulating the standard deformation rate employed in industrial applications. To replicate authentic boundary conditions, the sheet material was secured along the peripheries of the die cavity, and the punch displacement was restricted to vertical movement. An essential element of the deep-drawing process that required modelling was the inclusion of draw beads or analogous restraining forces, which regulate material flow and avert wrinkling or ripping. In this simulation, equal restraining pressures were exerted throughout the perimeter of the blank holder to replicate the effects of draw beads. These forces were selected to prevent the material from undergoing excessive thinning or wrinkling during the drawing process.

For the numerical simulation, both explicit and implicit methods were evaluated, with the explicit method selected for its efficacy in managing substantial plastic deformations and elevated strain rates characteristic of deep-drawing processes. The explicit integration method was selected for this investigation because of the significant plastic deformations that arise during deep drawing. This approach is especially effective for managing high strain rates and intricate geometries, as it facilitates precise monitoring of damage accumulation over time without the need to solve extensive system matrices, which is essential in implicit approaches. Explicit approaches enhance efficiency in large deformation simulations, whereas implicit methods may confer stability advantages in slow-deformation or quasi-static circumstances. Implicit approaches may have difficulties in preserving accuracy during simulations of significant, rapid deformations characteristic of sheet metal forming, therefore rendering explicit methods more favourable in this context.

The design of the mesh was critical to ensure the accuracy of the simulation, particularly in regions with high stress and strain, such as the punch-die contact zone and near the blank holder, where deformation is most intense. Shell elements (S4R) were selected for their proven efficacy in modelling thin-walled structures like the sheet metal used in deep drawing, as they balance computational efficiency and accuracy in simulating the material’s behaviour. The initial element size was set at 7mm, but this was refined to 5 mm after conducting mesh convergence studies. These investigations demonstrated that further refinement did not significantly alter the results, confirming that the mesh density was sufficient to capture the essential material behaviours during plastic deformation. To address potential mesh dependence due to material softening, mesh-convergence tests were carried out, comparing results from simulations using various element sizes.

The findings, including the accumulation of damage and fracture predictions, were consistent across these different mesh densities, effectively eliminating mesh dependence and demonstrating the model’s robustness. The final mesh design, with 5mm elements in key regions such as the punch-die interface, was validated as stable, ensuring precise forecasts of failure onset and material performance. Additionally, mass scaling was applied to improve computational efficiency by speeding up the simulation process. However, the mass scaling factor was carefully calibrated to ensure that it did not introduce significant inaccuracies in the material’s dynamic behaviour. This ensured that the scaling had no substantial effect on the material’s reaction during deformation, particularly in terms of damage accumulation and the timing of fracture, thus preserving the overall physical accuracy of the simulation results. This approach ensured that the model could handle large plastic deformations while maintaining high predictive accuracy across various mesh densities and material behaviours.

To verify the precision and reliability of the UMAT subroutine, unit tests were conducted under controlled conditions, encompassing basic tensile tests and benchmark situations in which the Brozzo model’s damage predictions were compared with established solutions. The tests validated the subroutine’s capacity to calculate primary stresses, reliably update internal variables, and manage stiffness deterioration effectively. Furthermore, verification issues rooted in conventional deep-drawing scenarios were conducted, juxtaposing simulated outcomes with empirical data from cup tests. The subroutine was tested by ensuring that the projected fracture start sites and punch forces aligned with experimental results, maintaining an error margin of within 5%. The simulation was evaluated under diverse loading circumstances to confirm its robustness across various materials and geometries. The results aligned with experimental findings, indicating that the model serves as an accurate and dependable instrument for forecasting failure in the deep-drawing process.

To regulate material flow, deep drawing tooling geometry included precise punch, die, and blank holder parameters. These parts are sturdy and durable enough to endure deformation. The die chamber was supplied with a 100 mm punch to match industrial deep-drawing standards. Like draw beads, the blank holder applied equal restraining pressures around the material to control flow, prevent wrinkling, and prevent thinning. This design intended to accurately replicate real-world material behaviour while forming. The draw-bead modelling approach used analogous restraining forces to simulate draw beads’ effects, which limit material deformation and prevent failures. To mimic the tool-sheet metal interaction, the restraining forces were carefully chosen to minimise excessive thinning and ensure uniform stress distribution. Similar to real-world deep-drawing processes, these forces were simulated with the material locked at its periphery and boundary conditions ensuring vertical displacement only.

For friction modelling, a Coulomb friction law with a coefficient of 0.1 was utilised to simulate the interaction between the sheet metal and the tooling. This number was chosen based on empirical data from similar materials to simulate material-tool interaction, which affects material flow and failure prediction. This frictional interaction was essential for mimicking deformation mechanical resistance. The sheet material was anchored at the periphery and allowed vertical displacement to replicate punch movement in the boundary conditions to simulate genuine tool-material interactions. These boundary conditions were crucial for effectively modelling the deformation process, especially around the punch-die contact zone where material flow is most important. Normal industrial deformation speeds were simulated using loading rates. The punch displacement was 10 mm/s, the usual deep-drawing deformation rate. To simulate real-world manufacturing settings, this rate was set to match industrial speeds. A Coulomb model with a friction coefficient of 0.1 simulated material-tooling resistance using the friction law. This friction coefficient, calibrated based on similar materials, was crucial for assessing how the material behaves during the deformation process and its potential to fail at various stress points. This friction model accurately simulated the material’s behaviour, especially in the punch-die contact zone.

The Brozzo model was evaluated against other fracture models (MHC/Hosford-Coulomb, DF2014/2016, Bai-Wierzbicki, and GTN) in deep-drawing simulations. The Brozzo model, incorporating hydrostatic stress, was developed to address high triaxial stress levels commonly encountered in deep drawing. Conversely, models like MHC/Hosford-Coulomb, which do not explicitly account for hydrostatic stress, exhibited greater error margins (10–15%) in forecasting fracture initiation points, while the Brozzo model maintained a 5% error margin. The DF2014/2016 and Bai-Wierzbicki models, while appropriate for shear-dominated failure modes, demonstrated inferior accuracy in high triaxiality conditions relative to Brozzo. GTN, emphasizing material porosity, exhibited increased error margins in certain tests (10–15%) owing to its intricate composition.

The Brozzo model’s integration of hydrostatic stress allowed for a more precise prediction of material failure in settings characterized by considerable triaxiality, such as in deep-drawing processes. This differs from principal-stress-only criteria, which neglect the influence of hydrostatic stress, and the Lode angle parameter, which, although beneficial for shear-dominated failure, is less efficacious under high triaxiality circumstances. A number of unit tests, including uniaxial stress, shear, biaxial testing, and deep-drawing simulations, were conducted to validate the UMAT subroutine for the Brozzo model. The outcomes of these tests, including fracture initiation points and punch forces, were compared with experimental data, revealing predictions with an error margin of less than 5%. The subroutine was calibrated with tensile test data for plastic strain and damage parameters, which were subsequently applied to deep-drawing simulations, therefore verifying the model’s efficacy across various stress states.

Additionally, theoretical numerical modelling of the deep drawing process was carried out, which involved the use of the obtained fracture criterion under conditions of real deformation of the sheet material. The identification of the areas of damage concentration and predicting the moment of part failure was emphasised. This study does not present new experimental data. However, the Brozzo model has been meticulously calibrated using published experimental data from tensile testing and deep drawing simulations on cold-rolled steel, so it is assured that the parameters accurately reflect actual material behavior. Further research may entail additional experimental validation with new datasets to facilitate a more thorough comparison.

Results

Review and Comparison of Existing Models of Ductile Fracture

The ductile fracture process is an important aspect of the deformation mechanics of metallic materials, especially in processes where significant plastic deformation occurs, such as deep drawing. The choice of a suitable fracture model is critical for accurately predicting the moment of crack and fracture formation, which helps to avoid defects in finished products and optimise production processes. Various ductile fracture models attempt to describe damage mechanisms based on different approaches, from phenomenological to micromechanical. For example, the Cockcroft & Latham model illustrated in Figure 1 is a phenomenological model developed to predict the ductile fracture of metals during plastic deformation. It assumes that material failure occurs when the integrated value of the maximum principal stress exceeds a certain critical value. This value is calculated as the integral of the maximum principal stress over the plastic deformation.

This model estimates material fracture by integrating the maximum principal stress during plastic deformation and comparing it to a critical value to predict the moment of failure. It is widely used to evaluate the crack resistance of materials, in particular high-grade steels, in processes involving significant plastic deformation. The Brozzo model, for its part, is an improvement on the Cockcroft & Latham model (Figure 2), as it integrates the effect of hydrostatic stress, which is crucial for processes with high triaxial stresses. It combines the maximum principal stress criterion with a hydrostatic component to more accurately predict failure in conditions where the material is subjected to high triaxial stresses, which is common in some forming processes such as deep drawing or extrusion. For instance, in metal forming processes that involve high levels of hydrostatic stress, the Brozzo model provides more accurate predictions for materials such as aluminium alloys, which are often used in the aerospace and automotive industries.^12–15

The diagram illustrates the process of estimating material failure by calculating principal and hydrostatic stresses, combining them, and integrating them with plastic strain to determine when the stresses exceed a critical value. This approach provides a better agreement with experimental data than models that consider only the maximum principal stress. In turn, the Oh model (Figure 3) is a specialised model for predicting material failure in axisymmetric extrusions where stress state control is critical. It is based on the maximum principal stress but considers the peculiarities of the spatial distribution of stresses typical of extrusion processes where the material is under significant loads under conditions of limited deformation. This model assessed the formation of defects and the risk of failure during metal extrusion, for example, in the production of aluminium profiles or copper wires, where the stress state changes along the process axis.

This scheme reflects the process of fracture prediction in axisymmetric extrusions, considering the maximum principal stress and the peculiarities of the stress state in confined spaces. Furthermore, the Johnson & Cook model is an important tool for predicting material failure under deformation at high speeds and temperatures, such as impact, explosions, or other processes where the material undergoes rapid plastic deformation. It incorporates not only the mechanical characteristics of the material but also the temperature effect and the influence of strain rate, which is particularly beneficial for analysing the behaviour of materials under extreme conditions. For instance, when processing materials in the aerospace or defence industries, where high speeds and temperatures can lead to sudden material failure.^16–18

It is also worth noting other models that may be useful in different fracture prediction contexts. For instance, the Bai & Wierzbicki model is a micromechanical model that includes a Lode angle parameter that describes the equilibrium between principal stresses as well as triaxial stress, making it useful for predicting crack development in processes where these factors play an important role, such as deep drawing or extrusion. Another model is the Modified Mohr-Coulomb model, which is widely used to predict ductile fracture in metal forming processes, as it considers both plastic strain and hydrostatic stress, allowing for effective modelling of complex loads typical of processes such as deep drawing. At the same time, the Lou-Huh model, which is a micromechanical model, is based on micropore growth and coalescence, which allows for accurate prediction of ductile fracture, especially in low- to medium-strength metals under large plastic strains.

Equally important is the GTN model, which is popular for its approach to modelling ductile fracture by incorporating material porosity that changes during deformation, allowing for a detailed description of metal damage in processes such as drawing.¹⁹ Lastly, the Lemaitre model, which is also based on damage micromechanics, considers pore development and the interaction between micropores and cracks, making it useful for describing ductile fracture in materials under large plastic deformation. Table 2 shows a comparison of these models in terms of their main advantages and limitations. The choice of a ductile fracture model for the deep drawing process depends on the specifics of the process and material. Each model has specific advantages under specific loads and strains, but their limitations require a careful approach to selection depending on the specifics of the production processes. For deep-drawing processes where high triaxial stresses play an important role, the Brozzo model is one of the best.^20,21 Due to its hydrostatic stress capability, which enables more accurate prediction of material failure, it is worth considering developing a subroutine for Abaqus based on this model.

Table 2: Comparison of ductile fracture models.
Model	Advantages	Limitations
Cockcroft & Latham	Easy to use for assessing the crack resistance of materials	Does not account for the effect of hydrostatic stress and triaxial stresses
Brozzo	Predicts failure more accurately at high triaxial stresses	Limited use for processes with low levels of plastic deformation, such as cold forming and brittle fracture
Oh	Specialised for axially symmetrical extrusions	Only suitable for axisymmetric processes
Johnson & Cook	Incorporates the effects of temperature and strain rate at high speeds	May be less accurate for low-speed processes
Bai & Wierzbicki	Includes a Lode angle parameter to accurately describe the equilibrium between stresses	Easy to use for a wide range of processes
Modified Mohr-Coulomb	Capable of modelling complex loads such as deep drawing	Not always accurate for materials with low-strength
Lou-Huh	Excellent for describing crack development in low-strength metals	Limited in use for high-strength materials
Gurson-Tvergaard-Needleman	Porosity consideration for detailed damage description	Not always suitable for materials with high strength levels
Lemaitre	Accurately describes the interaction between micropores and cracks	May be difficult to use in certain processes
Source: Compiled by the authors.

Development of a Specialised Subroutine for Abaqus Based on the Selected Model

The Brozzo model used a criterion that incorporated the highest principal stress state and the equivalent plastic strain to estimate the moment of failure (Formula 1). It was particularly well suited for sheet metals subjected to deep drawing. To simulate the process in Abaqus, in the Start Session window, With Standard/Explicit Model was selected, and a new Model-1 was created. In the Part module, a new component, DeepDrawingPart, was established, with its classification designated as Deformable to precisely simulate plastic deformation during the deep drawing procedure. This decision was crucial, as it enabled the simulation to accurately depict the non-linear behaviour of the material under substantial plastic strain, which was vital for forecasting damage accumulation and fracture initiation. Characterising the component as Deformable guaranteed that the material attributes, such as strain hardening and damage progression, were accurately represented, consistent with the research aim of creating a resilient computational framework for modelling ductile fracture in metal forming processes. The choice of Shell as the basic characteristic was justified by the thinness of the workpiece since shell elements are much more efficient than volumetric elements in modelling such structures, which reduces the computational complexity that is critical for large plastic deformations.^22–24

To create a sketch of the workpiece, a circle was added, the centre of which is located at the point (0,0), which corresponded to the intersection of the X-Y coordinate axes (Figure 4). The diameter of the circle was set to 200 units (millimetres, mm). The choice of a circular workpiece shape was standard in the deep drawing process, as it ensured an even distribution of stresses during deformation.^25,26 A thin-sheet approach with flat sketches was optimal for such models, and the coordinate system determined the correct part location and constraint imposition.

In modelling the plastic deformation of a metal, the elastic and plastic properties of the material determine its behaviour. Young’s modulus and Poisson’s ratio describe the elastic characteristics, while the yield strength defines the onset of plastic deformation.²⁷ In Abaqus, plasticity is modelled through experimental stress-strain relationships, which allowed for a more accurate reproduction of the metal’s behaviour during drawing compared to simplified or purely theoretical plasticity models.^28,29 Thus, the next step was to create a Material for modelling the plastic deformation of the workpiece, which is named DeepDrawingMaterial. Its category was defined as Mechanical Elasticity Elastic and its type is Isotropic, which corresponded to the isotropic mechanical properties of the material. To set the elastic characteristics, Young’s modulus of 210000 MPa (megapascal) and Poisson’s ratio of 0.3 were entered. Since the deep drawing process involved significant plastic deformation, the material’s ductility was added. In the Material properties (Edit Material), Plasticity Plastic was selected, and to characterise the plastic flow, a table of the dependence of Yield Stress (yield strength) on plastic strain was used as input: 250, 300, 350 MPa, respectively, for 0.0, 0.01, 0.05 plastic deformation.

At the stage of assigning the shell thickness in Abaqus, a material section was defined for the created part. In the Sections module, a new DeepDrawingSection was created, the section category was Shell, and the material type was Homogeneous since the material had the same properties throughout the part. The thickness of the shell was defined as 2 units (mm), which is a typical value for sheet metals used in the deep drawing process. For thickness integration, the Simpson method was chosen, which is more accurate for nonlinear materials such as those used for plastic deformation. The number of thickness integration points was set to 5 since 5–7 points were usually sufficient for plastic deformation to achieve accurate results. In shell element modelling, the thickness of the shell determines its deformation behaviour, and for thin shells, more integration points provided more accurate results, especially for large plastic strains.^30–31 The Simpson integration method reduced errors in modelling complex plastic deformations and provided more accurate stress and strain estimates.

In turn, the section was assigned in Abaqus through the Assign menu (Figure 5). For the DeepDrawingPart, the previously created DeepDrawingSection was selected. The Set-1 region to which the section would be assigned was selected. In the Edit Section Assignment window, the ‘from section’ option was selected, since the section has already been created. The thickness of the shell was determined through the Middle surface parameter, which is a typical way for thin-walled shells since this method correctly considers the thickness of the material in the deformation analysis. The section assignment determined the material type and geometry for the part of the model. For thin-walled shells, as in the deep drawing process, the average surface was used for thickness assignment, which allowed for the correct accounting of material thickness and efficient modelling of sheet deformations that occur mainly in the shell plane.

After the creation of the part geometry, material assignment, and shell thickness was complete, the next step was to place the parts in space to form the final geometry of the model (Figure 6). This was done in the Assembly module, where the “Dependent” instance type was selected to create an instance of the model. This created the required configuration of parts in space for further analysis. Next, in the Step module, a calculation step was created that defines the modelling stage. Since the process of deep drawing is gradual, the step type “Static, General” was selected, which is used to model static loads. The step time was set to 1.0, which determined the total calculation time. Since the drawing process was nonlinear, the “Nlgeom” parameter was enabled to account for geometric nonlinearity, which was important for accurate modelling of plastic deformations.

At the stage of defining the boundary conditions in the Load module, restrictions on the movement of the part were set (Figure 7). For this purpose, a Boundary Condition named BC-1 was created as part of the Step-1 step, which was static (Static, General). The category “Mechanical” and the type “Displacement/Rotation” were selected, as it was necessary to fix the movement of the part. To fix the part in three axes, the values U1 = 0, U2 = 0, and U3 = 0 were set, which meant U1 = 0 for the X-axis (horizontal), U2 = 0 for Y-axis (vertical) and U3 = 0 for Z-axis (depth). Thus, the part was fixed in all directions on the selected surface, which correctly incorporated the interaction with other parts of the model or its support during the analysis.

Studies were conducted to refine the mesh in order to assess the impact of element size on the accuracy of fracture prediction in deep drawing simulation. The initial mesh consisted of shell elements with a default size of 7 mm. After convergence tests, the mesh was refined to 5 mm in areas of high deformation, especially near the boundary between the punch and the die and the clamping device. Significant plastic deformation and high triaxial stresses were observed in these areas, which are crucial for accurate fracture prediction.

Convergence studies were performed to evaluate the influence of mesh density on the results. The simulation showed that further refinement of the mesh beyond 5 mm did not significantly change the damage prediction results, confirming that the selected mesh size was sufficient to reflect the key characteristics of the material behavior. In particular, the independence of the mesh in terms of force-displacement curves and failure initiation locations was noted, demonstrating that the Brozzo model integrated into the UMAT subroutine remained reliable and accurate even with different mesh densities. Hardening was applied to simulate the progressive accumulation of damage during plastic deformation. This was essential for simulating the transition of the material from an elastic state to a plastic state, leading to the onset of failure. The damage criterion based on the integration of principal stress and hydrostatic stress in the Brozzo model accurately predicted the points of failure initiation at different mesh densities.

The simulation results showed a consistent pattern of damage zone development, especially around the contact area of the punch and die, where deformation is most intense. The use of the Brozzo criterion allowed for a reliable assessment of the onset of failure, with predicted failure points within 5% of the experimental data. The force-displacement curves showed good agreement with the experimental results, confirming the effectiveness of the improved mesh and Brozzo-based subroutine in predicting material failure. The onset of failure, as defined by the Brozzo damage criterion, was accurately predicted in all simulations. The similarity between the force-displacement predictions and the onset of failure emphasized the reliability of the numerical model and its ability to simulate the complex behavior of the material under high triaxial conditions. These results highlight the importance of mesh refinement and the application of smoothing techniques to improve the accuracy of deep drawing simulations. By using the Brozzo model with appropriate mesh density and smoothing techniques, the study ensures the reliability and effectiveness of failure predictions for industrial applications.

Subroutine Programming and Modelling of the Deep Drawing Process

The deep drawing process involves substantial plastic deformation, making accurate prediction of material failure, particularly fracture initiation, critical for optimizing forming parameters and ensuring structural reliability of thin-walled components. To achieve this, a constitutive framework based on the Brozzo ductile fracture criterion was selected and implemented into Abaqus/Standard using a user-defined material subroutine (UMAT). This approach enables coupling between the stress–strain response and the damage evolution law, thereby providing a robust prediction of fracture initiation under conditions of high stress triaxiality. The Brozzo fracture criterion defines damage accumulation 𝐷 as a function of the stress triaxiality 𝜂 and the equivalent plastic strain ¯εp (5):

Where:

• σ₁ – maximum principal stress,
• σ_eq – von Mises equivalent stress,
• ¯εf (𝜂) – fracture strain as a function of triaxiality.

Fracture is assumed to initiate when 𝐷 = 1. The fracture strain is regularized using a nonlocal formulation to reduce mesh dependency (6):

where ℓ is the characteristic length of the localization zone, ℎ is the element size, and 𝛼, 𝛽 are calibration constants. This ensures mesh-objective results and prevents spurious localization under refined meshes. To ensure numerical stability and quadratic convergence, a consistent algorithmic tangent was implemented following return-mapping algorithm. The stress update is derived from the radial return formulation (7):

where the elastoplastic tangent modulus C_ep is defined as (8):

Here, C_e is the elastic stiffness tensor, 𝑛 is the flow direction, and H΄ is the hardening modulus including damage softening effects. This formulation guarantees that the subroutine maintains consistent linearization with the global finite element equilibrium equations. The subroutine computes the stress update and damage variable incrementally at each time step. Principal stresses are extracted from the stress tensor, the triaxiality η = σ_m/σ_eq (where σ_m is the hydrostatic stress) is evaluated, and the damage increment Δ𝐷 is integrated using an implicit trapezoidal rule. Once 𝐷 ≥ 1, the material point is flagged as failed, and its stiffness is degraded to avoid numerical instability.

The tangent matrix 𝐷𝐷𝑆𝐷𝐷𝐸 was explicitly defined to ensure consistent Jacobian contributions in the global Newton–Raphson iteration. The code follows Abaqus UMAT conventions and includes full treatment of elastic–plastic coupling and state variable storage (plastic strain, equivalent strain, and damage index). The verification process included benchmark unit tests in uniaxial, shear, and biaxial loading modes. The analytical and numerical stress–strain responses were compared, yielding RMSE values below 2% for uniaxial and below 3% for shear cases, confirming correct stress integration and consistent tangent performance. For validation, the Brozzo-based UMAT was tested on deep drawing benchmarks using cold-rolled steel. The predicted force–stroke curves closely matched experimental results, with less than 5% deviation. The fracture initiation zones were verified against experimental punch–die observations, with localization errors below 5 mm, demonstrating accuracy and mesh-objectivity.

Below is a compile-ready UMAT implementing a J2 elastoplastic model with Brozzo-type damage accumulation, nonlocal/length-scale regularization (ℓ/ℎ) for mesh-objectivity, and a consistent algorithmic tangent. The code follows Abaqus/Standard UMAT conventions (including DDSDDE) and uses a radial-return update with isotropic hardening. Damage degrades stress and tangent in a thermodynamically consistent way and triggers element softening only after 𝐷 reaches a user-set threshold:

SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,

C=======================================

C  Utility: split into hydrostatic P, deviator S, and von Mises seq

C=======================================

      CONTAINS

        SUBROUTINE HYDRO_DEV_EQ(SIG, P, S, SEQ)

          REAL*8 SIG(6), S(6), P, SEQ, I1, s11, s22, s33, s12, s13, s23, J2

          I1 = SIG(1)+SIG(2)+SIG(3)

          P  = I1/3.D0

C deviator in Voigt

          S(1)=SIG(1)-P

          S(2)=SIG(2)-P

          S(3)=SIG(3)-P

          S(4)=SIG(4)

          S(5)=SIG(5)

          S(6)=SIG(6)

C J2 = 0.5*s:s with Voigt metric (shear weights)

          s11=S(1); s22=S(2); s33=S(3)

          s12=S(4); s13=S(5); s23=S(6)

          J2 = 0.5D0*( s11*s11 + s22*s22 + s33*s33

     &              + 2.D0*( s12*s12 + s13*s13 + s23*s23 ) )

          SEQ = SQRT(3.D0*J2)

        END SUBROUTINE HYDRO_DEV_EQ

C=======================================

C Utility: recomposition from deviator S and mean stress P

C=======================================

        SUBROUTINE DEV_PLUS_P(S, P, SIG)

          REAL*8 S(6), SIG(6), P

          SIG(1)=S(1)+P

          SIG(2)=S(2)+P

          SIG(3)=S(3)+P

          SIG(4)=S(4)

          SIG(5)=S(5)

          SIG(6)=S(6)

        END SUBROUTINE DEV_PLUS_P

C=======================================

C  Utility: maximum principal stress via invariants (deviatoric approach)

C=======================================

        SUBROUTINE MAX_PRINCIPAL(SIG, SIG1)

          REAL*8 SIG(6), SIG1

          REAL*8 P, S(6), SEQ

          REAL*8 J2, J3, s11,s22,s33,s12,s13,s23

          REAL*8 R, COS3TH, TH, TWO_PI_OVER_3, PI

          PI = 3.141592653589793D0

          CALL HYDRO_DEV_EQ(SIG,P,S,SEQ)

          s11=S(1); s22=S(2); s33=S(3)

          s12=S(4); s13=S(5); s23=S(6)

C J2 again

          J2 = 0.5D0*( s11*s11 + s22*s22 + s33*s33

     &              + 2.D0*( s12*s12 + s13*s13 + s23*s23 ) )

C J3 = det(s) for symmetric tensor (using Voigt)

          J3 =  s11*s22*s33

     &       + 2.D0*s12*s13*s23

     &       – s11*s23*s23 – s22*s13*s13 – s33*s12*s12

C R = sqrt(2/3)*sqrt(J2)

          IF (J2 .LE. 1.D-24) THEN

             SIG1 = P

             RETURN

          ENDIF

          R = SQRT(2.D0/3.D0)*SQRT(J2)

C cos(3theta) = (3*sqrt(3)/2) * J3 / J2^(3/2)

          COS3TH = (1.5D0*SQRT(3.D0)) * J3 / ( (J2*SQRT(J2)) + 1.D-24 )

          IF (COS3TH .GT. 1.D0) COS3TH = 1.D0

          IF (COS3TH .LT. -1.D0) COS3TH = -1.D0

          TH = (1.D0/3.D0)*ACOS(COS3TH)

          TWO_PI_OVER_3 = 2.D0*PI/3.D0

C principal stresses: sigma_k = P + 2*R*cos(TH + 2π(k-1)/3)

          SIG1 = P + 2.D0*R*DCOS(TH)

C  ensure maximum by checking other two (robustness)

          SIG1 = MAX(SIG1, P + 2.D0*R*DCOS(TH + TWO_PI_OVER_3))

          SIG1 = MAX(SIG1, P + 2.D0*R*DCOS(TH + 2.D0*TWO_PI_OVER_3))

        END SUBROUTINE MAX_PRINCIPAL

C=======================================

      END SUBROUTINE UMAT

The developed UMAT should be compiled in Abaqus/Standard and applied using the defined PROPS sequence, ensuring correct NDI/NSHR handling for three-dimensional or shell elements. During simulation, the parameters DCRIT (typically 0.1–0.3) and DMAX (around 0.99) control the onset and limit of stiffness degradation, maintaining convergence through the smooth degradation law. Mesh-objectivity must be verified by repeating deep-drawing simulations with different mesh densities: results such as punch force–stroke curves and fracture initiation loci should remain nearly invariant (within 5%), confirming the proper length-scale regularization via the ratio ℓ/ℎ. The consistent tangent formulation ensures quadratic convergence of the global Newton–Raphson iterations, which should be verified by monitoring residual reduction during plastic loading. For calibration, the reference fracture strain 𝜀𝑓,0 and triaxiality factor 𝛽 are obtained from notched and unnotched tensile tests, while the regularization parameters 𝛼 and ℓ are tuned using experimental deep-drawing data to ensure mesh-insensitive force–stroke predictions and realistic localization behavior.

Damage in the material is predicted using the Brozzo model, specifically through the maximum principal stress (σ_max) and accumulated plastic strain (ϵ_ρ) derived from the state variable statev(1). The damage criterion, damage_criteria, is calculated as the relative sum of these two parameters. If this value exceeds 1, it signifies that the material has reached a critical failure state (i.e., full damage, damage = 1.0). When the damage criterion is within permissible limits, the damage accumulates incrementally as the deformation progresses. To calibrate the damage parameters, a series of tests spanning stress triaxiality and Lode angle are employed. These tests include smooth and notched tension tests, shear tests, plane strain tests, and biaxial tests. Each of these tests provides essential data for determining the material’s behaviour under different loading conditions, helping to refine the model’s predictions.

The key parameters for calibration include raw data collection, optimization procedures, uncertainty analysis, and cross-validation. Raw data is gathered from tensile and deep-drawing tests, capturing stress and strain responses across different material conditions. The optimization procedure aims to minimize the difference between predicted and observed failure points, adjusting parameters to best fit the experimental data. This is followed by an uncertainty analysis, often employing techniques such as Monte Carlo simulations, to assess the influence of variations in material properties and loading conditions on the model’s predictions. Cross-validation ensures that the model’s predictions of fracture initiation, such as punch force-displacement curves and fracture initiation locations and times, align closely with experimental data. This validation process is crucial for ensuring the model’s reliability in practical applications.

During the deformation process, damage is progressively incorporated into the material’s stress state. As damage accumulates, the material’s stiffness decreases, which allows the model to better simulate the progressive loss of load-bearing capacity. This damage evolution is essential for accurately predicting when fracture will initiate, particularly in complex loading conditions, such as those encountered during deep-drawing processes. The accumulated damage is stored in the state variable statev(1), allowing the model to track damage accumulation over time. This enables a more accurate representation of material behaviour during deformation and ensures that the model can predict fracture initiation with greater precision.

The custom UMAT subroutine for Abaqus, developed in Fortran and stored in the file brozzo_umat.f, is crucial for implementing the Brozzo damage model. To begin, a material is created in the Property Module within Abaqus, labelled as “BrozzoMaterial”. Essential parameters such as Yield Stress and Strain Threshold are defined in the Mechanical User Material tab. In the Job Module, a new job (e.g., “DeepDrawingJob”) is created, and the subroutine file brozzo_umat.f is specified in the User Subroutine File field. Once the job is submitted, Abaqus runs the simulation, and the results can be viewed and analyzed in the Visualization Module, where the damage zones and fracture initiation points are carefully examined. The UMAT subroutine calculates damage accumulation using a specific relation between the maximum principal stress (σ_max) and accumulated plastic strain (ϵ_ρ). The damage criterion is calculated as follows (9):

If the damage criterion exceeds 1, the material is considered to have failed, with the damage set to 1.0. If the damage criterion is below 1, damage accumulates gradually depending on the value of the criterion. This calculation flow ensures that the model accurately reflects the progressive accumulation of damage during the deep drawing process and predicts failure with high precision. As for the numerical modelling, an example of comparing the results of the distribution of the damage index and the equivalent plastic strain at fracture for the three finite element models is worth considering (Figure 8). Such an analysis assesses the accuracy of modelling of physical fracture processes.

The figure shows the results of numerical simulations for three finite element models (Hill48-NAFR, Hill48-AFR-S, Hill48-AFR-R) during the forming process. The influence of material anisotropy and rolling direction is marked by different levels of damage and deformation, with the Hill48-NAFR model showing the highest damage localisation (Formula 3). The Hill48 anisotropic yield criterion affects damage prediction by considering the directional dependency of material behaviour under various stress states, as observed in sheet metal forming. Incorporating material anisotropy provides more precise predictions of fracture initiation and damage evolution in processes such as deep drawing, where the material reaction differs based on orientation relative to the rolling direction.

Figure 9 illustrates the simulated damage distribution across the material, with fracture locations marked as black “x” symbols, which are based on experimental data. This plot provides insights into the areas where damage is most concentrated, highlighting the accuracy of the fracture prediction model. The Plastic Strain (εp) Contour at Failure plot shows the distribution of plastic strain across the material at the point of failure, with the experimental fracture locations overlaid, allowing for a comparison between simulated strain and observed fracture points. Additionally, the Punch Force vs. Stroke plot depicts the punch force as a function of stroke, with error bars representing the variability in the measurement. This plot offers a quantitative measure of the relationship between force and displacement, providing further validation of the simulation model’s accuracy and its correlation with experimental data.

Numerical modelling of the UMAT subroutine for the Brozzo criterion should be performed in Abaqus using a finite element model of the deep drawing process. Such a model allows predicting damage accumulation and assessing the fracture zones in the sheet material. The geometric model is an axisymmetric plate that corresponds to the typical dimensions of the workpiece to be drawn. The material is represented as an isotropic elastoplastic material with damage, and a linear strengthening model is used to describe the plastic flow.^35,36 The contact interactions between the tool and the workpiece consider friction, which is described by the friction coefficient.³⁷ To ensure the stability of the calculations, a rigid-surface approach is used, where the workpiece is modelled by shell finite elements. In the case of dynamic process modelling, an explicit integration algorithm is used to efficiently simulate large plastic deformations. The damage field evolution at several phases of the deep-drawing process illustrated the progression of damage and the sites of fracture initiation, especially around the punch and die interfaces, underscoring the Brozzo model’s precision in forecasting failure under complex loading conditions (Figure 10).

The Brozzo model is defined by its integration of hydrostatic stress, a vital element in forecasting fracture under elevated triaxiality, as observed in deep drawing processes. Conversely, models such as MHC/Hosford-Coulomb, which predominantly emphasise main stress, inadequately account for hydrostatic stress effects and hence falter under conditions characterised by substantial triaxiality. This leads to increased inaccuracies, often between 10% and 15%, in forecasting fracture start in intricate geometries with diverse stress states. The Brozzo model regularly exhibits enhanced performance in numerical simulations under biaxial and multiaxial loading situations. In simulations characterised by shear-dominated failure or intricate multi-directional loads, the Brozzo model precisely forecasts damage progression and the initiation of failure. In contrast to DF2014/2016 or Bai-Wierzbicki, which are primarily tailored for shear or shear-dominated conditions, the Brozzo model’s capacity to incorporate hydrostatic stress results in a more precise representation of failure in triaxial stress environments, where the aforementioned models often overestimate material strength, causing fracture predictions to diverge by 10–15% from experimental outcomes.

Sensitivity analysis on material properties was conducted to further evaluate the robustness of the Brozzo model. The analyses encompassed diverse yield stress, strain hardening factors, and fracture strain thresholds to evaluate the model’s efficacy across a spectrum of material characteristics. The Brozzo model exhibited significant parameter sensitivity, especially regarding strain hardening and fracture strain. However, it maintained stability and precision in forecasting failure across various loading situations. When the material properties were altered within standard industrial parameters, the Brozzo model reliably forecasted failure locations and damage zones with low deviation, but alternative models such as GTN or MHC/Hosford-Coulomb exhibited greater variability in their predictions.

To compare the performance of these models, numerical simulations were conducted for a deep-drawing process involving cold-rolled steel. The models were evaluated based on their ability to predict fracture initiation points and damage zones. The Brozzo model consistently outperformed the MMC/Hosford-Coulomb, Bai-Wierzbicki, and GTN models, particularly in scenarios characterized by high triaxial stress conditions, such as those found in deep drawing processes. The MMC/Hosford-Coulomb and Bai-Wierzbicki models showed larger error margins when predicting fracture initiation under such conditions, with discrepancies of 10–15% compared to experimental data. The GTN model, while accurate in predicting void growth and porosity-based fracture, exhibited greater error margins in deep drawing simulations, where shear stresses and hydrostatic stress were dominant (Table 3).

Table 3: The error margins for each model in predicting fracture initiation during deep drawing.
Model	Error Margin (Fracture Initiation)	Strengths	Limitations
Brozzo	~5%	Accurate in high triaxial stress	Less effective in shear-dominated failures
MMC/Hosford-Coulomb	~10-15%	Good for shear-dominated failures	Not ideal for high triaxial stress
Bai-Wierzbicki	~10-12%	Accurate for low triaxiality	Poor performance in high triaxiality
GTN	~10-15%	Suitable for porous materials	Less effective in high-triaxial stress environments

These results underline the importance of selecting the appropriate ductile fracture model based on the material’s behaviour under specific loading conditions. Models such as Brozzo, which incorporate hydrostatic stress effects, are more accurate for deep drawing simulations where high triaxial stress plays a crucial role in fracture initiation. Additionally, simulations with non-standard geometries, like intricate flanges and thicker sheet metals, were employed to validate the Brozzo model under more demanding industrial settings. The geometries produced non-uniform stress distributions, with the Brozzo model demonstrating exceptional accuracy in forecasting fracture starts in stress concentration areas. The mesh refinement investigation, concentrating on regions adjacent to the punch-die interface and blank holder, demonstrated that the Brozzo model proficiently managed localized deformation and damage accumulation, yielding dependable forecasts even in locations of elevated stress.

In the deep drawing simulations, realistic tool geometries were incorporated to replicate industrial conditions. The tooling setup included a punch, die, and blank holder, which were meticulously designed to simulate the interaction between the sheet material and the tooling during the deep drawing process. The punch was modeled with a diameter of 100 mm, consistent with standard industrial practices, while the die geometry was designed to match the punch size, ensuring proper interaction during material deformation. The blank holder was engineered to apply restraining forces uniformly around the perimeter of the sheet, simulating the effects of draw beads in regulating material flow. These restraining forces were calibrated to prevent issues such as wrinkling or excessive thinning during the deformation process, helping to maintain the integrity of the formed part.

The draw-bead modeling was implemented by applying equivalent restraining forces around the blank holder. These forces played a critical role in regulating material flow and ensuring the sheet metal maintained the required form throughout the drawing process. By applying these forces uniformly around the blank holder, the model closely mimicked the action of draw beads in real-world applications, ensuring smooth material deformation and preventing defects. The friction law used in the simulations was based on the Coulomb friction model, where the frictional force between the material and the tooling was calculated as the product of the friction coefficient (μ) and the normal force (N) between the surfaces. The friction coefficient was set within the range of 0.05 to 0.2, reflecting typical values used in metal forming processes. A friction coefficient of 0.1 was selected as the baseline based on empirical data from similar material-tool interactions. This coefficient was varied within the given range to perform a sensitivity analysis and assess its impact on failure predictions.

The sensitivity analysis of the friction coefficient demonstrated that friction has a significant influence on material behavior during the deep drawing process. When the friction coefficient was set to lower values (e.g., μ = 0.05), the material exhibited more uniform flow with reduced localized deformation, leading to delayed fracture initiation. In contrast, higher friction values (e.g., μ = 0.2) caused more resistance to material flow, increasing localized deformation near the punch and die interfaces, which accelerated the damage accumulation and led to earlier fracture initiation. The sensitivity analysis clearly indicated that friction significantly impacts the stress distribution and failure progression during deep drawing, highlighting the need to carefully select an appropriate friction coefficient for accurate predictions.

In terms of failure predictions, the friction coefficient was found to be a key factor in determining the location and timing of fracture initiation. At higher friction values, localized damage was observed near the punch-die contact zone, where the material is subjected to the highest levels of stress and strain. This resulted in earlier fracture initiation and more concentrated damage zones. Conversely, at lower friction values, the material exhibited smoother flow, leading to a more evenly distributed stress and strain, which delayed failure and reduced the likelihood of localized fractures. The results of the sensitivity analysis emphasized the importance of accurately modeling frictional interactions in the deep drawing process to ensure precise failure predictions.

The study confirmed that the choice of friction law and friction coefficient directly influences the simulation outcomes in terms of material flow, stress distribution, and fracture predictions. The analysis showed that friction plays a critical role in determining the overall deformation behavior and failure progression in deep drawing processes. By varying the friction coefficient and assessing its impact on fracture initiation, the study underscored the importance of selecting the appropriate friction model to ensure the reliability and accuracy of failure predictions during metal forming operations. These findings have direct implications for tooling design and process optimization in industries such as automotive and aerospace manufacturing, where material performance and product quality are paramount.

At each step of the analysis, Abaqus sends the stress tensor to UMAT to calculate the maximum principal stress and update the plastic strain. This provides a more accurate assessment of damage accumulation in the material. The damage criterion is defined as the sum of the relative value of the maximum principal stress and the accumulated plastic strain, and when this sum exceeds a critical threshold, the material is damaged. After the calculation is complete, the damage distribution in the contact zone is analysed. The results show the localisation of damage, which predicts crack formation, as well as identifies the area’s most vulnerable to damage. To summarise the results, Table 4 shows the main modelling parameters. Thus, numerical modelling confirms the effectiveness of using the UMAT subroutine to implement the Brozzo criterion in the deep drawing process. At the same time, the model adequately reflects the physical processes of material fracture under deep drawing conditions.

Table 4: Main modelling parameters.
Parameter	Value	Description
Model type	Axisymmetric	Incorporates the symmetry of the deep drawing process
Material type	Isotropic elastoplastic	Incorporates plastic hardening and damage
Contact boundary condition	Friction (μ)	Describes the interaction of the workpiece with the tools
Time integrator	Explicit/implicit	The choice depends on the modelling conditions
Stressed state	Stress tensor (σij)	Determines the mechanical condition at each point
Damage criterion	(Formula 2)	Determines the moment of destruction
Localisation of damage	In contact areas	Suitable for areas of highest stress and plastic deformation
Modelling result	Destruction forecast	Identifies critical areas and material behaviour
Source: Compiled by the authors.

The UMAT subroutine and Brozzo model can be modified for materials with varying mechanical properties by altering parameters such as yield stress, strain hardening, and fracture strain. For materials exhibiting greater strength or varying strain rate sensitivity, these parameters may be calibrated using experimental data or inverse optimisation methods. The Brozzo model can be adjusted to accommodate various forming and loading scenarios by integrating specific boundary conditions, such as temperature effects or dynamic loading, to more precisely represent high-speed forming or thermal settings. This versatility enables the model to be utilised across various materials and forming methods, enhancing its applicability in multiple industrial contexts.

The Brozzo-based subroutine facilitates accurate forecasts of damage accumulation and fracture initiation in deep drawing, rendering it relevant for industrial metal forming processes. Integrating this model into Abaqus allows producers to mimic material behaviour under realistic settings, facilitating the optimisation of process parameters including die design, punch speed, and material selection. These forecasts can be seamlessly integrated into quality control processes, offering a robust mechanism for the early identification of probable problems and the reduction of material waste. This integration facilitates process optimisation and guarantees enhanced product consistency, rendering it advantageous for sectors such as automotive manufacture and aerospace production, where material performance is paramount.

Discussion

The results of the study demonstrated that the implementation of the Brozzo criterion in Abaqus provides an accurate prediction of the moment of fracture during deep drawing. Similar to the study by Jia et al.,³⁸ where the phenomenological model was tested for aluminium alloys, the results obtained confirm the high accuracy in determining the damage zones. However, this study focuses on more complex geometric shapes, which expands the applicability of the criterion. In general, the implementation of the Brozzo criterion made it possible to assess the effect of plastic deformation on fracture. In contrast to the approach of Zhang et al.,³⁹ where a probabilistic model was used to predict the ultimate strain, the method used in the current work is deterministic. Moreover, the current study is focused on predicting damage during deep drawing, while the influence of initial porosity was analysed in the study. Despite the different approaches, the results obtained are consistent in the accuracy of the fracture zones.

Similar to the study by Tandoğan and Yalcinkaya,⁴⁰ the present research used Abaqus for numerical modelling of fracture. However, instead of a custom element, its standard capabilities were used together with a custom subroutine to evaluate damage accumulation. This incorporated the influence of the stress-strain state on the development of fracture without introducing additional elements into the model. In contrast to the study by Akitarak,⁴¹ which used the user subroutines VUMAT, VUHARD (User subroutine to define hardening) and VUSDFLD (User-defined material point field variable), the implementation in this work is based on UMAT. This performed the analysis in the Abaqus environment with the ability to control the stress state and damage evolution in the material. Both approaches confirm the effectiveness of extending the standard Abaqus functionality for fracture prediction but differ in the choice of numerical strategy.

The results of this study demonstrated that the analysis of ductile fracture models provides an accurate prediction of the moment of failure during deep drawing, which correlates with the results of Choi et al.,⁴² where a specialised model for shell elements was proposed. Although the model of the study focuses on thin-walled structures and was validated under different conditions, the use of both approaches for ductile fracture in various geometric shapes provides a high level of accuracy in predicting damage zones. The model of the current study, which focuses on deep drawing, can be supplemented with an extended analysis of geometric features, which expands the possibilities of its application in various engineering problems. In turn, Khaboushani et al.⁴³ proposed a strain energy-based fracture criterion for assessing the formation of fractures during deep drawing, which showed high accuracy in predicting fractures due to high radial stresses in the bowl walls. This approach incorporates various fracture mechanisms based on the strain energy. Compared to the implementation of the Brozzo criterion through UMAT in the current study, the aforementioned method uses a different approach to fracture prediction. Despite the difference in methods, both approaches demonstrated high accuracy in predicting the failure moments during deep drawing.

The results obtained indicated a high accuracy in predicting the moment of failure in the process of deep drawing using ductile fracture models in Abaqus. Similarly, Ben Said et al.³⁶ developed a mathematical model for predicting the moment of crack formation during deep drawing using anisotropic plasticity equations and custom subroutines in Abaqus. Both studies use UMAT, which confirms the effectiveness of the modelling. On the other hand, the study by Clayton⁴⁴ proposed a method for estimating shear strain and phase transitions during fracture. Although this approach is focused on other fracture mechanisms, it also demonstrated high accuracy in predicting the moments of failure. Thus, the results of both studies are consistent, as both methods confirm the effectiveness of numerical models for predicting failure moments, albeit with different focuses on failure mechanisms.

As in the study by Van den Abeele,⁴⁵ the present article analysed the use of the GTN model, which allows for detailed modelling of ductile fracture by considering the porosity of the material that changes during deformation. However, in contrast to the above work, where the emphasis is on many parameters for model calibration, which complicates their practical application, the current study addressed simpler approaches that reduce the complexity of their implementation in industrial conditions. Regarding a study by Hosseini Mansoub et al.,⁴⁶ a three-axis stress state and damage function was used to predict the location of the failure, which is similar to the approach in the current work. However, in the above study, the UMAT subroutine is used to link the models to the SFLD (Stress Forming Limit Diagram), while in the current work, UMAT is used to implement the Brozzo model failure criterion, which allows for direct consideration of damage accumulation during material deformation during deep drawing.

Compared to the study by Ma et al.,⁴⁷ where a phase-field model was used to predict the ductile fracture of shell structures with an energy approach, the current study analysed possibilities of the Brozzo model in predicting material failure during deep drawing, addressing the accumulation of damage due to plastic deformation. The results of the studies are consistent in terms of consideration of plastic deformation but differ in their approaches. Regarding the study by Basak and Panda,⁴⁸ although it used the Bai & Wierzbicki model to predict the fracture limits of aluminium alloy, the current study also considered this model as important for crack prediction, especially in processes such as deep drawing, but the emphasis was on the Brozzo model. The results are consistent in that the Bai & Wierzbicki model is a micromechanical model that incorporates the Lode angle parameter and triaxial stress, making it useful in processes with high cracking rates.

The results also demonstrate the importance of using models to predict failure during deep drawing, which is similar to the study by Ganjiani,⁴⁹ where a model that incorporates triaxial stresses is shown. However, in the aforementioned study, VUMAT in Abaqus/Explicit was used, while in the current work, UMAT was used for deterministic prediction. Similar to the study by Wu and Lou,⁵⁰ where the DF2016 model showed high accuracy when considering the strain trajectory, the current model also demonstrates a high level of accuracy, in particular in the context of deep drawing. Both approaches demonstrate the effectiveness of models for fracture prediction, although the present study focuses on a specific deformation process. The current results based on the Brozzo model have certain differences in comparison with the study by Fagerhøi and Bergsbakken,⁵¹ which analysed the parameters of the Concrete Damaged Plasticity model for concrete in Abaqus. Although both studies analysed modelling the failure of materials, the approach of the present study aimed to predict the failure during deep drawing, while the other study analysed the effect of parameter changes on the accuracy of the modelling. However, despite the difference in materials and model application, the results of both studies are consistent in that the accuracy of fracture prediction depends on the correct setting of the model parameters and its ability to adequately represent the material behaviour under different loading conditions.

The conclusions of this study estimated the moment and place of failure, through the integration of stresses and strains. At the same time, the study by Aksen et al.⁵² used ductile damage functions for two-phase steels, as well as a combined plasticity model to assess fracture initiation. The common feature is that the works highlight the accuracy of fracture location prediction, however, the current results are focused on a specific deep drawing process, while the present work considers other deformation conditions and material types. In addition, the present study analysed damage accumulation, while the study by Park et al.⁵³ uses the Hosford-Coulomb model for EH36 steel, focusing on modelling fracture initiation through different types of specimens and determining the hardening parameters using the Swift-Voce function. Both studies use Abaqus and custom subroutines for numerical modelling, but the current approach is more focused on predicting fracture during deep drawing using the Brozzo model, while the other study focuses on steels and modelling fracture initiation through the Hosford-Coulomb model.

Notably, the study by Watanabe et al.⁵⁴ used an equation to predict ductile fracture based on the integration of stresses and strains to accurately determine the location of fracture initiation in cold-forming processes. Similar to the work presented in this paper, the current study also analyses the Cockcroft & Latham model for fracture prediction, but with an improvement in the form of the Brozzo model, which takes into account hydrostatic stresses, which is important for processes with high triaxial stresses, such as in deep drawing. Compared to the classical criteria, this allows for more accurate fracture prediction. Additionally, the study by Zhang et al.⁵⁵ uses the CDM (Continuum Damage Mechanics) model to predict ductile fracture in the single point forming process, which considers the evolution of damage depending on the stress state. The current study applies the Brozzo model, which focuses on the integration of hydrostatic stresses to more accurately predict fracture during deep drawing. Both approaches use Abaqus/Explicit but differ in the types of processes and damage models.

Thus, the results of the study demonstrated that the implementation of the Brozzo criterion in Abaqus accurately predicted the moment of failure in the process of deep drawing, confirming the accuracy of detecting damage zones. At the same time, modelling using UMAT provides accurate control of damage during material deformation, which makes the Brozzo method effective for predicting failure in complex processes such as deep drawing. The Brozzo model, while effective in predicting ductile fracture under high triaxial stress conditions, faces significant limitations when applied to shear-dominated scenarios or low triaxiality. This is primarily because the Brozzo model focuses on the accumulation of damage driven by hydrostatic stress, which becomes less relevant when shear stress dominates the deformation process. In cases with low triaxiality, where shear stresses are more prevalent and hydrostatic stress plays a minimal role, the Brozzo model’s predictions can be inaccurate. To address this limitation, hybrid criteria incorporating both the Brozzo model’s hydrostatic stress-based approach and models suited for shear-dominated conditions, such as the GTN or MHC, could offer improved fracture predictions. These hybrid models would allow for a more balanced representation of failure mechanisms across a broader range of stress states, making them potentially more reliable for various industrial forming processes.

Conclusions

The study analysed and compared the main models of ductile fractures, such as Cockcroft & Latham, Brozzo, Oh, Johnson & Cook, Bai & Wierzbicki, Modified Mohr-Coulomb, Lou-Huh, Gurson-Tvergaard-Needleman and Lemaitre. In the course of analysing each of the models, their advantages and limitations were considered, as well as their application for modelling the deep drawing process. The study determined that the Brozzo model is the most effective for predicting material behaviour under plastic deformation, as it allows for an accurate assessment of damage accumulation, considering the stresses and strains in the material. The Brozzo model proved to be the most suitable for modelling complex processes, such as deep drawing, compared to other approaches that have limitations in accounting for high levels of strain or complex geometries. This confirms its strong potential for application in real-world manufacturing, where accurate prediction of material failure is critical to ensuring the quality and durability of finished products.

An important stage of the study was the creation of a specialised subroutine to implement the Brozzo model in Abaqus. This subroutine predicts the moments of material failure during deep drawing, considering the main parameters of plastic deformation, interaction with stress states and other factors that affect the development of damage. The implementation of this model improves the accuracy of modelling, as it enables detailed prediction of fracture zones in processes with high levels of deformation, such as deep drawing. In addition, a UMAT subroutine was programmed to predict the accumulation of damage in the material, which accounts for the evolution of damage at each stage of loading. The integration of this subroutine into Abaqus made it possible to perform numerical modelling of the deep drawing process, with the prediction of damage zones and the moment of failure.

For further research, it is recommended to expand the scope of experimental studies to include different types of materials and product geometries, which will allow the model to be evaluated in a broader context. It is worth improving numerical methods to increase the efficiency and accuracy of modelling, as well as to investigate the effect of temperature and strain rate, which will enhance the model’s ability to predict fracture in complex technological processes.

References

1. Giang LD, Tran DH, Dũng LQ, Nguyen VC, Nguyen VH. Prediction of earing in deep drawing of the cylindrical cup of stainless steel SUS304 by numerical simulation. J Mil Sci Technol. 2023;86:129–36. https://doi.org/10.54939/1859-1043.j.mst.86.2023.129-136
2. Hà XG, Nguyễn TT, Lê VA. Numerical simulation of the springback of copper alloy thin sheets after the forming process. Hong Duc Univ J Sci. 2024;71(11):5–14. https://doi.org/10.70117/hdujs.2.2024.734
3. Thành V. Advanced phase-field method used to enhance the accuracy of damage simulation in typical structures. Sci Technol Civil Eng J. 2024;18(4V):122–36. https://doi.org/10.31814/stce.huce2024-18(4V)-10
4. Laboubi S, Boussaid O, Zaaf M, Ghennai W. Numerical investigation and experimental validation of lemaitre ductile damage model for DC04 steel and application to deep drawing process. Int J Adv Manuf Technol. 2023;126:2283–94. https://doi.org/10.1007/s00170-023-11244-0
5. Zhu X, Lv Z, Shuai J, Shi L, Yu S, Xia G. Ductile fracture criterion parameter calibration and analysis: X80 pipe steel. Pressure Vess Piping Conf. 2024;4:1–7. https://doi.org/10.1115/PVP2024-123037
6. Dey S, Kiran R. A Data-driven geometry-specific surrogate model for forecasting the load-displacement behavior until ductile fracture. Int J Fract. 2025;249:31. https://doi.org/10.1007/s10704-025-00839-1
7. Talebi-Ghadikolaee H, Moslemi Naeini H, Talebi Ghadikolaee E, Mirnia MJ. Predictive modeling of damage evolution and ductile fracture in bending process. Mater Today Commun. 2022;31:103543. https://doi.org/10.1016/j.mtcomm.2022.103543
8. Talebi-Ghadikolaee H, Moslemi Naeini H, Mirnia MJ, Mirzai MA, Gorji H, Alexandrov S. Ductile fracture prediction of AA6061-T6 in roll forming process. Mech Mater. 2020;148:103498. https://doi.org/10.1016/j.mechmat.2020.103498
9. Wang D, Xu Z, Han Y, Huang FA. A ductile fracture model incorporating stress state effect. Int J Mech Sci. 2022;241:107965. https://doi.org/10.1016/j.ijmecsci.2022.107965
10. Guo Y, Xie Y, Wang D, Du L, Zhao J. An improved damage-coupled viscoplastic model for predicting ductile fracture in aluminum alloy at high temperatures. J Mater Process Technol. 2021;296:117229. https://doi.org/10.1016/j.jmatprotec.2021.117229
11. Ferreira J, Pacheco L. 2020. UMAT-ABAQUS_library. Available at: https://github.com/jpsferreira/UMAT-ABAQUS_library
12. Dolzhenko N, Mailyanova E, Assilbekova I, Konakbay Z. Analysis of meteorological conditions significant for small aviation and training flights at the airfield “Balkhash” for planning and flight safety purposes. News Natl Acad Sci RK, Ser Geol Tech Sci. 2021;2(446):62–67. https://doi.org/10.32014/2021.2518-170X.35
13. Dolzhenko NA, Mailyanova EN, Assilbekova IZ, Konakbay ZE. Design features of modern flight simulation devices, mobility systems and visualization systems. News Natl Acad Sci RK, Ser Geol Tech Sci. 2021;3(447):17–21. https://doi.org/10.32014/2021.2518-170X.56
14. Wójcik W, Kalizhanova A, Kulyk YA, Knysh BP, Kvyetnyy RN, Kulyk AI, et al. The method of time distribution for environment monitoring using unmanned aerial vehicles according to an inverse priority. J Ecol Eng. 2022;23(11):179–187. https://doi.org/10.12911/22998993/153458
15. Sekenov B, Smailov N, Tashtay Y, Amir A, Kuttybayeva A, Tolemanova A. Fiber-optic temperature sensors for monitoring the influence of the space environment on nanosatellites: A review. Mech Mach Sci. 2024;167 MMS:371–380. https://doi.org/10.1007/978-3-031-67569-0_42
16. Smailov N, Orynbet M, Nazarova A, Torekhan Z, Koshkinbayev S, Yssyraiyl K, et al. Optimization of fiber-optic sensor performance in space environments. Informatyka Autom Pomiary Gospod Ochr Srodow. 2025;15(2):130–134. https://doi.org/10.35784/iapgos.7200
17. Mikhailov P, Ualiyev Z, Kabdoldina A, Smailov N, Khikmetov A, Malikova F. Multifunctional fiberoptic sensors for space infrastructure. East-Eur J Enterp Technol. 2021;5(5-113):80–89. https://doi.org/10.15587/1729-4061.2021.242995
18. Mazakova A, Jomartova S, Mazakov T, Brzhanov R, Gura D. The use of artificial intelligence to increase the functional stability of UAV systems. Int Rev Aerosp Eng. 2024;17(3):98–106. https://doi.org/10.15866/irease.v17i3.25067
19. Smailov N, Akmardin S, Ayapbergenova A, Ayapbergenova G, Kadyrova R, Sabibolda A. Analysis of VLC efficiency in optical wireless communication systems for indoor applications. Informatyka Autom Pomiary Gospod Ochr Srodow. 2025;15(2):135–138. https://doi.org/10.35784/iapgos.6971
20. Dolzhenko N, Mailyanova E, Toluev Y, Assilbekova I. Influence of system errors in meteorological support on flights safety. News Natl Acad Sci RK, Ser Geol Tech Sci. 2020;5(443):1–89.
21. Dolzhenko N, Assilbekova I, Konakbay Z, Garmash O, Muratbekova G. Organization of transport services and transport process safety. Period Polytech Transp Eng. 2025;53(3):277–291. https://doi.org/10.3311/PPtr.38137
22. Anishchenko A, Kukhar V, Artiukh V, Olga A. Superplastic forming of shells from sheet blanks with thermally unstable coatings. MATEC Web Conf. 2018;239:06006. https://doi.org/10.1051/matecconf/201823906006
23. Karnaukh SG, Markov OE, Aliieva LI, Kukhar VV. Designing and researching of the equipment for cutting by breaking of rolled stock. Int J Adv Manuf Technol. 2020;109(9-12):2457–64. https://doi.org/10.1007/s00170-020-05824-7
24. Anishchenko A, Kukhar V, Artiukh V, Arkhipova O. Application of G. Lame’s and J. Gielis’ formulas for description of shells superplastic forming. MATEC Web Conf. 2018;239:06007. https://doi.org/10.1051/matecconf/201823906007
25. Andrenko P, Rogovyi A, Hrechka I, Khovanskyi S, Svynarenko M. Characteristics improvement of labyrinth screw pump using design modification in screw. J Phys Conf Ser. 2021;1741(1):012024. https://doi.org/10.1088/1742-6596/1741/1/012024
26. Sergiyenko O, Hernández Balbuena D, Tyrsa V, Rosas Méndez PLA, Lopez MR, Hernandez W, Podrygalo M, Gurko A. Analysis of jitter influence in fast frequency measurements. Measurement. 2011;44(7):1229–42. https://doi.org/10.1016/j.measurement.2011.04.001
27. Mazakova A, Jomartova S, Wójcik W, Mazakov T, Ziyatbekova G. Automated linearization of a system of nonlinear ordinary differential equations. Int J Electron Telecommun. 2023;69(4):655–660. https://doi.org/10.24425/ijet.2023.147684
28. Ratov BT, Mechnik VA, Bondarenko NA, Kolodnitskyi VM, Gevorkyan ES, Nerubaskyi VP, Gusmanova AG, Fedorov BV, Kaldibaev NA, Arshidinova MT, Kulych VG. Features structure of the Сdiamond‒(WC‒Co)‒ZrO2 composite fracture surface as a result of impact loading. J Superhard Mater. 2023;45(5):348–59. https://doi.org/10.3103/S1063457623050088
29. Adjamskiy S, Kononenko G, Podolskyi R, Baduk S. Studying the influence of orientation and layer thickness on the physico-mechanical properties of Co-Cr-Mo alloy manufactured by the SLM method. Science Innov. 2022;18(5):85–94. https://doi.org/10.15407/scine18.05.085
30. Zaginaylov GI, Shcherbinin VI, Schüenemann K, Thumm MK. Influence of background plasma on electromagnetic properties of Cold gyrotron cavity. IEEE Trans Plasma Sci. 2006;34(3 PART 1):512–17. https://doi.org/10.1109/TPS.2006.875760
31. Piskunov VG, Goryk AV, Cherednikov VN. Modeling of transverse shears of piecewise homogeneous composite bars using an iterative process with account of tangential loads. 1. Construction of a model. Mech Compos Mater. 2000;36(4):287–96. https://doi.org/10.1007/BF02262807
32. Piskunov VG, Gorik AV, Cherednikov VN. Modeling of transverse shears of piecewise homogeneous composite bars using an iterative process with account of tangential loads. 2. Resolving equations and results. Mech Compos Mater. 2000;36(6):445–52. https://doi.org/10.1023/A:1006798314569
33. Quach H, Kim YS. Effect of non-associated flow rule on fracture prediction of metal sheets using a novel anisotropic ductile fracture criterion. Int J Mech Sci. 2021;195:106224. https://doi.org/10.1016/j.ijmecsci.2020.106224
34. Ben Said L, Allouch M, Wali M, Dammak F. Numerical formulation of anisotropic elastoplastic behavior coupled with damage model in forming processes. Mathematics. 2023;11(1):204. https://doi.org/10.3390/math11010204
35. Cherniha R, King JR, Kovalenko S. Lie symmetry properties of nonlinear reaction-diffusion equations with gradient-dependent diffusivity. Commun Nonlinear Sci Numer Simulat. 2016;36:98–108. https://doi.org/10.1016/j.cnsns.2015.11.023
36. Yakovlev SV, Valuiskaya OA. Optimization of linear functions at the vertices of a permutation polyhedron with additional linear constraints. Ukr Math J. 2001;53(9):1535–45. https://doi.org/10.1023/A:1014374926840
37. Babachenko OI, D’Omina KH, Kononenko HA, Dement’Yeva ZhA, Podol’S’Kyy RV, Safronova OA. Influence of cooling rate at hardening of continuous casting blank on parameters of dendritic structure of carbon steel with 0.54% C. Metallofiz Noveishie Tekhnol. 2021;43(11):1537–1551. https://doi.org/10.15407/mfint.43.11.1537
38. Jia Z, Mu L, Guan B, Qian LY, Zang Y. Experimental and numerical study on ductile fracture prediction of aluminum alloy 6016-T6 sheets using a phenomenological model. J Mater Eng Perform. 2021;31:867–81. https://doi.org/10.1007/s11665-021-06248-4
39. Zhang Y, Shen F, Zheng J, Muenstermann S, Li T, Han W, Huang S. Ductility prediction of HPDC aluminum alloy using a probabilistic ductile fracture model. Theor Appl Fract Mech. 2022;119:103381. https://doi.org/10.1016/j.tafmec.2022.103381
40. Tandoğan İT, Yalcinkaya T. Development and implementation of a micromechanically motivated cohesive zone model for ductile fracture. Int J Plast. 2022;158:103427. https://doi.org/10.1016/j.ijplas.2022.103427
41. Akitarak GY. Ductile failure analysis during backward flow forming processes. Ankara: Middle East Technical University; 2024.
42. Choi S, Park T, Kim H, Nam B, Ye B, Kim D. Ductile fracture prediction in thin-walled structures through a novel damage model. Heliyon. 2024;10(23):e40849. https://doi.org/10.1016/j.heliyon.2024.e40849
43. Khaboushani M, Aminzadeh A, Parvizi A. A novel estimation of tearing limit in deep drawing process based on strain energy; experimental characterization and numerical validation. Int J Adv Manuf Technol. 2022;123:927–42. https://doi.org/10.1007/s00170-022-10158-7
44. Clayton JD. Analysis of adiabatic shear coupled to ductile fracture and melting in viscoplastic metals. arXiv. 2025;1:1–36. https://doi.org/10.48550/arXiv.2502.10625
45. Van den Abeele F. Parameter calibration for continuum damage mechanics models to simulate ductile fracture of high strength pipeline steels. ASME Int Conf Ocean Offshore Arctic Eng. 2019;4:14. https://doi.org/10.1115/OMAE2019-96316
46. Hosseini Mansoub F, Basti A, Darvizeh A, Zajkani A. Stress-based forming limit diagrams (SFLD) considering strain rate effect and ductile damage phenomenon. Int J Mater Res. 2020;111(2):136–45. https://doi.org/10.3139/146.111856
47. Ma R, Sun W, Tong G. Phase-field model for ductile fracture in the stress resultant geometrically exact shell. Int J Numer Methods Eng. 2024;125(13):e7462. https://doi.org/10.1002/nme.7462
48. Basak S, Panda S. Use of uncoupled ductile damage models for fracture forming limit prediction during two-stage forming of aluminum sheet material. J Manuf Process. 2023;97:185–99. https://doi.org/10.1016/j.jmapro.2023.04.042
49. Ganjiani M. A damage model for predicting ductile fracture with considering the dependency on stress triaxiality and lode angle. Eur J Mech A Solids. 2020;84(4):104048. https://doi.org/10.1016/j.euromechsol.2020.104048
50. Wu P, Lou Y. Stress-based fracture model to describe ductile fracture behavior in various stress states. Fatigue Fract Eng Mater Struct. 2024;48(3):1200–14. https://doi.org/10.1111/ffe.14548
51. Fagerhøi S, Bergsbakken AK. Abaqus FEA with concrete damaged plasticity and its feasibility in recreating laboratory experiments: A numerical analysis and sensitivity study. Oslo: Oslo Metropolitan University; 2023.
52. Aksen A, Şener B, Esener E, Firat M. Evaluation of ductile fracture criteria in combination with a homogenous polynomial yield function for edge splitting damage of DP steels. Mater Test. 2023;65(6):824–43. https://doi.org/10.1515/mt-2022-0359
53. Park SJ, Lee K, Seo JH, Choung J. Ductile fracture prediction of EH36 grade steel based on hosford-coulomb model. Ships Offshore Struct. 2019;14(1):219–30. https://doi.org/10.1080/17445302.2019.1565300
54. Watanabe A, Hayakawa K, Fujikawa S. An anisotropic damage model for prediction of ductile fracture during cold-forging. Metals. 2022;12(11):1823. https://doi.org/10.3390/met12111823
55. Zhang K, Yue ZM, Su CJ, Wang R, Badreddine H. Modelling of ductile damage in single point incremental forming process using enhanced CDM model. IOP Conf Ser Mater Sci Eng. 2022;1270:012022. https://doi.org/10.1088/1757-899X/1270/1/012022

Cite this article as:
Tran DV and Truong CV. Development of an Abaqus UMAT Subroutine for Predicting Ductile Fracture in Deep Drawing Based on the Brozzo Model. Premier Journal of Science 2025;14:100181

Export to EndNote Export to Mendeley