## Abstract

A micromechanics-based ductile fracture initiation theory is developed and applied for high-throughput assessment of ductile failure in plane stress. A key concept is that of inhomogeneous yielding such that microscopic failure occurs in bands with the driving force being a combination of band-resolved normal and shear tractions. The new criterion is similar to the phenomenological Mohr–Coulomb model, but the sensitivity of fracture initiation to the third stress invariant constitutes an emergent outcome of the formulation. Salient features of a fracture locus in plane stress are parametrically analyzed. In particular, it is shown that a finite shear ductility cannot be rationalized based on an isotropic theory that proceeds from first principles. Thus, the isotropic formulation is supplemented with an anisotropic model accounting for void rotation and shape change to complete the prediction of a fracture locus and compare with experiments. A wide body of experimental data from the literature is explored, and a simple procedure for calibrating the theory is outlined. Comparisons with experiments are discussed in some detail.

## 1 Introduction

Much has been written about low-triaxiality ductile fracture [1–4]; see Ref. [5] for synthesis and review. The basis of investigation is either deformation-induced anisotropy [6], initial anisotropy [6,7] or an inherent effect of the third stress invariant [8]. That voids change shape, and orientation at low stress triaxiality is observed in physical experiments [9,10] and can be modeled with various degrees of refinement [11]. However, implementing anisotropic constitutive relations is quite complex, as one needs to evolve tensorial internal variables, e.g., Refs. [12–15]. Hence, isotropic constitutive relations and failure models are generally preferred in practice [8,16,17] even if the physical basis is elusive.

The experimental database, from which facts—perceived or real—have been inferred, essentially involves thin-walled specimens. In all, a plane stress state prevails. A feature peculiar to plane stress is that when the biaxiality stress ratio is varied, the stress triaxiality *T*, which measures the hydrostatic stress in units of the deviatoric von Mises stress, and the Lode parameter *L*, a non-dimensional measure of the third stress invariant, do not vary independently from each other. Under such circumstances, when the strain to failure is plotted against *T*, the obtained locus often exhibits a cusp. Phenomenological failure models were proposed to account for this cusp by introducing an explicit dependence upon *L* of either the failure strain, e.g., Refs. [4,8,18], or an internal damage variable, e.g., Refs. [3,17]. However, both approaches fail to provide a rationale for key physical mechanisms.

In a recent work, Torki et al. [19] have provided a rationale for measured failure loci in plane stress. The physical mechanism is void-mediated failure, irrespective of stress state. The key concept is that microscopic failure occurs in bands with the driving force being a combination of normal and shear stresses, both resolved on the bands. Physically, the bands are defined by the spatial distribution of voids. For actual void distributions, which are clustered, a finite number of such bands may be available and the emergent behavior is anisotropic, regardless of the isotropy of prevailing behavior before band formation. For statistically isotropic void distributions, the band orientation may depend solely on the macroscopic stress state assuming small porosity levels. The process of determining the failure band orientation is similar to that for soils and rocks for which the Mohr–Coulomb model is most relevant. In ductile failure, the formation of bands is intimately connected with void coalescence, or better with inhomogeneous yielding, a terminology adopted to highlight the importance of early band formation under shear-dominated loading.

Here, the isotropic mechanism-based ductile failure theory developed in Ref. [19] is employed to reproduce fracture loci in sheet metal using experimental data from the literature. To this end, adjustable parameters are introduced in the porous material constitutive relations similar to the well-known Tvergaard parameters for Gurson-like constitutive relations [20,21]. An essential feature of the theory in Ref. [19], and arguably of any isotropic theory that proceeds from first principles, is prediction of infinite ductility in shear. To rationalize finite strains to failure under shear-dominated regimes, the aforementioned theory is supplemented with a simple anisotropic model [22]. A procedure for calibrating model parameters is outlined. The developed code is made available upon request.

## 2 Model Formulation

### 2.1 Isotropic Constitutive Relations.

The process of yielding in a material containing pores involves one of two modes illustrated in Fig. 1. Homogeneous yielding corresponds to diffuse plasticity and slow (Gurson-like) void growth, Fig. 1(a). In contrast, inhomogeneous yielding occurs when plasticity is localized in intervoid ligaments, Fig. 1(b). It corresponds to accelerated void growth leading to void coalescence. Note that the void growth does not necessarily mean growth of void volume fraction. Under intense shear loading, for instance, one void dimension may grow, relative to some initial size, with no significant volume change [22].

In a continuum description, the intervoid ligaments are represented by bands. When inhomogensous yielding is active, band orientation is solely determined by the macroscopic stress state for statistically isotropic and dilute void dispersions [19].

In Eq. (5), $\sigma eq=32\sigma \u2032:\sigma \u2032$ is the von Mises stress, *σ*_{m} is the mean normal stress, $\sigma \xaf$ is the matrix yield strength, *f* is the void volume fraction, *w* is the void aspect ratio (*w* > 1 for prolate voids), and $C\xaf$, *g*, and *κ* are functions of *f* and *w*; see Ref. [19] for $C\xaf$ and Ref. [23] for others. Also, *q*_{1} and *q*_{2} are Tvergaard parameters [21]. For spherical voids, $C\xaf=1,\kappa =3/2$, and *g* = 0, so that the Gurson–Tvergaard criterion [20,21] is recovered.

*τ*

_{n}and

*σ*

_{n}denote the band-resolved shear and normal tractions, respectively, $\tau \xaf=\sigma \xaf/3$ is the shear yield strength of the matrix,

*f*

_{b}and $w\xaf$ are, respectively, the band porosity and surrogate void aspect ratio, given in terms of

*f*and

*w*by

*K*denotes the complete elliptic integral of the first kind and $\beta (fb,w\xaf)$ is a scalar-valued function given in Appendix A. Also,

*p*

_{1}and

*p*

_{2}are Tvergaard-like parameters.

The formulation underlying Eqs. (3) and (4) is simplified from a multisurface representation [24] taking advantage of the fact that the two mechanisms of yielding are exclusive of each other and that no reversal is possible.

**e**

_{1}–

**e**

_{2}, at an angle

*φ*measured from

**e**

_{1}:

*σ*

_{1}>

*σ*

_{2}the corresponding principal stresses. Equation (10) is implicit since

*τ*

_{n}depends on band orientation $n$ (see Appendix B for a practical method of resolution).

*w*is given by

*α*

_{1}(

*w*) is a function of

*w*alone, while

*α*

_{2}(

*f*,

*w*) [23]. Also, $X\u02d9=Dkkp$ when $F=FH$ and $X\u02d9=Dkkp/c$ when $F=FI$, where $c=fbw$ denotes the ligament volume fraction.

*σ*

_{0},

*ϵ*

_{0}, and

*N*. The effective plastic strain $\u03f5\xaf$ is evolved using the identity:

### 2.2 Failure Criterion.

Under predominately tensile loading, the transition from homogeneous to inhomogeneous yielding corresponds to the onset of void coalescence and, putatively, to failure. When shearing is predominant, however, inhomogeneous yielding may set in much earlier than any definition of failure. It is thus appropriate, within the confines of an isotropic theory, to adopt as a common failure criterion that by which strain localization occurs. For the rate-independent case here, this corresponds to loss of ellipticity of the incremental problem [27,28]. A serious limitation of this approach under intense shear will be addressed in the following section.

*T*and

*L*, it is found that localization invariably occurs after the onset of inhomogeneous yielding within the range of void volume fractions considered. Also, the orientation of the shear band

**q**is, to first order, that of the plastic band of inhomogeneous yielding,

**q**=

**n**. It is sufficient, therefore, to compute the elasto-plastic tangent stiffness tensor $Ct$ based on the constitutive relation for inhomogeneous yielding. In doing so, one finds

**m**is a unit vector along the resolved shear traction on the plastic band with normal

**n**(determined as earlier for internal necking or void-sheet coalescence). Also,

### 2.3 Anisotropic Constitutive Relations.

In Ref. [19], Torki et al. have shown that failure criterion (20) predicts unlimited ductility in shear. This is all that an isotropic theory proceeding from first principles would predict under shear. In other words, any measured finite strain to failure in shear may not be due to an inherent effect of the third stress invariant, but rather to strong induced anisotropy. Constitutive relations that account for this have recently been developed by the authors [22,29]. For completeness, they are succinctly presented here.

The formulation proceeds formally from Eqs. (1)–(3), but anisotropic yield criteria are now used. The homogeneous yield criterion depends on void orientation **v**, $FH(\sigma ;f,w,v)=0$, while the inhomogeneous yield criterion depends on **v** and band orientation **n**, $FI(\sigma ;fb,w\xaf,v,n)=0$. Several criteria are now available in the literature, but evolution equations under inhomogeneous yielding are only found in Ref. [22].

Here, we use for $FH$ the criterion of Keralavarma and Benzerga [23], which for an isotropic matrix essentially reduces to the model in Ref. [30], and associated evolution equations. A description of the constitutive relations and their implementation may be found in Ref. [13].

*σ*

_{n}and

*τ*

_{n}is straightforward, unlike in the isotropic case (see Eq. (10) and Appendix B). Second, the surrogate parameters are defined in accordance with a periodic (hence anisotropic) void distribution. This is in contrast with the statistically isotropic distribution underlying Eqs. (7) and (8). Thus,

**n**and

**m**given. Also,

*C*=

**v**·

**n**and

*S*=

**v**·

**m**.

*λ*(ratio of out-of-band to in-band void spacing) through:

**v**=

**n**, the two would evolve differently under intense shear yet, available inhomogeneous criteria do not account for this; see Ref. [22].

*f*is governed by plastic matrix incompressibility, Eq. (11). That of

*w*follows from:

*c*is now given by

**v**is obtained from

**Ω**

^{v}is the deviation from the continuum spin due to the eigen-rotation of the void, calculated using Eshelby concentration tensors [32] after [13,33], and

**Ω**

^{l}is the shear-induced rotation that comes from mere distortion of void boundaries (dominant here), reading

*c*evolves according to

*λ*and

**n**:

Within this anisotropic framework, failure is taken to mean complete loss of stress carrying capacity, which occurs when *f*_{b} = 1 (or *f*_{b} = 1/*p*_{1} if *p*_{1} is used). This may occur while the overall void volume fraction is small, possibly unchanged from its initial value, which may be orders of magnitude smaller than *f*_{b}.

Thus, when the anisotropic theory is used, the localization-based failure criterion of Sec. 2.2 is not needed. This does not mean that strain localization is not possible, but the competition between failure by loss of load bearing capacity and strain localization is worth a separate study for the anisotropic theory.

## 3 Results

Before any detailed analysis, some anticipated features of the theory are worth emphasizing. In the isotropic version (Sec. 2.1), homogeneous yielding only depends on the mean normal stress and the von Mises stress. Thus, it is independent of the third stress invariant. Conversely, the constitutive relation for inhomogeneous yielding depends on band-resolved stresses (*σ*_{n}, *τ*_{n}), and this leads to dependence upon all stress invariants.

Therefore, it is void coalescence in bands that potentially leads to strong dependence upon both the stress triaxiality *T* and the Lode parameter *L*. This is illustrated in Fig. 2. The locus of strain to failure, *ϵ*_{f}, versus *T* and *L* is shown for an initial void volume fraction *f*_{0} = 10^{−3} and for three values of the initial void aspect ratio. Here and in all subsequent results, a single set of matrix hardening parameters is used with *σ*_{0} = 420 MPa, *ϵ*_{0} = 0.002, and *N* = 0.1. Also, *E* = 210 GPa and *ν* = 0.3 are the elastic constants.

Along proportional plane stress loading paths, *T* and *L* are related and a plane stress fracture locus may be inscribed onto the 3D locus, as illustrated in Fig. 2. A cusp generally appears in plane stress fracture loci although no cusp would appear in any cross section of the 3D locus (that is, at either constant *T* or constant *L*). These general features are similar to those of phenomenological failure models [4,8,18].

In the anisotropic version of the theory (Sec. 2.3), the material response depends on all stress components, but, in a program of loading specified by *T* and *L*, a strain to failure is obtained. It is emphasized in such case that the strain to failure may strongly depend on the orientation of the initial axes (both of the void **v**_{0} and the band **n**_{0}). The main feature of the anisotropic theory is that it *can* predict a finite value of *ϵ*_{f} in simple shear. More generally, the anisotropic theory predicts three possible regimes in shear: (i) infinite ductility without weakening, (ii) infinite ductility with weakening, and (iii) finite ductility with softening; see Ref. [22] for details. To our knowledge, no other theory is capable of such predictions at present.

### 3.1 Parametric Analysis

#### 3.1.1 Isotropic Theory.

There are two types of parameters in the isotropic theory, the initial values of internal variables, *f*_{0} and *w*_{0}, and the heuristic parameters *q*_{1}, *q*_{2}, *p*_{1}, and *p*_{2}. Unless otherwise noted, the default value of each of the latter is unity.

Figure 3 shows three sets of results, each corresponding to a specified value of *w*_{0} for a range of values of *f*_{0}. The larger the initial porosity, the lower the strain to failure across the range of triaxialities. However, the influence of *f*_{0} is much stronger near the cusp. The intensity of the cusp decreases with increasing *f*_{0} and may even disappear, as in the *f*_{0} = 0.05 case for *w*_{0} = 1/5, Fig. 3(a).

It is remarkable that in the biaxial tension regime (1/3 < *T* < 2/3) an order-of-magnitude variation in initial void volume fraction can have a negligible effect on ductility; for example, compare *ϵ*_{f} values for *f*_{0} = 0.005 and *f*_{0} = 0.05 in Fig. 3(a). In biaxial tension, neither triaxiality nor the initial void content seems to have a first-order effect on ductility. There are two exceptions to this. There is a sharp rise of *ϵ*_{f} for equibiaxial tension and, over a wider range, near uniaxial tension. These trends are due to internal necking being dominant under uniaxial tension and equibiaxial tension, and void sheeting prevailing in all other intermediate cases; for details, see Ref. [19].

In all cases, the ductility in shear is predicted to be infinite by the isotropic theory. If there is a cusp in the failure locus, then there is a minimum ductility in the generalized shear regime (0 < *T* < 1/3). This minimum lies closer to a pure shear state (*T* = 0) when the initial void volume fraction is small.

Figure 4 highlights the effect of initial void aspect ratio for *f*_{0} = 0.01. The main effect of *w*_{0} manifests around the cusp. In materials with no initial porosity, as expected in structural alloys, initially oblate voids are generally taken to represent a mechanism of nucleation by particle cracking, while spherical or prolate voids represent void nucleation by decohesion. In that sense, the nucleation mechanism may be viewed to indirectly influence the overall shape of the locus, in particular the emergence of a cusp and its magnitude.

Next, consider the effect of the Tvergaard parameters *q*_{1} and *q*_{2} in the criterion for homogeneous yielding, Eq. (5). This is illustrated in Fig. 5 for spherical voids. Introducing *q*_{1} amounts to multiplying the initial void volume fraction by *q*_{1} so that an effect similar to that of *f*_{0} (see Fig. 3(b)) would be expected. In effect, there are two differences between the influence of *f*_{0} and that of *q*_{1}. Not only does *q*_{1} multiply the effective initial porosity but also the instantaneous rate of growth of *f*, by way of normality to the yield criterion; see Eqs. (11) and (3). The main difference, however, is that the effect of *q*_{1} is confined to a relatively narrow range of triaxiality about the cusp, i.e., on either side of uniaxial loading, Fig. 5(a). Over much of the biaxial tension and generalized shear regimes, *q*_{1} has no effect on *ϵ*_{f}. This is due to the fact that the Lode parameter in these regions is close to *L* = 0 for which inhomogeneous yielding is active from the outset, which is independent of *q*_{1}.

Likewise, the second Tvergaard parameter *q*_{2} has no influence outside the triaxiality range about the cusp, Fig. 5(b). However, the effect of *q*_{2} appears to be more significant, both around the cusp and for equibiaxial tension. This is so because *q*_{2} appears inside the exponential term of the yield criterion and has thus a greater effect on porosity growth whenever homogeneous yielding is active.

Figure 6 shows the influence of the second pair of heuristic parameters, *p*_{1} and *p*_{2}, introduced in the inhomogeneous yield criterion, Eq. (6). Each parameter separately affects the entire fracture locus. When homogeneous yielding occurs first (|*L*| ∼ 1 or *T* about the cusp), failure occurs immediately after the transition to inhomogeneous yielding so that *p*_{1} and *p*_{2} have an effect. When inhomogeneous yielding prevails from the outset (|*L*| ∼ 0 or *T* away from the cusp), their influence is expected. The results in Fig. 6 show that the influence of parameter *p*_{1} is more gradual, which is a desirable feature in model calibration.

It is noteworthy that *p*_{1} and *p*_{2} have both an effect on ductility in biaxial tension. This is important since the Tvergaard parameters *q*_{1} and *q*_{2} had no effect on it, Fig. 5, while the microstructural parameters *f*_{0} and *w*_{0} had a small effect (except for very elongated voids and wide variations of initial void content, Fig. 3(c)).

Neither the microstructural parameters *f*_{0} and *w*_{0} nor the heuristic parameters *q*_{1}, *q*_{2}, *p*_{1}, and *p*_{2} have an effect on shear ductility, which is always predicted to be infinite. This appears as an unavoidable feature of the isotropic theory. The reason for this is that when stress triaxiality vanishes, there is no driving force for volumetric void growth, and void shape changes, if any, tend to make the voids spherical. As discussed in Ref. [22], an essential mechanism of failure in shear involves void rotation and elongation even if no volumetric growth occurs. This mechanism is inherently anisotropic.

#### 3.1.2 Anisotropic Theory.

In the anisotropic theory of Sec. 2.3, the initial state is defined by the triplet of scalars *f*_{0}, *w*_{0}, and *λ*_{0} in addition to directors **v**_{0} and **n**_{0}. Here, our focus is on analyzing the effect of microstructural parameters on the strain to failure in simple shear (*T* = 0 and *L* = 0), hereafter called shear ductility.

There can be more than one band of inhomogeneous yielding in the anisotropic formulation, but here only one band is used parallel to the shear direction such that **n**_{0} = **v**_{0}. This holds initially but upon subsequent deformation **n** ≠ **v**.

Figure 7 summarizes the key trends. The essential finding is twofold: (i) finite values of shear ductility are predicted; and (ii) the shear ductility can be as small as 0.01 and as large as 1.0 for initial void volume fractions in the usual range.

Figure 7(a) shows the effect of initial void volume fraction on shear ductility (i.e., value of *ϵ*_{f} in simple shear). Three sets of results are shown, which correspond to initially oblate, spherical, or prolate voids. In all cases, the initial void spacing is the same within the band or normal to it, so that *λ*_{0} = 1. The effect manifests strongly for small values of *f*_{0} and becomes small for porosities of order 0.01 and above. In the limit *f*_{0} → 0, *ϵ*_{f} → ∞ since no mechanism of failure other than void-mediated coalescence is accounted for.

Figure 7(b) depicts the effect of initial void aspect ratio on shear ductility for four values of *f*_{0} spanning five decades. In fact, what plays a key role in shear failure is not *f*_{0} per se, but the initial ligament ratio *χ*_{0}, defined as the in-plane void size to the in-plane void spacing. Basically, it is related to the band porosity via the relation $fb0=\chi 02$. Thus, the four listed values of *f*_{0} correspond to *χ*_{0} = 0.01, 0.05, 0.1, and 0.46 and initial band porosities of 10^{−4}, 2.5 × 10^{−3}, 0.01, and 0.2116, respectively. Figure 7(b) shows that oblate voids are more damaging with the effect of initial void aspect ratio decreasing gradually for elongated voids.

Clearly, the largest influence on shear ductility is that of the relative void spacing *λ*_{0}, Fig. 7(c). For given values of *f*_{0} and *w*_{0}, the shear ductility varies over a wide range (nominally from zero to infinity) depending on the initial anisotropy of void distribution. This is rooted in the fundamental influence of the relative ligament parameter *χ*_{0}, which is directly affected by void spacing. As illustrated in Fig. 8, the critical configuration for failure is when the rotating void has touched the cell boundaries such that some secondary linkup mechanism leads to total loss of load bearing capacity. When voids are initially more widely spaced along the shear direction than normal to it, it takes a lot more straining to achieve the critical condition; compare the *λ*_{0} < 1 and *λ*_{0} > 1 situations in Fig. 8.

#### 3.1.3 Synthesis.

The isotropic theory has three features: (i) infinite ductility in shear, (ii) cusp near uniaxial tension, and (iii) weakly varying ductility in biaxial tension. The parametric analysis reveals the following:

The initial void volume fraction and to a lesser extent the initial aspect ratio have a strong influence on the “intensity” of the cusp.

Among heuritsic parameters, the second Tvergaard parameter

*q*_{2}has the largest effect on the cusp (Fig. 5(b)).Microstructural parameters

*f*_{0}and*w*_{0}have a weak effect on ductility in biaxial tension; the Tvergaard parameters*q*_{1}and*q*_{2}have no effect on it.Heuristic parameters

*p*_{1}and*p*_{2}entering the inhomogeneous yielding criterion have a strong influence on biaxial tension ductility (Fig. 6).No parameter of the isotropic theory affects the infinite ductility in shear.

By way of contrast, the anisotropic theory has two essential features: (i) finite ductility in shear (unless *f*_{0} → 0) and (ii) infinite ductility under uniaxial loading (not shown). The parametric analysis reveals the following:

Initially oblate voids can be far more damaging than spherical or elongated ones.

The initial relative void spacing has a first-order effect on shear ductility.

These trends will be exploited later when dealing with model calibration on the basis of experimental data (Sec. 3.2.2).

### 3.2 Comparison With Experiments.

A large database on sheet metal fracture is now available. Here, we use data from Refs. [1,4,8] to demonstrate the applicability of the models. Focus is placed on overall trends concerning ductile crack initiation without any attempt to model precisely the plastic flow preceding failure or the nonproportional loading paths that may pertain to certain loading programs. One advantage of the proposed methodology is that it is high throughput in that one can quickly produce a fracture locus for sheet metal calibrated on experimental data of interest. Further refinement of such results is outside the scope of the analysis.

#### 3.2.1 Overview of Experiments.

Figure 9 depicts the specimens used in Refs. [1,4,8]. Specimens 1–5 were used by Bao and Wierzbicki [1]. Later, Bai and Wierzbicki [8] augmented their data set with additional specimens 11–14. Conversely, specimens 6–10 were used by Mohr and Marcadet [4]. The studies in Refs. [1,8] used Al 2024-T351, while in Ref. [4], two dual-phase (DP) steels and one TRIP steel were investigated. A TRIP steel was also used in Ref. [8], but the trends are similar, and so the data are not used here. Not all specimens produce a plane stress state (e.g., 4, 5, 10, and 13). When postnecking deformation is small, as in Al alloys at the ambient, the uniaxial stress state that prevails in specimens 4 and 13 may be viewed as plane stress. In punching (specimen 10), the disk thickness is small so that the deformation is stretch dominated, thus resulting in a quasi-equibiaxial tensile stress state.

The results of experiments are typically produced in terms of a fracture locus giving a measure of fracture strain, $\u03f5f*$, versus a measure of stress state, often a strain-weighted average triaxiality, $T\xaf$. If the location of crack initiation is known, then an estimate of $\u03f5f*$ may be obtained either by direct measurement, e.g., using digital image correlation [35] if crack initiation is superficial, or using a hybrid method. The latter is most common and was indeed used in Refs. [1,4,8]. In the hybrid method, the specimen is meshed and simulated in finite elements using a high-fidelity plasticity model such that the computed global (specimen-level) response matches the measured one. Up to the global displacement where it is *believed* that a crack initiates, local quantities, such as effective plastic strain, $\u03f5\xaf$, triaxiality, and Lode parameter, are extracted from the calculations and used as basis for plotting experimental data with $\u03f5f*$ being identified with the value of $\u03f5\xaf$ at global failure. If some triaxiality variations occur at the presumed location of failure, $T\xaf$ is used for plotting purposes. In Ref. [4], a variant method is used to account for nonporportional loadings. The procedure has limitations but is useful since local quantities of interest are generally inaccessible in experiments.

In initially crack-free specimens, Fig. 9, a damage to cracking transition occurs. An estimate of the fracture strain may be obtained by simply integrating the constitutive relations of Sec. 2 for an elementary volume element under two conditions: (i) the material volume where cracking first occurs is subject to nominally uniform deformation up to crack initiation; (ii) stress triaxiality variations do not lead to strong deviations from the fracture locus under proportional loading. Under such conditions, the application of the proposed theory presents immense advantages. It requires knowledge of initial microstructural parameters *f*_{0} and *w*_{0} as well as (if needed) *λ*_{0}, **v**_{0}, and **n**_{0} and the heuristic parameters *q*_{i} and *p*_{i}.

#### 3.2.2 Parameter Calibration Procedure.

Two methods were used for model identification using experimental data. The first employs an in-house optimizer, while the second follows a simple procedure. Plane stress fracture data can be dispersed in a rather erratic way. In addition, the shape of fracture loci depicted in Sec. 3.1 presents features that are not well suited for optimization (due to sharp gradients near the cusp and the shear limit). In all cases, the optimizer produced results that were at best comparable to the simple procedure. Thus, in what follows, we shall present results of the simple procedure, which is streamlined below:

Measure the initial void volume fraction,

*f*_{0}, or the initial volume fraction of void nucleation sites,*f*_{p}.If neither measurement is available, assume one of two things:

Void nucleation by particle cracking, in which case take

*f*_{0}≤ 10^{−3}and*w*_{0}∼ 0.1–0.2, orVoid nucleation by decohesion, in which case take

*f*_{0}=*f*_{p}and*w*_{0}=*w*_{p}, the particle aspect ratio. Estimates of*f*_{p}and*w*_{p}are available in the literature for most structural alloys.

Use heuristic parameter

*p*_{1}to calibrate the biaxial tension “plateau”; see Fig. 6(a).Use heuristic parameter

*q*_{2}to calibrate the intensity of the cusp; see Fig. 5(b).Use the relative void spacing,

*λ*_{0}, to calibrate the shear ductility using the anisotropic theory; see Fig. 7(c).Obtain the final fracture locus as the minimum of the isotropic locus (which generally works for biaxial tension and the cusp) and the anisotropic locus (which general works near the shear limit).

Refine as needed by small manipulations of the values of

*w*_{0}and*f*_{0}.

The availability of microscopic measurements lowers the uncertainty on some microstructural parameters. For example, *λ*_{0} is proposed above as an adjustable parameter. In principle, it could be inferred from measurements. However, these are tedious [6] and typically not reported by most experimentalists. Conversely, precise knowledge of *f*_{0} and *w*_{0} may not be needed given the small effects they have, say on the biaxial plateau; see Figs. 3 and 4. Unless otherwise noted, *q*_{i} and *p*_{i} parameters are restricted to lie between 1 and 2. The so-estimated fracture strain would constitute a lower bound given that some strain is usually required before nucleation. The reader would also note than no mention was made of plasticity model identification. This will be discussed further in Sec. 4, as it is inconsequential to the main issue at hand.

#### 3.2.3 Fracture Loci.

Experimental data reported in Ref. [1] are replotted in Fig. 10. For ease of reference, the data points are numbered following the scheme of Fig. 9. We emphasize that the “experimental locus” is in fact the result of a hybrid procedure; see Sec. 3.2.1.

Figure 10 shows simulation results for *f*_{0} = 0.001 assuming either nucleation by particle cracking (*w*_{0} = 0.1) or by decohesion from equiaxed particles (*w*_{0} = 1). The biaxial tension regime is not well discriminated by the data.

The only available data point 5 is for a shallow round notch (Fig. 9), which does not truly correspond to a plane stress state. This affects step 3 of the procedure for determining *p*_{1}. The height of the cusp is lower with *w*_{0} = 0.1 using *q*_{2} = 1 and is found to be too high for *w*_{0} = 1 using the maximum value of *q*_{2} = 2. The shear ductility is simulated using *λ*_{0} = 0.38 for *w*_{0} = 0.1 and *λ*_{0} = 0.9 for *w*_{0} = 1. Since the *w*_{0} = 0.1 locus is closer to experimental data, one may hypothesize that nucleation by particle cracking is more likely in this material and that, at least in the severely sheared region of specimen 1, the voids were more widely spaced in the shear direction (*λ*_{0} < 1).

Both fracture loci in Fig. 10 exhibit a double cusp. The first cusp about *T* = 1/3 is inherent to the locus. The second one is the artefact of two intersecting loci. Indeed, the branch of the locus closest to shear (say 0 < *T* < 0.2) corresponds to the anisotropic theory.

Next, data on a 1.4mm sheet of TRIP steel are reproduced in Fig. 11. The data are of better quality in the biaxial tension regime thanks to using notched flat specimens 8 and 9 (Fig. 9), which lead to *p*_{1} = 1.4 for both nucleation scenarios. In this material, the cusp seems to be better captured when assuming *w*_{0} = 1, which suggests more equiaxed nucleation sites. A value of *q*_{2} = 1 was used for the *w*_{0} = 0.1 case. Using a larger value of *q*_{2} would improve the cusp prediction but would lower the prediction of data point 7. Both loci pass through point 6, which is captured thanks to *λ*_{0}.

Data for the two DP steels in Ref. [4] are replotted in Fig. 12 along with our predictions, using as above two nucleation scenarios. For the DP780 steel, the *w*_{0} = 0.1 locus is closer to the experimental one, although the data do not discriminate the cusp. Either nucleation scenario seems reasonable for the DP590 steel.

## 4 Discussion

Until low-triaxiality fracture was investigated in some detail [1], engineering failure predictions were mostly based on void growth models [16,20,36] with a few notable exceptions reviewed in Ref. [37]. Such models do not account for any effects of the third stress invariant, *J*_{3}, nor do they account for anisotropies induced by low-triaxiality loading. The effect of *J*_{3} on void growth has since been analyzed by means of cell model analyses, e.g., Ref. [38]. It was eventually incorporated in porous material plasticity equations from first principles, e.g., see Ref. [39] and references therein. When examined, the effect of *J*_{3} on void growth is modest and by no means can explain the trends seen in sheet metal fracture, particularly in the plane stress limit.

Constitutive formulations were therefore developed to account for inherent *J*_{3} effects (isotropic theory in Ref. [19] and Sec. 2.1) or apparent Lode effects (anisotropic theory in Ref. [22] and Sec. 2.3). Both sets of constitutive relations contain idealizations, as discussed in Refs. [19,22]. However, the analyses in Refs. [19,22] and their application here demonstrate that first-order *J*_{3} effects in ductile fracture are associated with void coalescence, not void growth. In doing so, a connection is naturally drawn to widely adopted phenomenological failure models. Notable among these is the so-called modified Mohr–Coulomb model [4,8]. The connection is made through the concept of failure in microscopic bands such that an effect of (band-resolved) shear and normal tractions emerges.

There are fundamental differences between our formulation and that of phenomenological uncoupled models. The latter suffer from three drawbacks:

They lack physical basis.

They assume a locally hardening, hence invertible, response.

They do not refer to any internal state variable.

Such drawbacks have many consequences, which will be further analyzed in a forthcoming publication. Here, it suffices to draw the reader’s attention that the range of predictions enabled by uncoupled models is within the reach of the present formulation. In addition, the new theory has the immense advantage of connecting to microscopic aspects of damage accumulation to fracture.

The link to microscopic aspects of fracture initiation enables analyzing propagating uncertainties in macroscopic measurements. To illustrate this, consider the case of Al 2024-T351 initially used in Ref. [1], later augmented with additional data in Ref. [8] then revisited by Papasidero et al. [40]. Figure 13 replots data from Ref. [40] where a quantitative discrepancy was noted between experimental data for the same material (Al 2024) in the same heat condition (T351). Given that the same hybrid procedure was used in both Refs. [1,40], uncertainties related to the procedure itself should be minimal. One possible source of discrepancy is therefore microstructural variations. This is illustrated in Fig. 13(a).

The microstructural parameters used for either set in Fig. 13(a) are quite close. It is possible to discriminate the two sets using an order of magnitude variation of *f*_{0} and assuming much elongated particles in Ref. [40] as per the trends in Fig. 3(c). However, such variations are unlikely. Also, as shown in Fig. 13(b), assuming a different nucleation scenario only affects the secondary cusp and not the biaxial tension plateau. Therefore, we conclude that microstructural variations are unlikely to explain the discrepancy among the two datasets in Refs. [1,40] for the biaxial tension regime. We hypothesize that the key difference is rooted in using a different set of specimens among the two studies. Indeed, tubular specimens were used in Ref. [40], whereas sheet specimens were mostly used in Ref. [1], except for the only datapoint used to calibrate the biaxial tension regime (shallow round notch).

For shear failure, our anisotropic theory predicts a strong influence of microstructural variations that have to do with the effect of *λ*_{0}; see Fig. 7(c).

Another issue that is often raised is whether shear ductility is lower than tensile ductility. First, we note that phenomenological models are ill-equipped to address the question. In a Mohr–Coulomb formulation, for instance, it is posited ab initio that failure occurs in shear. The results from our isotropic theory show that a finite ductility in shear must account for induced anisotropies. When these are included, local variations of void spacing parallel to the shear direction can lead to significant scatter, far more significant than what may be found for tension.

The availability of data for *T* ∼ 0.33 seems to justify the cusp predicted by the theory. For example, Bai and Wierzbicki have used specimens 11–16 (see Fig. 9), and this enriches the data about the cusp; Fig. 13(c). Also, the biaxial tension specimen 14 is better than a round notched bar (specimen 5) for estimating the plateau although the state of plane stress induced by the groove is not in the plane of the sheet. Given the potential anisotropy of sheet metal, this induces additional uncertainty in the data.

The calculations carried out here assumed uniform deformation and proportional loading. Several issues encountered in matching experimental data would involve deviations from both as well as the possibility of failure due to shear band formation; see Refs. [41,42]. All such situations would require full boundary value problem solutions. For materials capable of much straining beyond necking, not only does the plane stress assumption break down but also the assumption of proportional loading. Also, when the anisotropic theory is used, it is possible that failure by instability precedes the complete loss of stress carrying capacity implied by void impingement. The choice of the failure criterion was motivated by simplicity and is not expected to affect the overall trends for the composite fracture loci. When instabilities are not critical, one advantage of the methodology proposed here is that it enables a quick assessment of trends over the full stress state regime along with uncertainty analysis due to microstructural variations.

Another aspect that is much emphasized in the literature concerned with phenomenological uncoupled models, e.g. Refs. [4,8,18], is the attention given to high-fidelity simulation of plasticity. There is no particular difficulty in doing the same with the present theory. We leave it to potential users to refine the plasticity part when making direct comparisons. However, we emphasize that a good prediction of plastic flow does not guarantee a good prediction of crack initiation. This issue is illustrated in Appendix C.

## 5 Conclusions

Analyses of plane stress ductile fracture have both practical and fundamental importance. Theoretical considerations based on recent developments of porous material plasticity lead to a rich body of results.

To the extent that the fracture process may be viewed as isotropic, a cusp is generally predicted in the locus giving the fracture strain versus stress triaxiality.

The cusp is rooted in the dependence of void coalescence upon the third stress invariant.

The initial void volume fraction or void shape has minimal effect on ductility in biaxial tension except near uniaxial tension or equi-biaxial tension.

The isotropic theory cannot rationalize finite values of shear ductility.

The anisotropic theory predicts a strong effect of relative void spacing on shear ductility

Heuristic parameters were introduced in the constitutive relations to enable better matching with experimental data.

## Acknowledgment

A. A. B. acknowledges support from the National Science Foundation (Grant No. CMMI-1932975). We thank Vignesh Radhakrishnan for supplying the cell model results of Fig. 14.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. The authors attest that all data for this study are included in the paper.

### Appendix A: Relations for Inhomogeneous Yielding

*β*appearing in Eq. (6) is given by

### Appendix B: Orientation of Void Sheeting

*τ*

_{n}obeys the biquadratic equation:

*τ*

_{n}must be in the range 0 <

*τ*

_{n}< |

*σ*

_{1}−

*σ*

_{2}|/2. Once

*τ*

_{n}is obtained numerically, the band orientation

**n**(i.e., the angle

*φ*) is obtained from Eq. (10).

### Appendix C: Effect of Yield Strength on Ductility

Ductile failure is not brittle fracture. An error on stress prediction does not necessarily translate into meaningful errors on fracture strain prediction. To illustrate this for void-mediated ductile failure, consider single-void cell model simulations. Axisymmetric loading is assumed such that a quarter cell is modeled and the ratio of overall axial to lateral stress is kept fixed. Under these conditions, *L* = −1 and *T* is a constant. The void attributes (volume and shape) are taken to represent average properties of the material (void volume fraction and void aspect ratio). The matrix is taken to obey an elastoplastic behavior, Eqs. (1) and (2), with the plasticity represented by an assoicated von Mises criterion and power-law hardening as in Eq. (13). The calculations were performed using abaqus. For illustration, two simulations were carried out using *T* = 0.66, *f*_{0} = 0.01, *w*_{0} = 1, *λ*_{0} = 1, *N* = 0.1, and *ϵ*_{0} = 0.002. The values of elastic constants *E* and *ν* were as in the text. In one simulation, the initial yield stress was *σ*_{0} = 420 MPa, in the other half this value.

Figure 14(a) shows the overall cell effective stress, *σ*_{eq}, versus effective strain, *E*_{eq}, responses. Even for an initial void volume fraction of 0.01, each response is initially close to that of the void-free matrix. As expected, the cell with a hard matrix exhibits a harder effective response. When the effective stress is normalized with *σ*_{0}, the two curves are barely distinguishable, Fig. 14(b). The evolution of the void volume fraction in Fig. 14(c) is also insensitive to the yield stress. By way of consequence, the strain to failure, which is defined by the transition from homogeneous yielding to inhomogeneous yielding (i.e., the abrupt load drop in Fig. 14) is insensitive to the yield stress.

The implication of this simple illustration is as follows. One stress strain curve in Fig. 14(a) could be thought of as corresponding to a measured response, while the other to a simulated response suffering from a 200% error on stress prediction. Yet, the simulated response predicts exactly the correct ductility. This is due to the latter depending essentially on the triplet {*f*_{0}, *w*_{0}, *λ*_{0}}, not on *σ*_{0}.