Crack Face Contact Modeling is Essential to Predict Crack-Parallel Stresses Open Access

Phase-field fracture models provide a powerful approach to modeling fracture, potentially enabling the unguided prediction of the crack growth in complex patterns. To ensure that only tensile stresses and not compressive stresses drive crack growth, several models have been proposed that aim to distinguish between the compressive and tensile loads. However, these models have a critical shortcoming: they do not account for the crack direction, and hence they cannot distinguish between crack-normal tensile stresses that drive crack growth and crack-parallel stresses that do not. In this study, we apply a phase-field fracture model, developed in our earlier work, that uses the crack direction to distinguish crack-parallel stresses from crack-normal stresses. This provides a transparent energetic formulation that drives cracks to grow in when crack faces open or slide past each other, while the cracks respond like the intact solid when the crack faces contact under normal compressive loads. We compare our approach against two widely used approaches, Spectral splitting and the volumetric-deviatoric splitting, and find that these predict unphysical crack growth and unphysical stress concentrations under loading conditions in which these should not occur. Specifically, we show that the splitting models predict spurious crack growth and stress concentration under pure crack-parallel normal stresses. However, our formulation, which resolves the crack-parallel stresses from the crack-normal stresses, predicts these correctly.

Phase-field fracture models introduce a new field, in addition to the deformation, that indicates if a material element is intact or fractured, and then regularize this phase-field using higher-order derivatives. The higher-order derivatives endow the phase-field with sufficient smoothness to enable the use of standard numerical methods to predict the crack growth without specifying crack paths beforehand. In other words, rather than sharp singular cracks as in the classical approach to fracture, cracks are now regularized to be nonsingular variations of the phase-field, making them amenable to standard numerical approaches. The elastic response is coupled to the phase-field by setting that the regions that are fractured cannot carry any load, i.e., the elastic energy density goes to zero in fractured regions.

However, the combination of crack regularization to a finite volume and setting the fractured elastic energy to zero causes unphysical behavior, such as crack growth under compressive loads. To address these issues, Miehe et al. [1] proposed a partition of the elastic energy into tensile and compressive parts (denoted the “spectral” split), and allowed the damaged region to sustain only the compressive part. An alternative class of approaches, following Amor et al. [2], proposed a partition of the energy into compressive hydrostatic, tensile hydrostatic, and deviatoric parts (denoted the “volumetric-deviatoric” or “VolDev” split), and allowed the damaged region to sustain only the compressive hydrostatic part. This has been extended to the finite deformation setting in Refs. [3–6]; however, none of these works consider the orientation of the crack. Other modifications of the energy decompositions include Ref. [7] that proposed the star-convex model as a modification of the volumetric-deviatoric decomposition; Ref. [8] that extended the spectral decomposition to account for mixed mode fracture; and Ref. [9] to account for anisotropy. These approaches are reviewed and compared in Refs. [10–13] and elsewhere.

However, all of the spectral and VolDev splitting approaches have a key shortcoming: they do not account for the orientation of the crack. Compressive stresses along the crack face and across the crack face are both treated in exactly the same way. For instance, tension across the crack faces that drives the crack growth is not distinguished from tension along the length of the crack, that does not directly drive crack growth but it is important in other ways. For instance, Bazant et al. [14] highlighted that phase-field models do not accurately capture the influence of crack-parallel stresses, which significantly impact the fracture energy and the fracture process zone width; they argue that the inability of phase-field models to incorporate the tensorial nature of stress interactions leads to discrepancies when compared to experimental results from the gap test. Similarly, Nguyen et al. [15] emphasized that phase-field models lack the ability to explicitly account for the history-dependent effects of T-stress, which are crucial for realistic simulations of mixed-mode fracture in concrete and quasi-brittle materials.

An important departure from the idea of splitting approaches is the work by Steinke and Kaliske [16], introduced the idea that the crack orientation plays an essential role in determining the response. Based on this idea, Refs. [16–19] proposed modifications of the stress tensor to account for the difference in crack-normal and crack-parallel tractions. While our approach in Refs. [20,21] builds on the ideas of Ref. [16], an important distinction is that we modify the energy rather than the stress, and hence our model is hyperelastic. In contrast, as shown in Ref. [20], the stress-based approach is not generally hyperelastic. Our approach enables us to distinguish between the crack-parallel stresses and crack-closing stresses as well as allow for crack faces sliding past each other. Our approach applies a QR (upper triangular) decomposition of the deformation gradient in the basis associated with the local crack normal, enabling us to distinguish between the crack faces sliding past each other—which is formulated to cost no energy; crack opening—which is also formulated to cost no energy; and crack-closing and crack face contact—which are formulated to have an elastic response that is identical to the intact material. This is shown schematically in Fig. 1.

In this article, we compare the phase-field fracture model developed by Hakimzadeh et al. [20,21], who distinguished between the crack-parallel stresses and crack-closing stresses, against the spectral [1] and VolDev [2] splitting approaches. Specifically, in Sec. 3, we examine the effect of tensile and compressive crack-parallel stresses and show that the crack growth and stresses predicted by the spectral and VolDev approaches are unphysical, in contrast to the model from Ref. [20].

A key limitation of the standard phase-field fracture model is its inability to distinguish between tensile, shear, and compressive loadings across or along the crack face, often leading to unphysical crack propagation under compression. Various modifications have been proposed to address this issue, but they fail to incorporate the crack orientation explicitly [1,2,6–10,27]. As noted in Ref. [16] and further demonstrated in Refs. [20,21] and elsewhere, correctly capturing the crack behavior under mixed-mode loading requires a formulation that accounts for the crack face normal.

To overcome this limitation, we adopt the phase-field formulation introduced in Ref. [20], which incorporates an effective crack energy density $Wd$ that depends on the crack orientation. Using the QR (or Gram–Schmidt) decomposition of the deformation gradient $∇y$⁠, e.g., Ref. [28], in the basis of the crack enables us to distinguish the various modes shown in Fig. 1. We then assigned zero energy to modes in which the crack faces slide past each other or open up, and assign the intact response to modes in which the crack faces contact. This ensures a physically consistent response, particularly under compressive loading, by preventing unphysical interpenetration of crack faces, as demonstrated through several examples in Refs. [20,21].

The key feature of this model is the introduction of the effective crack energy density $Wd$ that is active in the fractured volumes. It imparts the properties of an idealized sharp crack to the regularized (finite volume) phase-field crack. The function $Wd$ takes as argument the deformation gradient, as well as a vector field $d$ that corresponds to the normal to the crack face which is essential to be able to decompose the deformation into the modes shown in Fig. 1. However, $d$ is not an independent field; rather it is a function of $ϕ$ as described in Ref. [21]. For brevity, we do not present the detailed expressions for $Wd$ but refer instead to Refs. [20,21]. In this article, we use an isotropic linear elastic response, with Lamé moduli $λ$ and $μ$⁠.

A simplistic approach to modeling the energy density of damaged volumes is to assume that the damaged region carries zero elastic energy. However, this assumption often leads to unphysical behaviors, such as the interpenetration of crack faces, incorrect mechanical response when the crack closes, and crack growth under compression. To address these shortcomings, alternative formulations have been proposed, notably by the spectral split by Miehe et al. [1] and the VolDev split by Amor et al. [2], which decompose the energy and selectively associate specific energy components with the damaged region. While these methods offer improvements over the zero-energy assumption, they do not account for crack orientation in the energetic decomposition, leading to unphysical predictions under complex stress states.

In the Appendix of Ref. [20], we identified significant challenges associated with spectral splitting. Specifically, this approach does not accurately capture the crack response under complex loading conditions, such as uniaxial tension parallel to the crack. Despite these limitations, many studies on the fracture modeling continue to rely on this approach or variants like those proposed by Lo et al. [29], which builds on the work of Ref. [1]. Recent studies in studying the crack behavior under complex loadings, including Ref. [30], utilize this extended method by Lo et al. [29]. Other examples of fracture modeling include recent work by Clayton et al. [31], which adopts the formulation introduced by Miehe et al. [32] and Stocek et al. [33], which is built on the frameworks of Refs. [34] and [1].

As discussed in the Appendix of our previous work [20], the regularized approach proposed by Amor et al. [2] encounters limitations in accurately capturing the crack response under complex stress states. These include uniaxial tension parallel to the crack plane and uniaxial compression normal to the crack. Despite these limitations, this approach remains widely used in several studies that model crack behavior in complex loadings.

We employ the FEM for numerical simulations, using the open-source FEniCS library, which has been extensively utilized in mechanics applications [35,36]. The details are provided in Ref. [37] for the spectral split approach, Ref. [38] for the VolDev split approach, and Ref. [20] for our approach.

Here we analyze the crack growth and the stress distribution when the crack is subject to crack-parallel (tensile and compressive) stresses using spectral splitting, VolDev splitting, and our approach.

Consider a specimen containing a crack with normal $e2$⁠, subjected to a far-field stress of $σ0e1⊗e1$⁠, with $σ0>0$⁠. Assuming isotropic elasticity, this corresponds equivalently to the following far-field strain:

We highlight two key unphysical aspects of this result. First, when the Poisson ratio $ν$ is positive, we find an unphysical compressive stress in the $e2$ direction. Second, the stress along the $e1$ direction is zero, whereas under a far-field stress parallel to the crack, the material is expected to behave as an intact material, sustaining the far-field stress as it is. This incorrect stress distribution leads to spurious crack growth. The artificial lack of stiffness in tension parallel to the crack creates spurious stress concentrations near the crack tip, as the intact material ahead of the tip sustains the far-field stress in the $e1$ direction, while the cracked region does not. Moreover, this discrepancy can significantly affect the calculation of the T-stress [14].

In the VolDev model, we obtain $σsplit=0$⁠. However, the material should fully sustain the far-field tensile stress parallel to the crack. This unphysical response leads to spurious stress concentrations at the crack tip, which influence both the crack driving force and the calculation of the T-stress.

We next numerically evaluate the crack path and stress distribution in a finite domain. For a clear comparison, we use the same domain, mesh, and material properties across all three approaches, and first set $ν=0$ to eliminate crack faces separating or contacting due to the crack-parallel loading. We then examine having $ν>0$ and $ν<0$⁠.

In our calculations, we set the vertical displacement and the horizontal tractions to zero on the top and bottom faces; the horizontal displacement is controlled on the left and right faces; and the vertical tractions are set to zero on the left and right faces.

We first examine the effect of tensile loading in Fig. 2, which shows the initial crack, the subsequent crack growth, and the stress distribution at the first load step, for all three models. For simplicity of interpretation, we set $ν=0$⁠. We observe that our approach predicts essentially a uniform stress distribution, that corresponds to the cracked specimen having the same response as an intact specimen, and there is no crack growth. Both of these features correspond to physical expectations.

In both the spectral and VolDev splits, the stresses are significantly non-uniform with spurious stress concentrations at the crack tips that were anticipated by the simple analysis above. The spurious stress distribution, in turn, drives unphysical crack growth.

We next examine the effect of compressive loading in Fig. 3, which shows the initial crack, the subsequent crack growth, and the stress distribution at the first load step, for all three models; we set $ν=0$⁠. We observe again that our approach predicts essentially a uniform stress distribution, that corresponds to the cracked specimen having the same response as an intact specimen, and there is no crack growth. Both of these features correspond to physical expectations.

Further, in both the spectral and VolDev splits, the stresses are significantly non-uniform with spurious stress concentrations at the crack tips that were anticipated by the simple analysis above. The spurious stress distribution, in turn, drives unphysical crack growth that should simply not occur given that the specimen is under compression.

In the calculations above, we have set $ν=0$ to aid the interpretation of the stress fields. Here we compare the models with crack-parallel tension and compression, with positive and negative $ν$⁠. We examine a horizontal crack in a large rectangular specimen, with the top and bottom traction-free, and the displacement fixed—to be zero in the vertical direction and nonzero in the horizontal direction—on the sides.

Figure 4 compares the stresses and deformations under crack-parallel compression with $ν>0$⁠. In our model, the stress field is uniform and the deformation corresponds to uniform barreling due to the lateral expansion, as expected physically. However, in both the spectral splitting and VolDev models, the stress is clearly inhomogeneous when it should not be, and the deformation around the crack shows unphysical features.

Figure 5 compares the models under crack-parallel compression with $ν<0$⁠. In our model, as expected, there is lateral contraction that causes the crack faces to lose contact and move apart; the stress field is therefore inhomogeneous with stress concentrations at the crack tips, as is typical. However, in both the spectral splitting and VolDev models, while they appear to capture the lateral contraction, the stress distribution is not at all as expected around crack tips, and the deformation shows unphysical features near the crack tips.

Figure 6 compares the models under crack-parallel tension with $ν>0$⁠. In our model, as expected, there is lateral contraction that causes the crack faces to lose contact and move apart, and the stress field is therefore inhomogeneous with stress concentrations at the crack tips. However, in both the spectral splitting and VolDev models, while they capture the lateral contraction, the stress distribution is not at all as expected around crack tips, and the deformation shows unphysical features near the crack tips.

Finally, Fig. 7 compares the models under crack-parallel tension with $ν<0$⁠. In our model, as expected, there is lateral expansion and the overall deformation is uniform as the crack faces remain in contact; the stress field is similarly homogeneous around the crack. However, in both the spectral splitting and VolDev models, the deformation is completely unphysical: the material within the damaged region shrinks laterally, corresponding to the crack faces inter-penetrating. Further, the stress distribution is not at all as expected around crack tips.

We have shown that the widely-used spectral and VolDev splitting approaches lead to unphysical crack growth and spurious stress concentrations in the important settings of crack-parallel stresses. We have further shown that accounting for the crack face contact through the crack-response model proposed in Ref. [20] provides the correct predictions in these settings, suggesting the use of this model in situations where crack-parallel stresses play an important role [14].

We have idealized the crack as having no friction between the crack faces for simplicity. For future studies aiming to infuse more realistic physics into the model, it is important to further incorporate friction in the model [17,18]. Another important area of application is to couple the model to poromechanics, e.g., by adapting Refs. [28,30,39–42]. A further important area is to combine the crack-response model with crack nucleation models, following Refs. [43–47]. For all of these future directions, it is essential to predict the correct stress fields that in turn provide the correct crack growth driving force.

Our focus in this article has been on the small-strain regime for simplicity. However, the crack-response model of Ref. [20] is valid in the large deformation setting and can be readily incorporated in such models, e.g., Refs. [13,48–51].

We thank Meenu Krishnan for useful discussions. We acknowledge financial support from NSF (2108784, 2012259, and 2342349), ARO (MURI W911NF-24-2-0184), AEI (PID2021-124195NB-C32, PCI2024-155023-2, CEX2023-001347-S funded by MICIU/AEI/10.13039/501100011033), and ERC (834728); and NSF for XSEDE computing resources provided by Pittsburgh Supercomputing Center.

There are no conflicts of interest.

The code implementing all three approaches for the calculations in this article is available online.²