Phase field modeling of ductile fracture and application in metal forming

. A phase field model for ductile fracture is coupled with the Modified Mohr Coulomb (MMC) model for plastic damage evolution and subsequent crack growth. An energy-based damage threshold is applied to control degradation due to ductile damage. The model is implemented through a user subroutine. MMC parameters from the literature are utilized and found to be compatible with the model, accurately reproducing material response curves in a variety of loading conditions for 6016-T4 aluminum alloy. The influence of model parameters is demonstrated and additionally the Nakazima test is simulated to demonstrate the capability of the model predicting the formability of the material through a failure locus. The model is found capable of reproducing experimentally observed crack paths and quantitative material behavior.


Introduction
The phase field methodology has been applied to an increasing variety of fracture cases over the last few years.Based on Griffith's energy balance [1], the variational approach to fracture was developed by Francfort and Marigo [2], providing solutions to several limitations of the original theory.Bourdin et al. [3] regularized this approach through the introduction of a functional to describe the crack surface and this became the basis for the phase field model for brittle fracture, later formalized by Miehe et al. [4].The framework is expanded to encompass more complex mechanisms such as dynamic crack growth (e.g.[5]), anisotropic crack growth (e.g.[6]) and fatigue crack growth (e.g.[7]) among many.
Adapting the phase field method to a ductile fracture framework involves the introduction of a plasticity related component in addition to the elastic energy density and fracture surface energy terms.The defining characteristic of the various models introduced in the literature is the coupling mechanism for plasticity.Duda et al. [8] opted to not couple the plastic damage with phase field evolution, leading to a model described as a brittle fracture model applied to an elastoplastic solid, with limited plastic deformation at the crack tip before failure.Ambati et al. [9] coupled their model by introducing a coupling term that acts through the degradation function.Models in [10,11] act through the introduction of a separate degradation function for the plastic dissipation term.It should be noted that the latter incorporates non-local plastic effects as well in their formulation.In a framework inspired by the phase field approach to fatigue, Yin et al. [12] degrade the fracture toughness through the accumulated plastic strain.To control the degradation over the material loading history, some models introduce a plastic damage threshold, where plastic damage influences the material response once this value is passed.Both plastic strain thresholds (e.g.[13]) and energy thresholds (e.g.[10]) have been employed in the literature.
The micromechanical mechanisms of ductile fracture demonstrate a high dependency on the stress state characterized by stress triaxiality and Lode parameter (see [14,15]).Ductile fracture models incorporating both parameters have been applied to numerical studies of metal forming (e.g.[16,17]).Recently, attempts have been made to incorporate these effects in ductile phase field models.Borden et al. [10] incorporate stress triaxiality through the addition of a section of the Johnson-Cook damage model to the plastic damage accumulation.Li et al. [18] introduced a damage threshold that evolves with the modified Mohr Coulomb model (MMC).Vajari et al. [19] introduced a model that degrades the fracture toughness through the Stress-Weighted Ductile Fracture Model.
While most studies are focused on demonstrating the ability of various ductile models to recreate known crack paths from literature, quantitative studies accurately capturing the loaddisplacement response of the material remain limited.This study aims to simulate a ductile phase field model, capable of reproducing crack paths for various boundary value problems as well as the material response as accurately as possible.In the scope of this study, we introduce the MMC model to the plastic damage evolution, potentially allowing the model to accurately capture the material response in a broad range of loading conditions.
In this work, a non-linear J2 plasticity model is coupled with the phase field paradigm through coupled-temperature displacement elements, with temperature acting as the phase parameter.This is done by employed by exploiting a similarity between the phase field driving force and steady state heat transfer equation (see [20]).The in-built finite strain framework in the commercial finite element (FE) solver Abaqus is employed in conjunction with user subroutines.The framework is assessed through a comparative study with the experimental and numerical data from [21] while utilizing their hardening and MMC parameters.Load displacement curves are corroborated for various specimens.Furthermore, to demonstrate the utility of the phase field framework in metal forming processes, we utilize the Nakazima test [22], which is employed to determine a material's formability for deep drawing operations.Sheet specimens of various dimensions are subjected to a punch test, with the results subsequently used to plot a failure limit curve for the studied material.
The paper is organized as follows.Section 2 addresses the theory behind the ductile fracture phase field model employed in this paper and its various features.Section 3 covers several benchmark examples, as well as their load-displacement curves recreated through simulations.Section 4 includes the Nakazima test simulations and the achieved failure limit curve.The paper is concluded with a summary.

Theory and Methodology
Phase field fracture relies on the introduction of an auxiliary variable known as the phase field parameter ϕ ∈ [0,1] where a value of 0 represents undamaged material and 1 represents a complete loss of stress carrying capacity.The incorporation of plasticity in the phase field paradigm is based on an additional plastic term introduced to the variational form of Griffith's energy balance, as where E s is the fracture surface energy, while E e and E p refer to the elastic and plastic energy components of the total energy functional, respectively.Both the fracture and elastic contributions may be studied in more detail from [4].The plastic contribution is thus defined by where W p is the energy corresponding to the plastic damage and g is the degradation function taken identical for both elastic and plastic contributions of energy.This brings certain computational advantages, allowing the retention of the original yield function without any phase field components as the yield surface and equivalent stress degrade identically.A very simple alteration to the brittle strong form becomes adequate to incorporate plastic damage leading to where G c represents the fracture toughness and  0 is the length scale parameter governing the diffusiveness of the crack.H is the local history variable given by the maximum elastic energy over the deformation history.The definition of W p is the dominating factor in the evolution of damage as it quantitatively dominates H. Commonly represented by plastic energy dissipation, it is given as follows: W p = ∫ σ eq εė q p dt (4) where εė q p is the equivalent plastic strain rate and σ eq is the equivalent stress.This particular formulation has its drawbacks.In the case of brittle fracture, the elastic damage accumulates at an accelerating pace due to the linear elastic nature of the problem.Plastic contribution to damage, however, accumulates significantly slower, leading to a gradual degradation over a large strain range rather than a dramatic drop past a critical strain level, which is a far more realistic representation.The addition of a plastic dissipation threshold W c p solves this problem, with significant changes made to the role of existing parameters.G c no longer represents the fracture toughness of the material, with this role passed on to W c p .Furthermore, the length scale is demoted back to a regularization parameter, no longer influencing the material fracture toughness.Instead, G c only influences the rate of degradation past the threshold, controlling the reduction in stress following critical damage, where it may be inflated to ease convergence issues.W p is now replaced by 〈W p − W c p 〉 where 〈. 〉 is the Macaulay operator.An alternative family of degradation functions proposed by Alessi et al. [23] was chosen over the more common quadratic function. g= where k=1 retrieves the original quadratic function.Lower values of k delay premature degradation, though lowering the value has consequences in terms of convergence issues.A value of 0.5 was found to sufficiently reduce premature degradation while presenting no serious convergence problems.Furthermore, the equivalent plastic strain at failure, ε f , for the MMC model (see [21]) is given as follows: where K, n and C � i are material specific constants.L is the Lode parameter and T is the stress triaxiality with both defined as where σ 1 , σ 2 and σ 3 are principal stresses while σ h and σ eq are the hydrostatic and equivalent stresses, respectively.The incorporation of the MMC model opens up the range of specimens and deformation cases that may be studied, where the model can accurately capture damage at evolving stress states.
While in the classic approach, damage evolves through equivalent plastic strain and failure strain, as D= ∫ ε̇e q p ε f dt (8) where D=1 represents complete fracture and ε f is the failure strain.In the context of the phase field paradigm, damage evolves through plastic dissipation that drives the phase field parameter to retain an energy-based framework.We employed the following definition for W p : In the current implementation, the addition of the MMC model helps scale the plastic contribution in regard to the stress state.While classical uncoupled fracture models, employ a strain-based incremental measure of damage where the material is considered fully broken at a damage value of 1, this framework uses W c p as the fracture indicator.J2 plasticity model with nonlinear isotropic hardening is coupled with the phase field paradigm through a staggered scheme, where plastic and elastic damage contributions from the previous increment influence the evolution of the phase field parameter in the current one.An extended Voce hardening rule given below is utilized to define isotropic hardening.
(1-exp(-C i ε eq p )) (10) where σ y is the yield stress, σ 0 is the initial yield stress value, and Q i and C i are hardening parameters.The plasticity framework and the phase field evolution are coupled through the degradation function g(ϕ).J2 plasticity model is solved with classical radial return algorithm in the user subroutine (UMAT) and at the end of each increment stress and material stiffness are degraded through g(ϕ) retrieved from the previous converged increment, as σ=g(ϕ)σ 0, ℂ=g(ϕ)ℂ 0 (11) where σ 0, and ℂ 0 are the undegraded stress tensor and material tangent stiffness tensor, respectively, retrieved from the J2 plasticity solution.The phase field evolution is solved by using the heat generation subroutine (HETVAL) exploiting the similarities between the phase field strong form (Eq. 3) and steady-state heat transfer equation (see [20] for the details of the implementation).Elastic and plastic energy contributions from the previous converged increment are used to drive the phase field evolution in the current increment.

Ductile Failure Simulations
Several benchmark examples are simulated to gauge the proposed model's ability to capture ductile failure crack paths as well as the material response.Firstly, 4 specimen geometries, namely notched tension (NT3, NT10), plane strain tension (PST) and in-plane shear (ISS), are solved using the proposed phase field formulation for 6016-T4 aluminum alloy.The details of the specimen dimensions and experimental data are given in [21].The MMC model variables, calibrated and tested in the aforementioned paper, are used in the ductile fracture phase field model.The parameter set given in Table 1 is employed in the presented phase field framework.
Table 1.Material and model parameters.All simulations are performed with the implicit finite element solver Abaqus.Models are meshed with a fully integrated 3D (hexahedral) temperature-displacement coupled element (C3D8T in Abaqus) as shown in Fig. 1.NT3, NT10 and ISS are pulled from the center of the top pin and held from the center of the bottom pin while allowing for rotation around the pins.Pin centers are connected to the specimen with the MPC beam constraint.For the PST specimen, the clamped regions are assumed to be rigid and modeled as rigid bodies in FE as shown in Fig. 1.The PST specimen is pulled upwards from the top rigid section and held fixed from the bottom rigid section.Only the middle portion of the ISS specimen is modeled for FE simulations to save computational time.In order to apply boundary conditions, top and bottom pin centers are connected to the FE mesh with the MPC beam constraint as shown in Fig. 1 with black lines.The details of the specimen dimensions and experimental procedure can be found in [21].The average element size in the gauge section is 0.05 mm with 5 elements in the thickness direction in all FE models.Solutions are performed using automatic step size control with a maximum step size of 0.001 s and a total time of 1 s which is found to be an efficient value in terms of solution time and convergence.[21] and the present phase field results.Load vs.

Fig. 2. Comparison of experimental
Displacement curves for 4 specimen geometries.
In Fig. 2, the comparison of force vs. displacement is depicted for 4 specimens.For the given MMC parameters, the material response performed well for each of the four specimens compared to experimental data.The fracture strain for the ISS specimen is exaggerated slightly, but referring to [21] confirms that this is a feature of the MMC fit, not the phase field model.It can be concluded that the developed phase field framework is compatible with classical uncoupled ductile failure models such as Johnson-Cook and MMC.The simulated crack paths are shown in Fig. 3.A straight horizontal crack is developed in the tension specimens resulting in a cup-cup fracture surface.For the ISS test, a slanted crack is developed in the gauge region.All crack paths are compatible with the experimental results presented in [21].In order to illustrate the effect of length scale in the current model, two additional simulation results are shown in Fig. 4 with the PST specimen.It is clear that the length scale has an insignificant effect on the failure strain but the diffusiveness of the phase field is increased with increasing length scale.Using a smaller length scale is desirable to have relatively precise crack paths; however, this requires an equivalently finer mesh as comparable scales for both can lead to the mesh distribution influencing the crack path.In all simulations, we set the length scale as 1 mm.Furthermore, in Fig. 4, the effect of critical plastic dissipation, W c p , is depicted.This is the main calibration parameter of the proposed model that defines the ductility limits of the material where larger critical dissipation delays the fracture strain value.
Fig. 4. Effect of length scale and critical plastic dissipation, W c p , on fracture.

Nakazima Forming Test Simulations
In order to demonstrate the capability of the presented framework in predicting material formability, the Nakazima test [22] is simulated to reproduce a failure limit curve which assesses the formability of a material in the context of the deep drawing process.The dimensions of the Nakazima setup are adopted from [24].Specimen thickness is taken as 1mm.The finite element model of the test is given in Fig. 5.Both the punch and the die are taken as rigid bodies and the test sample is meshed with fully integrated 3D (hexahedral) temperature-displacement coupled elements.Outer edge of the Nakazima sample and the die is fixed in all three directions.Displacement boundary condition is applied to the punch while fixing the displacement in transverse directions.Incrementation is retained from the previous section.Frictionless contact is assumed between all interacting surfaces.The punch is pushed vertically until the specimen experiences crack initiation and subsequent crack growth.In Fig. 6, the crack pattern developed for the various specimens may be observed.For narrow specimens, the crack initiates horizontally at the center followed by slanted crack growth which is a pattern that can be observed experimentally (see [24,25]) but most damage models fail to capture it.A horizontal crack development is obtained for wider specimens in agreement with the experimental observations.
The failure limit curve is plotted by taking the major and minor strain histories at the center point of each specimen.The curves are depicted in Fig. 7.A critical facet of plotting a forming limit curve (FLC) is choosing the failure point for the specimen.In this study, we choose the point right before the crack initiation predicted by the phase field model.Hence, we preferred to refer to the curve as a failure limit curve rather than a forming limit curve.The shape of the curve obtained matches the prevalent trend in the literature.A possible source of inaccuracies the current study may be the assumption of isotropic hardening behavior as sheet metals tend to be anisotropic.The study in [25] compared the performance of an anisotropic yield function vs. an isotropic one and concluded that the anisotropic yield function is found to be more successful at reproducing majorminor strain history curves.a route for the non-local coupled implementation of classically uncoupled damage models.Finite element simulation results demonstrated good correspondence with experimental data in terms of both crack paths and material force-displacement response.The Nakazima test is analyzed in FE to show a possible forming application of the developed phase field model, and the results are found to be promising.The framework is a simplistic implementation of a non-local coupled ductile failure model functioning with an implicit solver.The versatility of this methodology could allow future works to simulate more complex ductile fracture behavior such as anisotropic plasticity and damage models.

Fig. 7 .
Fig. 7. Failure loci plot for Nakazima specimens.Summary A ductile phase field fracture model is developed to simulate crack paths and accurate material responses for various benchmark cases.It demonstrates compatibility with the MMC model allowing the use of existing calibrated parameters for various materials.The framework provides