Phase-field Modeling of Brittle Fracture#

The variational formulation of the Griffith principle of (brittle) fracture [FM98] introduces a phase-field regularization of sharp cracks (\(\Gamma\)) via the crack surface density function \(\gamma(\phi,\nabla\phi)\) where the (damage) phase field \(0\leq\phi\leq1\) describes a smeared crack over a distance modulated by the length scale parameter \(l_0\). A general expression for \(\gamma(\phi,\nabla\phi)\) is:

\[ \gamma(\phi,\nabla\phi)=G_c \frac{1}{c_0l_0}[\alpha(\phi)+l_0^2|\nabla\phi|^2] \]

where \(\alpha(\phi)\) is the geometric crack function and \(c_0= 4\int^{1}_{0}\sqrt{\alpha(\phi)}d\phi\) a scaling factor. The crack surface energy \(\Psi_{\Gamma}\) can thus be calculated by integrating over the domain \(\Omega\):

\[ \Psi_{\Gamma}=\int_{\Gamma}G_c dA \approx \int_{\Omega}G_c \gamma(\phi;\nabla\phi)dV \]

where \(G_c\) [\(J/m^2\)] is the critical energy release rate, often related to the material fracture toughness. With the above regularization, the strain energy stored in the solid can be expressed as the (elastic) strain energy \(\Psi(\epsilon(\bm{u}))\) multiplied by an appropriate degradation function \(g(\phi)\). The total energy of the solid is then given by

\[ W(\bm{u},\phi)=\int_{\Omega}g(\phi)\Psi(\epsilon(\bm{u}))dV + \int_{\Omega}G_c \frac{1}{c_0l_0}[\alpha(\phi)+l_0^2|\nabla\phi|^2]dV \]

Minimization of \(W(\bm{u},\phi)\) with respect to \(\bm{u}\) and \(\phi\) provides solution to Griffith’s (brittle) fracture problem, providing appropriate boundary conditions are met and the crack irreversibility constraint is resolved.

\(W(\bm{u},\phi)\) can be regarded as the free energy of a gradient damage model if \(l_0\) is linked to the material tensile strength, \(\sigma_c\). It can be shown, in particular, that \(\sigma_c^2= \beta \frac{EGc}{l_0}\), where \(\beta\) depends on the specific model (e.g. \(\beta=3/8\) for AT1 and \(\beta=27/256\) for AT2)[TanneLB+18].

Damage driving force and crack irreversibility#

Voids growth and coalescence is inhibited by hydrostatic pressure. In other words, crack propagation does not release strain energy associated with negative volumetric strains. Therefore, suitable expressions for \(\Psi\) must prevent hydrostatic stress components from contributing to damage evolution. For linear elasticity, we consider one by [AMM09] based on the following deviatoric and volumetric split:

\[ \Psi^{+}= \begin{cases} \Psi_{dev} + \Psi_{vol} = \mu(\epsilon_{dev} : \epsilon_{dev}) + 0.5 \bulk tr^2(e) \quad \text{if} \quad \trace \epsilon \geq 0 \\ \Psi_{dev} = \mu(\epsilon_{dev} : \epsilon_{dev}) \qquad \qquad \qquad \qquad \qquad \text{if} \quad \trace \epsilon < 0 \end{cases} \]

where \(\mu\) and \(\bulk\) are the shear and bulk moduli, respectively.

To avoid crack healing during unloading, \(\Psi^{+}\) can be handled as a history variable evaluated at quadrature points [MHW10]. The irreversibility condition, formally expressed as \(\dot\phi\geq0\), is then fulfilled by setting

\[ \Psi^{+}=max(\Psi^{+}_{n},\Psi^{+}_{n+1}) \]

Weak form#

The strong form of the phase-field problem of brittle fracture can be generalized for linear elasticity (in the absence of body forces) as:

\[ \begin{cases} \nabla \cdot \sigma_{degr} = \bm0 \qquad \text{in} \;\Omega \\ \nabla \cdot \bm h(\phi) +H(\phi) = 0 \qquad \text{in} \;\Omega \end{cases} \]

with

(159)#\[ \bm{h}(\phi)= \frac{2l_0}{c_0}G_c\nabla\phi \qquad H(\phi)=-g'(\phi)\Psi^+ - \frac{G_c}{c_0l_0}\alpha'(\phi) \]

where \(\sigma_{degr}\) is the degraded Cauchy stress. Alongside boundary conditions for displacement (\(\bm{u} = \bar{u} \; \text{on} \; \partial\Omega_D\) and \(\bm \sigma_{degr} \cdot \bm{n} =\bar{\bm{t}} \quad \text{on} \; \partial\Omega_N\)) we have,

\[ \phi = 1 \quad \text{on} \; \Gamma \quad , \quad \nabla\phi \cdot \bm n=0 \quad \text{on} \; \partial\Omega_N \]

The weak form of (159) can be derived and stated as finding \(\bm u \in \mathcal{V} \) and \(\phi \in \mathcal{Q}\) such that:

(160)#\[ \begin{cases} \int_{\Omega}(\sigma_{degr}:\nabla v)dV = 0 \qquad \forall v \in \mathcal{V} \\ \int_{\Omega}(\bm h \cdot\nabla q - Hq)dV = 0 \qquad \forall q \in \mathcal{Q} \end{cases} \]

Note

The current implementation uses a single field of components \(u_x\), \(u_y\), \(u_z\), and \(\phi\).

Note

To further prevent interpenetrating cracking surfaces, damage should only affect the deviatoric part of the stress under compressive strains, that is:

\[ \begin{cases} \sigma_{degr} = g(\phi)\sigma \quad \trace \epsilon \geq 0 \\ \sigma_{degr} = \sigma_v + g(\phi)\sigma_{dev} \quad \trace \epsilon < 0 \end{cases} \]

However, the so-called (non-variationally consistent) hybrid formulation in [AGDL15], whereby \(\sigma_{degr} = g(\phi)\sigma\), has gained popularity and shown to perform well in tensile loading conditions.

AT1 and AT2 models#

Quadratic, cubic, and quartic expressions for \(g(\phi)\) have been proposed in the literature, as reviewed in [KSchluterMuller15]. The selected function should decrease monotonically between \(g(0)=1\) and \(g(1)=0\), and be \(g'(1)=0\). Here, we use the widely adopted expression

(161)#\[ g(\phi) = (1 - \eta)(1 - \phi)^2 + \eta \]

where \(\eta<<1\) prevents numerical instabilities associated with the otherwise complete loss of stiffness for \(\phi=1\).

Unlike that of \(g(\phi)\), the choice of \(\alpha(\phi)\) significantly influences predictions of damage evolution and overall mechanical response. Following Ambrosio and Tortorelli, the two standard models have been termed AT1, in which \(\alpha(\phi)=\phi\), and AT2 model, in which \(\alpha(\phi)=\phi^2\). The two yields \(c_0=2\) and \(c_0=8/3\), respectively.

In AT2, damage starts accumulating as load is applied and the initial response is non-linear for any material. In contrast, AT1 introduces a threshold for \(\Psi^{+}\) below which damage does not occur [TanneLB+18]. The presence of an elastic stage makes AT1 appear most suitable to model the mechanical behaviour of materials in which damage is delayed by inelastic phenomena, such as in elastoplastic and viscoelastic materials.

Note

Adopting AT2 with degradation function (161) guarantees \(0\leq\phi\leq1\) [MWH10]. For other combinations, SNES’s semi-smooth solver for variational inequalities can be used to bound \(\phi\), e.g. [AKBP21].

AT1: Weak form and linearization (Not yet implemented)#

Substituting AT1 characteristic functions into (160) gives the weak form of the AT1 model with linear elasticity. Find \(\bm u \in \mathcal{V} \) and \(\phi \in \mathcal{Q}\) such that

\[ \begin{cases} \int_{\Omega}(\sigma_{degr}:\nabla v)dV = \bm0 \qquad \forall v \in \mathcal{V} \\ \int_{\Omega}[ \frac{3}{4}l_0G_c\nabla p \cdot \nabla q - (2(1 - \eta)(1 - \phi)\Psi^{+} - \frac{3G_c}{8l_0})q]dV = 0 \qquad \forall q \in \mathcal{Q} \end{cases} \]

Straightforward linearization of the relevant terms gives:

\[ d\sigma_{degr} = \frac{\partial\sigma_{degr}}{\partial\epsilon}+\frac{\partial{\sigma}}{\partial\phi}=[(1-\eta)(1 - \phi)^2 +\eta]d\sigma - 2(1-\eta)(1-\phi)d\phi\sigma \]

\[ dh=\frac{\partial h}{\partial\epsilon}+\frac{\partial h}{\partial\phi}= \frac{3}{4}l_0G_c\nabla d\phi \]

\[ dH=\frac{\partial H}{\partial\epsilon}+\frac{\partial H}{\partial\phi}= 2(1-\eta)(1-\phi)d\Psi^+ - 2(1-\eta)\Psi^ {+}d\phi \]

The term \(d\Psi^+\) is linked to crack irreversibility condition. If \(d \Psi^{+}_{n+1} < \Psi^{+}_{n}\) , \(d\Psi^{+}=0\), else:

\[ d\Psi^{+}= \begin{cases} d\Psi_{dev} + d\Psi_{vol} = d\sigma : \epsilon \quad \text{if} \quad \trace \epsilon \geq 0 \\ d\Psi_{dev} = d\sigma_{dev} : \epsilon_{dev} \qquad \text{if} \quad \trace \epsilon < 0 \end{cases} \]

AT2: Weak form and linearization#

Substituting AT2 characteristic functions into (160) gives the weak form of the AT2 model with linear elasticity. Find \(\bm u \in \mathcal{V} \) and \(\phi \in \mathcal{Q}\) such that

\[ \begin{cases} \int_{\Omega}(\sigma_{degr}:\nabla v)dV = \bm0 \qquad \forall v \in \mathcal{V} \\ \int_{\Omega}[l_0G_c\nabla\phi \cdot\nabla q - (2(1-\eta)(1 - \phi)\Psi^{+} - \frac{G_c}{l_0}\phi)q]dV = 0 \qquad \forall q \in \mathcal{Q} \end{cases} \]

Straightforward linearization of the relevant terms gives:

\[ d\sigma_{degr} = \frac{\partial\sigma_{degr}}{\partial\epsilon}+\frac{\partial{\sigma}}{\partial\phi}=[(1-\eta)(1 - \phi)^2 +\eta]d\sigma - 2(1-\eta)(1-\phi)d\phi\sigma \]

\[ dh=\frac{\partial h}{\partial\epsilon}+\frac{\partial h}{\partial\phi}= l_0G_c\nabla d\phi \]

\[ dH=\frac{\partial H}{\partial\epsilon}+\frac{\partial H}{\partial\phi}= 2(1-\eta)(1-\phi)d\Psi^+ - (2(1-\eta)\Psi^ {+} + \frac{G_c}{l_0})d\phi \]

As for AT1 above, if \(d \Psi^{+}_{n+1} < \Psi^{+}_{n}\) , \(d\Psi^{+}=0\), else:

Viscous regularization#

To prevent convergence loss at solution jumps associated with an abrupt crack propagation, a viscous regularization [MHW10] has been implemented by complementing the damage residual with the term \(\xi\dot{\phi}\) with (damage) viscosity \(\xi\).

The term thus appears in the expressions for \(H(\phi)\) and \(dH(\phi)\) above. For instance, in the AT2 case,

\[ H(\phi)= 2(1-\eta)(1-\phi)\Psi^+ - \frac{G_c}{l_0}\alpha'(\phi) - \xi\dot{\phi} \]

and

\[ dH= 2(1-\eta)(1-\phi)d\Psi^+ - (2(1-\eta)\Psi^ {+} + \frac{G_c}{l_0})d\phi -\xi * \mathrm{shift_v} \]

with \(\mathrm{shift_v}\) evaluated at time \(\mathrm{shift_v} = \frac{\partial\dot{\phi}}{\partial\phi}|_{\phi_n}\).