Micromechanical modelling of ductile fracture in pipeline steel using a bifurcation-enriched porous plasticity model

In this paper, we investigate the possibility of predicting ductile fracture of pipeline steel by using the Gurson–Tvergaard–Needleman (GTN) model where the onset of void coalescence is determined based on in situ bifurcation analyses. To this end, three variants of the GTN model, one of which includes in situ bifurcation, are calibrated for a pipeline steel grade X65 using uniaxial and notch tension tests. Then plane-strain tension tests and Kahn tear tests of the same material are used for assessment of the credibility of the three models. Explicit finite element simulations are carried out for all tests using the three variants of the GTN model, and the results are compared to the experimental data. The capability of the simulation models to capture onset of fracture and crack propagation in the pipeline steel is evaluated. It is found that the use of in situ bifurcation as a criterion for onset of void coalescence in each element makes the GTN model easier to calibrate with less free parameters, all the while obtaining similar or even better predictions as other widely used formulations of the GTN model over a wide range of different stress states.


Introduction
It is well established that ductile fracture results from the nucleation, growth and coalescence of voids and flaws driven by plastic flow of the matrix material (Anderson 2005). Based on this knowledge, Gurson (1977) established a micromechanical model to represent the yield condition for a porous material with spherical microvoids based on an upper bound analysis. The Gurson model is an extensively used porous plasticity model, where material softening is incorporated due to the evolution of microscopical voids. There is only one microstructural variable in the model associated with the porosity of the material, namely the void volume fraction f , which can be treated as an internal damage variable. Some noteworthy extensions of the Gurson model are the q i -parameters proposed by Tvergaard (1981Tvergaard ( , 1982, void nucleation by Chu and Needleman (1980), void coalescence by Tvergaard and Needleman (1984) and Needleman and Tvergaard (1984), which combined yield the wellknown Gurson-Tvergaard-Needleman (GTN) model.
Large plastic deformation is frequently accompanied by strain localization [see e.g. the work of Cox and Low (1974) and Costin (1979)] into bands of intense straining. Experiments show that such localization of the plastic strain is a frequent precursor to failure (Rice 1976). Once strain localization takes place, large strains typically accumulate inside the band and quickly lead to failure, without substantially affecting the strains in the surrounding material. In numerical models, a softening mechanism must be present in the constitutive equations of the material in order to trigger strain localization (Rudnicki and Rice 1975). It follows that strain localization can be captured by porous plasticity models, like the GTN model, or coupled damage models since such models are able to describe strain softening as a result of damage evolution. The onset of strain localization may occur as either a bifurcation from a homogeneous deformation state or it may be triggered by some assumed initial imperfection, as formulated in a quite general context by Rice (1976). Accordingly, strain localization can be studied using either an imperfection analysis or a bifurcation analysis. Some notable studies adopting the imperfection analysis to investigate strain localization are Needleman and Rice (1978), Hutchinson and Tvergaard (1981), Saje et al. (1982), Mear and Hutchinson (1985), Gruben et al. (2017) and Morin et al. (2018), while the bifurcation approach has been studied by e.g. Besson et al. (2001), Haddag et al. (2009), Chalal andAbed-Meraim (2015) and Erice et al. (2020). A benefit of the bifurcation approach is that it does not require any fitting parameter, such as the initial imperfection size needed in the imperfection band approach. The mathematical formulations of these two approaches are quite analogous, and if the initial imperfection is set to be zero in the imperfection analysis, the problem is reduced to a bifurcation analysis. The two approaches have been reviewed and summarized by e.g. Rice (1976), Needleman and Rice (1978), Yamamoto (1978) and Morin et al. (2018). Doghri and Billardon (1995) adopted localization analysis in simulations of actual structures. In real structures, the ideal situation of an infinite medium in which the band-like discontinuity appears is no longer a valid assumption, so it was proposed to compute Rice's condition for localization during the finite element (FE) calculation. Macro-crack initiation was assumed to occur when the condition for localization was met, and the simulation was then stopped as stability and uniqueness of the solution were no longer ensured. Becker (2002) incorporated a failure criterion in the Gurson model based on Drucker's condition for material stability (Drucker 1957) and a simplified version of the bifurcation analysis by Rice (1976). When instability or localization occurred, the critical void volume fraction at failure was assumed reached and fracture of the element was imposed by setting the stress to zero. Becker used this model in FE simulations of fracture and fragmentation of an expanding ring experiment with good agreement. Besson et al. (2003) discussed the possibility of using the result of a localization analysis to affect the constitutive equations, instead of simply deleting elements as strain localization occurs. They compared simulations with the porous plasticity models by Gurson (1977) and Rousselier (1987), and evaluated the localization criterion proposed by Rice (1976) in the post-processing of the numerical results. By indicating when strain localization occurred on the global response, they showed that Rice's localization criterion tends to slightly underestimate the occurrence of macroscopic failure. As the use of localization as a fracture indicator might be somewhat conservative, they proposed to instead modify the constitutive equations by for instance determining the critical value f C of the void volume fraction for onset of coalescence by the bifurcation analysis.
In a numerical study using finite element unit cell simulations for a wide range of proportional stress states, Tekoglu et al. (2015) found that macroscopic localization into a normal band or a shear band occurs either simultaneously or prior to void coalescence depending on the stress triaxiality T. It was concluded that at low and moderate stress triaxiality, T\1, macroscopic localization occurs simultaneously with void coalescence, while at higher triaxiality macroscopic localization occurs before void coalescence. Thus, the localization phenomenon can be considered in many instances to be a precursor to and a severe warning against initiation of ductile fracture. Furthermore, the importance of the Lode parameter L on ductile fracture has been investigated extensively in the literature both experimentally, e.g. Barsoum and Faleskog (2011), and numerically, e.g. Dunand and Mohr (2014). As shown by Tekoglu et al. (2015) and Morin et al. (2018) amongst others, strain localization is a phenomenon reacting strongly to the Lode parameter L, where states of generalized shear (L ¼ 0) are prone to earlier localization than states of generalized tension (L ¼ À1) and generalized compression (L ¼ þ1). These two findings combined suggest that strain localization could be a powerful indicator for incipient ductile fracture in many instances. Erice et al. (2020) adopted in situ bifurcation as a fracture initiation indicator in FE simulations of a plane-strain material block deformed in various degrees of compression before reverse loading in tension was imposed. The occurrence of strain localization agreed with crease formation observed in the experiments, which triggered brittle fracture during the reverse loading in tension.
Following the proposal by Besson et al. (2003) and the work by Erice et al. (2020), we combine the GTN model with in situ bifurcation analysis in this study, adopting strain localization into a normal band or a shear band as an indicator for incipient coalescence. Thus, f C is never calibrated but determined during the simulation with unique values in each material point. This way, onset of coalescence becomes dependent on the stress state expressed in terms of the stress triaxiality T and the Lode parameter L. The bifurcation-enriched GTN model is compared with two other more standard versions of the GTN model to evaluate the usefulness of the bifurcation analysis. To calibrate and evaluate the credibility of the GTN models, mechanical tests of X65 grade pipeline steel are used. The parameters of the three GTN models are identified from uniaxial and notch tension tests, while their predictive capability with regards to fracture initiation and crack propagation is evaluated by use of planestrain tension and Kahn tear tests. Explicit finite element simulations are carried out for all tests using the three variants of the GTN model and the results are compared to the experimental data.

Porous plasticity model
The GTN model is adopted to describe the porous plastic material using a corotational stress approach under the assumption of small elastic strains. Thus, linear hypoelasticity is applied to describe the elastic behavior of the isotropic material, which is defined by Young's modulus E and Poisson's ratio m.
Following the work of Gurson (1977), Tvergaard (1981Tvergaard ( , 1982, Tvergaard and Needleman (1984) and Needleman and Tvergaard (1984), the yield condition of the GTN model is defined as: where r eq ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3r 0 : r 0 =2 p is the von Mises equivalent stress, r 0 being the deviatoric part of the Cauchy stress tensor r, f Ã is the effective void volume fraction, and q 1 ; q 2 ; q 3 are the parameters introduced by Tvergaard (1981). We also adopt the typical convention that q 3 q 2 1 . The flow stress of the matrix material r M is defined by the extended Voce rule: where r 0 is the initial yield stress, Q i and C i are hardening parameters, and p is the equivalent plastic strain.
The associated flow rule is adopted to define the plastic rate-of-deformation tensor D p , which can be written as: where _ k is the plastic multiplier, determined from the consistency condition. The equivalent plastic strain rate _ p is defined by relating the plastic power of the voided material to the plastic dissipation within the matrix material: so that the equivalent plastic strain is found as p ¼ R t 0 _ pdt À . The evolution of the void volume fraction f is controlled by the growth of existing voids and nucleation of new voids: where A N governs the nucleation of voids. We have adopted here the continuous nucleation rule with a calibration procedure based on tensile tests of smooth and notched cylindrical specimens, as proposed by Zhang et al. (2000). An advantage of the continuous nucleation rule is that it has only one parameter which can be uniquely determined from experiments. An alternative would be to apply the more sophisticated nucleation rule proposed by Chu and Needleman (1980), but at the cost of introducing two additional parameters. As the nucleation rule is not related to the underlying physical mechanisms for void nucleation, the void volume fraction f is here considered to be a damage parameter. The effective void volume fraction f Ã ðf Þ accounts for the effects of rapid void coalescence at failure and is taken to have the form proposed by Tvergaard and Needleman (1984) and Needleman and Tvergaard (1984): where f C is the critical void volume fraction. When the void volume fraction reaches this critical value, void coalescence occurs and the void growth accelerates until f F is reached, which is the void volume fraction at macroscopic failure. Finally, f Ã U ¼ f Ã f F ð Þ is the effective void volume at macroscopic failure, which is typically set to f Ã U ¼ 1=q 1 assuming that q 3 q 2 1 . In the numerical implementation of the material model, element deletion is set to occur as f 0:98f F to avoid numerical problems.
The general conditions for a bifurcation in an elastic-plastic solid, corresponding to the localization of the deformation into a planar band, were derived by Rudnicki and Rice (1975) and Rice (1976). Assuming plastic loading both outside and inside of the band, the condition for bifurcation, or loss of ellipticity, is given by: The acoustic tensor A t is defined as: where C t is the tangent modulus tensor, n is the unit normal to the localization band, and the tensor R is defined by: As suggested by Rudnicki and Rice (1975), we assume the planar band normal to be contained in the plane defined by the maximum and minimum principal stress directions e I and e III , where the indices can be correlated with the principal stresses r I ! r II ! r III . This allows to define the normal with a single angle /, defined between the band normal n and the maximum principal stress direction e I in the plane defined by e I and e III , which can take values from 0 / p=2. To find the band normal with the condition min det A t n ð Þ ð Þ ½ ¼0, a total of 50 bands linearly distributed in between e I and e III is used. To limit the computational cost linked to the bifurcation analyses, the tangent modulus tensor C t and the determinant of the acoustic tensor A t are evaluated just before the beginning of material softening [as defined in Morin et al. (2018)] until bifurcation is reached. Due to the numerical evaluation of the bifurcation condition, the strict definition of min det A t n ð Þ ð Þ ½ ¼0 is relaxed and evaluated as min det A t n ð Þ ð Þ ½ 0. The bifurcation criterion has been implemented together with the GTN model in a user-material subroutine in the ABAQUS non-linear explicit time integration FE solver. The numerical evaluation of the bifurcation condition leads to an average increase of CPU cost by 20% for the simulations presented in this study compared with the standard GTN model.

Material
In the present study, a pipeline steel grade X65 was investigated. The nominal chemical composition of the X65 steel is given in Table 1. The material was delivered as a pipe made from a hot-rolled plate. An extensive experimental program was carried out to establish the stress-strain behavior and ductility of the material under quasi-static loading conditions at room temperature. The specimens were extracted from the side of the pipe opposite to the longitudinal weld according to the standard DIN EN 10002, which is the side of the pipe exposed to the largest deformations during the manufacturing of the pipe from the plate. The longitudinal direction of all specimens coincided with the longitudinal direction of the pipe, which means that the fracture path is always along the circumference of the pipe. Thus, the influence from the anisotropy induced by the pipe rolling on the results should be reduced as much as possible.

Uniaxial and notch tension tests
The geometry and dimensions of the tensile specimens are provided in Fig. 1. The left side of the subfigures in Fig. 1 shows the cross-sections of the specimens, while a side-view of the specimens is shown to the right. The four specimens will henceforth be denoted Uniaxial, R2, R0.8 and R0.2, sorted in order of decreasing notch radius. The side faces of the notch in the R0.8 and R0.2 specimens were for manufacturing reasons inclined with an angle of 17:5 and 45 , respectively. After manufacturing, a control of the resulting geometry of the uniaxial and notched specimens using edge tracing on digital pictures was performed. No significant deviations from the nominal geometry were found for the uniaxial, R2 and R0.2 specimens. In contrast, large deviations from the nominal geometry were discovered for the R0.8 specimens, as illustrated in Fig. 2 where the red contour indicates the nominal geometry. As the intention of the uniaxial and notch tension tests was to calibrate the porous plasticity model, the tests were performed with the received R0.8 specimens, but the finite element model of this specimen was defined entirely based on the measured geometry obtained by edge tracing.
Quasi-static tension tests at room temperature were performed with a 100 kN universal testing machine. The loading speed was constant and equal to 0:15 mm/ min for the uniaxial tension specimens. Since the notched specimens had a shorter gauge length, the loading speed was decreased to 0:10 mm/min in order to obtain approximately the same initial strain rate for the uniaxial and notched specimens. The testing machine registered the applied force and displacement, while a laser extensometer was used to measure the current diameters in the thickness and circumferential directions of the pipe at the smallest crosssectional area of the specimen (Fourmeau et al. 2013). The height of the laser beams was continuously adjusted during the tests to make sure that the minimal diameter was always detected. Two repeat tests were carried out for each test configuration. Due to the anisotropy induced by the pipe rolling, the cross-sectional area A of the specimens became elliptic with deformation and is given by is the geometric mean of the diameters D x and D y in the thickness and circumferential directions of the pipe, respectively, which coincided with the transverse directions of the specimens. The diameter reduction at the minimum crosssection is calculated as where the effect of anisotropy is embedded in D. The true stress r and logarithmic strain e over the minimum cross-section are calculated as where A 0 ¼ pD 2 0 =4 is the initial cross-sectional area, D 0 being the initial minimum diameter of the gauge region. Plastic incompressibility and small elastic strains are assumed in these formulas. It is noted that r and e represent average values over the minimum cross-section for the notched specimens and for the uniaxial tension specimen after the occurrence of necking.

Plane-strain tension test
The geometry of the plane-strain tension specimen is given in Fig. 3. The plane-strain tension tests were performed with the same 100 kN universal testing machine as for the uniaxial and notch tension tests. The speed was set to 0:30 mm/min for the plane-strain tension specimen to obtain approximately the same initial strain rate as in the uniaxial tension tests, since the gauge region was twice as long. The applied force and crosshead displacement were continuously registered by the testing machine. However, as the crosshead displacement is affected by the machine stiffness, digital image correlation (DIC) was used to measure the displacement instead. The in-house DIC software eCorr developed by Fagerholt et al. (2013) was used for this purpose. In order to use DIC, each specimen was painted with a speckle-pattern and pictures were taken by a high-resolution camera with a frequency of 1 Hz. By defining a virtual vector on the specimen surface, the DIC code provides its elongation in a similar way as by using an extensometer. Due to a slightly skewed deformation of the plane-strain tension specimens, the vector elongation was determined by specifying seven vectors of length 22.64 mm placed symmetrically across the whole gauge length perpendicularly to the propagating crack and distributed evenly to span the entire width of the specimen. The average elongation of these seven vectors was applied when reporting the experimental results and when these results were compared with the finite element simulations.

Kahn tear test
The geometry of the Kahn tear test specimen is given in Fig. 4. The Kahn tear tests were performed using a 250 kN universal testing machine. As for the planestrain tension tests, the specimens were painted with a speckle-pattern and DIC was used to record the displacement field, applying the same camera and frame rate as in the plane-strain tension tests. All the Kahn tear tests were run at a speed of 1 mm/min. Force and crosshead displacement were measured continuously, and each of the five Kahn specimens tested were stopped at different levels of deformation. Computer Tomography (CT) scanning of the five partly fractured specimens was performed in a Nikon XT H225 ST MicroCT machine to get a better Tunneling refers to the scenario where the throughthickness crack front grows faster in the center section of the specimen than closer to the surface. This phenomenon is often observed in ductile fracture experiments.
As for the plane-strain tension tests, the crosshead displacement measured during the Kahn tear tests was not used, but rather the elongation of a vector determined by the use of DIC. The vector elongation was found by specifying a vector placed symmetrically across the V-notch in the direction normal to the crack path. The length of the vector was 28.53 mm, and it was placed 2 mm inside of the notch tip. The elongation of this vector was applied when reporting the experimental results and in comparisons with the finite element simulations.

Finite element procedures
Finite element simulations of all the experimental tests were carried out to calibrate the parameters of the GTN model and to evaluate its predictive capabilities.  We use the uniaxial and notch tension tests for parameter calibration and the plane-strain tension and Kahn tear tests for evaluation of the model's credibility. The main aim of the study is to evaluate the use of bifurcation analysis to determine the onset of coalescence in the GTN model. To this end, three variants of the GTN model are applied: • GTN-1-One free parameter: A N . The simplest possible GTN model, included in order to assess predictions obtained using a single parameter description of damage evolution. This is essentially the GTN model without void coalescence, or equivalently the Gurson model with void nucleation. • GTN-2-Two free parameters: A N and f F . In this model, f C is predicted by in situ bifurcation analyses, whereas A N and f F are calibrated from experimental data. It is thus assumed that void coalescence occurs simultaneously to bifurcation. • GTN-3-Three free parameters: A N ; f C and f F .
In this model, all three parameters describing the damage evolution are calibrated from experimental data.
In all three cases, the yield stress and work hardening were calibrated based on the uniaxial tensile test and the Bridgman-Le Roy correction [consisting of the correction equation for round tensile specimens proposed by Bridgman (1952) and the empirical expression for the neck geometry parameter proposed by Le Roy et al. (1981)] to estimate the equivalent stress after necking. Nominal values for steel were used for the elastic constants. The forcediameter reduction curves to fracture from the uniaxial and notch tension tests were used to calibrate the parameters of the GTN model related to void nucleation and coalescence, i.e., A N for GTN-1, A N and f F for GTN-2 and A N , f C and f F for GTN-3. The initial porosity was set to zero in all cases. Genetic optimization was used to determine the actual damage parameters by inverse modelling of the uniaxial and notch tension tests. With all parameters duly calibrated, the credibility of the three variants of the GTN model was evaluated by simulation of crack initiation in the plane-strain tension tests and crack initiation and propagation in the Kahn tear tests.
ABAQUS/Explicit was used as the numerical solver in all finite element (FE) simulations to allow for ductile fracture by element erosion. The loading in the simulations was applied as a prescribed velocity, which was ramped up with a smooth transition function over the first 5% of the step duration. Uniform mass scaling was applied to increase the stable time step in the simulations; thus, it was ensured that the kinetic energy was negligible compared to the internal energy in all simulations to appropriately simulate the quasi-static loading process. All numerical models were discretized using 8-node linear brick volume elements with reduced integration (type C3D8R in ABAQUS).
It is well known that the GTN model shows a strong mesh sensitivity in the post-localization stage, as does any coupled damage model without the use of a regularization technique [see e.g. Tvergaard and Needleman (1995) and Kenik et al. (2012)]. Therefore, the mesh-size effects were investigated to determine a fixed element size suitable for simulation of the uniaxial and notch tension, plane-strain tension and Kahn tear tests. The rather small axisymmetric tension test specimens applied in this study limited the element size, and thus a characteristic element size in the radial direction of 0:068 mm was used in the critical gauge region of all the FE models.

Uniaxial and notch tension tests
The primary purpose of the uniaxial and notch tension tests was to calibrate the parameters of the three variants of the GTN model. With one exception, oneeighth of the tensile specimens was modelled, by utilizing three orthogonal symmetry planes. For the R0.8 specimen, where edge tracing was used to capture the exact geometry of the specimen, only two orthogonal symmetry planes were applied, as the full length of the specimen had to be modelled. The finite element models of the uniaxial and notch tension test specimens are shown in Fig. 5. The mesh was refined in the longitudinal direction of the tensile specimens to avoid excessive element distortion in the neck or notch region upon deformation. The smallest elements in the neck or notch region have initial dimensions of 0:068 mm in the radial direction and approximately 0:015 mm in the longitudinal direction. Upon large plastic deformations, these elements elongated significantly in the longitudinal direction and contracted in the radial direction, resulting in a good aspect ratio inside the neck or notch region at fracture.
As already mentioned, damage parameters of the three variants of the GTN model were determined through inverse calibration using a genetic optimization algorithm in Python. The algorithm proposed by Storn and Price (1997) was applied in this study. Genetic optimization is a metaheuristic method inspired by the process of natural selection, commonly used to generate global optima to optimization problems by relying on bio-inspired operators such as mutation, crossover and selection. The optimization process started from a population of generated individuals with evenly distributed values for A N and f C , where each population in the iterative optimization process is called a generation. The new generation of candidate solutions are then used in the next iteration of the algorithm. The L 1 loss function employed in the optimization is given as: where DD f;i and D b D f;i are the experimental and numerical diameter reductions at incipient fracture of the specimens, respectively, and n s is the number of specimens. Using the genetic algorithm, the loss function L 1 was minimized to find A N for GTN-1, A N and f F for GTN-2 and A N , f C and f F for GTN-3. In the optimization process, the yield strength and work hardening parameters were kept constants.

Plane-strain tension test
In the simulations of the plane-strain tension tests, three orthogonal symmetry planes were applied to model only one-eighth of the specimen. The finite element model of the plane-strain tension specimen is shown in Fig. 6. In the center of the specimen where crack initiation was expected to occur, a circular part with finer mesh was partitioned to enable use of the same characteristic element size as applied for the axisymmetric tensile specimens. The partitioned part  was enforced to deform together with the rest of the specimen using a tie constraint. This modelling strategy was chosen because experimental data after crack initiation were not available. The reason for this was that the crack propagated so fast that no data points could be captured by the camera between crack initiation and a fully ruptured specimen. It follows that the finite element model is valid only until crack initiation, whereas crack propagation would have been affected by the symmetry planes and different mesh size along the crack path. We are thus only addressing crack initiation in the simulations of the plane-strain tension test.

Kahn tear test
In the simulations of the Kahn tear tests, the finite element model was divided into different parts as well, coupled together using a tie constraint. Thus, elements in the vicinity of the crack path were assigned the characteristic size used in the simulations of the uniaxial and notch tension tests, while elements further away only undergoing elastic deformation had larger dimensions. A single symmetry plane with normal vector in the thickness direction was utilized in the finite element model of the Kahn tear test specimen, which is shown in Fig. 7. It follows that slant fracture across the thickness of the specimen could not be captured with this modelling strategy.

Uniaxial and notch tension tests
Simulations of the uniaxial and notch tension tests were carried out to calibrate the damage parameters of the three variants of the GTN model using genetic optimization as described in Sect. 4.1. The force versus diameter reduction curves from the experiments and the simulations are shown in Fig. 8. The curves for both repeated tests are plotted to display the scatter in the experimental results. The markers designate the last registered datapoint from the experiments still with load-carrying capacity, as well as the point of first element deletion in the numerical simulations. For the uniaxial and R0.2 specimens fracture occurs abruptly, and no data points are captured between crack initiation and full fracture.
In contrast, the crack grows more slowly for the R2 and R0.8 specimens, leading to a tail in the response curve with steeper slope. When calibrating the GTN models, the predicted sudden slope change in the response curve was taken as the point of crack initiation for the R2 and R0.8 specimens.  The global response exhibited by the experiments is captured very well by both models incorporating void coalescence, i.e., the GTN-2 and GTN-3 models, which give essentially the same prediction accuracy. The GTN-1 model, on the other hand, overestimates the ductility for the uniaxial and R2 specimens (lower stress triaxiality) and slightly underestimates the ductility for the R0.8 and R0.2 specimens (higher stress triaxiality).
It is noted that the force level in the tests with the R2 and R0.8 specimens is similar, due to the deviations from nominal geometry of the latter specimen. Nevertheless, the ductility is significantly lower for the R0.8 specimen, which is captured by all three variants of the GTN model. It is further noted that the ductility of the R0.8 specimen is even less than the ductility of the R0.2 specimen, which has a much sharper notch and thus a higher force level. Also, this effect is well captured by the various models.
The stress-strain curves from the uniaxial and notch tension tests and the numerical simulations are shown in Fig. 9 in terms of the average true stress r and the average logarithmic strain e over the minimum cross-section area of the specimens. As in Fig. 8, the markers designate the last registered datapoint from the experiments still with load-carrying capacity as well as the point of first element deletion in the numerical simulations. It is seen that the introduction of a notch significantly increases the flow stress due to the triaxial tensile stress field inside the notch region, and as a result, the ductility is reduced. It is interesting to note that due to the machining error, the true stress and logarithmic strain for the R0.8 specimen hardly differs from the R2 specimen, although the ductility is significantly different.
In Fig. 10, the result of the calibration procedure for the GTN-3 model with A N and f C as free parameters is shown for illustration purposes. It was found by trial and error that f F ¼ 2f C gave good results for both the uniaxial and notched specimens and this relation was therefore fixed throughout. The lighter markers converge towards solutions with lower value of the loss function L 1 , i.e., better solutions. The dotted lines correspond to the best individual solution after 13 epochs, which was deemed adequate enough to be used as a converged set of material parameters.
The yield stress and work hardening parameters obtained from the uniaxial tension test after Bridgman-LeRoy correction and the damage parameters for the three variants of the GTN model obtained by genetic optimization using the uniaxial and notch tension test data are reported in Table 2. Recall that in the GTN-2 model, the critical value of the void volume fraction for incipient coalescence, f C , was determined by in situ bifurcation analysis and thus varied with the problem at hand. In contrast to e.g. Nonn and Kalwa (2010), who determined the initial void volume fraction f 0 in the GTN model for the simulation of an X100 pipeline steel based on the volume fraction of large inclusions in the microstructure, we have herein treated the damage parameters in the GTN model in a more empirical manner.   Fig. 11, the simulated void volume fraction in the critical element is shown for the uniaxial and notch tension tests as a function of the diameter reduction. The critical element was always the centermost element. The void volume fraction when void coalescence occurs (f f C ) is marked with a square, while element failure is marked with a diamond. As seen from the figure, the higher nucleation rate A N of the GTN-1 model makes the rapid void growth occur significantly earlier than for the other two models. The high porosity at element deletion for the GTN-1 model (f F ¼ 0:98 Â q À1 1 ) results in relatively large global deformations before the first element is deleted.
For higher triaxialities, the rapid growth of the void volume fraction implies that the exact value of f C hardly affects the point of crack initiation. Considering the evolution of the void volume fraction for the R0.2 specimen, it is apparent that the void growth occurs essentially for a constant diameter reduction after a critical value of f is reached, which is roughly equal to 0:04. It should although be emphasized that the diameter reduction is a global quantity, and the exponential void growth in the critical element is accompanied by a strong increase in the local plastic strain. Based on the global response in Fig. 8 and the void growth in Fig. 11 obtained by the two GTN models incorporating coalescence, namely GTN-2 and GTN-3, it can be concluded that they give essentially overlapping results even though the value of f C is very different. It thus seems like the nucleation rate A N and the porosity at failure f F are significantly more important material parameters for the prediction of ductile fracture in the uniaxial and notch tension tests. Note that the in situ bifurcation analysis gives different values of f C for the four axisymmetric specimens, see Fig. 11.
Another interesting finding is that as the stress triaxiality increases in the critical element when the notch radius gets smaller, bifurcation occurs at significantly lower void volume fraction f as shown by the blue square markers in Fig. 11. For the uniaxial tension specimen, the calibrated value of f C for the GTN-3 model almost perfectly coincides with the value of f C predicted by the in situ bifurcation analysis of the GTN-2 model. But as the notch radius decreases and the stress triaxiality increases in the central element, the differences in f C become significant, although this does not result in a visible change in the predicted global response in Figs. 8 and 9. Significant  It should also be noted how crack initiation in the R0.2 specimen occurs at a larger diameter reduction than the R0.8 specimen, and how this behavior is captured by all three variants of the GTN model. As seen in Fig. 11, the void volume fraction at bifurcation is decreasing significantly as the stress triaxiality increases, occurring at extremely low porosity for the R0.2 specimen, but still this is not the specimen with the lowest ductility. The reason for the higher ductility of the R0.2 specimen than of the R0.8 specimen is that the equivalent plastic strain primarily evolves in the notch root instead of in the center of the specimen. Accordingly, the critical element in the center of the R0.2 specimen has a markedly slower evolution of the void volume fraction than the central element of the R0.8 specimen, and thus the ductility of the R0.2 specimen is higher. Figure 12 displays the equivalent plastic strain in the critical element versus the diameter reduction for all tensile specimens. The rapid increase of the equivalent plastic strain at a critical value of the diameter reduction is evident, especially at higher stress triaxiality, and coincides with the occurrence of bifurcation, or strain localization, in the central element for the GTN-2 model. Whereas the GTN-2 and GTN-3 models predict about the same local plastic strain at fracture, a markedly larger value is predicted by the GTN-1 model as void coalescence is not accounted for. In the uniaxial, R2 and R0.8 specimens, the strain localization in the critical element occurs for decreasing values of the equivalent plastic strain as the notch radius gets smaller and the local stress triaxiality higher. The R0.2 specimen breaks this trend because the plastic deformation mainly occurs in the vicinity of the notch root instead of in the specimen center. Figure 13 shows the equivalent plastic strain field at failure for the four specimen geometries. It is evident that as the stress triaxiality increases with decreasing notch radius, the equivalent plastic strain at fracture decreases significantly, and that the plastic strain is highly concentrated around the notch root in the R0.2 specimen. Even though the stress triaxiality in the center of the specimen increases as the notch gets sharper, the stress triaxiality closer to the surface of the specimen is hardly affected by the notch geometry and attains values close to T ¼ 2=3. The low stress triaxiality close to the surface explains why the significant plastic strain in the vicinity of the notch root of the R0.2 specimen is not enough to trigger crack initiation. Figure 14 compares the experimental and simulated force versus vector elongation curves for the planestrain tension tests plotted until fracture initiation. The vector elongation is defined as L=L 0 À 1, where L 0 is the initial vector length, and L the current vector length. All variants of the GTN model predicts the force level with good accuracy, and the deviation between the experimental and numerical response curves is about the same as the deviation between the Fig. 11 Simulated void volume fraction versus diameter reduction curves for the uniaxial and notch tensile tests, where void coalescence or in situ bifurcation is marked with a square and crack initiation by a diamond Fig. 12 Equivalent plastic strain versus diameter reduction curves from simulations of the uniaxial and notch tension tests, where void coalescence or in situ bifurcation is marked with a square and crack initiation by a diamond three repeated tests that were deformed to fracture. The average vector elongation at fracture initiation in the experiments is 0.068, whereas the numerical values are 0.088, 0.066 and 0.072 for the GTN-1, GTN-2 and GTN-3 models, respectively. The fracture initiated in the center of the specimen in a normal mode which justifies the use the symmetry planes when only the initiation of fracture is targeted.

Plane-strain tension test
The ductility of the material in plane-strain tension is adequately predicted by the GTN-2 and GTN-3 models that incorporate accelerated void growth at coalescence. The difference in ductility predicted with these two models is about the same as the experimental scatter. In contrast, the GTN-1 model markedly overestimates the experimentally obtained ductility. It is interesting that even though the GTN models were calibrated based on force versus diameter reduction curves from axisymmetric specimens (i.e., with Lode parameter L about equal to À1 in the critical element), both the GTN-2 and GTN-3 models were able to accurately predict crack initiation in the flat planestrain tension specimen that is deformed predominantly in generalized shear (i.e., with Lode parameter L about equal to 0 in the critical element). Even if the GTN-2 model has lower nucleation rate A N compared to the two other GTN models, it predicts failure at a lower elongation. The reason for this result is that bifurcation occurs already at f C % 0:02 owing to the shear-dominated deformation mode.
Experimental and simulated equivalent plastic strain fields on the specimen surface are compared in Fig. 15, where the simulation results are obtained by the GTN-2 model. The DIC image shown is the last image before specimen fracture in the experiment, whereas in the simulation with the GTN-2 model the strain field is plotted at the instance of first element failure. While the experimental and simulated fields are qualitatively, and to some extent quantitatively, similar, the plastic strain is larger in the center of the specimen in the simulation than in the experiment. The grey area in the central part of the predicted strain field shows that the strain exceeds the maximum level of the scale bar. The reason for this is most likely that the resolution of the mesh used in the DIC analysis is limited by the size of the speckle pattern and the resolution of the pictures, as well as the imaging frequency. The mesh size in the DIC analysis was chosen as small as possible while still giving consistent results, but it is still significantly coarser (6.4 times coarser) than the FE mesh in the critical center region of the plane strain specimen.
In the experiments the strain field was always slightly non-symmetric, as seen in Fig. 15, while in the simulation it is symmetric due to the use of a perfect specimen geometry and the three orthogonal  Possible explanation for the slightly non-symmetrical strain field could be a small offset in the test setup, small imperfections in the specimen geometry or spatial variations in the material characteristics due to the pipe rolling. However, except for the lower maximum strain and the slightly nonsymmetric deformation mode in the experiments, the simulation gives results in good agreement with the experiments. Figure 16 depicts how the void volume fraction in the critical element evolves as a function of the vector elongation in the simulations of the plane-strain tension test with the three GTN models. It should be noted how early bifurcation (or strain localization) occurs in the GTN-2 model under predominant shear deformation, as marked by the blue square. Even though strain localization occurs at a significantly lower value of the void volume fraction than the calibrated f C used for the GTN-3 model, crack initiation as predicted by the GTN-2 and GTN-3 models occurs at only slightly different global deformations. It is also interesting to note how the GTN-2 model has the lowest nucleation rate but ultimately fails first out of the three GTN models, due to the Lode parameter sensitivity of the bifurcation criterion. Thus, the in situ bifurcation analysis incorporates Lode parameter dependence in the GTN model without introducing additional free parameters and evidently predicts results in good agreement with the experiments. Also note the rather significant additional elongation of the specimen from bifurcation occurs (blue square) until material point failure occurs (blue diamond) in Fig. 16. This indicates that strain localization as a failure criterion (i.e., element erosion occurring at bifurcation) would be overly conservative for shear-dominated loading conditions. Figure 17 shows how the equivalent plastic strain evolves as a function of the vector elongation as predicted by the three GTN models. The rather excessive equivalent plastic strain before element erosion occurs should be noted particularly for the GTN-1 model without accelerated void growth at coalescence. At higher stress triaxialities, the equivalent plastic strain grows rapidly at bifurcation, as illustrated by the uniaxial and notch tension tests in Fig. 11. In Fig. 17, the rather gradual growth of the Fig. 15 Comparison of the equivalent plastic strain field as measured by DIC in the experiment (left) and predicted (right) by the GTN-2 model immediately before crack initiation. The grey area in the central part of the predicted strain field shows elements with strains that exceed the maximum level of the scale bar Fig. 16 Void volume fraction f versus vector elongation from simulations with the three GTN models, where void coalescence or in situ bifurcation is marked with a square and crack initiation by a diamond equivalent plastic strain predicted by all GTN models is a result of the moderate stress triaxiality in the plane-strain tension test, which explains the significant amount of global deformation needed from the occurrence of bifurcation until the critical element fails. Figure 18 illustrates the evolution of the stress triaxiality T and the Lode parameter L as a function of the equivalent plastic strain in the critical element. The stress triaxiality starts at about 0.577 and increases to about 1.5 at crack initiation, whereas the Lode parameter starts out with a value close to zero but decreases rapidly towards À1 as bifurcation takes place. Thus, the plane-strain tension specimen deforms under generalized shear until strain localization and then the deformation state changes rapidly towards generalized tension. This type of behavior has also been reported by e.g. Morin et al. (2018). The critical element in the simulation with the GTN-2 and GTN-3 models fails very close to L ¼ À1, while in the simulation with the GTN-1 model L starts to drift back towards a generalized shear stress state before the element ultimately fails.

Kahn tear test
The experimental and simulated global response of the Kahn tear test specimen is shown in Fig. 19. The diamond markers in the plots denote where each of the five Kahn specimens were stopped, in order to study the tunneling effect and the fracture mode at different global deformations using CT scans.
It is evident that all three variants of the GTN model can capture the global response throughout the Kahn tearing test, and there are only some minor deviations in the predictions and experimental results towards the end of the test. Even if all three GTN models give satisfactory results, the GTN-2 model with in situ bifurcation analysis gives the most accurate results. The mean squared error (MSE) in the force over the entire range of the vector elongation of the most deformed test specimen for the GTN-1, GTN-2 and GTN-3 models are 0.815, 0.158 and 0.539 kN 2 , respectively. The results shown in Fig. 19 illustrate how the GTN-1 model predicts the highest force and thus highest ductility, which is due to the lack of a coalescence criterion. The GTN-2 model predicts the lowest ductility, which is as expected from the results obtained for the plane-strain tension test. A lower force level at the same vector elongation indicates that the crack tip has propagated the furthest, due to faster Fig. 17 Equivalent plastic strain versus vector elongation curves from the FE simulations with the GTN models, where void coalescence or in situ bifurcation is marked with a square and crack initiation by a diamond Fig. 18 Stress triaxiality T (left) and Lode parameter L (right) plotted against the equivalent plastic strain in the critical element as predicted by the GTN models, where void coalescence or in situ bifurcation is marked with a square and crack initiation by a diamond element erosion and thus lower ductility. However, all three models capture both crack initiation and crack propagation with satisfactory accuracy. Figure 20 compares strain fields obtained from an experiment by means of DIC analysis with predictions by the GTN-2 model at the three vector elongations marked with a circle in Fig. 19. The first row shows the effective strain field on the surface of the specimen obtained by the DIC analysis, while the second and third rows present the equivalent plastic strain field predicted at the surface and center of the specimen by the GTN-2 model. It is evident that the measured and predicted plastic strain fields are qualitatively similar but predicted values on the surface are in general larger. The gray region close to the crack path in the simulation has values exceeding the maximum value calculated in the DIC analysis, which is the maximum value of the common scale bar. As discussed for the plane-strain tension specimen, this is primarily the result of the significantly coarser mesh applied in the DIC analysis. The simulations further seem to predict a slightly larger compression zone compared to what was calculated in the DIC analysis.
The equivalent plastic strain field in the center of the Kahn tear specimen is included in the third row in order to display the tunneling phenomenon, i.e., that the crack propagates faster in the interior of the specimen than at the surface owing to the higher stress triaxiality. The edge of each level of deleted elements during the tunneling phenomenon is shown in Fig. 20, making it possible to evaluate the tunneling depth. By comparing with the CT scans of the specimens deformed to different degrees, it was found that the tunneling phenomenon was accurately predicted throughout the entire deformation history. The CT Fig. 19 Experimental and simulated force versus vector elongation curves for the Kahn tear test Fig. 20 Comparison of plastic strain fields as calculated on the surface in the DIC analysis (upper row) and predicted by the GTN-2 model on the surface (middle row) and center (lower row) of the specimen scans show that the tunneling depth, i.e., the difference between the crack length in the center and at the surface, decreases as the specimen deforms, which was also captured in all the numerical simulations (as is visible in the third row of Fig. 20).
Comparisons of the predicted and experimental fracture modes along the crack path are shown in Fig. 21. In the FE model, through-thickness symmetry was applied to reduce the computational time, and the only possible fracture mode is the symmetric mode illustrated to the left in Fig. 21. An FE model without through-thickness symmetry was also simulated, and the fracture mode was unchanged. The three different fracture modes observed in the experiments are shown to the right in Fig. 21. In the five Kahn tear tests, the crack propagation switched between these three failure modes, without any discovered pattern in which one was preferred at a given stage of the experiment. The only trend found was that the most common failure mode was the one shown in the center, i.e., the cup-and-cone fracture mode, occurring slightly more often than the other two observed modes. The third mode observed in the experiments, i.e., the cup-and-cup fracture mode, was the same as the one predicted in the FE simulations. This failure mode is governed by void growth and eventually coalescence in the center of the specimen first, before the crack grows outward to the surface, and the result is a doubly symmetric mode. Figure 22 presents the void volume fraction in every element that fails along the center of the specimen in the FE simulation using the GTN-2 model with in situ bifurcation. The step time is just the time normalized to be unity at the end of the simulation. In the implementation of the GTN models, elements are deleted when the void volume fraction f reaches 0:98f F , as illustrated by the dashed line in Fig. 22. It is apparent that the elements that fail can be split into two different groups with distinctly different behavior, namely elements that fail during crack initiation and elements that fail during crack propagation.
During crack initiation, strain localization occurs at very high void volume fractions. The first element to fail does not experience strain localization before reaching a void volume fraction of 0:98f F . The fact that the critical element does not experience bifurcation before failure is likely the reason why all the three GTN models predict similar global response up to crack initiation, see Fig. 19. The crack initiates in an element slightly inside the notch surface, in the center of the specimen. The crack then propagates backwards towards the notch surface, as well as inwards resulting in the tunneling effect. What can be described as a steady state crack propagation phase occurs as fracture  propagates from the core of the specimen and reaches the surface of the specimen, which is at a step time approximately equal to 0.1. Figure 22 illustrates this phenomenon, at which strain localization occurs at remarkable low void volume fractions for the rest of the simulation. The void volume fraction as a function of step time in every element that fails along the surface of the specimen is shown in Fig. 23 for the GTN-2 model with in situ bifurcation. Figure 23 illustrates that the first element to erode on the surface of the specimen occurs as discussed at a step time approximately equal to 0.1, corresponding to the beginning of what we herein will describe as the steady state crack propagation phase. All elements along the surface experiences bifurcation at very low void volume fractions, and there is no marked difference between the crack initiation and propagation phases, even if the three first elements to experience failure have a higher void volume fraction at bifurcation than the succeeding ones. However, the trend is that the void volume fraction at bifurcation increases steadily for the elements failing along the surface as the crack progresses (in contrast to the elements along the center).
The evolution of stress triaxiality and Lode parameter in the critical element where crack initiation occurs is shown in Fig. 24. As mentioned, bifurcation does not occur before the critical element fails. At crack initiation, both the stress triaxiality T and the Lode parameter L are monotonically decreasing.
For comparison, Fig. 25 illustrates the evolution of stress triaxiality and Lode parameter for all elements along the center of the specimen during what could be described as the steady-state crack propagation phase, which corresponds to all elements along the center of the specimen experiencing fracture after a step time equal to 0.1. The lines plotted with black color correspond to the end of the simulation, while the colored curves (red or blue) correspond to the beginning, and as the crack propagates, the colors shift gradually from red or blue to black. The stress state in the first * 40 elements that fail, i.e., those not included in Fig. 25, drifts from the stress state experienced at crack initiation, as illustrated in Fig. 24, towards the typical stress state experienced during the steady-state crack propagation, as shown in Fig. 25. These elements are not included to illustrate the consistency of the behavior in the steady-state crack propagation phase. By comparison of Figs. 24 and 25, the significant difference in the evolution of the stress state with plastic straining between crack initiation and the steady-state crack propagation phase is evident, and further the significantly larger plastic strain at fracture at crack initiation is noted.
The elements on the opposite side of the specimen relative to the notch where the crack initiates are subjected to compression (L ¼ þ1; T % À0:46) before the stress state changes into plane strain (L ¼ 0) and further into generalized tension (L ¼ À1) after strain localization occurs. The initial evolution of the stress state in these elements is represented by the horizontal lines at the bottom of the plots in Fig. 25. It is interesting to note that the elements first deformed plastically in compression, are the ones for which bifurcation occurs at the lowest equivalent plastic strain, due to an increase in the stress triaxiality. Figure 22 illustrated that bifurcation occurred at increasingly lower void volume fractions as the crack propagated, which fits well with the slight increase in stress triaxiality shown in Fig. 25.
Another interesting observation is the significant change in the stress state that occurs immediately after the bifurcation criterion is fulfilled, which is readily seen by considering the evolution of the Lode parameter. Bifurcation is preceded by a sudden shift towards generalized tension as well as an increase in stress triaxiality. Up until bifurcation occurs, the elements are deformed under generalized shear and a stress triaxiality comparable to that of the plane-strain tension specimens. In the present work, in situ bifurcation was incorporated into the GTN model as a void coalescence criterion. Comparisons were made between the bifurcation-enriched GTN model and two standard variants of the GTN model, as well as with extensive experimental results. The material parameters were calibrated based on uniaxial and notch tension tests, while plane-strain tension tests and Kahn tear tests were used for evaluation of the credibility of the three GTN models.
For the uniaxial and notch tension tests used for calibration, it was demonstrated that as long as void coalescence was incorporated in the GTN model, accurate predictions of crack initiation could be obtained. For the plane-strain tension test, the GTN-2 model, adopting bifurcation as the coalescence criterion, gave a slightly conservative estimate of crack initiation, while the GTN-3 model with a fixed critical void volume fraction for coalescence, resulted Stress triaxiality T and Lode parameter L as a function of the equivalent plastic strain in the center of the specimen throughout the steady-state crack propagation phase. The colored curves (red or blue) correspond to the beginning of the phase, and as the crack propagates, the colors shift gradually from red or blue to black in a slightly non-conservative result. The difference is due to the bifurcation criterion being dependent on both the stress triaxiality and the Lode parameter. For the Kahn tear test, all three GTN models gave satisfactory results for crack initiation and propagation, but the GTN-2 model gave results in closest agreement with the experiments.
By employing in situ bifurcation, the idea of letting f C depend on the stress triaxiality as proposed by e.g. Steglich and Brocks (1998) is taken one step further, incorporating a dependency on both the stress triaxiality and the Lode parameter. It has been shown herein that by using in situ bifurcation to determine f C uniquely in each element, we are able to capture the effect of significant variations in both stress triaxiality and Lode parameter. At the same time the need for calibration of f C from experimental data is removed, while obtaining similar or even better predictions than other widely used formulations of the GTN model over a wide range of different stress states. The results thus indicate that the benefit of incorporating in situ bifurcation analysis in combination with the GTN model is twofold: less free parameters in need of calibration as well as possibly better prediction of crack initiation and propagation (provided strain localization occurs prior to ductile fracture).