Calibration of the modified Mohr-Coulomb fracture model by use of localization analyses for three tempers of an AA6016 aluminium alloy

This paper presents a novel calibration procedure of the modified Mohr-Coulomb (MMC) fracture model by use of localization analyses and applies it for three tempers of an AA6016 aluminium alloy. The localization analyses employ the imperfection band approach, where metal plasticity is assigned outside the band and porous plasticity is assigned inside the band. Ductile failure is thus assumed to occur when the deformation localizes into a narrow band. The metal plasticity model is calibrated from notch tension tests using inverse finite element modelling. The porous plasticity model is calibrated by use of localization analyses where the deformation histories from finite element simulations of notch and plane-strain tension tests are prescribed as boundary conditions. Subsequently, localization analyses are used to establish the failure locus in stress space for proportional loading conditions and thus to determine the parameters of the MMC fracture model. Finite element simulations of notch tension and in-plane simple shear tests as well as two load cases of the modified Arcan test are used to validate the calibrated fracture model. The predictions by the simulations are in good agreement with the experiments, even though some deviations are seen for each temper. The results demonstrate that localization analyses are a cost-effective and reliable tool for predicting ductile failure, reducing the number of mechanical tests required to calibrate the MMC fracture model compared to the hybrid experimental-numerical approach usually applied.


Introduction
Modelling and simulation of ductile fracture in metallic materials is an active research field where significant progress has been made over the last decades. This research is important since industries like the automotive industry want to utilize the materials to the brink of failure. Thus, the demand for accurate predictions of fracture by numerical simulations is increasing. Reliable design of structural components against ductile fracture requires a robust numerical framework able to accurately describe the damage and fracture properties of the material. In many lightweight metals, which have received special attention by the automotive industry in recent years, strength and ductility are inversely proportional properties. As strength is often favoured in this case, the ductility imposes a great challenge in design of safety components of such materials.
Nucleation, growth and coalescence of microscopic voids at various length scales is known to be the physical mechanism governing ductile failure. Studies have agreed that the stress state affects the ductility of a metallic material [1][2][3] . The influence of the hydrostatic stress state was discovered early and has since been included in several fracture models.
non-proportional loading paths [5] . However, the extension of this type of model to non-proportional loading is found by Benzerga et al. [6] to be in contradiction with micromechanical observations, and therefore questionable.
Fracture models that incorporate stress triaxiality and Lode parameter dependence usually have the disadvantage of a comprehensive calibration scheme, requiring multiple model parameters to be calibrated from different types of mechanical tests. Often a hybrid experimentalnumerical approach is employed, where the stress and strain histories are extracted from the critical element in a finite element simulation of the test. The optimal set of parameters is then found by comparing the curves from the simulations to those of the experiments. The accuracy of this approach relies on experiments that cover a range of stress states and exhibit close to proportional loading paths all the way to fracture. The latter requirement is difficult to fulfil for most stress states. As an alternative to this approach, localization analyses may be used to predict ductile failure. An underlying assumption here is that strain localization is a precursor to failure, and thus may be regarded as the onset of fracture. Mechanical tests are then only needed to calibrate the constitutive equations used in the localization analyses and there is no requirement for proportional loading paths. Morin et al. [7] combined unit cell simulations and localization analyses to predict failure of a steel under non-proportional loading. The numerical results were validated against experimental results reported by Basu and Benzerga [8] and found to be in good agreement. The versatility and effectiveness of the localization analyses were demonstrated by Morin et al. [9] where failure loci of metals were generated from localization analyses and applied to an advanced high-strength steel subjected to proportional loading paths. The results were evaluated against 3D unit cell analyses by Dunand and Mohr [10] and proven to give comparable results in a fraction of the computational time. Gruben et al. [11] applied an experimental-numerical approach to determine the strain localization and ductile fracture of two dual-phase steels. Four tests that covered stress states from simple shear to equi-biaxial tension were conducted. Numerical simulations of the tests were performed, and the failure strains were estimated by comparison to the experimental data. Localization analyses by use of the imperfection band approach were conducted to predict the onset of localization. The results indicated that the localization analyses provided conservative values of the failure strains and that the stress state inside the band tends to move towards a generalized shear state prior to localization. Bergo et al. [12] used unit cell simulations and localization analyses to calibrate failure loci for three different steels. The study was confined to generalized tension stress states and thus only the dependence on stress triaxiality was included. The uncoupled plasticity and fracture models were calibrated based on a single uniaxial tension test and micromechanical simulations using unit cells, metal and porous plasticity and localization theory. The predicted ductility was somewhat conservative for Weldox 460E and non-conservative for Weldox 900E, but accurate for Weldox 700E. It was emphasized that the accuracy of the localization analyses relies heavily on an accurate calibration of the porous plasticity model. However, it was noted that this approach is well suited to reduce the experimental programme required to calibrate fracture models.
Wierzbicki et al. [13] , Li et al. [14] , Gruben et al. [15] and Bai et al. [16] have all presented studies comparing the predictive capabilities of different uncoupled fracture criteria for various steels and aluminium alloys. Several of the criteria evaluated are heuristic extensions of wellknown criteria like the Mohr-Coulomb (MC), Cockcroft-Latham (CL), Rice-Tracey (RT) and Wilkins criteria to name a few. As Wierzbicki et al. [13] pointed out, the quality of a fracture criterion intended for industrial application may be roughly measured by the performance and the cost. Here, performance is defined by the accuracy of simulations compared to experimental tests, while the cost is related to the number of mechanical tests needed to calibrate the model and the complexity related to this. Among the criteria with only one mechanical test required for calibration is the CL criterion. It has been used with success in many studies and is favoured by many for its simplicity and versatility [17,18] . However, it is evident that the criterion has its limitations as it postulates lower ductility in generalized tension than in generalized shear for low and moderate stress triaxiality. The modified Mohr-Coulomb (MMC) fracture model was proposed by Bai and Wierzbicki [19] , where parameters used to control both the dependence of the stress triaxiality and the Lode parameter were introduced. The model provides a monotonic decrease in ductility for increasing stress triaxiality, which is coupled with an asymmetric function of the Lode parameter. The model has been proposed in various versions and is extensively used in the literature. Bai and Wierzbicki [19] calibrated the MMC fracture model for an AA2024 T351 aluminium alloy and a TRIP690 steel, and validated it against various mechanical tests. Accurate predictions of fracture initiation were obtained, but it was noted that the predictions were less accurate in generalized tension. Dunand and Mohr [20] investigated the predictive capability of the MMC fracture model by comparing predictions to fracture experiments on TRIP780 steel. Nine different experiments were used in the comparison covering wide ranges of stress triaxiality and Lode parameter. Fracture initiation was correctly predicted in all simulations of the tests. It was suggested that the underlying physics of the fracture model is of less importance than its mathematical flexibility, implying that even though phenomenological fracture models are motivated by micromechanical observations, their ability to fit experimental data is a superior characteristic. However, a fracture model with high flexibility allows erroneous calibration and requires detailed knowledge by the user.
A modification of the MMC fracture model denoted the Hosford-Coulomb (HC) fracture model was proposed by Mohr and Marcadet [5] , where the von Mises equivalent stress was replaced by the Hosford equivalent stress in combination with the normal stress acting on the plane of maximum shear. The HC fracture model is based on the extensive study on 3D unit cells by Dunand and Mohr [10] , thus it has a micromechanical foundation in contrast to the MMC fracture model. The HC fracture model was presented on locus form, and damage accumulation was taken care of by an integral-based approach. Fracture experiments on three different steels were conducted and three experiments were used to calibrate the fracture criterion. When compared against the experiments, the simulations showed good agreement for the three materials and fracture was accurately predicted in all cases. Gorji and Mohr [21] and Zhang et al. [22] investigated ductile fracture in the aluminium alloy AA6016. In Gorji and Mohr [21] , the HC fracture model was employed in combination with an anisotropic plasticity model to predict shear fracture in deep drawing tests. Eight cup drawing experiments were used in the calibration process to increase the robustness of the fracture model. The results show that the plasticity and fracture models can predict the location and the onset of fracture with good accuracy. In Zhang et al. [22] , an anisotropic Drucker yield function and a fracture criterion proposed by Lou et al. [23] were employed, where a simple shear test and two notch tension tests with different radii were used in the calibration of the fracture criterion. By comparing the experiments to the simulations it was found that the onset of fracture was accurately predicted in tests ranging from simple shear to uniaxial tension.
In the present study, we apply an isotropic plasticity and fracture model to predict ductile fracture in various experiments, where specimens are taken from AA6016 aluminium alloy sheets in three different tempers. The MMC fracture model was selected and calibrated by use of localization analyses based on two mechanical tests. The study is a natural extension to the work by Bergo et al. [12] where a stress triaxiality dependant fracture model was investigated together with an isotropic plasticity model. The aim of the study is to assess the accuracy of the calibrated MMC fracture model by comparison against mechanical tests, where the model's ability to predict fracture initiation and crack propagation is evaluated for a range of stress states. The results demonstrate that the use of localization analyses to calibrate a fracture model has

Materials
Experiments were conducted on three different tempers of the aluminium alloy AA6016. The materials were delivered as 1.5 mm thick plates with in-plane dimensions 625 mm × 625 mm in tempers T4, T6 and T7 by Hydro Aluminium Rolled Products in Bonn. This alloy is mainly used in the automotive industry as outer body panels due to its excellent surface quality, good formability, and high strength. To obtain the various tempers, all plates were first solution heat-treated at 530 °C before being air quenched to reach temper T4. Tempers T6 and T7 were then obtained for some of the plates by artificial ageing for 5 h at 185 °C and for 24 h at 205 °C, respectively. The chemical composition of the alloy is given in Table 1 . The yield strength of the tempers ranges from about 135 MPa for T4 to 245 MPa for T6, and the ultimate tensile strength ranges from roughly 200 MPa for T7 to just below 300 MPa for T6. All mechanical tests were carried out with the longitudinal axis along the rolling direction, unless specified otherwise. The initial thickness of all specimens was measured and found to be similar to the nominal plate thickness of 1.5 mm. An Instron 5985 series universal testing machine was used in all tests, where the force was measured by a 30 kN load cell attached to the actuator. A Prosilica GC2450 camera orientated perpendicular to the in-plane axes of the specimen captured images from all tests. All specimens were spray-painted with a speckle pattern to enable 2D digital image correlation (2D-DIC) by use of the in-house software eCorr [24] .

Uniaxial tension tests
In Granum et al. [25] , uniaxial tension tests in three different directions with respect to the rolling direction (0°, 45°and 90°) of the plates were conducted. Additional uniaxial tension tests in the rolling direction of the plate were conducted only for temper T4 in conjunction with the material test programme presented in this study. This was done to monitor the natural ageing that occurs in temper T4 under prolonged room temperature storage, resulting in solute clustering. This effect is known to slightly strengthen the alloy. The specimen had a gauge length of 70 mm and a width of 12.5 mm, and is depicted in Fig. 1 a). The tests were conducted with a crosshead velocity of 2.1 mm/min, resulting in an initial strain rate of 5 × 10 − 4 s − 1 in the gauge region. A virtual extensometer with an initial length of 60 mm was used to extract displacements by use of 2D-DIC.

Notch tension tests
Notch tension tests in the rolling direction with two different notch radii were conducted, with geometry inspired by the specimens used in Erice et al. [26] . The two specimens are denoted NT10 and NT3, and the geometries are depicted in Fig. 1 b) and Fig. 1 c), respectively. The NT10 specimen had a notch radius of 10 mm, while the NT3 specimen had a notch radius of 3.35 mm. The minimum width of the notch was 5 mm for both geometries. The specimens were tested with a crosshead velocity of 0.6 mm/min. Force measurement from the load cell and images from the camera were recorded at 2 Hz. Two sets of virtual extensometers available in eCorr were used in the post-processing of the experiments, one global and one local. The initial length of the global and local virtual extensometers was 15 mm and 2 mm, respectively, for both test geometries and the virtual extensometers were placed centric to the notch radius.

Plane-strain tension tests
The geometry of the plane-strain tension (PST) specimen is depicted in Fig. 1 d). It had a length of 100 mm and a width of 40 mm. The opening of the notch was 10 mm and the width at the narrowest point inside the notch was measured to 17.33 mm. The force was measured by the load cell and images were taken by the camera at 2 Hz. The tests were conducted with mechanical clamps, where the clamped region on each side of the gauge was approximately 40 mm × 35 mm. Prior to testing, the clamped regions of the specimen were sanded down to ensure good grip between the clamps and the specimen. The tests were conducted with a crosshead velocity of 0.4 mm/min. Two global virtual extensometers with an initial length of 10 mm positioned approximately 16 mm from the centre of the notch were used to extract the displacements in the DIC calculations. The displacements from two virtual extensometers in eCorr were used to ensure that no rotations were enforced during loading. A local virtual extensometer with a gauge length of 2 mm placed centric across the notch was used to obtain local measurements from the tests.

In-plane simple shear tests
The in-plane simple shear (ISS) specimen had a gauge length of 5 mm and the geometry is depicted in Fig. 1 e). The tests were conducted with a crosshead velocity of 0.15 mm/min in an attempt to obtain an initial strain rate in the gauge region of 5 × 10 − 4 s -1 . A virtual extensometer in eCorr spanning across the gauge region with an initial length of 10.05 mm was used to extract displacements. A camera aimed perpendicular to the in-plane surface captured images and force measurements were recorded by the load cell, both at 1 Hz.

Modified Arcan tests
Six modified Arcan specimens were cut from a plate of each temper. The geometry of the specimen is given in Fig. 1 f). The specimen was clamped by four loading brackets using 12 M6-bolts as shown in Fig. 2 . Two different load cases were applied by altering the loading direction ; three with = 90°and three with = 45°as shown in Fig. 2 a) and Fig. 2 b), respectively. All tests were conducted with a crosshead velocity of 1 mm/min, similar to those carried out in [27] .
The tests are abbreviated "Arcan " to distinguish between the two load cases. In the Arcan90 tests, the specimen is loaded in a tension mode where the load axis coincides with the specimen's longitudinal axis. The pin connecting the mounting brackets to the test machine allows the mounting brackets and specimen to rotate when loaded, enabling the specimen to tear open. In the Arcan45 tests, the specimen is subjected to a mixed-mode loading. The width at the narrowest point in the gauge section was measured to 32 mm with negligible variations between the specimens. A virtual extensometer of initial length 10.5 mm spanning across the notch along the longitudinal axis of the specimen in eCorr was used to extract displacements by use of 2D-DIC.

Uniaxial tension tests
Engineering stress-strain curves of representative tests from Granum et al. [25] are plotted in Fig. 3 a). A slight variation in elongation at fracture between the different tensile directions is seen for each of the tempers, while the flow stress showed a maximum deviation of 3%. The Lankford coefficients R were calculated for all tension tests and are listed in Table 2 , where denotes the angle with respect to the rolling direction. All coefficients are below unity and somewhat higher in the rolling direction. The similarity in flow stress between the three directions suggests that the alloy exhibits isotropic properties with respect to strength and strain hardening. However, the Lankford coefficients suggest that the material is prone to thinning and exhibits moderate  These values represent the strain where fracture initiates at the surface of the specimen. The engineering stress-strain curves for three sets of temper T4 tests are presented in Fig. 3 b). The T4-1 curve is a representative test in the rolling direction from Fig. 3 a). T4-2 is a test conducted in conjunction with the modified Arcan tests and T4-3 is a test conducted in conjunction with the plane-strain tension and in-plane simple shear tests. The numbering of the tests indicates the order they were performed in, i.e., the T4-1 test was conducted at an earlier point in time than the T4-2 test, which was performed prior to the T4-3 test. It is seen that natural ageing increases the strength as the T4-2 and T4-3 curves are shifted upwards compared to the T4-1 curve. However, the natural ageing seems to have reached saturation when the Arcan tests were done, as the difference between the T4-2 and T4-3 curves is negligible. This effect was also investigated by Engler et al. [28] for the same material and the reader is referred to this work for a thorough discussion of this phenomenon. Due to the negligible differences between the T4-2 and T4-3 curves, the need for multiple calibrations of the plasticity model for this temper was deemed unnecessary.

Notch tension tests
The experimental results from the NT10 and NT3 tests are shown in Fig. 4 . The repeatability of the notch tension tests was excellent, and thus only one test per configuration is plotted. The displacement was extracted from the global virtual extensometer while the logarithmic strain is calculated from the displacement measured by the local virtual extensometer. The force levels between the two test geometries are similar while the displacements are slightly larger for the largest notch radius. The local strain is also seen to be higher in the tests with the largest notch radius. This is expected as the sharper notch radius confines the gauge section more than the larger notch radius, resulting in higher stress triaxialities within this region.   Fractured specimens from the six different notch tension tests were examined by visual inspection and are shown in Fig. 5 . Only one test per configuration is shown as negligible differences were seen between fracture surfaces of repeated tests. A slant fracture surface was found for all tests, even though some NT3 tests displayed rough shear lips. In general, the fracture surfaces were rougher for temper T7 than for tempers T4 and T6, and shear lips were more prominent in the NT3 tests than in the NT10 tests.

Plane-strain tension tests
The force-displacement and logarithmic strain-displacement curves from representative plane-strain tension tests are shown in Fig. 6 . Duplicate tests are not shown due to the excellent repeatability. The results are in accordance with the trends seen for the notch tension tests, where the most prominent difference being the similarity in elongation at fracture between tempers T6 and T7. The fracture initiated in the centre of the specimen for all tests and propagated in a straight line towards the free edges, perpendicular to the loading direction as seen in Fig. 7 . For the tests in temper T6, the crack propagated instantly resulting in a sudden loss of load-carrying capacity, while the tests in temper T4 exhibited slightly slower crack propagation. Only the T6-3 test experienced complete fracture, where the specimen was pulled apart. In the rest of the tests, the force level dropped below a threshold limit at which the test was stopped automatically. For the tests in temper T7, the crack propagated slowly and it took approximately 40 s from initiation to completion. By inspection of the fractured T6-3 specimen shown in Fig. 7 , a slant fracture surface was observed, where the crack was seen to flip to the other admissible shear band.

In-plane simple shear tests
The force-displacement curves from the shear tests are shown to the left in Fig. 8 where the numbered markers coincide with the numbered images on the right-hand side. Owing to slight scatter, results from all repeat tests are presented. The strain fields obtained by 2D-DIC reveal that strains localized in a thin band across the gauge section in all tests (not shown for brevity). By inspection of the images, fracture seemingly occurred simultaneously within this band, as the origin of fracture initiation was difficult to pinpoint exactly. In-plane rotations were observed in all tests, resulting in an angle of the localized deformation band compared to the loading direction, as evident from the images in Fig. 8 . The angle of the band with respect to the loading direction was similar between repetitions and tempers. The drop in the force-displacement curves, particularly for temper T7, is presumed to occur due to the combined effect of material softening and area reduction in the localized deformation band. The shear tests indicated that temper T7 has superior ductility compared to tempers T4 and T6, and the high ductility makes it difficult to pinpoint the onset of fracture in the temper T7 tests. By inspection of the images, the onset of fracture is anticipated to occur somewhere between point 2 and 3 in the representative forcedisplacement curve for temper T7 in Fig. 8 . All specimens displayed a smooth and flat fracture surface through the thickness, and there were negligible differences among repetitions and tempers.

Modified Arcan tests
The force-displacement curves from the modified Arcan45 and Ar-can90 tests are shown in Fig. 9 a) and Fig. 9 b), respectively. Overall, the trends are in accordance with the other mechanical tests presented showing that temper T6 gives the highest peak force followed in turn by tempers T4 and T7. The peak forces are consistently higher in the Arcan90 tests than in the Arcan45 tests, and the ratio between the peak forces in the two load cases is almost identical amongst the tempers. The displacement at peak force is smaller in the Arcan45 tests than in the Arcan90 tests, but the final displacement is larger in the former. This is linked to the crack paths being longer in the Arcan45 tests than in the Arcan90 tests, especially for tempers T4 and T6. As Fig. 9 b) indicates, the Arcan90-T6 test experienced a sudden loss of load-carrying capacity as the crack propagated instantly across the specimen between two images of the test. Even though temper T7 is the most ductile alloy condition, the tests in temper T4 exhibits the largest displacements. This is linked to the combination of adequate strength, work-hardening and ductility in temper T4, which seems to be more favourable than the high ductility and low work-hardening seen for temper T7 in these tests.   Fractured specimens from the modified Arcan tests are shown in Fig. 10 . In the upper part of the figure, the three different fracture modes observed in the tests are displayed. A curved crack path was observed in both the Arcan45-T4 and the Arcan45-T6 tests. In the Arcan45-T6 tests, the crack was arrested approximately 10 mm from the edge and the tests were stopped as the force level dropped below a lower threshold. In the Arcan45-T4 tests, the crack was arrested approximately 5 mm from the edge when stopped, but during the dismantling of the specimens from the loading brackets, the specimens were pulled apart. By inspection of the fracture surfaces from the Arcan45-T4 tests, all three tests exhibited the alternating slant fracture phenomenon. The crack initiated and propagated in a slant fracture mode until arrested. By inspection of the Arcan45-T6 specimen, the fracture initiated in a V-mode and propagated in this mode for approximately 3 mm before it transitioned to a slant fracture mode. The alternating slant fracture phenomenon was not observed in any of these tests. The V-mode was also seen in Gruben et al. [29] for Arcan45 specimens made of Docol 600DL steel. In the Arcan45-T7 tests, fracture initiated in a V-mode within the same area as for the Arcan45-T4 and Arcan45-T6 tests, and the crack initially followed a similar curved path. However, the crack path deflected abruptly after a few millimetres and propagated perpendicularly to the longitudinal axis of the specimen as shown in Fig. 10 . The fracture surfaces were rough with evident shear lips. All the Arcan90 tests exhibited a similar crack path for all three tempers, but differences in the fracture surfaces were seen. In the Arcan90-T4 tests, the fracture initiated and propagated in a slant fracture mode where the crack was seen to flip abruptly throughout the deformation in all tests. One of the Arcan90-T6 tests initiated in a Vmode before a transition to slant fracture was seen, while the two others initiated in a slant fracture mode. Both the Arcan90-T4 and Arcan90-T6 tests showed a smooth fracture surface, while the Arcan90-T7 tests had a rougher fracture surface, similar to what was seen in the Arcan45-T7 tests.

Stress state parameters
Any admissible stress state can be expressed by the three ordered principal stresses 1 ≥ 2 ≥ 3 given by where 0 ≤ ≤ π 6 is the deviatoric angle, vM = √ 3 2 is the von Mises equivalent stress, and h = I 1 /3 is the hydrostatic stress. Here, J 2 and I 1 are the second principal deviatoric stress invariant and the first principal stress invariant, respectively. The stress state may be conveniently described by the stress triaxiality T and the Lode parameter L . The stress triaxiality is defined as the ratio between the hydrostatic stress h and the von Mises equivalent stress vM , viz. The Lode parameter describes the deviatoric stress state, and is defined either in terms of the deviatoric angle or the ordered principal stresses ( 1 , 2 , 3 ) as The range of the Lode parameter is L ∈ [ − 1, + 1], where L = − 1, 0 and + 1 represent states of generalized tension, generalized shear and generalized compression, respectively.

Plasticity model
The constitutive relation is given by the high-exponent yield surface proposed by Hershey [30] and Hosford [31] , the associated flow rule and an extended Voce hardening rule. Even though the material exhibits moderate plastic anisotropy according to the Lankford coefficient, an isotropic plasticity model is applied. The equivalent stress is given by the principal stresses as where a is a parameter controlling the curvature of the yield surface. This parameter is set to a = 8 in this study as suggested by Hosford [32] based on polycrystal plasticity calculations. The yield function is expressed as where M is the matrix material flow stress, 0 is the yield stress, R is the hardening variable and p is the equivalent plastic strain. The hardening variable is defined by an extended Voce hardening rule on the form where R i are hardening terms that saturate at different levels of plastic strain. The hardening parameters Q i and C i work in pair controlling the strain hardening of the material.

MMC fracture model
Fracture in the simulations is governed by a modified Mohr-Coulomb fracture model. The version of the Mohr-Coulomb model used in this study was inspired by the one proposed by Bai and Wierzbicki [19] . They transformed the original model into locus form where the failure strain ̄ f , i.e., the equivalent plastic strain at failure, was defined in terms of the stress triaxiality T and the Lode angle parameter ̄. The latter is a normalized parameter of the Lode angle and is within the range ̄∈ [ −1 , +1 ] . In this study, the Lode parameter L is used instead of the Lode angle parameter ̄, and the expression for the failure strain ̄ f is then given as [19] where The model has six parameters that must be calibrated: ̂ 1 governs the pressure dependence; ̂ 2 and K control the overall ductility; ̂ 3 determines the Lode parameter dependence; ̂ 4 governs the asymmetry of the failure locus between states of generalized tension and compression; and increasing values of n shift the ductility upwards and decrease the stress triaxiality and Lode parameter dependence [19] .
Damage is introduced by the damage variable D which is set to evolve with increments of the equivalent plastic strain p , given as The material is undamaged in its initial configuration, i.e., D = 0, and fracture initiates when D = 1. Whereas the failure locus is valid for proportional loading only, the damage accumulation rule is assumed to account for non-proportional load paths in an approximate way.

Finite element modelling
The finite element (FE) simulations of the mechanical tests used in the calibration process were conducted with the implicit solver of Abaqus [33] with displacement-controlled loading. All simulations with the MMC fracture model were conducted using the explicit solver of Abaqus with velocity-controlled loading. The specimens were discretized by use of 8-node linear brick elements with selective reduced integration, denoted C3D8 in Abaqus. Fracture is modelled by element erosion, where elements are removed when D in Eq. (9) reaches unity. In the simulations of the NT3, NT10 and PST tests, three symmetry planes were utilized, and the reduced models were optimized with respect to the number of elements. The validity of the reduced models with optimized mesh design was verified against selected simulations of the full specimen with a dense, uniform mesh. The differences in the predicted crack initiation and propagation between simulations with the optimized and uniform mesh designs were negligible. The maximum time step in the implicit simulations was selected so that each simulation had around 200 increments. All simulations were conducted with 5 elements over the half-thickness and an in-plane discretization in the gauge region with a characteristic element size of 0.15 mm. This resulted in initially cubic-shaped elements in the gauge region.
An in-plane view of the FE models with the mesh depicted on the initial geometry of the specimens is shown in Fig. 11 . In all simulations, except for the simulations of the PST tests, the load was assigned to a reference node positioned in the centre of the pinhole. An MPC beam constraint was used to connect the reference node to the element set on the boundary of the specimen, to mimic the effect of a pin pulling the specimen. This is visualized by the red lines and the blue reference nodes in Fig. 11 . This approach limits the number of elements in the FE models significantly and presumes that the omitted part moves as a rigid body. Also, the rotations induced by the pinned link is recognized in this modelling approach by allowing the reference node to rotate inplane. In the simulation of the PST test, the clamped area was assumed to behave as a rigid body and thus omitted in the model. The boundary conditions were thus assigned to the edges bordering to the clamped regions of the model. The ISS and modified Arcan test specimens were modelled according to Fig. 11 d)-e), where only through-thickness symmetry was utilized. In the ISS model, two reference nodes connected to the edges of the specimen were applied to allow for in-plane rotations.
In the modified Arcan model, the part of the specimen gripped by the clamping frame was presumed to move as a rigid body and was thus omitted from the model. In-plane rotations were accounted for by inserting reference nodes coinciding with the centre of the two loading pins shown in Fig. 2 . The reference nodes were connected to the edges of the specimen by an MPC beam constraint. This approach simplifies the model substantially and enables feasible computational times with the desired discretization. In all explicit simulations, the velocity was ramped up over the first 10% of the simulation time, and the energy balance was carefully checked to ensure quasi-static loading conditions.

Localization analyses
The localization analyses were conducted using the imperfection band approach, following the work by Rice [34] . A detailed description of the approach used in this study can be found in Morin et al. [9] , and only a brief overview is given herein. A solid, homogenous body that is homogeneously deformed is considered. Within this body, a thin planar band is assumed to exist where stress and strain rates are allowed to be different from their values outside the band. However, continuing equilibrium and compatibility across the imperfection band are enforced. The normal to the band is denoted n and a sketch of a solid body with a planar band is depicted in Fig. 13 . Localization is set to occur when the strain rate inside the band becomes infinite. The critical orientation of the band is not known on beforehand and localization analyses covering a range of band orientations must be conducted to find the one  exhibiting the lowest ductility. The failure strain is taken as the equivalent plastic strain outside of the imperfection band at localization inside the band.
The solid body is governed by the plasticity model described in Section 4.2 . Inside the band, a porous plasticity model is used to represent the material behaviour, which enables a simple approach to introduce an imperfection by nucleation and growth of voids. The heuristic extension of the Gurson-Tvergaard model [35,36] proposed by Daehli et al. [37] is adopted, where the yield condition is given by where eq is the Hershey-Hosford equivalent stress defined in Eq. (4) , M is the flow stress of the matrix material according to Eq. (5) , q 1 , q 2 , q 3 are the Tvergaard parameters, h is the hydrostatic stress, and f is the void volume fraction. The evolution of void volume fraction is defined as ̇ = ̇ g + ̇ n + ̇ s = ( 1 − ) tr p + n ̇ + s ( ′ ) ′ ∶ p eq (11) where ̇ g is the void growth rate, ̇ n is the void nucleation rate and ̇ s represents the contribution from void softening in shear to the porosity evolution [38] . The parameters A n and k s govern void nucleation and voids softening in shear, respectively. Furthermore, ′ is the deviatoric stress tensor, D p is the plastic rate-of-deformation tensor defined by the associated flow rule and ( ′ ) is a function of the second and third invariant of the deviatoric stress tensor, J 2 and J 3 , respectively, viz.
By including the term for void softening in shear in the evolution equation, the physical meaning of the void volume fraction f is lost and it should be interpreted as a damage parameter, as suggested by Nahshon and Hutchinson [38] . The reader is referred to Daehli et al. [37] for a detailed account of the porous plasticity model used in the localization analyses.

Calibration of hardening parameters
The calibration procedure follows a similar approach as employed in e.g. Mohr and Marcadet [5] . The uniaxial tension tests were used to obtain an initial estimate of the hardening parameters. A spreadsheet was used to fit two of the three hardening terms to the true stress-strain curve up to necking. Inverse modelling of the NT10 tests by use of the   optimization tool LS-OPT [39] was then employed to calibrate the hardening parameters. The initial estimate was used as a starting point in the optimization, where the first hardening term was kept fixed and the second and third hardening terms could change. Sequential FE simulations were then conducted with different hardening parameters where LS-OPT employed a genetic optimization algorithm to find the optimal set of parameters. The force-displacement curves from the NT10 tests were used as targets in the optimizations. The finite element model employed is presented in Section 4.4 and the calibrated hardening parameters are displayed in Table 3 . The force-displacement and logarithmic strain-displacement curves from the simulations are plotted as solid lines together with the experiments as crosses for the two notch tension tests in Fig. 12 . The good agreement between the simulations and experiments for both tests illustrates the validity of the calibrated plasticity model.

Calibration of the MMC fracture model
The predictive capability of the localization analyses relies on an accurate description of the material behaviour inside and outside the imperfection band. The plasticity model used outside the band was calibrated as described in the previous section. For the porous plasticity model used inside the band, the parameters introduced by Tvergaard [36] were given standard values, i.e., q 1 = 1.5, q 2 = 1.0, 3 = 2 1 = 2 . 25 . The nucleation rate A n and the void shearing parameter k s were calibrated based on localization analyses following the process depicted in Fig. 13 , whereas the initial void volume fraction f 0 was set to zero. The deformation gradient F ( t ) was extracted from the critical element in the through-thickness centre of an NT10 and a PST simulation and assigned as boundary conditions in localization analyses. A series of localization analyses with varying nucleation rate A n and void shearing parameter k s was conducted. According to Nahshon and Hutchinson [38] , the void shearing parameter is suggested to be in the range 1 ≤ k s ≤ 3 for structural alloys, thus three values within this range were investigated: k s = {1.0, 2.0, 3.0}. The optimal value of A n was found when localization occurred for a strain outside the band similar to the experimental failure strain for a given k s , as depicted in the plot in Fig. 13 . The optimal value of k s was not searched for and the calibrated value of k s was chosen from the three values investigated. The experimental failure strain was taken from the critical element in a simulation based on both global and local measurements from DIC results. In general, failure was found to initiate just before the final dip in the force levels, as it was assumed that strain localization initiated at this point. The plot in Fig. 13 shows the results from the localization analyses for temper T6 where the opti-  mal value of the nucleation rate A n is depicted. The calibrated values of the parameters in the porous plasticity model are given in Table 4 .
With the metal and porous plasticity models calibrated, localization analyses with proportional loading conditions were conducted, i.e., assigning a deformation where the stress triaxiality T and the Lode parameter L are constant during the entire simulation. From these analyses, failure strains covering a considerable region of the stress space were obtained, even if localization was not achieved for all the applied stress states. The parameters of the MMC fracture model were finally optimized against this cloud of points in a Python script. Approximately 100 points were used in each optimization to ensure a solid basis for the identification. The basin-hopping algorithm [40] available in the SciPy package [41] was employed to determine the optimal set of model parameters. The algorithm aims to find the global minimum of a cost function, here defined as the difference between the failure strains calculated with MMC fracture model and localization analyses for a wide range of stress states. Within the global stepping algorithm, a local minimization is carried out using a sequential least squares programming (SLSQP) method [42] . Approximately 60 basin-hopping iterations were performed to find the optimal set of parameters using the default parameters of the algorithm (the reader is referred to [41] for more details upon these numerical aspects). The calibrated parameters are given in Table 5 .
The calibrated MMC fracture model and the target points calculated by the localization analyses (given the abbreviation SLM) are shown in Fig. 14 for selected values of the stress triaxiality. From the figure it is evident that the main trends are captured in the calibration of the fracture  model, even though the fit around generalized shear ( L = 0) is not always good. The dependence on the stress triaxiality around L = 0 is small according to the localization analyses, resulting in the points for L = ± 0.2 in some cases overlapping each other. This behaviour is not accurately captured by the calibrated fracture model, where an evident dependence on the stress triaxiality is seen around L = 0. The fit is accurate for generalized tension and compression for all three tempers. The somewhat flat failure locus predicted by the localization analyses is amongst other linked to the use of a Hershey yield surface with exponent a = 8. The influence of the yield surface curvature on ductile failure was investigated by Daehli et al. [37] , where Hershey yield surfaces with exponent a = 2 (i.e., equal to the von Mises yield surface) and a = 8 were investigated by use of localization analyses. The results suggest that a yield surface with a = 8 displays a flatter failure locus compared to the yield surface with a = 2. The reader is referred to Daehli et al. [37] for further details on the influence of the yield surface curvature on ductile failure. The results from the calibration of the MMC fracture models are shown in Fig. 15 . The overall shape of the three fracture surfaces in the left column of the figure is similar, where both an evident stress triaxiality and Lode parameter dependence is seen. The monotonic decrease in ductility for increasing stress triaxiality is shown in the middle column of the figure for generalized compression ( L = 1), generalized tension ( L = − 1) and generalized shear ( L = 0). The rate of decrease in ductility for increasing stress triaxiality is similar for generalized compression and generalized tension for all tempers. However, the rate of decrease in ductility for increasing stress triaxiality is much lower for generalized shear, especially for tempers T4 and T6. Generalized shear exhibits the lowest ductility followed by generalized tension and generalized compression, where the difference between the two latter is small. Temper T7 exhibits a much stronger dependence on the stress triaxiality for generalized shear compared to tempers T4 and T6. This is clearly visualized in the right column where the failure loci are plotted as a function of the Lode parameter for selected values of the stress triaxiality. It should be noted that this strong stress triaxiality dependence was not predicted by the localization analyses and is a result of the calibration of the MMC fracture model. The ductility in simple shear is considerably higher for temper T7 than for temper T6, while temper T4 is somewhere in between. The asymmetry of the failure loci is evident for all tempers where the lowest ductility for selected values of the stress triaxiality is found for slightly negative values of the Lode parameter. The high ductility for T = 0 and L = ± 1 for temper T6, higher than for both tempers T4 and T7, is somewhat peculiar. The localization analyses did not provide results for these stress states and this part of the failure locus is obtained by extrapolation. However, for the stress states achievable by experiments, the failure locus for temper T6 exhibits lower ductility than the failure loci for tempers T4 and T7.
The highest and lowest ductility on the plane stress failure locus is found for the stress states representing uniaxial tension ( L = − 1, T = 0.33) and plane-strain tension ( = 0 , = 1∕ √ 3 ), respectively. It is worth noting that the ductility is higher in uniaxial tension than in equibiaxial tension ( L = 1, T = 0.67) for all tempers. This difference would not be possible to express with the Hosford-Coulomb fracture model, where the ductility is forced to be equal in uniaxial and equi-biaxial tension. The cusp on the plane-stress failure locus for uniaxial tension and the valley towards simple shear ( L = 0, T = 0) for all tempers categorizes these material's fracture behaviour as Lode parameter dominated. A stress triaxiality dominated fracture behaviour would exhibit higher ductility in simple shear compared to uniaxial tension, and thus no cusp would appear in the plane-strain failure locus for uniaxial tension.

Material tests
The results from the simulations of the material tests on the notch tension (NT), plane-strain tension (PST) and in-plane shear (ISS) specimens with the MMC fracture model are shown in Figs. [16][17][18][19] . The experimental results are presented as discrete crosses, while the solid lines represent the simulations. By comparison of the response curves for temper T4 shown in Fig. 16 , the predictions by the MMC fracture model are in general found to be good. When comparing the force-displacement curves for the four tests, the agreement for the notch tension tests is excellent, while there are some deviations for the PST and ISS tests. These deviations are expected when modelling a moderately anisotropic material with an isotropic yield surface. The onset of fracture is accurately predicted for the NT10 test, while it is slightly conservative for the NT3 and PST tests. For the ISS test, the response curves deviate already at yielding and the onset of fracture is predicted for a larger displacement than in the experiment. The reason for this deviation is linked to the texture of the alloy, requiring an anisotropic yield criterion to capture the behaviour as discussed in Achani et al. [43] . Engler et al. [28] investigated the microstructure and texture of an AA6016 sheet in temper T4 where a characteristic cube recrystallization texture was found. This texture leads to a relatively high yield stress in shear compared to uniaxial tension [43] . Thus, the yield stress in a shear test is not expected to be accurately predicted with an isotropic yield surface. Considering that the texture is not altered by the heat-treatment, this behaviour is expected for all tempers. The conflicting fracture prediction from the PST and ISS tests illustrates the difficulties of finding a set of parameters that accurately describes the onset of fracture in both tests.
The results from simulations with the MMC fracture model for temper T6 are shown in Fig. 17 . The prediction of fracture in the PST test is slightly conservative, while the fracture predictions for the NT10, NT3 and ISS tests are slightly non-conservative. The impressive accuracy in the predictions of the temper T6 tests is evident as the least accurate prediction is obtained for the NT10 test, which was used in the calibration. As expected, the deviations between the force-displacement curves for the ISS test initiated already at yield. Despite this, accurate prediction on the onset of fracture is also obtained for this test.
The results from the simulations with the MMC fracture model for temper T7 are shown in Fig. 18 . The agreement between the experimental and numerical response curves and the predicted onset of fracture is excellent for the NT10, NT3 and PST tests. The only notable deviation amongst these tests is the shift in the local strain for the PST tests, resulting in a slightly higher strain at the onset of fracture in the simulation compared to the experiment. As mentioned earlier, the exact onset of fracture in the ISS test is difficult to determine based on the forcedisplacement curves. The displacement at which the onset of fracture is predicted in the simulation appears to be a reasonable estimate when inspecting images from the test at a similar displacement.
The predictive capability of the MMC fracture model has been demonstrated in terms of global and local response parameters for four different material tests. In Fig. 19 , the stress state histories extracted from simulations of the same experiments are presented. The solid lines are taken from simulations where the critical damage parameter is set artificially high and the end of the lines indicate the onset of fracture in the experiments. Fracture in the experiments was determined based on both local and global measurements for NT10, NT3 and PST tests, while global measurements were used for ISS. The dots indicate the onset of fracture predicted by the MMC fracture model. The stress states covered by the experiments include mainly negative values of the Lode parameter and positive values of the stress triaxiality. The stress state histories are taken from the through-thickness centre element for the NT10, NT3 and PST tests, which corresponds to the element subjected to the largest equivalent plastic strain. In the simulations of the ISS tests, the element   subjected to the largest value of the damage parameter D at the displacement of fracture in the experiment is taken as the critical element. The damage parameter was used to determine the critical element since the largest value of the equivalent plastic strain was found on the throughthickness surface within the notch. This region is heavily affected by the in-plane rotations that occur in the tests, which makes the element subjected to the largest value of equivalent plastic strain an unsuitable choice for the ISS tests. The chosen critical element in all ISS tests is located on the in-plane surface within the gauge region where strains localize. Among the eight integration points within the critical element, the one subjected to the largest value of equivalent plastic strain is chosen in all tests. When comparing the predictions by the MMC fracture model with the experimental values in Fig. 19 , i.e., comparing the dots to the end point of the solid lines, the trends are similar to the ones in Figs. [16][17][18] . By inspection of Fig. 19 , it is evident that the stress state histories are quite similar among the different tempers apart from in the ISS test. The reason for this is the varying position of the critical element among the simulations of the ISS test. The difference in strength, work-hardening and ductility between the three tempers results in different deformation processes which affect the position of the critical element, emphasizing the difficulties faced with an in-plane simple shear test. Roth and Mohr [44] investigated the challenges related to determining the strain to failure for simple shear for a wide range of sheet metals. Amongst the study's conclusions, it was stated that the shape of the specimen plays a significant role and that different material properties such as strength, work-hardening and ductility require different shapes of the specimen. By inspection of the stress state histories for the ISS tests, temper T6 is closest to exhibit a proportional loading path, suggesting that the geometry of the ISS test specimen is suitable for the material properties of this temper. Ideally, both the stress triaxiality and the Lode parameter should be equal to zero all the way to fracture in a shear test. Especially the ISS test for temper T7 exhibits a loading path where T and L vary markedly throughout the deformation, making it less suitable to use as target in a calibration process. This is one of the reasons why the PST test was chosen in order to calibrate the void shearing parameter k s in the porous plasticity model and not the ISS test. Considering the rather simple plasticity model chosen and that only two material tests were used in the calibration, the predictions by the MMC fracture model are deemed satisfactory.

Modified Arcan tests
To further assess the predictive capabilities of the MMC fracture model, the modified Arcan tests with = 45°and = 90°were simulated. Here, the model's ability to predict both crack initiation and propagation is tested. In Fig. 20 the force-displacement curves from experiments and simulations of the Arcan45 tests are shown to the left, with corresponding crack paths on the undeformed configuration to the right. Despite the small inaccuracies seen in the predictions for the material tests, excellent agreement between the experimental and numerical results is seen for the Arcan45 tests, both in terms of force-displacement curves and crack paths. The onset of fracture is initiated at the correct displacement and position on the specimen for all three tempers. Additionally, the simulated crack propagation occurs mostly along the correct paths at similar velocities as the experimental ones. Even the somewhat surprising straight crack path seen in the Arcan45-T7 test was predicted accurately. The slanted fracture surface observed in the experiments was not predicted in any of the simulations. To predict slanted fracture, the through-thickness symmetry must be abandoned and a much denser mesh accompanied by a coupled damage model is most likely required [45] . However, the crack in the temper T7 simulation propagated in a tunnelling mode from initiation to complete fracture. The stress triaxiality inside the notch where fracture initiated was between 0.3 and 0.4, while the Lode parameter was approximately equal to − 1 for all tempers. Just in front of the propagating crack, a region with stress triaxiality between 0.6 and 0.7 and a Lode parameter close to zero was present for all tempers. Fig. 21 shows the experimental and numerical force-displacement curves for the Arcan90 tests. The onset of fracture was accurately predicted for tempers T4 and T7, while for temper T6 fracture occurred slightly later in the simulation than in the experiment. Additionally, there was a slight deviation in the force level in parts of the response curve before the onset of fracture for temper T6 of unknown reasons. The crack propagation was accurately predicted for tempers T4 and T7, where the agreement between the experimental and numerical response curves was good throughout the whole deformation process. The velocity of the propagating crack for temper T6 was not accurately captured in the simulations, where a lower velocity than in the tests was predicted.
The onset of fracture in the simulation was found to occur a few millimetres within the notch and not at the free surface. This occurred even though the largest value of the equivalent plastic strain was found on the free surface. However, by inspection of the stress state at the onset of fracture, the region inside the notch was found to be subjected to a higher stress triaxiality and a Lode parameter closer to zero than on the free surface. Fracture initiated in this region before propagating perpendicularly to the loading direction. The agreement between the crack pattern in the experiment and simulation was excellent for all tempers and is not shown for brevity.
The strain fields of an Arcan45-T6 and Arcan90-T7 test are shown in Fig. 22 from both DIC and FE simulations. The strain fields were taken when the crack had propagated approximately halfway through the specimen in both tests. The magnitude of the strains is consistently higher in the FE simulations than in the DIC simulations owing to the denser mesh used in the former. However, the qualitative trends are similar between the two sets of simulations for both tests. A narrow zone with localized strains in front of the propagating crack is correctly predicted in both cases. In the Arcan45-T6 test, there is a band with slightly higher strains across the specimen which is not fully developed in the FE simulations and thus only partially predicted.

Conclusions
This paper has presented a novel calibration procedure of the modified Mohr-Coulomb (MMC) fracture model by use of localization analyses of the imperfection band type and applied it for three tempers of the aluminium alloy AA6016. By this approach, the onset of ductile fracture is set to coincide with strain localization. Notch tension (NT10) and plane-strain tension (PST) tests were used in the calibration of the metal and porous plasticity models assigned to the material outside and inside the imperfection band, respectively. The fracture model was validated against notch tension (NT3) and in-plane simple shear (ISS) tests in addition to two loading cases of the modified Arcan test. The goal was to assess the predictive capabilities of the proposed modelling approach, where localization analyses were used as basis to calibrate the MMC fracture model. Finite element simulations with the Hershey-Hosford plasticity model combined with the MMC fracture model were in most cases able to accurately predict the onset of fracture. For both loading cases of the modified Arcan tests, initiation and crack propagation were predicted with good accuracy in all but one simulation. Considering that only two experiments have been used in the calibration of the plasticity and fracture models, the potential of this approach is emphasized. The use of localization analyses to predict ductile fracture is an effective and well-suited tool for design of components and structures of metallic materials, where the need for an extensive test programme could be reduced.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.