Unlike micromechanics failure models that have a well-defined crack path, phase-field fracture models are capable of predicting the crack path in arbitrary geometries and dimensions by utilizing a diffuse representation of cracks. However, such models rely on the calibration of a fracture energy (Gc) and a regularization length-scale (lc) parameter, which do not have a strong micromechanical basis. Here, we construct the equivalent crack-tip cohesive zone laws representing a phase-field fracture model, to elucidate the effects of Gc and lc on the fracture resistance and crack growth mechanics under mode I K-field loading. Our results show that the cohesive zone law scales with increasing Gc while maintaining the same functional form. In contrast, increasing lc broadens the process zone and results in a flattened traction-separation profile with a decreased but sustained peak cohesive traction over longer separation distances. While Gc quantitatively captures the fracture initiation toughness, increasing Gc coupled with decreasing lc contributes to a rising fracture resistance curve and a higher steady-state toughness—both these effects cumulate in an evolving cohesive zone law with crack progression. We discuss the relationship between these phase-field parameters and process zone characteristics in the material.
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, lc, 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, lc 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 . For all these degradation functions, the phase-field solution should converge to a discrete crack solution in the limit as lc approaches zero . Others have argued that lc 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 Gc in the phase-field model only approximately relates to the classical Griffith critical energy release rate in the limit of lc → 0 .
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 KI remote displacement loading. In Sec. 3, we systematically study the influence of Gc and lc 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.
2 Problem Formulation
2.1 Phase-Field Fracture Modeling.
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.  and is conducted within a user element subroutine in the commercial finite element software, abaqus.
2.2 Boundary Value Problem.
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, ν, Gc/ED, lc/D). In this study, we set ν = 0.3 and direct attention to the phase-field related parameters (Gc/ED, lc/D). Note that lc 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 Gc approaches the fracture toughness only in the limit of lc → 0. However, convergence studies show that larger (finite) values of lc are often required for the phase-field fracture models to match experimental results [24–28,56–58], suggesting that (lc, Gc) can be regarded as phenomenological material properties. In the following, we will parametrically vary Gc/(ED) from 0.005 to 0.03 and lc/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 x2 = 0 where ϕ = 1.
3.1 Fracture Resistance and Crack Growth Process.
Figure 2 shows the effects of Gc and lc on the fracture resistance (R-) curves. For all cases, Gc (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 lc/D and Gc/(ED). As shown in Fig. 2(b), increasing lc/D from 4 to 14 at a prescribed Gc/(ED) = 0.015 decreases the rate of increase of Γ by ∼50%. Comparatively, reducing Gc/(ED), and hence Γini, from 0.03 to 0.005 at a fixed lc/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, Gc and lc, we show in Fig. 3 contours of the von Mises stress, σe, and the phase-field damage variable, ϕ, for three combinations of Gc/(ED) and lc/D, and at three crack growth instances, . Comparison between Gc/(ED) = 0.03 versus 0.015 at lc/D = 4 shows that the stresses at the crack-tip are notably larger at higher Gc/(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 lc/D from 4 to 14 but at a fixed Gc/(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, Gc/(ED) is a quantitative measure of the fracture initiation toughness. In addition, the higher von Mises stress-fields surrounding the propagating crack at higher Gc in Fig. 3 suggests that increasing Gc 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, lc/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 lc/D corresponds to a larger and more diffused process zone, while a smaller lc/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 lc/D results in lower fracture resistance. These effects of Gc and lc are qualitatively in agreement with the stress–strain response of a homogeneously deformed 1D material with , which shows that the failure stress scales with [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 Gc/(ED) = 0.03 and 0.015 with lc/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 Gc/(ED) = 0.015 with lc/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.
3.2 Crack-Tip Cohesive Zone Laws.
Unlike local micromechanics fracture models where damage is confined to a narrow process zone ahead of the crack , 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 Gc/(ED) and lc/D on the distributions of the crack-tip cohesive tractions, t2, cohesive separations, δ2, and the phase-field damage parameter, ϕ, at Δa/D = 25. Both t2(x1) and ϕ(x1) are obtained directly from the phase-field finite element calculations along x2 = 0. To obtain the equivalent cohesive separation distributions in a linear elastic background material, we impose the measured t2(x1) along x2 = 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 δ2(x1) = 2u2(x1) along x2 = 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 lc/D with the same Gc/(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 Gc/(ED) with the same lc/D (Fig. 5(a)) proportionally increases both the peak and post-peak traction distributions, while maintaining the same traction distribution profile. Increasing Gc/(ED) causes a pronounced increase in the crack-tip separations versus the effects of lc/D (compare Fig. 5(b) versus Fig. 6(b)). However, increasing lc/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 Gc/(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 (t2, δ2) envelope) are two important variables which control the fracture resistance. While Gc/(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 Gc/(ED) increases both the cohesive strength and energy, but the functional form of the cohesive traction-separation relationship remains the same. Increasing lc/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 lc/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 lc/D = 14. These changes in shape of the cohesive zone law with lc/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., , relative to the separation energy Γ. Thus, Γs/Γ = 0 at δ2 = 0 represents a fully intact material, while Γs/Γ = 1 when the cohesive zone is fully developed at δ2 = δ0 (inset in Fig. 5(b)) denotes a fully separated material. A comparison of the evolution of ϕ and Γs/Γ with δ2 across various (lc, Gc) 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 lc and Gc combinations. Note that the cohesive energy of each of these (lc, Gc, Δ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, δ0, particularly at high Gc/(ED) of 0.025 where a 16% increase in the cohesive strength is observed as the crack propagates from 10 to 115. At higher lc/D, the cohesive zone laws rapidly converge with Δa/D, indicating a transition to steady-state crack growth.
4 Discussions and Conclusion
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 Gc and/or small lc, 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 Gc is a direct quantitative measure of Γini, in agreement with prior studies [42,53], while both Gc and lc 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 Gc 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 Gc. In contrast, lc 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 lc. This change in the shape of the traction-separation relationship is associated with the increasing size of the diffuse damage zone with lc—a larger lc creates both a thicker and longer fracture process zone with reduced cohesive strength. By comparison, increasing Gc has negligible influence on the process zone size.
Our studies, therefore, demonstrate that lc quantitatively relates to the size of the fracture process zone. In the limit of lc → 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 Γ = Gc, since there would be no energy dissipation in the background material. However, convergence studies often lead to the adoption of relatively large values of lc/D > 10 [24,25,29,62], which introduces a finite thickness fracture process zone. Differing from these numerical studies, our results suggest that lc should be calibrated and selected based on the physical thickness of the fracture process zone. For example, a small value of lc/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 lc/D will better represent the more diffused process zone associated with micro-crazing in polymers [66,67]. This would provide lc/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 lc/D.
The cohesive strength within the process zone depends on both Gc and lc. Uniform deformation studies show that the critical strength in the stress–strain response of a phase-field element scales with [23,42,59], and the cohesive strength of our equivalent crack-tip cohesive zone laws appear to follow similar trends. Conceivably, once lc is calibrated to represent the size of the fracture process zone, Gc 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, Gc and lc, in phase-field fracture models, by constructing the equivalent crack-tip cohesive zone laws representative of the diffuse damage process. We demonstrate that Gc can be perceived as the fracture initiation toughness, and quantitatively controls both the cohesive strength and energy. While lc was initially introduced as a mathematical construct, we show that lc should be physically interpreted as a measure of the thickness of the diffuse damage (fracture process) zone. By adopting finite values of lc, together with appropriately calibrated values of Gc, 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.
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.