Advancements in Phase-Field Modeling for Fracture in Nonlinear Elastic Solids under Finite Deformations

Abstract

The prediction of failure mechanisms in nonlinear elastic materials holds significant importance in engineering applications. In recent years, the phase-field model has emerged as an effective approach for addressing fracture problems. Compared with other discontinuous fracture methods, the phase-field method allows for the easy simulation of complex fracture paths, including crack initiation, propagation, coalescence, and branching phenomena, through a scalar field known as the phase field. This method offers distinct advantages in tackling complex fracture problems in nonlinear elastic materials and exhibits substantial potential in material design and manufacturing. The current research has indicated that the energy distribution method employed in phase-field approaches significantly influences the simulated results of material fracture, such as crack initiation load, crack propagation path, crack branching, and so forth. This impact is particularly pronounced when simulating the fracture of nonlinear materials under finite deformation. Therefore, this review outlines various strain energy decomposition methods proposed by researchers for phase-field models of fracture in tension–compression symmetric nonlinear elastic materials. Additionally, the energy decomposition model for tension–compression asymmetric nonlinear elastic materials is also presented. Moreover, the fracture behavior of hydrogels is investigated through the application of the phase-field model with energy decomposition. In addition to summarizing the research on these types of nonlinear elastic body fractures, this review presents numerical benchmark examples from relevant studies to assess and validate the accuracy and effectiveness of the methods presented.

1. Introduction

With advancements in science and technology, new materials, including nonlinear elastic materials, have been developed. Typical nonlinear elastic materials, including rubber polymer materials, composite materials, hydrogels, and biological materials, are extensively utilized in various applications. Consequently, accurate prediction of failure mechanisms in nonlinear elastic materials is of paramount importance in engineering applications.

In recent years, phase-field simulation methods for crack propagation have received wide attention. This approach, also referred to as the fracture variational approach, enables the simulation of complex crack patterns, including crack intersections, without the need for direct modeling of discontinuities. The initial proposal of this method, which minimizes potential energy to solve the fracture problem, can be attributed to Francfort and Marigo et al. [1]. They drew inspiration from the works of Mumford and Shah et al. [2], as well as Ambrosio and Tortorrelli et al. [3], on the Griffith brittle fracture theory. Additionally, Miehe et al. [4] introduced an alternative method based on the principles of continuum mechanics and thermodynamics for small deformations. This method incorporated a crucial mechanism to differentiate the effects of tension and compression on crack propagation.

This review focuses on the utilization of phase-field methods for simulating fracture in rubber-like polymer materials under finite deformation. This method has been further developed to simulate dynamic fracture and has shown good agreement with multiple benchmark problems [5,6,7,8,9]. Miehe [10], Ambati [11], and McAuliffe [12] have also studied the application of the phase-field method in elastoplastic materials. They found that the accumulation of plastic deformation led to a reduction in elastic strength yet had no effect on the plastic response. Furthermore, Ambati et al. [13] extended the phase-field model for kinematically linear elastic fractures to three-dimensional finite strain environments. Subsequently, Borden et al. [14] made further advancements to the model by deriving micro-force control equations for energy potential under finite deformation. Although previous studies by Miehe and Schänzel et al. [15], as well as Chen et al. [16], also employed phase-field methods to model fracture in hyper-elastic solids under finite deformation, their approach included an artificial viscosity-based regularization formula that lacked the ability to differentiate between tensile and compressive effects during crack growth, which impeded the accuracy of numerical simulations.

To overcome the constraints of prior methodologies, Borden et al. [14], as well as Arriaga and Waisman et al. [17] split the elastic deformation energy into volumetric and deviatoric parts within the phase-field framework to distinguish the contributions of tensile and compressive stresses to fracture of elastoplastic materials under small deformations. Although the spectral decomposition of the logarithmic strain tensor into tensile and compressive parts [18] emerged more reasonable according to the studies by Miehe et al. [4] and Borden et al. [8], Borden et al. pointed out that this method may lead to reduced numerical accuracy and may make convergence hard [14]. A novel strain energy decomposition method for phase-field modeling of fracture in nonlinear elastic materials was proposed by Tang and Zhang et al. [19]. In contrast to the methods utilized by Miehe and Hofacker et al. [9], as well as Borden et al. [8], this approach involved the decomposition of energy into two parts: one that pertained to the principal stretch and the other exclusively comprising nonlinear terms of volumetric deformation. By decomposing the energy in this manner, the distinct functions of tension and compression could be identified, facilitating the direct simulation of actual fracture and leading to a more realistic depiction of the physical process.

This paper reviews existing phase-field models utilized for studying fracture in nonlinear elastic materials under finite deformations. Section 2 provides an overview of the available phase-field models for studying nonlinear elastic fracture under finite deformation, presenting the fundamental equations for different types of nonlinear elastic fracture modeling. In Section 3, simulation methods for nonlinear elastic fracture of tension–compression asymmetry materials are discussed. Section 4 investigates the fracture behavior of hydrogels. Lastly, Section 5 summarizes the current research landscape on nonlinear elastic fracture, identifies unresolved issues, suggests future research directions, and concludes the review.

2. Fracture Study of Nonlinear Elastic Materials with Tension–Compression Symmetry

2.1. Relevant Work on Finite Deformation Phase-Field Fracture Models

This section presents a chronological overview of recent advancements in simulating fracture propagation using the phase-field method under finite deformation conditions within the last five years. Dittmann et al. [20] proposed a unified computational framework that can simulate large deformation, thermodynamics, fracture, and contact problems in a unified and thermodynamically consistent manner. Fang et al. [21] studied the effect of yield function and material hardening on the behavior of elastic–plastic solid materials by introducing a unified yield criterion and isotropic hardening in the phase-field model. Tang et al. [19] developed a new energy decomposition method within a finite deformation framework, considering both the stretching and volume parts. This method incorporates the distinct contributions of stretching and compression into the phase-field model for simulating crack initiation and propagation. Yin et al. [22] established an anisotropic phase-field model under finite strain based on a second-order phase-field theory. They derived an equivalent fracture surface energy density function by evaluating the fracture state of both the homogeneous matrix and fiber materials using a single phase-field variable.

Tarafder et al. [23] devised a phase-field model under finite deformation for multi-phase composite materials, enabling the simultaneous prediction of cracks in individual phases and interface debonding. Numerical examples illustrate the notable distinctions in fracture behavior between finite deformation and small deformation solutions. Tian et al. [24] proposed a new dynamic phase-field model capable of simulating crack propagation in nonlinear deformation while achieving crack propagation speeds close to asymptotic limits. In addition, they developed an adaptive distorted mesh deletion scheme to tackle the issue of FE (finite element) mesh distortion in fractures involving large deformations. Peng et al. [25] presented a phase-field method, employing ES-FEM, to simulate the fracture behavior of hyper-elastic materials under finite deformation. Thomas et al. [26] proposed a phase-field model under finite strain that accounts for damage anisotropy by decomposing the right Cauchy–Green strain tensor into anisotropic components. The model introduces a suitable weak solution concept, allowing for spatial and temporal discretization. Eldahshan et al. [27] developed a plastic fracture phase-field model within a large plastic strain framework by integrating an adaptive isotropic mesh re-partitioning strategy with a plastic fracture phase-field model.

Swamynathan et al. [28] proposed a methodology to split the energy density into tension and compression components in phase-field analysis of finite deformation fracture. The remarkable feature of this model is that the combination of tension and compression precisely equals total reference energy of the undamaged system, making it applicable to different finite strain energy densities. Hu et al. [29] presented a novel plastic fracture phase-field model in a consistent variational framework under finite deformation kinematics. Peng et al. [30] devised a fourth-order phase-field model that considers the decomposition of strain energy during finite deformation to simulate the fracture behavior of hyper-elastic materials under finite deformation. Dan et al. [31] formulated a finite deformation finite element model to analyze the electric–elastic coupling response and crack evolution in heterogeneous piezo-electric composite materials. Phase-field calculation models necessitate a high grid resolution in the vicinity of cracks. To address the high gradient problem of two phase-field order parameters, an adaptive wavelet-enhanced hierarchical finite element framework was employed, optimizing the selection of wavelet-based hierarchical enrichment functions. Koutromanos et al. [32] investigated the toughness fracture of structural steel under cyclic loading by combining the phase-field method with large deformation mechanics calculations, accounting for the materials’ elastic–plastic behavior. Han et al. [33] introduced plastic driving force and fracture toughness degradation into crack calculations, proposing a new toughness fracture phase-field model. This model allows for the generation of various cracks by adjustment of the activation conditions of the plastic driving force and reducing the fracture toughness. Zhao et al. [34] employed phase-field variables of crack surface density to simulate crack propagation and internal damage variables, thereby bridging the micro-damage of multiple network elastic bodies and macroscopic fractures in the finite deformation model.

The advancements in these studies demonstrate the continuous development and improvement of the phase-field method in simulating fracture propagation under finite deformation conditions, providing a robust tool and theoretical foundation for investigating complex material fracture behavior.

2.2. Material Model and Phase-Field Formulation

The deformation gradient tensor serves as the description of the deformation measure under finite deformation. As the object undergoes deformation, point $X$ in the unaltered reference configuration is transformed to a point $x$ in the present deformed configuration through the deformation gradient (see Figure 1), then the deformation gradient and its Jacobian are represented.

$$F=\frac{\partial x}{\partial X}$$

(1)

$$J=\det \left(F\right)$$

(2)

The discontinuity surfaces $\mathsf{\Gamma }$ within represent an aggregation of discrete cracks. Following Griffith’s fracture theory, the energy required to produce a unit area of fracture surface is equivalent to the critical fracture energy density $G_c$. To avoid the issue of numerically handling the discontinuity that represents a crack, the fracture surface $\mathsf{\Gamma }$ is approximated using a phase field $p\left(X,t\right)$:

$${\displaystyle {\int }_{\mathsf{\Gamma }}{G}_c}dS={\displaystyle {\int }_B{G}_c\left[\frac{{\left(p-1\right)}^2}{4l}+\frac{{\left(p-1\right)}^2}{4l}+l\frac{\partial p}{\partial X_1}\frac{\partial p}{\partial X_1}\right]}dX$$

(3)

where $l$ is a model parameter that controls the width of the smooth approximation of the crack. All computations in this review are conducted in Lagrangian coordinates and all integrals are computed within the unaltered reference configuration. The phase field $p$ takes a value of 1 outside the crack and 0 inside the crack (refer to Figure 1).

Equations (1)–(3) presented above describe a nonlinear elastic model for tension–compression symmetry materials. It should be noted that this model is applicable to both tension–compression asymmetry materials and hydrogels, which will be discussed later in this review. Hence, further elaboration is deemed unnecessary.

The simplest neo-Hookean model is used as an archetype. Currently, this method can be extended to incorporate more complex hyper-elastic material models, including the Mooney–Rivlin, Odgen, and other advanced hyper-elastic models proposed by Davidson et al. [35] and Li et al. [36].

When not considering phase fields, the free energy density W of the neo-Hookean model is as follows:

$$W=\frac{1}{2}\mu \left(I_1-3-2\ln J\right)+\frac{1}{2}{K}_B{\left(\ln J\right)}^2$$

(4)

where $I_1=tr\left(F{F}^T\right),J=\det \left(F\right)$, $\mu$ is initial shear modulus, and $K_B$ is the initial bulk modulus.

Taking into account the differences in fracture behavior under tension and compression, as well as the variations in the decomposition methods of the free energy under finite and small deformations, the free energy $W$ is as follows:

$$W=W_1+W_2=\frac{1}{2}\mu {\displaystyle \sum _{i=1}^3\left({\lambda }_i^2-1-2\ln {\lambda }_i\right)}+\frac{1}{2}{K}_B{\left(\ln J\right)}^2$$

(5)

where ${\lambda }_i$ is the $F$ principal stretch, $J={\lambda }_1{\lambda }_2{\lambda }_3$, and $I_1={\displaystyle \sum _{i=1}^3{\lambda }_i^2}$. Here, $W$ is decomposed into $W_1=\frac{1}{2}\mu {\displaystyle \sum _{i=1}^3\left({\lambda }_i^2-1-2\ln {\lambda }_i\right)}$ and $W_1$ is a linear function; $W_2=\frac{1}{2}{K}_B{\left(\ln J\right)}^2$ and $W_2$ is a nonlinear function.

2.2.1. Decomposing Energy into Tension and Compression Components

Inspired by Miehe et al. [4] and Borden et al. [8], Tang and Zhang et al. [19] decomposed the energy into tension and compression components. The form of the free energy $W$ coupled with the phase field $p$ for simulating fracture is given by:

$$G_{rub}=\left[\left(1-\kappa \right)p^2+\kappa \right]W^{+}+W^{-}$$

(6)

$\kappa$ is introduced for improving numerical stability; its value is much smaller than 1.

$W^{+}$ and $W^{-}$ are defined by:

$$W^{+}=W\left({\lambda }_i^{+}, J^{+}\right), {\lambda }_i^{+}=\{\begin{array}{l}{\lambda }_i, {\lambda }_i>1\\ 1, {\lambda }_i\le 1\end{array}, J^{+}=\{\begin{array}{l}J, J>1\\ 1, J\le 1\end{array}$$

(7)

$$W^{-}=W\left({\lambda }_i^{-}, J^{-}\right), {\lambda }_i^{-}=\{\begin{array}{l}1, {\lambda }_i>1\\ {\lambda }_i, {\lambda }_i\le 1\end{array}, J^{-}=\{\begin{array}{l}1, J>1\\ J, J\le 1\end{array}$$

(8)

2.2.2. Decomposing Energy into Isochoric and Volumetric Components

Swamynathan et al. [28] developed a model that permits the evolution of the phase-field in the presence of local tensile conditions while maintaining the compressive rigidity of the material. Building upon the invariant-based work of Hesch et al. [37], they split the energy function into isochoric and volumetric components.

$$W_0\left(\overline{I_1},\overline{I_2},J\right)=W_{iso}\left(\overline{I_1},\overline{I_2}\right)+W_{vol}\left(J\right)$$

(9)

Here, $\overline{I_1}$ and $\overline{I_2}$ are the first and the second invariant of the isochoric part of the right Cauchy–Green deformation tensor given by

$$\overline{I_1}=J^{-2/3}{I}_1, \overline{I_2}=J^{-4/3}{I}_2$$

(10)

where $I_1$ and $I_2$ are the first and second invariant of the right Cauchy–Green tensor. The undamaged energy density was divided into contributions of tension and compression

$$W\left(\overline{I_1},\overline{I_2}, J\right)=W_{iso}^{+}\left({\overline{I_1}}^{+}, {\overline{I_2}}^{+}\right)+W_{vol}^{+}\left(J^{+}\right)+W_{iso}^{-}\left({\overline{I_1}}^{-}, {\overline{I_2}}^{-}\right)+W_{vol}^{-}\left(J^{-}\right)$$

(11)

where positive and negative components of energies are determined using the Heaviside functions as

$$W_{iso}^{+}={\displaystyle \sum _{i}H\left(I_i-3\right)W_{iso}}\left(\overline{I_i}\right)$$

(12)

$$W_{iso}^{+}={\displaystyle \sum _i\left[1-H\left(I_i-3\right)\right]W_{iso}}\left(\overline{I_i}\right)$$

(13)

The ultimate expression of the energy function is presented as

$$G\left(F,p\right)=g\left(p\right)\left[W_{iso}^{+}\left({\overline{I_1}}^{+}, {\overline{I_2}}^{+}\right)+W_{vol}^{+}\left(J^{+}\right)\right]+W_{iso}^{-}\left({\overline{I_1}}^{-}, {\overline{I_2}}^{-}\right)+W_{vol}^{-}\left(J^{-}\right)$$

(14)

where the degradation function is defined as

$$g\left(p\right)={\left(1-p\right)}^2$$

(15)

2.3. Numerical Verification

To investigate the performance of the aforementioned model, Swamynathan et al. [28] conducted a study on the classical three-point bending test under finite deformation and plane strain conditions. Figure 2a illustrates the geometry and boundary conditions of the test. A monotonic increasing displacement load was applied at the center of the specimen’s upper boundary, while a vertical slit of appropriate length was introduced at the lower boundary. The mesh was divided using quadrilateral elements, as shown in Figure 2b. During the crack initiation stage of the test, the displacement $\mathsf{\Delta }u=1\times {10}^{-4}\mathrm{mm}$. As the crack propagation stage commenced, the displacement $\mathsf{\Delta }u$ decreased to 1 × 10⁻⁵ mm.

The evolution of the phase field during different stages of displacement loading is illustrated in Figure 3; the crack path results obtained from all three models were consistent. However, in the model proposed by Amor et al. [38], compressive damage was observed in the region where the specimen was loaded, while in the models proposed by Tang and Zhang et al. [19], as well as by Swamynathan et al. [28], this phenomenon was alleviated. It can be concluded that the latter two models accurately distinguished the local states of tension and compression. In this demonstration, a model derived from the research of Amor et al. is presented as a comparison. Note that, for clarity in illustrating the crack propagation, the regions within the continuum where $p>0.98$ have been removed.

2.4. Summary

This chapter begins by summarizing the research achievements of phase-field fracture models for nonlinear elastic materials under finite deformation conditions in the past five years. Two typical phase-field models are then presented: one decomposes the free energy into tension and compression components, while the other decomposes the free energy into isochoric and volumetric components. Subsequently, the simulation results of these two models are compared using a classical three-point bending test. The advantage of the approach that decomposes the free energy into tension and compression components lies in its ability to be extended to tension–compression asymmetric materials, while the approach that decomposes the free energy into isochoric and volumetric components exhibits a higher level of fit with the existing literature.

Currently, for complex geometric shapes and material models such as visco-plasticity and elasto-plasticity, further improvement in mesh sensitivity and parametric studies is required. The ES-FEM method, which combines the finite element method with the idea of meshless methods, offers advantages such as precise results and a lack of susceptibility to element distortion. Combining the phase-field method with ES-FEM for finite deformation fracture modeling shows promising potential for development.

3. Fracture Study of Nonlinear Elastic Materials with Tension–Compression Asymmetry

Many natural or artificially synthesized materials exhibit tension–compression asymmetry, showing different behaviors under tensile and compressive loads. These tension–compression asymmetric materials encompass natural materials such as brain tissue and artificially synthesized materials such as polyester (polyethylene terephthalate), exhibiting a significant disparity between their tensile and compressive elastic moduli.

Previous research has established theoretical models to describe the mechanical properties of tension–compression asymmetric materials. In the framework of small deformations, Ambartsumyan’s pioneering work [39,40] on the theory of bimodulus elasticity is summarized in monographs. Du and Guo [41], along with Du et al. [42], formulated a series of variational principles and related definitions for bimodulus materials based on the constitutive relationship proposed by Ambartsumyan [39,40].

3.1. Research Progress

The contributions of the past five years in the field of phase-field fracture of nonlinear elastic tension–compression asymmetric materials are shown in chronological order. Li et al. [43] investigated the application of a dynamic gradient damage model as a phase-field model in dynamic fracture problems. The phase-field model offers a distinct advantage by providing a unified framework for addressing 2D and 3D crack evolution problems. A comparative analysis was conducted between two distinct damage constitutive laws and tension–compression asymmetry formulas. However, compared with traditional methods based on crack sharpness descriptions, it may incur higher computational cost. Tang et al. [19] constructed a Lagrangian equation and derived a coupled equation group describing deformation and phase-field evolution through crack separation energy and fracture energy in the phase-field approximation. The model exhibits good robustness and efficiency and its validity is confirmed through multiple examples, aligning well with experimental observations. This model offers the benefit of effectively managing tension and compression asymmetry within finite deformation, rendering it suitable for simulating fractures in rubber-like materials characterized by such asymmetry.

Shahba et al. [44] developed a FES phase-field simulation framework for modeling the fracture processes in anisotropic elastic materials exhibiting tension–compression asymmetry. The article also discussed the issue of numerical instability that adversely affects computational fracture simulations. To address this, a solution based on viscous stabilization was proposed to effectively overcome the convergence problem in the nonlinear finite element solver. Zhang et al. [45] presented a method for modeling nonlinear elastic material with tension–compression asymmetry. The method avoids the non-zero stress issue by decomposing the energy into stretch and compression parts and simulating the asymmetry by directly setting the ratio of shear and bulk modulus under tension and compression. Cheng et al. [46] introduced a wavelet-enriched adaptive FE model for solving a coupled crystal plasticity-phase field model, allowing for the simulation of crack propagation within poly-crystalline micro-structures. The advantage of this model is no need to assume a prior crack path. Their research indicates that the stored elastic strain energy due to material anisotropy and tension–compression asymmetry, along with the defect energy from the slip systems, are the main driving forces for crack propagation under finite deformation.

A model coupling crystal plasticity and phase-field modeling was developed by Tu et al. [47] for analyzing crack initiation and propagation in poly-phase–poly-crystalline micro-structures of 7000 aluminum alloy. You et al. [48] proposed a plastic damage-coupled phase-field framework that can describe various fracture modes of quasi-brittle materials. The framework allows for the incorporation of different tension crack initiation criteria and constitutive relations, enabling an accurate description of tension–compression asymmetry in mechanical response. However, the comprehensive yield surface smoothness of the model presents some issues and capturing the integration of tensile and compressive failure remains a challenging theoretical problem. Li et al. [49] demonstrated the use of two common tension–compression asymmetry phase-field fracture models to solve plane stress problems. The correctness of the proposed plane stress formula was verified through benchmark testing examples. Pathrikar et al. [50] proposed a gauge field theory model for quantifying and evolving micro-crack defects in solid deformation. The model considers kinematic aspects of deformation and damage, describing various characteristics of brittle damage, such as tension–compression asymmetry, stiffness degradation, and energy functionals, including crack contributions.

A coupled phase-field model was proposed by Hao et al. [51] to simulate fracture under high-speed impact. The model incorporates tension–compression asymmetry and introduces a new historical variable to enforce the irreversible condition of crack propagation. Ziaei-Rad et al. [52] proposed a method that decomposes the constitutive relation into crack driving and persistent parts, specifically designed for materials exhibiting anisotropic/orthotropic behavior in phase-field fracture to explain the tension–compression asymmetry. This model is applicable to any anisotropic elastic material in a 3D setting. Building upon this foundation, two existing tension–compression asymmetric models, namely the volume-deviator model and the no-tension model, were expanded to incorporate anisotropic properties. The extended models accurately capture anisotropic constitutive behavior and tension–compression asymmetry in crack response, exhibiting qualitative agreement with the anticipated behavior of orthogonal materials.

Various researchers have developed phase-field simulation frameworks and models to study fracture processes in nonlinear elastic materials with tension–compression asymmetry. These studies contribute to a better understanding of crack propagation and fracture behavior in complex materials.

3.2. Tension–Compression Asymmetry Phase-Field Model

3.2.1. Decomposing Energy into Tension and Compression Components

Zhang et al. [45] extended the approach of energy decomposition into tension and compression components to tension–compression asymmetric materials. Consistent with Section 2.2 above, body $B$ is surrounded by a curved surface $\partial B$ (see Figure 1). As the object undergoes deformation, point $X$ in the unaltered reference configuration is transformed to a point $x$ in the present deformed configuration through the deformation gradient. The right Cauchy–Green tensor can be defined as

$$C=F^{T}F=Q\mathsf{\Lambda }{Q}^T$$

(16)

where $Q$ consists of the orthogonal eigenvectors of $C$ and $\mathsf{\Lambda }=diag\left({\lambda }_1^2,{\lambda }_2^2,{\lambda }_3^2\right)$ is a diagonal matrix of principal right Cauchy–Green tensor and ${\lambda }_i,i=1\dots 3$ are the principal stretches. Suppose $a_m{\mathsf{\mid }}_{m=1\dots 3}$ are eigenvalues of matrix $C$, then the following equation can be proven

$$a_{m}dC={\mathit{a}}_m\otimes {\mathit{a}}_m$$

(17)

where ${\mathit{a}}_m$ is the eigenvector corresponding to the eigenvalue $a_m$. As we know, $a_m={\lambda }_m^2$ is the eigenvalue of matrix $C$.

The neo-Hookean model serves as a fundamental basis and is extended to the Ogden model. To differentiate between tension and compression, the corresponding shear moduli are defined as

$$\mu \left({\lambda }_i\right)=\{\begin{array}{l}{\mu }^{+}, {\lambda }_i>1\\ {\mu }^{-}, {\lambda }_i\le 1\end{array}$$

(18)

the corresponding bulk modulus is then defined as

$$K\left(J\right)=\{\begin{array}{l}{K}^{+}, J>1\\ K^{-}, J\le 1\end{array}$$

(19)

where $J=\det (F)={\lambda }_1{\lambda }_2{\lambda }_3$. Through the adoption of these shear and bulk moduli formulations, the strain energy of the initial neo-Hookean model, which incorporates tension–compression asymmetry, can be reformulated as follows:

$$W={\displaystyle \sum _{i=1}^3\mu \left({\lambda }_i\right)}\left[\frac{1}{2}\left({\lambda }_i^2-1\right)-\ln {\lambda }_i\right]+\frac{1}{2}K\left(J\right){\left(\ln J\right)}^2$$

(20)

3.2.2. Decomposing Energy into Elasticity, Plasticity, and Surface Energy Parts

With a focus on the influence of tension and shear, You et al. [48] introduced the definition of crack surface density per unit volume within the fracture process zone as follows:

$${\gamma }_1\left(p,\nabla p\right)=\frac{1}{2l_{c1}}{p}_1^2+\frac{l_{c1}}{2}{\left(\nabla p_1\right)}^2$$

(21)

$${\gamma }_2\left(p,\nabla p\right)=\frac{1}{2l_{c2}}{p}_2^2+\frac{l_{c2}}{2}{\left(\nabla p_2\right)}^2$$

(22)

where $l_{c1}$ and $p_1$ are the length scale parameter and the phase-field damage variable of regularization under tension, while $l_{c2}$ and $p_2$ are the length scale parameter and the phase-field damage variable of regularization under shear.

Considering the Griffith-type critical energy release rates $g_{c1}$ for tension and $g_{c2}$ for shear, the system’s Helmholtz free energy density $G_H$ is introduced within the thermo-dynamic framework. $G_H$ is broken down into four parts: elasticity, plasticity, and two surface energies:

$$G_H\left({\epsilon }^e,{\epsilon }^p,\alpha ,p_1,p_2\right)=G^e\left({\epsilon }^e,p_1\right)+G^p\left({\epsilon }^p,\alpha ,p_2\right)+G_1\left(p_1,\nabla p_1\right)+G_2\left(p_2,\nabla p_2\right)$$

(23)

where

$$G_1\left(p_1,\nabla p_1\right)=g_{c1}{\gamma }_1\left(p_1,\nabla p_1\right)$$

(24)

$$G_2\left(p_2,\nabla p_2\right)=g_{c2}{\gamma }_2\left(p_2,\nabla p_2\right)$$

(25)

$$G^e\left({\epsilon }^e,p_1\right)=\frac{1}{2}\left[{\left(1-p_1\right)}^2+\kappa \right]{\epsilon }^e:{ℂ}^0:{\epsilon }^e$$

(26)

$$G^p\left({\epsilon }^p,\alpha ,p_2\right)=g\left(p_2\right)\left(\frac{1}{2}{\epsilon }^p:\mathscr{H}:{\epsilon }^p+\frac{1}{2}h{\alpha }^2\right)$$

(27)

In Equation (22), ${\epsilon }^e$ and ${\epsilon }^p$ are elastic and plastic strains, respectively. ${ℂ}^0$ represents the elasticity tensor of the material in its undamaged state. To ensure the well-posedness of the system in regions with partially broken domains, a small positive parameter $\kappa$ is introduced. A fourth-order kinematic hardening modulus tensor, denoted as $\mathscr{H}$, is positively definite, and $h$ is a non-negative parameter related to the isotropic hardening of the material. In this case, the researchers did not explore various degradation functions and the degradation function follows a standard quadratic form.

$$g\left(p_2\right)={\left(1-p_2\right)}^2$$

(28)

3.3. Numerical Verification

Numerical examples are employed in this section to demonstrate the fracture modeling capability of the tension–compression asymmetrical phase-field model for nonlinear elastic materials. Initially, this section presents the research findings of Zhang et al. [45]. In the simulations, ${\mu }^{+}$ and $K^{+}$ are kept unchanged, while ${\mu }^{-}$ is varied; the ratio of bulk moduli $K^{+}/K^{-}$ and the ratio of shear moduli ${\mu }^{+}/{\mu }^{-}$ are kept equal.

To investigate the three-point bending behavior of a solid specimen with a vertical notch under plane strain conditions, the model depicted in Figure 4 is established. The figure illustrates the boundary conditions and the discrete FE mesh employed in the model. A finite element model is established, where the solid is discretized by triangular elements.

Figure 5 illustrates the contours of the phase-field $p$ on the deformed configuration with two modulus ratios ${\mu }^{+}/{\mu }^{-}=5.0,0.5$. The nucleation of new cracks occurs at imposed displacements of 0.97 and 0.61, respectively, leading to the subsequent propagation of the pre-existing crack.

You et al. [48] also conducted a three-point bending test; Figure 6 shows the model and the finite element mesh partitioning results of the simulation experiment. Figure 7 compares the predicted crack evolution processes for $g_{c1}=0.045\mathrm{N}/\mathrm{mm}$ and $g_{c1}=0.029\mathrm{N}/\mathrm{mm}$, respectively. The crack initiates from the notch tip and propagates towards the upper border of the block. However, when $g_{c1}=0.029\mathrm{N}/\mathrm{mm}$, the crack propagation speed is faster.

3.4. Summary

This section presents a summary of the accomplishments of previous researchers in constructing models for nonlinear elastic tension–compression asymmetric materials. Particularly, the methods of Zhang et al. [45] and You et al. [48] are demonstrated. Zhang et al.’s strategy demonstrates excellent numerical stability through the direct setting of shear modulus and bulk modulus ratios, enabling the simulation of tension–compression asymmetry. Conversely, You et al.’s model effectively captures pressure sensitivity and unilateral effects. Nevertheless, a trade-off has been made regarding the smoothness of the combined yield surfaces. Researchers have proposed numerous damage constitutive laws and tension–compression asymmetry formulas, along with various methods to simulate the fracture behavior of materials displaying such asymmetric characteristics. However, certain theoretical issues remain, such as ensuring the smoothness of the integrated yield surface and capturing the integrated failure behavior during both tensile and compressive loading. Furthermore, the high computational cost associated with the phase-field model may impede its application in large-scale problems. Currently, researchers are actively addressing these challenges to enhance the accuracy and efficiency of the phase-field model in simulating fracture propagation under finite deformation in tension–compression asymmetric materials.

4. Fracture Behavior Research of Hydrogels

Water-containing soft solids extensively exist in synthetic materials such as hydrogels, as well as in natural materials such as elastic proteins and materials designed to mimic elastic proteins. With advancements in biomedical and chemical engineering, these materials are extensively used in various applications, including the development of scaffolds and flexible electronic devices. Notably, these materials possess the ability to maintain their structural integrity and withstand macroscopic cracks. It is noteworthy that these water-containing soft solids exhibit softer mechanical properties compared with their dehydrated counterparts, which are typically harder. This distinction in mechanical behavior can be utilized to achieve complex functionalities, such as camouflage.

4.1. Research Progress

From a micro-mechanical perspective, water-containing soft solids can be viewed as polymer networks with a high degree of hydrophilicity and water-swelling, where hydration plays a critical role in determining their elastic properties. When the notched specimen is subjected to subcritical stretching and maintained at a constant state, ‘delayed fracture’ [53] occurs, where the crack propagates following an incubation period of adequate duration, attributed to the diffusion of water within the polymer network. The exceptional toughness of double-network hydrogels has gained considerable attention and qualitative crack models have been developed to elucidate its underlying mechanisms. These models aim to explain the superior toughness exhibited by double-network hydrogels compared with conventional hydrogels [54,55,56,57].

From the perspective of macroscopic fracture mechanics, extensive research has been conducted on the fracture behavior of these hydrated materials. In this context, many numerical methods [58] have been proposed to investigate how crack propagation initiates by assuming the existence of pre-existing cracks and analyzing the crack tip region. This approach also offers a means of measuring material toughness, allowing for the evaluation of the material’s brittleness or ductility [59,60,61,62,63]. Classical fracture mechanics theories have been extended to hydrogels and even utilized to examine fatigue fracture [64,65]. Zhang et al. [66] developed a model that couples the cohesive zone and the Mullins effect to quantitatively predict the fracture energy of soft and tough hydrogels, as well as the fracture field near the crack tip. Long et al. [67] presented an analytical solution describing the deformation and stress field in the vicinity of a static crack tip in a soft elastic material. Guo et al. [68] studied the stress and deformation fields in the vicinity of a crack tip within a self-healing hydrogel specimen featuring a unilateral notch. They employed a model that incorporates physical crosslink fracture dynamics to describe the behavior of finite deformations in nonlinear visco-elastic materials. Yu et al. [69] developed asymptotic and finite element analysis methods to analyze the crack tip field of stable crack propagation in polymer gels based on linear poroelastic theory.

However, the aforementioned methods are unable to predict the initiation time and location of cracks and they may face challenges when predicting crack initiation and propagation. To partially overcome these challenges, a potential solution is to introduce a diffuse crack model that relies on the fracture phase field.

Böger et al. [70] proposed and implemented a variational framework, employing the finite element method, to establish a phase-field fracture model for gels. Their numerical simulations unveiled the occurrence of crack initiation and subsequent propagation throughout the dehydration process in gels. Mao et al. [71] formulated a phase-field model by adopting Gurtin’s [72,73] original virtual method. The constitutive equation for gel fracture contains two novel physical components: (1) the entropy and internal energy caused by the stretching of the Kuhn segments and the polymer cross-linked network structure and (2) the damage and failure of the polymer network attributed to the alteration in internal energy. Mao et al. [71] have illustrated the successful application of the phase-field approach in simulating the fracture behavior of double-notched specimens subjected to uniaxial tension. Although a fully coupled diffusion–deformation fracture framework holds promise for simulating real-world problems, it introduces more parameters that are not easily calibrated through experiments, particularly those related to water diffusion. When the deformation time scale significantly outweighs the diffusion time scale, neglecting water diffusion becomes plausible [54,55,66]. This assumption facilitates the calibration of material parameters involved in experimental studies.

Böger et al. [74] introduced a phase-field model for crack propagation in water gels, integrated it into a variational framework, and implemented it using an operator splitting algorithm. They conducted creep tests driven by diffusion to investigate the slow mass transport phenomena in water gels, specifically crack initiation and evolution induced by drying. Zhang et al. [75] proposed a graph finite element phase-field hybrid model to simulate the fracture behavior of a water gel-based curved shell, abstracting the water gel in the biomedical field as a shell and effectively improving the computational efficiency. Zheng et al. [76] proposed a diffusion fracture-based rate-independent variational principle phase-field method, which can effectively simulate the fracture process of super elastic materials and water gels under different boundary conditions, and verified the robustness of the method. Subsequently, they integrated the phase-field evolution equation with the analogy of heat transfer and diffusion laws to establish a theoretical model of the temperature-sensitive water gel fracture simulation phase-field method that considers diffusion-coupled large deformation. This model eliminates the need for pre-defined cracks and can simulate cracks of various types of temperature-sensitive water gels under different boundary conditions. Liu et al. [77] established an anisotropic model for porous-visco-hyper-elastic damage based on the theory of porous media (TPM). This model aimed to analyze the relationship between visco-hyper-elasticity, time-dependent behavior of fluid transport coupling, and fracture behavior of water gel composites. During the numerical verification process, the aging mechanism of the anisotropic water gel composite fracture behavior was discussed.

Researchers have made significant advancements in the field of phase-field fracture models for gels. These innovative approaches hold great promise for understanding and analyzing the fracture behavior of gels in real-world applications.

4.2. The Material and Phase-Field Models for Hydrogels

Numerous models have been developed to analyze the finite deformation behavior of hydrogels [73,78,79,80,81]. The theoretical formulations typically reference the free energy density of hydrogels to the state of the network that is dry and devoid of solvent. However, in reality, soft solids such as hydrogels contain a significant number of water molecules. Therefore, the reference state is defined as the initially hydrated state rather than the dehydrated state and is expressed as the reference configuration in its undeformed state.

4.2.1. Decomposing Energy into Tension and Compression Components

Zhang et al. [82] employed the approach of decomposing energy into tension and compression components to conduct fracture simulations of hydrogel materials. The initial polymer volume fraction is denoted by ${\varphi }_0$ and its corresponding free energy density is defined by the following equation:

$$\frac{G}{kT/\upsilon }=\frac{1}{2}N\upsilon {\varphi }_0\left(\overline{I}-3\right)+\frac{k_B}{kT/\upsilon }\left[J^2-1+2{\varphi }_0\left(1-J\right)\right]/{\varphi }_0+\left(J-{\varphi }_0\right)\ln \left(1-\frac{{\varphi }_0}{J}\right)+\chi {\varphi }_0^2\left(1-\frac{1}{J}\right)-\frac{{\mu }_0}{kT}\left(J-1\right)-\frac{\mathsf{\Delta }\mu }{kT}\left(J-{\varphi }_0\right)$$

(29)

with the first invariant $\overline{I}=J^{-2/3}I\left(I=tr\left(b\right),J=det\left(F\right)\right)$.

In the aforementioned equation, $N$ represents the number of chains per unit volume, $\upsilon$ is the volume per solvent molecule, $\chi$ is the Flory–Huggins parameter that characterizes the enthalpy of mixing between polymer and solvent molecules, $k$ is the Boltzmann constant, $T$ is the temperature, and $k_B$ is the pseudo bulk modulus. Given that the swollen state is adopted as the undeformed reference configuration, the chemical potential at this state, denoted as ${\mu }_0$, can be represented by the following expression:

$$\frac{{\mu }_0}{kT}=\frac{2k_B}{kT/\upsilon }\left(1-{\varphi }_0\right)/{\varphi }_0+\chi {\varphi }_0^2+{\varphi }_0+\ln \left(1-{\varphi }_0\right)$$

(30)

$\mathsf{\Delta }\mu /kT$ represents the change in chemical potential with respect to the reference configuration. The initial (undeformed) state is denoted by $\mathsf{\Delta }\mu /{\mu }_0=0$. Constrained deswelling of the gel is represented by $\mathsf{\Delta }\mu /{\mu }_0>0$, while $\mathsf{\Delta }\mu /{\mu }_0<0$ represents swelling.

To ensure resistance in the compressed state, the free energy density $G$ is subsequently decomposed into isochoric and volumetric components.

$$G\left(I,J\right)=G_d\left(\overline{I}\right)+G_{\nu }\left(J\right)$$

(31)

where

$$\frac{G_d}{kT/\upsilon }=\frac{1}{2}N\upsilon {\varphi }_0\left(\overline{I}-3\right)$$

(32)

$$\frac{G_{\nu }}{kT/\upsilon }=\frac{k_B}{kT/\upsilon }\left[J^2-1+2{\varphi }_0\left(1-J\right)\right]/{\varphi }_0+\left(J-{\varphi }_0\right)\ln \left(1-\frac{{\varphi }_0}{J}\right)+\chi {\varphi }_0^2\left(1-\frac{1}{J}\right)-\frac{{\mu }_0}{kT}\left(J-1\right)-\frac{\mathsf{\Delta }\mu }{kT}\left(J-{\varphi }_0\right)$$

(33)

Subsequently, a phase-field methodology is proposed to characterize the fracture process within the critical region in hydrogels. This method aims to quantify the reduction in stiffness of the material within the failure zone. Inspiration for the ensuing expression for the hydrogel’s free energy, which is coupled with the phase-field $p$, was drawn from the works of Miehe et al. [4] and Borden et al. [8], in the realm of linear elasticity for small deformations.

$$G_{gel}=\left(\overline{I},J,p\right)=\left[\left(1-K\right)p^2+K\right]\left[G_d\left(\overline{I}\right)+G_{\nu }^{+}\left(J\right)\right]+G_{\nu }^{-}\left(J\right)$$

(34)

where

$$G_{\nu }^{+}\left(J\right)=G_{\nu }{H}^{+}$$

(35)

$$H^{+}=\{\begin{array}{l}1, \frac{\partial G}{\partial J}\ge 0\\ 0, \frac{\partial G}{\partial J}<0\end{array}, H^{-}=\{\begin{array}{l}0, \frac{\partial G}{\partial J}\ge 0\\ 1, \frac{\partial G}{\partial J}<0\end{array}$$

(36)

The sign of $H^{\pm }$ is dependent on whether the material is subjected to hydrostatic pressure or compression. Essentially, the superscript of $H^{\pm }$ is indicative of the sign of the hydrostatic (or mean) Cauchy stress:

$${\sigma }_m=\frac{1}{3}tr\left(\sigma \right)=\frac{\partial G}{\partial J}$$

(37)

4.2.2. Decomposing Energy into Stretch and Mixing Components

Compared with the method of Zhang et al., Zheng et al. [83] decomposed undamaged free energy density $G$ into the part of stretch and mixing. Note that the model is established based on temperature-sensitive hydrogel.

$$G\left(F,p\right)=g\left(p\right)G_0\left(F\right)$$

(38)

where

$$G_0\left(F\right)=\frac{1}{2}N{k}_{B}T\left[F^{T}F-3-2\ln \left(J\right)\right]+k_{B}T\left[D\ln \left(\frac{\nu D}{1+\nu D}\right)+\frac{\chi D}{1+\nu D}\right]-\mu D$$

(39)

where $N$ is the referential chain density, $k_B$ is the Boltzman constant, $T$ is the absolute temperature, $\nu$ is the nominal volume of a solvent molecule, $D$ is the concentration of the solvent in the hydrogel, $\chi$ is the interaction parameter, and $\mu$ is the chemical potential in the solvent.

To incorporate the phase field model with the large deformation behavior of hydrogels, the potential energy of a hydrogel body is formulated as

$$W_{\mathrm{int}}=E\left(u,p\right)+W\left(p\right)$$

(40)

where

$$E\left(u,p\right)={\displaystyle {\int }_{\mathsf{\Omega }}G\left(F,p\right)}dV$$

(41)

$$W\left(p\right)={\displaystyle {\int }_{\mathsf{\Omega }}{g}_c\gamma \left(p,\nabla p\right)}dV$$

(42)

where $W$ (fracture energy) represents the cumulative fracture surface area multiplied by $g_c$, the critical fracture energy, and

$$\gamma \left(p,\nabla p\right)=\frac{1}{2l_c}{p}^2+\frac{l_c}{2}{\left(\nabla p\right)}^2$$

(43)

4.3. Numerical Verification

4.3.1. Decomposing Energy into Tension and Compression Components

A finite element model was constructed, as illustrated in Figure 8a, where the ${\varphi }_0=0.15$. Figure 8b,c displays the deformation patterns at loading point displacements of $d/H=1.0$ and $d/H=1.5$, respectively.

As the loading point displacement $d$ increases, the crack phase field undergoes evolution. At a loading displacement of $d/H=1.5$, Figure 8d exhibits an isosurface map of the crack phase field in the initial configuration. The contour analysis of the crack phase field demonstrates the initiation of the crack at the center of the bottom surface, followed by its linear propagation towards the top surface. The simulation results were experimentally validated, showing consistency with the experimental results, which indicates the accuracy of the simulation. However, due to mesh deformation, the crack propagation was limited to approximately half the height of the specimen.

4.3.2. Decomposing Energy into Stretch and Mixing Components

The notched hydrogel was investigated by Zheng et al. [83]. Figure 9 depicts the geometry of the sample, where the initial crack length measures 5 mm. In the experimental configuration known as pure shear, the bottom side of the rectangular specimen is fixed, while the top side experiences displacement while retaining its original length.

The deformation patterns of hydrogel specimens in different crack states are illustrated in Figure 10, demonstrating the evolution of cracks until final fracture. Figure 10a depicts the region surrounding the initial notch tip of the hydrogel specimens, showcasing the localized crack initiation. In Figure 10b, the fracture propagation is illustrated, where the load reaches its peak magnitude. Subsequently, there is a rapid decline in load until the complete fracture of the specimen. Figure 10c depicts the subsequent crack extension until the specimen separates into two parts.

4.4. Summary

Consistent with their previous work on tension–compression symmetric [19] and asymmetric materials [45], Zhang et al.’s [82] approach decomposed the energy into tensile and compressive parts, ensuring both fracture capture and numerical stability. Zheng et al. [83] established a phase field fracture model for temperature-sensitive hydrogels, which can be extended to other hyper-elastic models.

In recent years, the phase-field method has emerged as a prominent numerical approach for studying fracture models of hydrogels. Based on the continuity description of energy functional and phase-field variable, the phase-field method enables the simulation of mechanical and thermo-dynamic properties of hydrogels, including initial damage, fracture process, and fracture morphology. Previous studies have demonstrated the applicability of the phase-field method in analyzing various fracture behaviors of hydrogels, including those of temperature-sensitive and chemically reactive hydrogels. Moreover, the phase-field method can be combined with other models and experimental data for validation and application, such as in conjunction with continuum mechanics models and finite element methods to predict the mechanical response and fracture behavior of hydrogels.

While progress has been made in studying fracture models of hydrogels using the phase-field method, several challenges still need to be addressed. First, existing phase-field models often demand significant computing resources and time, necessitating further optimization and enhancement of their computational efficiency. Second, the phase-field method requires more accurate physical and material parameters to describe the mechanical and thermo-dynamic properties of hydrogels. Therefore, additional experimental research and parameter optimization is necessary. Lastly, the phase-field method struggles to effectively simulate the multi-scale and multi-physical process characteristics of hydrogels, underscoring the necessity for the advancement of phase-field models capable of capturing such multi-scale and multi-physical phenomena. Resolving these challenges will contribute to a deeper understanding of the fracture behavior of hydrogels and offer valuable insights for the design and application of hydrogels in biomedical and engineering fields.

5. Summary and Future Directions

Over the past decade, the phase-field method has gained widespread recognition in studying the fracture of nonlinear elastic materials under finite deformation. This method demonstrates its capability to capture complex fracture patterns across various material models and structures and provides insights into different scenarios under diverse loading conditions. Consequently, it offers a viable solution for addressing challenges beyond the scope of commonly used models. However, it necessitates substantial efforts in mathematical formulas, numerical implementation, and model validation. This comprehensive review paper aims to elucidate critical aspects of fracture in nonlinear elastic materials under finite deformation and highlight its notable contributions.

In addition, this review introduces several formulations for the phase-field evolution equation of nonlinear elastic materials and some methods for analyzing such materials. Among these, a novel energy decomposition method that incorporates strain energy is highlighted, particularly emphasizing considering stretching and compression within the framework of finite deformation. This approach aims to predict crack initiation and propagation in the phase-field model, providing an invaluable resource for studying fracture behavior in nonlinear elastic materials. Furthermore, the review covers several numerical implementation methods, algorithms, and experimental studies.

Given the significant challenge associated with predicting the fracture of nonlinear elastic materials under finite deformations, utilizing the phase-field method in addressing this issue is currently limited. Therefore, there remains a pressing need for further research to comprehensively explore this area.

• The research presented in this review was conducted within the quasi-static phase field framework. However, there is currently a lack of comprehensive research using the methodology of phase field to investigate the physical mechanisms of dynamic crack propagation, especially the conditions leading to crack branching. In addition, accurately determining the location of the crack tip to calculate crack propagation velocity remains a challenge. Therefore, it is necessary to establish a phase-field model that incorporates dynamic loading with rate effects and rapid crack propagation.
• Fatigue crack formation and propagation contribute significantly to the failure of many engineering structures; predicting these phenomena continues to pose challenges in modeling and simulation. Some researchers have begun employing phase-field methods to study the fatigue crack in brittle or quasi-brittle materials. Hence, the utilization of the phase-field approach to explore fatigue crack problems in nonlinear elastic materials is a novel idea.
• While the phase-field models discussed in this review have demonstrated their capability to predict fracture strains of specimens under quasi-static loading using the finite element mesh, they face limitations in capturing crack propagation during large deformation, especially in the case of three-point bending. This phenomenon may be attributed to severe deformation of the meshes.