Cohesive Zone Interpretations of Phase-Field Fracture Models

The ability to model and predict crack initiation and propagation in both brittle and ductile solids is pivotal to achieving optimal structural design for fracture resistance. This is especially apparent with the advancement of additive manufacturing technologies which are now capable of rapidly producing material structures of varying complexities [1–3]. Micromechanical failure models typically rely on a local approach, where the fracture event is localized within a well-defined fracture process zone, which is embedded within a continuum constitutive model representing the background material [4,5]. The micromechanisms for fracture, such as void growth and coalescence within the fracture process zone, are modeled either discretely or through damage constitutive models, such as the Gurson model or cohesive zone laws [6–8]. Such micromechanical models have enabled fundamental studies to elucidate the coupling between complex mechanisms of fracture within the process zone and plastic deformation of the background material [9–11]. A major limitation, however, is that the process zone (and hence the crack path) has to be established a priori [12–17], which complicates the modeling of convoluted crack patterns, including crack branching and merging.

The phase-field approach to fracture departs from the discontinuous description of failure in local micromechanics-based models. The formulation is based on the variational approach of the classical Griffith energy balance for brittle fracture and regularizes the topology of the sharp crack as a diffuse damage zone instead of a discontinuity [18–23]. As such, the model is able to handle topologically complex fractures and has been widely adopted to solve challenging fracture problems, including hydrogen-assisted cracking [24,25], fracture in viscoelastic materials, biomaterials, or anisotropic solids with different material symmetry [26–30], as well as dynamic brittle fracture with complex branching [31,32].

To-date, there is no clear micromechanical basis for the parameters used in the phase-field fracture models. In particular, the amount of crack regularization is controlled through a prescribed length-scale parameter, l_c, which some have perceived as purely a mathematical construct to allow the Griffin crack (Fig. 1(a), left) to be smeared over a diffused continuum zone (Fig. 1(a), right) [33,34]. In the phase-field formulation, l_c arises through a degradation function to describe the material behavior as it transistions between fully intact and damage states. Classical approach for phase field assumes a stifffness-based degradation function, usually in the form of a hyperbole function, as it allows for an accurate reproduction of linear elastic fracture mechanics response [35–37]. In recent years, an energy-based degradation function has been proposed where certain mechanisms of damage can be incorporated directly into the phase-field model [38]. For all these degradation functions, the phase-field solution should converge to a discrete crack solution in the limit as l_c approaches zero [33]. Others have argued that l_c represents a specific material property that is closely connected to the critical stress for crack nucleation [39–42]. Because of this diffuse crack representation, the fracture energy density term G_c in the phase-field model only approximately relates to the classical Griffith critical energy release rate in the limit of l_c → 0 [43].

To provide mechanistic insights into the above phase-field parameters, one approach is to extract a micromechanics representation of an equivalent fracture process zone of the phase-field fracture model through a cohesive zone law. This cohesive zone law constitutes the relationship between interfacial tractions in equilibrium with the surrounding body and the cohesive separations compatible with the deformation fields of the surrounding body [9,10]. A general view is that the cohesive strength (peak traction) and the cohesive energy (area under the traction-separation relation) are the two primary material parameters governing the macroscopic fracture behavior [13–16]. However, the functional form of the cohesive zone law has been reported to reflect the fracture micromechanisms [44,45]. This has led to the development of inverse techniques to systematically uncover the exact functional form of the cohesive zone laws governed by different failure processes under both monotonic fracture and cyclic fatigue [7,46–49].

In this paper, we explore the relationship between parameters of a phase-field fracture model and the fracture resistance and crack growth mechanics under mode I loading by investigating the equivalent crack-tip cohesive zone laws representing the phase-field fracture model based on a hyperbole stiffness degradation function. Section 2 describes the finite element implementation of the phase-field fracture model, along with the formulation of a small-scale yielding, modified boundary layer model with imposed monotonic K_I remote displacement loading. In Sec. 3, we systematically study the influence of G_c and l_c on both the macroscopic fracture resistance and microscopic fracture processes as quantified through an equivalent crack-tip cohesive zone law. Section 4 discusses the relationship between these phase-field parameters and key process zone characteristics in the material and concludes with a summary.

By expressing (5) in a weak form, we discretize (u, ϕ) in a standard finite element scheme and formulate the residuals and the stiffness matrices. This numerical implementation is described in detail in Ref. [24] and is conducted within a user element subroutine in the commercial finite element software, abaqus.

From dimensional analysis, both the macroscopic fracture resistance, Γ, in (8), and the spatial distribution of microscopic field quantities, σ_ij/E and ϕ, depend on the following dimensionless geometric-material parameters (Γ/ED, ν, G_c/ED, l_c/D). In this study, we set ν = 0.3 and direct attention to the phase-field related parameters (G_c/ED, l_c/D). Note that l_c was originally introduced as a mathematical construct to transform a discrete crack surface into smooth continuum gradient parameters to represent a smeared (diffuse) crack (Fig. 1(a)), implying that (2) → (1) and G_c approaches the fracture toughness only in the limit of l_c → 0. However, convergence studies show that larger (finite) values of l_c are often required for the phase-field fracture models to match experimental results [24–28,56–58], suggesting that (l_c, G_c) can be regarded as phenomenological material properties. In the following, we will parametrically vary G_c/(ED) from 0.005 to 0.03 and l_c/D from 4 to 14 to elucidate the effects of these phase-field parameters on both the microscopic crack growth processes and the macroscopic fracture resistance. For all our calculations, we define the current crack-tip location as the furthest distance along x₂ = 0 where ϕ = 1.

Figure 2 shows the effects of G_c and l_c on the fracture resistance (R-) curves. For all cases, G_c (value denoted by square symbols) quantitatively captures the fracture initiation toughness, Γ_ini, as defined by the energy release rate or J-integral (Γ in (8)) at the first instance of crack growth (i.e., at an elemental crack length of Δa = D). The initial crack growth is characterized by the rapid rise in Γ. The rate of increase of Γ slows down significantly for crack growth beyond about 10D, but the continued rise of Γ with crack length Δa beyond this point suggests the continued toughening of the background material, albeit at rates that depend on both l_c/D and G_c/(ED). As shown in Fig. 2(b), increasing l_c/D from 4 to 14 at a prescribed G_c/(ED) = 0.015 decreases the rate of increase of Γ by ∼50%. Comparatively, reducing G_c/(ED), and hence Γ_ini, from 0.03 to 0.005 at a fixed l_c/D = 4 in Fig. 2(a) decreases the rate of increase of Γ by several-folds.

To elucidate the toughening mechanisms associated with the phase-field parameters, G_c and l_c, we show in Fig. 3 contours of the von Mises stress, σ_e, and the phase-field damage variable, ϕ, for three combinations of G_c/(ED) and l_c/D, and at three crack growth instances, $Δa/D=5,40,and80$⁠. Comparison between G_c/(ED) = 0.03 versus 0.015 at l_c/D = 4 shows that the stresses at the crack-tip are notably larger at higher G_c/(ED) across all Δa/D, with the maximum von Mises stress increasing by nearly two-fold. Interestingly, the contours for ϕ remain almost unchanged. By contrast, increasing l_c/D from 4 to 14 but at a fixed G_c/(ED) = 0.015 dramatically increases the width of the diffuse damage zone ϕ. This smearing of damage over a wider region also corresponds with an almost two-fold decrease in the maximum von Mises stress.

As previously shown in Fig. 2, G_c/(ED) is a quantitative measure of the fracture initiation toughness. In addition, the higher von Mises stress-fields surrounding the propagating crack at higher G_c in Fig. 3 suggests that increasing G_c also increases the stress-carrying capacity ahead of the growing crack, which ultimately leads to improved fracture resistance over the entire transient crack growth regime. On the other hand, l_c/D can be interpreted as a measure of the thickness of the fracture process zone defined by ϕ > 0. As shown in Fig. 3, a larger l_c/D corresponds to a larger and more diffused process zone, while a smaller l_c/D results in a narrower process zone. Since the material stiffness monotonically degrades with g(ϕ) in (3) for ϕ > 0, the reduced stress-carrying capacity in the process zone ahead of the crack-tip at higher l_c/D results in lower fracture resistance. These effects of G_c and l_c are qualitatively in agreement with the stress–strain response of a homogeneously deformed 1D material with $∇ϕ=0$⁠, which shows that the failure stress scales with $Gc/lc$ [23,42,59].

Both the contours of σ_e and ϕ in Fig. 3 are evolving with Δa/D, inferring that the process zone is still developing with crack propagation. This change in near-tip condition with the transition from crack initiation to crack growth is responsible for the rising R-curves in Fig. 2. We show in Fig. 4 the evolution of the near-tip process zone by centering the contours of ϕ about the current crack-tip. Similar contours of ϕ (process zone sizes) are observed between G_c/(ED) = 0.03 and 0.015 with l_c/D = 4. In both cases, the damage contours are continuously expanding at a near constant rate as the crack propagates from Δa = 5D to 80D. In contrast, the contours of ϕ for G_c/(ED) = 0.015 with l_c/D = 14 are nearly two- to three-fold larger. A comparatively smaller increase in process zone size is also observed as the crack grows from Δa/D = 5D to 40D versus Δa/D = 40D to 80D, indicating that the crack growth is reaching its steady-state.

Unlike local micromechanics fracture models where damage is confined to a narrow process zone ahead of the crack [5], diffuse nonlocal approaches such as the phase-field fracture model dissipate damage over a wider region and across many elements as shown in Figs. 3 and 4 [12,60,61]. To provide a homogenized view of the diffuse crack-tip process in the phase-field fracture model, we construct the equivalent local traction-separation relationship constituting the cohesive zone law embedded within a linear elastic background material. Figures 5 and 6 show the effects of G_c/(ED) and l_c/D on the distributions of the crack-tip cohesive tractions, t₂, cohesive separations, δ₂, and the phase-field damage parameter, ϕ, at Δa/D = 25. Both t₂(x₁) and ϕ(x₁) are obtained directly from the phase-field finite element calculations along x₂ = 0. To obtain the equivalent cohesive separation distributions in a linear elastic background material, we impose the measured t₂(x₁) along x₂ = 0 of the same finite element mesh with the same K-field displacements as boundary conditions, albeit with a linear elastic (i.e., non-phase-field) constitutive relationship of the same E and ν, and compute δ₂(x₁) = 2u₂(x₁) along x₂ = 0. In this fashion, we effectively project the diffuse damage ϕ within the finite process zone onto a zero-thickness cohesive zone embedded within a linear elastic body, as shown schematically in the inset of Fig. 5(b).

Increasing l_c/D with the same G_c/(ED) (Fig. 6(a)) is found to reduce the peak cohesive tractions. All post-peak traction distributions, however, converge to the same path further ahead of the crack-tip. In contrast, increasing G_c/(ED) with the same l_c/D (Fig. 5(a)) proportionally increases both the peak and post-peak traction distributions, while maintaining the same traction distribution profile. Increasing G_c/(ED) causes a pronounced increase in the crack-tip separations versus the effects of l_c/D (compare Fig. 5(b) versus Fig. 6(b)). However, increasing l_c/D also increases the length of the crack-tip cohesive zone, as evidenced by both the larger non-zero separations (Fig. 6(b)) and the slower decay of damage (Fig. 6(c)) further ahead of the crack-tip. In comparison, increasing G_c/(ED) has negligble influence on the cohesive zone size—the separation and damage distributions (Figs. 5(b) and 5(c)) both decay to zero at nearly identical distances ahead of the crack-tip.

Together, the above cohesive traction (Figs. 5(a) and 6(a)) and separation (Figs. 5(b) and 6(b)) distributions are used to construct the traction-separation relationships in Fig. 7 constituting the crack-tip cohesive zone laws in a linear elastic background material. Generally, the cohesive strength (peak cohesive traction) and the cohesive energy (area under the (t₂, δ₂) envelope) are two important variables which control the fracture resistance. While G_c/(ED) corresponds to the fracture initiation toughness, Γ_ini, the cohesive energy here quantitatively corresponds to the fracture resistance Γ at the specific crack instant where the cohesive zone law was constructed (Δa/D = 25). As shown in Fig. 7(a), increasing G_c/(ED) increases both the cohesive strength and energy, but the functional form of the cohesive traction-separation relationship remains the same. Increasing l_c/D in Fig. 7(b) decreases the cohesive strength, and to a lesser extent, the cohesive energy. More importantly, a significant change in the functional form of the cohesive zone law is observed—from a rapid increase and subsequent decrease in cohesive tractions at short separation distances for l_c/D = 4 to the flattening of the cohesive zone law resulting in a sustained, but lower, peak traction over longer separation distances (i.e., longer cohesive zone) for l_c/D = 14. These changes in shape of the cohesive zone law with l_c/D can be used to represent changes in the fracture micromechanisms.

In the phase-field model, the damage parameter ϕ provides a measure of the transition from an undamaged (ϕ = 0) to a fully cracked (ϕ = 1) material. An alternative “cohesive zone” approach to quantify crack growth is to consider the cohesive energy dissipated as the separation progresses, as quantified by the evolving area encompassed by the traction-separation response, i.e., $Γs(δ2)=∫0δ2t2(δ2′)dδ2′$⁠, relative to the separation energy Γ. Thus, Γ_s/Γ = 0 at δ₂ = 0 represents a fully intact material, while Γ_s/Γ = 1 when the cohesive zone is fully developed at δ₂ = δ₀ (inset in Fig. 5(b)) denotes a fully separated material. A comparison of the evolution of ϕ and Γ_s/Γ with δ₂ across various (l_c, G_c) combinations in Fig. 8 demonstrates that the damage or crack growth assessment from both approaches are in good quantitative agreement.

The rising R-curves in Fig. 2 suggest that the cohesive zone law is evolving as the crack propagates. Accordingly, we show in Fig. 9 the evolving cohesive zone laws reconstructed from the traction and separation distributions at increasing Δa/D across various l_c and G_c combinations. Note that the cohesive energy of each of these (l_c, G_c, Δa) combinations are quantitatively in perfect agreement with the calculated energy release rate, Γ(Δa/D), in Fig. 2, since all of the dissipation energy in the process zone is projected onto an equivalent zero-thickness cohesive zone. While the shape of the cohesive zone law generally remains unchanged, increasing Δa/D increases both the peak cohesive traction as well as the total separation, δ₀, particularly at high G_c/(ED) of 0.025 where a 16% increase in the cohesive strength is observed as the crack propagates from $Δa/D=$ 10 to 115. At higher l_c/D, the cohesive zone laws rapidly converge with Δa/D, indicating a transition to steady-state crack growth.

In the micromechanics modeling of fracture, a brittle or ductile fracture response often correlates with the extent of plastic dissipation in the background material. In the case of brittle fracture, damage is often confined to the thin process zone ahead of the crack, and the limited plastic dissipation results in a flat R-curve representing rapid and unstable crack propagation. In the case of ductile fracture, the development of significant background plastic dissipation results in a rising R-curve. The phase-field fracture model we have adopted is widely considered to be “brittle” as we have assumed an elastic background material. We remark that “ductile” phase-field fracture models based on elasto-plastic background materials have also been proposed in the open literature [36,56,57,62,63]. However, our crack growth simulations for this “brittle” phase-field fracture model demonstrate a rising R-curve in cases with large G_c and/or small l_c, suggesting that this seemingly “brittle” model can effectively simulate the fracture response of ductile materials.

Akin to crack growth in an elasto-plastic material, the macroscopic fracture resistance Γ in our phase-field model can be delineated into two primary contributions: the fracture initiation toughness, Γ_ini, and phase-field dissipation energy, Γ_p. Our simulation results in Fig. 2 show that G_c is a direct quantitative measure of Γ_ini, in agreement with prior studies [42,53], while both G_c and l_c have profound effects on Γ_p. Similar to the growth of the plastic zone size during transient crack growth, the transition from crack initiation to steady-state crack growth in our phase-field fracture model is marked by the development of a diffuse damage zone which evolves with Δa (Figs. 3 and 4). This diffuse zone can be treated as an evolving fracture process zone in the transient crack growth regime and is distinct from micromechanics fracture models which almost always assume a fixed process zone size.

One approach to quantify the phase-field fracture process is to project the diffuse damage into an equivalent crack-tip cohesive zone law which evolves with crack growth. In doing so, the cohesive energy quantitatively equates to the total fracture resistance, Γ(Δa), which encompasses both the contributions of Γ_ini and Γ_p(Δa). Our results in Figs. 5–7 show that G_c has a strong influence on both the cohesive strength and cohesive energy but has negligible effects on the shape of the cohesive zone law—this suggests that the underlying crack growth mechanics remain unchanged with G_c. In contrast, l_c significantly changes the functional form of the cohesive zone law, and a transition from a sharp traction-separation profile resembling a bilinear cohesive zone model to a flattened traction-separation profile resembling a trapezoidal cohesive zone model is observed with increasing l_c. This change in the shape of the traction-separation relationship is associated with the increasing size of the diffuse damage zone with l_c—a larger l_c creates both a thicker and longer fracture process zone with reduced cohesive strength. By comparison, increasing G_c has negligible influence on the process zone size.

Our studies, therefore, demonstrate that l_c quantitatively relates to the size of the fracture process zone. In the limit of l_c → 0, the process zone physically collapses to a zero-thickness cohesive zone, and the phase-field fracture model would indeed be equivalent to a sharp crack propagating within an elastic medium. The fracture process would always be brittle with Γ = G_c, since there would be no energy dissipation in the background material. However, convergence studies often lead to the adoption of relatively large values of l_c/D > 10 [24,25,29,62], which introduces a finite thickness fracture process zone. Differing from these numerical studies, our results suggest that l_c should be calibrated and selected based on the physical thickness of the fracture process zone. For example, a small value of l_c/D should be used to represent the narrow process zone for cracking in a brittle rock-like material such as shale or concrete [64,65], while a larger l_c/D will better represent the more diffused process zone associated with micro-crazing in polymers [66,67]. This would provide l_c/D with a stronger physical basis. The “brittle” phase-field fracture model we adopt does not account for background plasticity. Without loss of generality, our simulations suggest that the background plastic dissipation in elasto-plastic materials constituting Γ_p can heuristically be treated as a diffuse damage zone represented by appropriately calibrated values of l_c/D.

The cohesive strength within the process zone depends on both G_c and l_c. Uniform deformation studies show that the critical strength in the stress–strain response of a phase-field element scales with $Gc/lc$ [23,42,59], and the cohesive strength of our equivalent crack-tip cohesive zone laws appear to follow similar trends. Conceivably, once l_c is calibrated to represent the size of the fracture process zone, G_c can in turn be calibrated to fit the cohesive strength representing the appropriate crack growth mechanisms, such as void growth and coalescence, fiber pull-out, and phase transformation.

In conclusion, we have obtained new physical insights into the role of the energy and length-scale parameters, G_c and l_c, in phase-field fracture models, by constructing the equivalent crack-tip cohesive zone laws representative of the diffuse damage process. We demonstrate that G_c can be perceived as the fracture initiation toughness, and quantitatively controls both the cohesive strength and energy. While l_c was initially introduced as a mathematical construct, we show that l_c should be physically interpreted as a measure of the thickness of the diffuse damage (fracture process) zone. By adopting finite values of l_c, together with appropriately calibrated values of G_c, the “brittle” phase-field fracture model can be used to simulate a wide variety of material systems, from brittle rock materials with narrow process zones to ductile metals with diffuse background plastic dissipation. We remark that a commonly used stiffness-based hyperbole degradation function is considered in this study. Extracting the equivalent cohesive zone laws for other degradation functions is a subject of future work.

The authors acknowledge the support provided by National Science Foundation under Grant Nos. NSF-CMMI-2009684 and NSF-DMR-18-09696.

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.