Phase-field Modeling of Fracture

Small strain

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(l_0, \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(l_0, \phi)\) is:

\[ \gamma(l_0, \phi) = \frac{1}{c_0 \, l_0} \left(\alpha(\phi)+l_0^2 \lvert \nabla\phi \rvert^2 \right), \]

where \(\alpha(\phi)\) is the geometric crack function and \(c_0= 4\int^{1}_{0}\sqrt{\alpha(\phi)}\diff \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 d \Gamma \approx \int_{\Omega} G_c \gamma(l_0, \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(\bm\varepsilon(\bm{u}))\) multiplied by an appropriate degradation function \(g(\phi)\). The total energy of the solid is then given by

(243)\[ W(\bm{u},\phi)=\int_{\Omega}g(\phi)\psi(\bm\varepsilon(\bm{u}))dv + \int_{\Omega} \frac{G_c}{c_0l_0} \left( \alpha(\phi)+l_0^2 \lvert \nabla\phi \rvert^2 \right) dv \]

Minimization of \(W(\bm{u},\phi)\) with respect to \(\bm{u}\) and \(\phi\) solves the Griffith’s (brittle) fracture problem but does not address crack irreversibility nor the effect of loading direction on damage evolution observed experimentally.

\(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, \(\bm{\sigma}_c\). It can be shown, in particular, that \(\bm{\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, consistent with cracking preserving the part of strain energy associated with compressive strain components. Expression for \(W(\bm{u},\phi)\) in (243) can then be recast as:

(244)\[ \bar{W}(\bm{u},\phi)=\int_{\Omega}(g(\phi)\psi^+(\bm\varepsilon(\bm{u})) + \psi^-(\bm\varepsilon(\bm{u})))dv + \int_{\Omega} \frac{G_c}{c_0l_0} \left( \alpha(\phi)+l_0^2 \lvert \nabla\phi \rvert^2 \right) dv \]

Several expressions for \(\psi^+\) can be found in the literature, including some function of tensile critical stresses or positive principal stretches [MHW10]. For linear elasticity, we consider one by [AMM09] based on the following deviatoric and volumetric split:

(245)\[ \psi^{+}= \begin{cases} \psi_{\text{vol}} + \psi_{\text{dev}} = \frac{\bulk}{2} (\trace \bm\varepsilon)^2 + \mu(\bm\varepsilon_{\text{dev}} : \bm\varepsilon_{\text{dev}}) \quad \text{if} \quad \trace \bm\varepsilon \geq 0 \\ \psi_{\text{dev}} = \mu(\bm\varepsilon_{\text{dev}} : \bm\varepsilon_{\text{dev}}) \qquad \qquad \qquad \qquad \qquad \text{if} \quad \trace \bm\varepsilon < 0 \end{cases} \]

where \(\mu\) and \(\bulk\) are the shear and bulk moduli, respectively. For \(\psi^-\), this gives

(246)\[ \psi^{-}= \begin{cases} 0 \qquad \quad \text{if} \quad \trace \bm\varepsilon \geq 0 \\ \psi_{\text{vol}} \qquad \text{if} \quad \trace \bm\varepsilon < 0 \end{cases} \]

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

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

Weak form

For the linear damage model we write the functional as

(247)\[ \begin{aligned} \Pi (\bm u, \phi) &= \bar{W}(\bm u, \phi) - \Pi_{\text{ext}} (\bm u), \\ \Pi_{\text{ext}} (\bm u) &= \int_{\Omega}{\bm{u} \cdot \rho \bm{g}} \, dv + \int_{\partial \Omega}{\bm{u} \cdot \bar{\bm{t}}} \, da \end{aligned} \]

and by invoking the stationarity of \(\Pi(\bm u, \phi)\) with respect to \(\bm u\) and \(\phi\), we obtain the weak form as: Find \((\bm u, \phi) \in \mathcal{V} \times \mathcal{Q} \subset H^1 \left( \Omega \right) \times H^1 \left( \Omega \right)\), such that

(248)\[ \begin{aligned} R_u(\bm u, \phi) &\coloneqq \int_{\Omega} \nabla \bm v \tcolon \bm{\sigma}_{\text{degr}} \, dv - \bm L_{\text{ext}}(\bm v) =0, \qquad \forall v \in \mathcal{V} \\ R_{\phi}(\bm u, \phi) &\coloneqq \int_{\Omega} \left( q \, H + \nabla q \cdot \frac{2G_c l_0}{c_0}\nabla\phi \right) \, dv = 0, \qquad \forall q \in \mathcal{Q} \end{aligned} \]

where \(\bm L_{\text{ext}}(\bm v) = \int_{\Omega}{\bm{v} \cdot \rho \bm{g}} \, dv + \int_{\partial \Omega}{\bm{v} \cdot \bar{\bm{t}}} \, da\) is the external body and traction forces with \(\bm{\sigma}_{\text{degr}} \cdot \bm{n} =\bar{\bm{t}} \quad \text{on} \; \partial\Omega_N\), and we replaced the virtual displacement and phase field \(\delta \bm u\) and \(\delta \phi\) with the test functions \(\bm v, q\), i.e., \(\delta \bm u = \bm v, \, \delta \phi = q\).

The degraded Cauchy stress \(\bm{\sigma}_\text{degr}\) and \(H\) in (248) are derived as

(249)\[ \begin{aligned} \bm{\sigma}_{\text{degr}} &= \frac{\partial \left(g(\phi) \psi^{+} +\psi^{-} \right)}{\partial \bm u} =\begin{cases} g(\phi) \, \bm{\sigma}, \qquad \qquad \quad \text{if} \quad \trace \bm\varepsilon \geq 0 \\ \bm{\sigma}_{\text{vol}} + g(\phi) \, \bm{\sigma}_{\text{dev}}, \quad \text{if} \quad \trace \bm\varepsilon < 0 \end{cases}, \\ \bm{\sigma} &= \bm{\sigma}_{\text{vol}} + \bm{\sigma}_{\text{dev}} = \bulk \trace(\bm\varepsilon) \bm{I} + 2 \mu \, \bm\varepsilon_{\text{dev}} \end{aligned} \]

(250)\[ H = g'(\phi) \psi^{+} + \frac{G_c}{c_0l_0}\alpha'(\phi). \]

Alongside boundary conditions for displacement (\(\bm{u} = \bar{u} \; \text{on} \; \partial\Omega_D\) and \(\bm{\sigma}_{\text{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 \]

Note

The so-called (non-variationally consistent) hybrid formulation in [AGDL15], whereby \(\bm{\sigma}_{\text{degr}} = g(\phi)\bm{\sigma}\), has shown to perform well in tensile loading conditions when staggered (i.e. alternate minimization) solution schemes are adopted.

Newton linearization

As explained in Newton linearization, we need the Jacobian form of the (248) which is: Find \((\diff \bm{u}, \diff \phi) \in \mathcal{V} \times \mathcal{Q}\) such that

(251)\[ \begin{aligned} \int_{\Omega} \nabla \bm{v} \tcolon \diff \bm{\sigma}_{\text{degr}} dv &= -R_u, \quad \forall \bm{v} \in \mathcal{V}, \\ \int_{\Omega} \left( q \diff H + \nabla q \cdot \frac{2G_c l_0}{c_0}\nabla \diff \phi \right) \, dv &= -R_{\phi}, \quad \forall q \in \mathcal{Q} \end{aligned} \]

with

(252)\[ \diff H = \frac{\partial H}{\partial\bm u} \diff \bm u +\frac{\partial H}{\partial\phi} \diff \phi = g'(\phi) \diff \psi^{+} + g^{''}(\phi) \diff \phi \, \psi^{+} + \frac{G_c}{c_0l_0} \alpha^{''}(\phi) \, \diff \phi, \]

(253)\[ \begin{aligned} \diff \bm{\sigma}_{\text{degr}} &= \begin{cases} g(\phi) \, \diff \bm{\sigma} + g' \diff \phi \, \bm{\sigma}, \qquad \qquad \qquad \text{if} \quad \trace \bm\varepsilon \geq 0 \\ \diff \bm{\sigma}_{\text{vol}} + g(\phi) \, \diff \bm{\sigma}_{\text{dev}} + g' \diff \phi \, \bm{\sigma}_{\text{dev}}, \quad \text{if} \quad \trace \bm\varepsilon < 0 \end{cases}, \\ \diff \bm{\sigma} &= \diff\bm{\sigma}_{\text{vol}} + \diff\bm{\sigma}_{\text{dev}} = \bm{\sigma}_{\text{vol}}(\diff \bm u) + \bm{\sigma}_{\text{dev}}(\diff \bm u). \\ \end{aligned} \]

where \(\bm{\sigma}\) is defined in (249). The term \(\diff \psi^+\) is linked to crack irreversibility condition. If \(\diff \psi^{+}_{n+1} < \psi^{+}_{n}\), \(\diff \psi^{+}=0\), else

(254)\[ \diff \psi^{+} = \begin{cases} \bm{\sigma} \tcolon \diff \bm\varepsilon , \qquad \text{if} \quad \trace \bm\varepsilon \geq 0 \\ \bm{\sigma}_{\text{dev}} \tcolon \diff \bm\varepsilon, \quad \text{if} \quad \trace \bm\varepsilon < 0 \end{cases} \]

Note

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

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

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

where \(\eta \ll 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=8/3\) and \(c_0=2\), 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 (255) 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].

To complete the linearization given in (251) we need the first and second derivatives of \(g(\phi)\)

\[ g'(\phi) = -2(1 - \eta)(1 - \phi), \qquad g^{''}(\phi) = 2(1 - \eta), \]

and \(\alpha(\phi)\) which for AT1 are

\[ \alpha'(\phi) = 1, \qquad \alpha^{''}(\phi) = 0, \]

and for AT2 are

\[ \alpha'(\phi) = 2 \phi, \qquad \alpha^{''}(\phi) = 2. \]

Viscous regularization

To prevent convergence loss at solution jumps associated with an abrupt crack propagation, a viscous regularization is implemented to introduce a damping force opposing rapid changes in damage [MHW10] with (damage) viscosity \(\xi\). Expressions for \(H(\phi)\) and \(\diff H(\phi)\) above are updated to include the viscous term as

\[ H(\phi)= g'(\phi)\psi^{+} + \frac{G_c}{c_0l_0}\alpha'(\phi) + \xi \frac{\partial \phi}{\partial t} \]

and

\[ \diff H= g^{''}(\phi) \psi^{+} \diff \phi + g'(\phi) \diff \psi^{+} + \frac{G_c}{c_0l_0} \alpha^{''}(\phi) \diff \phi + \xi * \mathrm{shift_v} \diff \phi \]

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

Hyperelasticity, initial configuration

To extend the small strain model described in previous section, we need to use hyperelastic strain energy. For example, using volumetric-isochoric split for Neo-Hookean material defined in (69), we can write \(\psi^{+}\) and \(\psi^{-}\) as

(256)\[ \psi^{+}= \begin{cases} \psi_{\text{vol}} + \psi_{\text{iso}} = \bulk \, V(J) + \frac \secondlame 2 \left( \mathbb{\bar{I}_1} - 3 \right) \quad \text{if} \quad J - 1 \geq 0 \\ \psi_{\text{iso}} = \frac \secondlame 2 \left( \mathbb{\bar{I}_1} - 3 \right) \qquad \qquad \qquad \qquad \text{if} \quad J - 1 < 0 \end{cases} \]

(257)\[ \psi^{-}= \begin{cases} 0 \qquad \quad \text{if} \quad J - 1 \geq 0 \\ \psi_{\text{vol}} \qquad \text{if} \quad J - 1 < 0 \end{cases} \]

Substituting above energies into (244) and taking derivative with respect to \(\bm u\) and \(\phi\) gives the weak form for hyperelastic model as

(258)\[ \begin{aligned} R_u(\bm u, \phi) &\coloneqq \int_{\Omega_0} \nabla_X \bm v \tcolon \bm{F}\bm{S}_{\text{degr}} \, dV - \bm L_{\text{ext}}(\bm v) =0, \qquad \forall v \in \mathcal{V} \\ R_{\phi}(\bm u, \phi) &\coloneqq \int_{\Omega_0} \left( q \, H + \nabla q \cdot \frac{2G_c l_0}{c_0}\nabla\phi \right) \, dV = 0, \qquad \forall q \in \mathcal{Q} \end{aligned} \]

where the degraded stress is defined by

(259)\[ \begin{aligned} \bm{S}_{\text{degr}} &=\begin{cases} g(\phi) \, \bm{S}, \qquad \qquad \qquad \qquad \text{if} \quad J - 1 \geq 0 \\ \bm{S}_{\text{vol}} + g(\phi) \, \bm{S}_{\text{iso}}, \qquad \qquad \quad \text{if} \quad J - 1 < 0 \end{cases}, \\ \bm{S} &= \bm{S}_{\text{vol}} + \bm{S}_{\text{iso}} = \bulk J \, V' \, \bm{C}^{-1} + \secondlame J^{-2/3}\left( \bm{I} - \frac{1}{3}\mathbb{{I}_1}\bm{C}^{-1} \right) \end{aligned} \]

where we have used (76). The linearization of (258) can be written as

(260)\[ \begin{aligned} \int_{\Omega_0} \nabla_X \bm{v} \tcolon \left(\diff \bm{F} \bm{S}_{\text{degr}} + \bm{F} \diff \bm{S}_{\text{degr}} \right) dV &= -R_u, \quad \forall \bm{v} \in \mathcal{V}, \\ \int_{\Omega_0} \left( q \diff H + \nabla q \cdot \frac{2G_c l_0}{c_0}\nabla \diff \phi \right) \, dV &= -R_{\phi}, \quad \forall q \in \mathcal{Q} \end{aligned} \]

with

(261)\[ \begin{aligned} \diff \bm{S}_{\text{degr}} &= \begin{cases} g(\phi) \, \diff \bm{S} + g' \diff \phi \, \bm{S}, \qquad \qquad \qquad \text{if} \quad J - 1 \geq 0 \\ \diff \bm{S}_{\text{vol}} + g(\phi) \, \diff \bm{S}_{\text{iso}} + g' \diff \phi \, \bm{S}_{\text{iso}}, \quad \text{if} \quad J - 1 < 0 \end{cases}, \\ \diff \bm{S} &= \diff\bm{S}_{\text{vol}} + \diff\bm{S}_{\text{iso}}. \\ \end{aligned} \]

where \(\diff\bm{S}_{\text{iso}}\) and \(\diff\bm{S}_{\text{vol}}\) are derived in (156) and (157), respectively. Note that adding viscosity for hyperelastic will be similar to small strain case, but in \(H\) and its linearization \(\diff H\) we need to use the Neo-Hookean energy given in (256) and its linearization

(262)\[ \diff \psi^{+} = \begin{cases} \bm{S} \tcolon \diff \bm{E} , \qquad \text{if} \quad J - 1 \geq 0 \\ \bm{S}_{\text{iso}} \tcolon \diff \bm{E}, \quad \text{if} \quad J - 1 < 0 \end{cases} \]

Hyperelasticity, current configuration

To write hyperelasticity damage model in current configuration we need to compute the invariants of left Cauchy-Green tensor \(\bm b\) in (256) and write the residual as

(263)\[ \begin{aligned} r_u(\bm u, \phi) &\coloneqq \int_{\Omega_0} \nabla_x \bm v \tcolon \bm{\tau}_{\text{degr}} \, dV - \bm L_{\text{ext}}(\bm v) =0, \qquad \forall v \in \mathcal{V} \\ r_{\phi}(\bm u, \phi) &\coloneqq \int_{\Omega_0} \left( q \, H + \nabla q \cdot \frac{2G_c l_0}{c_0}\nabla\phi \right) \, dV = 0, \qquad \forall q \in \mathcal{Q} \end{aligned} \]

where the degraded stress \(\bm{\tau}_{\text{degr}}\) is defined by

(264)\[ \begin{aligned} \bm{\tau}_{\text{degr}} &=\begin{cases} g(\phi) \, \bm{\tau}, \qquad \qquad \qquad \qquad \text{if} \quad J - 1 \geq 0 \\ \bm{\tau}_{\text{vol}} + g(\phi) \, \bm{\tau}_{\text{iso}}, \qquad \qquad \quad \text{if} \quad J - 1 < 0 \end{cases}, \\ \bm{\tau} &= \bm{\tau}_{\text{vol}} + \bm{\tau}_{\text{iso}} = \bulk J \, V' \, \bm{I} + \secondlame J^{-2/3}\left( \bm{b} - \frac{1}{3}\mathbb{{I}_1}\bm{I} \right) \end{aligned} \]

where we have used (77). The linearization of (263) can be written as

(265)\[ \begin{aligned} \int_{\Omega_0} \nabla_x \bm{v} \tcolon \left( \diff \bm{\tau}_{\text{degr}} - \bm{\tau}_{\text{degr}} \left( \nabla_x \diff \bm{u} \right)^T \right) dV &= -r_u, \quad \forall \bm{v} \in \mathcal{V}, \\ \int_{\Omega_0} \left( q \diff H + \nabla q \cdot \frac{2G_c l_0}{c_0}\nabla \diff \phi \right) \, dV &= -r_{\phi}, \quad \forall q \in \mathcal{Q} \end{aligned} \]

where \(\diff \bm{\tau}_{\text{degr}} - \bm{\tau}_{\text{degr}} \left( \nabla_x \diff \bm{u} \right)^T = \nabla_x \diff \bm{u} \ \bm{\tau}_{\text{degr}} + \bm{F} \diff \bm{S}_{\text{degr}} \bm{F}^T\) and

(266)\[ \begin{aligned} \bm{F}\diff \bm{S}_{\text{degr}} \bm{F}^T &= \begin{cases} g(\phi) \, \bm{F}\diff \bm{S} \bm{F}^T + g' \diff \phi \, \bm{\tau}, \qquad \qquad \qquad \text{if} \quad J - 1 \geq 0 \\ \bm{F} \diff \bm{S}_{\text{vol}}\bm{F}^T + g(\phi) \, \bm{F} \diff \bm{S}_{\text{iso}} \bm{F}^T + g' \diff \phi \, \bm{\tau}_{\text{iso}}, \quad \text{if} \quad J - 1 < 0 \end{cases}, \\ \end{aligned} \]

where \(\bm{F}\diff\bm{S} \bm{F}^T\) is derived in (129) and finally

(267)\[ \diff \psi^{+} = \begin{cases} \bm{\tau} \tcolon \diff\bm\epsilon , \qquad \text{if} \quad J - 1 \geq 0 \\ \bm{\tau}_{\text{iso}} \tcolon \diff\bm\epsilon, \quad \text{if} \quad J - 1 < 0 \end{cases} \]