Phase-Field Modeling of Fracture under Compression and Confinement in Anisotropic Geomaterials

TL;DR

The paper develops a fully anisotropic phase-field fracture framework capable of handling crack-face contact under compression and confinement in layered geomaterials. It combines a transversely isotropic homogenized elastic response with an anisotropic, orientation-sensitive fracture energy and a variationally defined crack-normal, solved via a mixed FEM scheme. Validation shows good agreement with fully-resolved layered simulations and with confinement experiments, and the model can predict wing cracks without prescribing external cracks. This approach enables efficient, general analyses of anisotropic fracture in geostructural materials and lays groundwork for future extensions to friction, plasticity, and poromechanics.

Abstract

Strongly anisotropic geomaterials undergo fracture under compressive loading. This paper applies a phase-field fracture model to study this fracture process. While phase-field fracture models have several advantages, they provide unphysical predictions when the stress state is complex and includes compression that can cause crack faces to contact. Building on a phase-field model that accounts for compressive traction across the crack face, this paper extends the model to anisotropic fracture. The key features include: (1) a homogenized anisotropic elastic response and strongly-anisotropic model for the work to fracture; (2) an effective damage response that accounts consistently for compressive traction across the crack face, that is derived from the anisotropic elastic response; (3) a regularized crack normal field that overcomes the shortcomings of the isotropic setting, and enables the correct crack response, both across and transverse to the crack face. To test the model, we first compare the predictions to phase-field fracture evolution calculations in a fully-resolved layered specimen with spatial inhomogeneity, and show that it captures the overall patterns of crack growth. We then apply the model to previously-reported experimental observations of fracture evolution in laboratory specimens of shales under compression with confinement, and find that it predicts well the observed crack patterns in a broad range of loading conditions. We further apply the model to predict the growth of wing cracks under compression and confinement. The effective crack response model enables us to treat the initial crack simply as a non-singular damaged zone within the computational domain, thereby allowing for easy and general computations.

Figures (10)

• Figure 1: The top row shows different loadings, and the middle and lower rows show the idealized deformation for intact and cracked specimens, respectively. Based on this idealization, we assign zero energy to modes (a) and (b). The dotted lines in the second and third row show the undeformed configuration. This decomposition follows steinke2019phase.
• Figure 2: Polar plot of $\left(\frac{\epsilon}{4} \nabla \phi \cdot {\mathbold A} \nabla \phi + \frac{\epsilon ^3}{32} \nabla\nabla\phi : \mathbb{A} : \nabla\nabla\phi\right)$ as a function of $\theta$, where we set $\nabla\phi \sim \cos\theta\sin \theta$ and $\nabla\nabla\phi \sim \cos^2\theta\cos\theta\sin\theta\cos\theta\sin\theta\sin^2 \theta$, following li2015phase. We use $\epsilon = 0.04$, and (a) $\alpha = \alpha_1 = 0$, (b) $\alpha = \alpha_1 = 10$, and (c) $\alpha = \alpha_1 = 100$. We use $\alpha = \alpha_1 = 100$ throughout our calculations. We note that $\alpha$ and $\alpha_1$ must lie in the open range $(-1,\infty)$ for ${\mathbold A}$ and $\mathbb{A}$ to be positive-definite teichtmeister2017phase. The structural director is set to be ${\mathbold a} = {\mathbold e}_2$ (i.e., ${\mathbold a}_{\perp} = {\mathbold e}_1$).
• Figure 3: Left: the vector field ${\mathbold d}$ is not normal to the crack faces when we incorporate an anisotropic penalization of ${\mathbold d}$ in a model that treats ${\mathbold d}$ as a primary variable, generalizing the isotropic model of hakimzadeh2022phase. Right: the model proposed in this paper where the scalar field $\phi$ is a primary variable and ${\mathbold d}$ is computed from the kinematic constraints \ref{['eq:bulk1']}. It is clear that (a) does not work and (b) does.
• Figure 4: (a) The geometry of the fully-resolved specimen showing the inhomogeneous layers. Red indicates Material A (the "strong" material) and blue indicates Material B ("weak" material). The layers have flat interfaces between them, but they appear rugged due to the meshing. (b) A typical finite element mesh used in the calculations. (c) A close-up of the mesh around the center hole.
• Figure 5: Comparison of crack growth in the fully-resolved and homogenized settings at various layer orientations under uniaxial compression. The horizontal faces are subject to displacement-controlled compression, while the vertical faces are traction-free, following tien2006experimental.