Ductile Fracture in Plane Stress

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.

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), $σeq=32σ′:σ′$ is the von Mises stress, σ_m is the mean normal stress, $σ¯$ is the matrix yield strength, f is the void volume fraction, w is the void aspect ratio (w > 1 for prolate voids), and $C¯$⁠, g, and κ are functions of f and w; see Ref. [19] for $C¯$ and Ref. [23] for others. Also, q₁ and q₂ are Tvergaard parameters [21]. For spherical voids, $C¯=1,κ=3/2$⁠, and g = 0, so that the Gurson–Tvergaard criterion [20,21] is recovered.

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.

What is of particular importance here is that Eq. (5) depends on the first and second stress invariants (via σ_eq and σ_m), whereas Eq. (6) also depends on the third invariant via the resolved shear stress and resolved normal stress.

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.

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(σ;f,w,v)=0$⁠, while the inhomogeneous yield criterion depends on v and band orientation n, $FI(σ;fb,w¯,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].

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₁ if p₁ 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.

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₀ = 10⁻³ 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 σ₀ = 420 MPa, ϵ₀ = 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₀ and the band n₀). 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.

There are two types of parameters in the isotropic theory, the initial values of internal variables, f₀ and w₀, and the heuristic parameters q₁, q₂, p₁, and p₂. 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₀ for a range of values of f₀. The larger the initial porosity, the lower the strain to failure across the range of triaxialities. However, the influence of f₀ is much stronger near the cusp. The intensity of the cusp decreases with increasing f₀ and may even disappear, as in the f₀ = 0.05 case for w₀ = 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.005 and f₀ = 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.01. The main effect of w₀ 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₁ and q₂ in the criterion for homogeneous yielding, Eq. (5). This is illustrated in Fig. 5 for spherical voids. Introducing q₁ amounts to multiplying the initial void volume fraction by q₁ so that an effect similar to that of f₀ (see Fig. 3(b)) would be expected. In effect, there are two differences between the influence of f₀ and that of q₁. Not only does q₁ 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₁ 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₁ 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₁.

Likewise, the second Tvergaard parameter q₂ has no influence outside the triaxiality range about the cusp, Fig. 5(b). However, the effect of q₂ appears to be more significant, both around the cusp and for equibiaxial tension. This is so because q₂ 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₁ and p₂, 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₁ and p₂ 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₁ is more gradual, which is a desirable feature in model calibration.

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

Neither the microstructural parameters f₀ and w₀ nor the heuristic parameters q₁, q₂, p₁, and p₂ 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.

In the anisotropic theory of Sec. 2.3, the initial state is defined by the triplet of scalars f₀, w₀, and λ₀ in addition to directors v₀ and n₀. 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₀ = v₀. 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 λ₀ = 1. The effect manifests strongly for small values of f₀ and becomes small for porosities of order 0.01 and above. In the limit f₀ → 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₀ spanning five decades. In fact, what plays a key role in shear failure is not f₀ per se, but the initial ligament ratio χ₀, 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=χ02$⁠. Thus, the four listed values of f₀ correspond to χ₀ = 0.01, 0.05, 0.1, and 0.46 and initial band porosities of 10⁻⁴, 2.5 × 10⁻³, 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 λ₀, Fig. 7(c). For given values of f₀ and w₀, 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 χ₀, 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 λ₀ < 1 and λ₀ > 1 situations in Fig. 8.

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:

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

2. Among heuritsic parameters, the second Tvergaard parameter q₂ has the largest effect on the cusp (Fig. 5(b)).

3. Microstructural parameters f₀ and w₀ have a weak effect on ductility in biaxial tension; the Tvergaard parameters q₁ and q₂ have no effect on it.

4. Heuristic parameters p₁ and p₂ entering the inhomogeneous yielding criterion have a strong influence on biaxial tension ductility (Fig. 6).

5. 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) and (ii) infinite ductility under uniaxial loading (not shown). The parametric analysis reveals the following:

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

2. 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).

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.

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, $ϵf*$⁠, versus a measure of stress state, often a strain-weighted average triaxiality, $T¯$⁠. If the location of crack initiation is known, then an estimate of $ϵf*$ 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, $ϵ¯$⁠, triaxiality, and Lode parameter, are extracted from the calculations and used as basis for plotting experimental data with $ϵf*$ being identified with the value of $ϵ¯$ at global failure. If some triaxiality variations occur at the presumed location of failure, $T¯$ 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₀ and w₀ as well as (if needed) λ₀, v₀, and n₀ and the heuristic parameters q_i and p_i.

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:

1. Measure the initial void volume fraction, f₀, or the initial volume fraction of void nucleation sites, f_p.

2. If neither measurement is available, assume one of two things:

   • Void nucleation by particle cracking, in which case take f₀ ≤ 10⁻³ and w₀ ∼ 0.1–0.2, or

   • Void nucleation by decohesion, in which case take f₀ = f_p and w₀ = w_p, the particle aspect ratio. Estimates of f_p and w_p are available in the literature for most structural alloys.

3. Use heuristic parameter p₁ to calibrate the biaxial tension “plateau”; see Fig. 6(a).

4. Use heuristic parameter q₂ to calibrate the intensity of the cusp; see Fig. 5(b).

5. Use the relative void spacing, λ₀, to calibrate the shear ductility using the anisotropic theory; see Fig. 7(c).

6. 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).

7. Refine as needed by small manipulations of the values of w₀ and f₀.

The availability of microscopic measurements lowers the uncertainty on some microstructural parameters. For example, λ₀ 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₀ and w₀ 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.

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.001 assuming either nucleation by particle cracking (w₀ = 0.1) or by decohesion from equiaxed particles (w₀ = 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₁. The height of the cusp is lower with w₀ = 0.1 using q₂ = 1 and is found to be too high for w₀ = 1 using the maximum value of q₂ = 2. The shear ductility is simulated using λ₀ = 0.38 for w₀ = 0.1 and λ₀ = 0.9 for w₀ = 1. Since the w₀ = 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 (λ₀ < 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.4 for both nucleation scenarios. In this material, the cusp seems to be better captured when assuming w₀ = 1, which suggests more equiaxed nucleation sites. A value of q₂ = 1 was used for the w₀ = 0.1 case. Using a larger value of q₂ 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 λ₀.

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.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.

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₃, nor do they account for anisotropies induced by low-triaxiality loading. The effect of J₃ 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₃ 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₃ 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₃ 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:

1. They lack physical basis.

2. They assume a locally hardening, hence invertible, response.

3. 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₀ 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 λ₀; 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.

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.

1. 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.

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

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

4. The isotropic theory cannot rationalize finite values of shear ductility.

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

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

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.

There are no conflicts of interest.

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.

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.01, w₀ = 1, λ₀ = 1, N = 0.1, and ϵ₀ = 0.002. The values of elastic constants E and ν were as in the text. In one simulation, the initial yield stress was σ₀ = 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 σ₀, 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₀, w₀, λ₀}, not on σ₀.

Furthermore, when a fracture locus is determined using the so-called hybrid method [4,8,18], the problem of error estimation does not seem to have received the attention it deserves. A comprehensive analysis may be found in Ref. [43].