¹¹institutetext: Christopher J. Larsen ²²institutetext: Department of Mathematical Sciences, Worcester Polytechnic Institute, Worcester, MA 01609-2280, USA
²²email: [email protected] ³³institutetext: John E. Dolbow ⁴⁴institutetext: Department of Mechanical Engineering, Duke University, Durham, NC 27708, USA
⁴⁴email: [email protected] ⁵⁵institutetext: Oscar Lopez-Pamies ⁶⁶institutetext: Department of Civil and Environmental Engineering, University of Illinois, Urbana–Champaign, IL 61801-2352, USA
⁶⁶email: [email protected]

A Variational Formulation of Griffith Phase-Field Fracture
with Material Strength

C. J. Larsen J. E. Dolbow O. Lopez-Pamies

Abstract

In this expository Note, it is shown that the Griffith phase-field theory of fracture accounting for material strength originally introduced by Kumar, Francfort, and Lopez-Pamies (J Mech Phys Solids 112, 523–551, 2018) in the form of PDEs can be recast as a variational theory. In particular, the solution pair $({\bf u},v)$ defined by the PDEs for the displacement field ${\bf u}$ and the phase field $v$ is shown to correspond to the fields that minimize separately two different functionals, much like the solution pair $({\bf u},v)$ defined by the original phase-field theory of fracture without material strength implemented in terms of alternating minimization. The merits of formulating a complete theory of fracture nucleation and propagation via such a variational approach — in terms of the minimization of two different functionals — are discussed.

Keywords:

Brittle materials; Crack nucleation; Crack propagation; Strength; Fracture energy

1 Introduction

Kumar et al. (2018a; 2020) and Kumar and Lopez-Pamies (2020) have recently established that any attempt at the formulation of a complete macroscopic theory of fracture nucleation and propagation in nominally elastic brittle solids must account for three intrinsic material properties: ( $i$ ) the elasticity of the solid, ( $ii$ ) its strength, and ( $iii$ ) its critical energy release rate.

In the most basic setting, that of an isotropic linearly elastic brittle solid occupying an open bounded domain $\mathrm{\Omega}\subseteq\mathbb{R}^{3}$ that is subjected to monotonic and quasistatic loading conditions¹¹1We will come back to the more general case of non-monotonic and dynamic loading conditions in Section 4 below. in time $t\in[0,T]$ , these properties are characterized by:

i.

the stored-energy function

$W({\bf E})=\mu\,{\rm tr}\,{\bf E}^{2}+\dfrac{\lambda}{2}\left({\rm tr}\,{\bf E}\right)^{2}$ (1)

describing elasticity of the solid, that is, the elastic energy (per unit undeformed volume) stored by the solid when deformed a strain ${\bf E}$ ,
ii.

the strength surface

$\mathcal{F}(\mbox{\boldmath$\sigma$})=0$ (2)

describing the strength of the solid, that is, the set of critical stresses at which the solid fractures when subjected to monotonically increasing uniform stress $\sigma$ , and
iii.

the critical energy release rate

$G_{c}$ (3)

describing the intrinsic fracture energy of the solid, that is, the amount of energy (per unit undeformed area) required to create new surface in the solid from an existing crack.

In contrast to the specific functional forms (quadratic in ${\bf E}$ and constant) of the stored-energy function (1) and critical energy release rate (3), there is no one-size-fits-all functional form for the strength surface (2), since different solids can feature very different strength surfaces. In this work, for definiteness, we shall restrict attention to strength surfaces of the Drucker-Prager form

\mathcal{F}(\mbox{\boldmath$\sigma$})\equiv\sqrt{\mathcal{J}_{2}}+\dfrac{\sigma_{\texttt{ts}}}{\sqrt{3}\left(3\sigma_{\texttt{hs}}-\sigma_{\texttt{ts}}\right)}\mathcal{I}_{1}-\dfrac{\sqrt{3}\sigma_{\texttt{hs}}\sigma_{\texttt{ts}}}{3\sigma_{\texttt{hs}}-\sigma_{\texttt{ts}}}=0,

(4)

where $\mathcal{I}_{1}={\rm tr}\,\mbox{\boldmath$\sigma$}$ and $\mathcal{J}_{2}=\frac{1}{2}\,{\rm tr}\,\mbox{\boldmath$\sigma$}^{2}_{D}$ stand for the first and second principal invariants of the stress $\sigma$ and the deviatoric stress $\mbox{\boldmath$\sigma$}_{D}=\mbox{\boldmath$\sigma$}-\frac{1}{3}({\rm tr}\,\mbox{\boldmath$\sigma$}){\bf I}$ tensors, while the material constants $\sigma_{\texttt{ts}}>0$ and $\sigma_{\texttt{hs}}>0$ denote the uniaxial tensile and hydrostatic strengths of the solid, that is, they denote the critical stress values at which fracture nucleates under states of monotonically increased uniaxial tension $\mbox{\boldmath$\sigma$}={\rm diag}(\sigma>0,0,0)$ and tensile hydrostatic stress $\mbox{\boldmath$\sigma$}={\rm diag}(\sigma>0,\sigma>0,\sigma>0)$ , respectively.

Remark 1

According to our choice of signs in (4), any stress state such that

\left(3\sigma_{\texttt{hs}}-\sigma_{\texttt{ts}}\right)\mathcal{F}(\mbox{\boldmath$\sigma$})\geq 0

is in violation of the strength of the solid. The strength of the vast majority of solids, especially hard solids (e.g., glass), is such that $3\sigma_{\texttt{hs}}-\sigma_{\texttt{ts}}>0$ . Nonetheless, $3\sigma_{\texttt{hs}}-\sigma_{\texttt{ts}}<0$ for some soft solids (e.g., natural rubber).

Remark 2

The two-material-parameter strength surface (4), originally introduced by Drucker and Prager (1952) to model the yielding of soils, is arguably the simplest model that has proven capable of describing reasonably well the strength of many nominally brittle solids (Lopez-Pamies, 2023; Kumar et al., 2020; 2022; 2024; Kamarei et al., 2024), thus its use here as a representative template.

Remark 3

For given uniaxial tensile and hydrostatic strengths $\sigma_{\texttt{ts}}$ and $\sigma_{\texttt{hs}}$ , the strength surface (4) predicts the shear, biaxial tensile, and uniaxial compressive strengths

\sigma_{\texttt{ss}}=\dfrac{1}{\sqrt{3}}\left(-\dfrac{1}{3}+\dfrac{\sigma_{\texttt{hs}}}{\sigma_{\texttt{ts}}}\right)^{-1}\sigma_{\texttt{hs}},\quad\sigma_{\texttt{bs}}=\left(\dfrac{1}{3}+\dfrac{\sigma_{\texttt{hs}}}{\sigma_{\texttt{ts}}}\right)^{-1}\sigma_{\texttt{hs}},\quad{\rm and}\quad\sigma_{\texttt{cs}}=\left(-\dfrac{2}{3}+\dfrac{\sigma_{\texttt{hs}}}{\sigma_{\texttt{ts}}}\right)^{-1}\sigma_{\texttt{hs}},

which are defined as the critical stress values at which fracture nucleates under uniform states of monotonically increased shear stress $\mbox{\boldmath$\sigma$}={\rm diag}(\sigma>0,-\sigma,0)$ , biaxial tension $\mbox{\boldmath$\sigma$}={\rm diag}(\sigma>0,\sigma>0,0)$ , and uniaxial compression $\mbox{\boldmath$\sigma$}={\rm diag}(\sigma<0,0,0)$ . Direct use of these relations (or analogous ones for any other multi-axial stress state of interest) allows to rewrite (4) in terms of different pairs of critical strength constants. For hard solids, it is customary to use $\sigma_{\texttt{ts}}$ and $\sigma_{\texttt{cs}}$ , while for soft solids it is more convenient to use $\sigma_{\texttt{ts}}$ and $\sigma_{\texttt{hs}}$ . In this work, as indicated by (4), we favor the latter parametrization.

1.1 The problem

Consider hence a body made of an isotropic linearly elastic brittle solid with stored-energy function (1), Drucker-Prager strength surface (4), and critical energy release rate (3) that, initially, at time $t=0$ , occupies the open bounded domain $\mathrm{\Omega}$ . We denote the boundary of the body by $\partial\mathrm{\Omega}$ , its outward unit normal by ${\bf N}$ , and identify material points in the body by their initial position vector

{\bf X}\in\mathrm{\Omega}.

The body is subjected to a body force (per unit undeformed volume) $\textbf{b}({\bf X},t)$ , a displacement $\overline{{\bf u}}({\bf X},t)$ on a part $\partial\mathrm{\Omega}_{\mathcal{D}}$ of the boundary, and a surface force (per unit undeformed area) $\overline{{\bf s}}({\bf X},t)$ on the complementary part $\partial\mathrm{\Omega}_{\mathcal{N}}=\partial\mathrm{\Omega}\setminus\partial\mathrm{\Omega}_{\mathcal{D}}$ . In response to these stimuli — all of which, again, are assumed to be applied monotonically and quasistatically in time — the position vector ${\bf X}$ of a material point in the body will move to a new position specified by

{\bf x}={\bf X}+{\bf u}({\bf X},t),

(5)

where ${\bf u}({\bf X},t)$ is the displacement field. We write the associated strain at ${\bf X}$ and $t$ as

{\bf E}({\bf u})=\dfrac{1}{2}\left(\nabla{\bf u}+\nabla{\bf u}^{T}\right).

In addition to the deformation (5), the applied body force and boundary conditions may result in the nucleation and subsequent propagation of cracks in the body. We describe such cracks in a regularized fashion via the phase field

v=v({\bf X},t)

taking values in the range $[0,1]$ .

1.2 The phase-field fracture theory of Kumar, Francfort and Lopez-Pamies (2018)

According to the phase-field fracture formulation put forth by Kumar, Francfort and Lopez-Pamies (2018), the displacement field ${\bf u}_{k}({\bf X})={\bf u}({\bf X},t_{k})$ and phase field $v_{k}({\bf X})=v({\bf X},t_{k})$ at any material point ${\bf X}\in\overline{\mathrm{\Omega}}=\mathrm{\Omega}\cup\partial\mathrm{\Omega}$ and at any given discrete time $t_{k}\in\{0=t_{0},t_{1},...,t_{m},t_{m+1},...,t_{M}=T\}$ are determined by the system of coupled partial differential equations (PDEs)

\left\{\begin{array}[]{ll}{\rm Div}\left[(v_{k}^{2}+\eta_{\varepsilon})\dfrac{\partial W}{\partial{\bf E}}({\bf E}({\bf u}_{k}))\right]+\textbf{b}({\bf X},t_{k})=\textbf{0},&\,{\bf X}\in\mathrm{\Omega}\vspace{0.2cm}\\ {\bf u}_{k}({\bf X})=\overline{{\bf u}}({\bf X},t_{k}),&\,{\bf X}\in\partial\mathrm{\Omega}_{\mathcal{D}}\vspace{0.2cm}\\ \left[(v_{k}^{2}+\eta_{\varepsilon})\dfrac{\partial W}{\partial{\bf E}}({\bf E}({\bf u}_{k}))\right]{\bf N}=\overline{\textbf{s}}({\bf X},t_{k}),&\,{\bf X}\in\partial\mathrm{\Omega}_{\mathcal{N}}\end{array}\right.

(6)

and

\left\{\begin{array}[]{lll}\varepsilon\,\delta^{\varepsilon}G_{c}\bigtriangleup v_{k}=\dfrac{8}{3}v_{k}W({\bf E}({\bf u}_{k}))+\dfrac{4}{3}c_{\texttt{e}}({\bf X},t_{k})-\dfrac{\delta^{\varepsilon}G_{c}}{2\varepsilon},&\,\mbox{ if }v_{k}({\bf X})<v_{k-1}({\bf X}),&{\bf X}\in\mathrm{\Omega}\vspace{0.2cm}\\ \varepsilon\,\delta^{\varepsilon}G_{c}\bigtriangleup v_{k}\geq\dfrac{8}{3}v_{k}W({\bf E}({\bf u}_{k}))+\dfrac{4}{3}c_{\texttt{e}}({\bf X},t_{k})-\dfrac{\delta^{\varepsilon}G_{c}}{2\varepsilon},&\,\mbox{ if }v_{k}({\bf X})=1\;\mbox{ or }\;v_{k}({\bf X})=v_{k-1}({\bf X})>0,&{\bf X}\in\mathrm{\Omega}\vspace{0.2cm}\\ v_{k}({\bf X})=0,&\,\mbox{ if }v_{k-1}({\bf X})=0,&{\bf X}\in\mathrm{\Omega}\vspace{0.2cm}\\ \nabla v_{k}\cdot{\bf N}=0,&&\,{\bf X}\in\partial\mathrm{\Omega}\end{array}\right.

(7)

with, e.g., initial conditions ${\bf u}({\bf X},0)\equiv\textbf{0}$ and $v({\bf X},0)\equiv 1$ , where $\nabla v_{k}=\nabla v({\bf X},t_{k})$ , $\bigtriangleup v_{k}=\bigtriangleup v({\bf X},t_{k})$ , $\eta_{\varepsilon}$ is a positive constant of $o(\varepsilon)$ , and where, as elaborated below, $c_{\texttt{e}}({\bf X},t)$ is a driving force whose specific constitutive prescription depends on the particular form of strength surface $\mathcal{F}(\boldsymbol{\sigma})=0$ , while $\delta^{\varepsilon}$ is a non-negative coefficient whose specific constitutive prescription depends in turn on the particular form of $c_{\texttt{e}}({\bf X},t)$ .

Remark 4

The inequalities in (7) stem from the facts that, by definition, the phase field is bounded according to $0\leq v\leq 1$ and, by constitutive assumption, fracture is an irreversible process, in other words, healing is not allowed.

Remark 5

The parameter $\varepsilon$ in (7), with units of $length$ , regularizes sharp cracks. Accordingly, by definition, it can be arbitrarily small. In practice, $\varepsilon$ should be selected to be smaller than the smallest material length scale built in (6)-(7), which comes about because of the different units of the elastic stored-energy function $W({\bf E})$ ( $force/length^{2}$ ), the strength function $\mathcal{F}(\mbox{\boldmath$\sigma$})$ ( $force/length^{2}$ ), and the critical energy release rate $G_{c}$ ( $force/length$ ); see, e.g., Appendix C in (Kumar et al., 2024) and Appendix B in (Kamarei et al., 2024). As a rule of thumb, it typically suffices to set $\varepsilon<3G_{c}/(16\mathcal{W}_{\texttt{ts}})$ , where $\mathcal{W}_{\texttt{ts}}$ is given by expression (9)₁ below.

Remark 6

With regard to the preceding remark, it is important to emphasize that models of fracture in which $\varepsilon$ — or any such type of a length scale parameter — cannot be taken arbitrarily small are not approximations of any sharp fracture model. Interestingly, by virtue of the finite value of $\varepsilon$ that may be chosen in their implementation, those models may feature an apparent strength. That strength is an artifact of such an $\varepsilon$ , one that disappears as $\varepsilon\searrow 0$ , and not actual material strength.

Remark 7

On their own, the governing equations (6) and (7) are standard second-order PDEs for the displacement field ${\bf u}_{k}({\bf X})$ and the phase field $v_{k}({\bf X})$ , the latter of which is additionally subjected to standard variational inequalities. Accordingly, their numerical solution is amenable to a finite element (FE) staggered scheme in which (6) and (7) are discretized with finite elements and solved iteratively one after the other at every time step $t_{k}$ until convergence is reached. In the next section, as the main result of this Note, we show that the solution pair $({\bf u}^{\varepsilon}_{k},v^{\varepsilon}_{k})$ computed in such a staggered approach corresponds in fact to the fields that minimize separately two different functionals.

The driving force $c_{\emph{{e}}}({\bf X},t)$ and coefficient $\delta^{\varepsilon}$ .

For a solid whose strength is characterized by the Drucker-Prager strength surface (4), Kumar et al. (2020; 2022) and Kamarei et al. (2024) have worked out different constitutive prescriptions for the driving force $c_{\texttt{e}}({\bf X},t)$ and coefficient $\delta^{\varepsilon}$ that are equivalent in the limit as $\varepsilon\searrow 0$ , but contain different corrections of $O(\varepsilon^{0})$ . Here, for definiteness, we consider the constitutive prescription proposed by Kamarei et al. (2024). It reads

\displaystyle c_{\texttt{e}}({\bf X},t)=\alpha_{2}\sqrt{\mathcal{J}_{2}}+\alpha_{1}\mathcal{I}_{1}\qquad{\rm and}\qquad\delta^{\varepsilon}=f\left(\dfrac{\sigma_{\texttt{hs}}}{\sigma_{\texttt{ts}}}\right)\dfrac{3G_{c}}{16\mathcal{W}_{\texttt{ts}}\varepsilon}+g\left(\dfrac{\sigma_{\texttt{hs}}}{\sigma_{\texttt{ts}}}\right),

(8)

where

\displaystyle\alpha_{1}=\dfrac{1}{\sigma_{\texttt{hs}}}\delta^{\varepsilon}\dfrac{G_{c}}{8\varepsilon}-\dfrac{2\mathcal{W}_{\texttt{hs}}}{3\sigma_{\texttt{hs}}},\qquad\alpha_{2}=\dfrac{\sqrt{3}(3\sigma_{\texttt{hs}}-\sigma_{\texttt{ts}})}{\sigma_{\texttt{hs}}\sigma_{\texttt{ts}}}\delta^{\varepsilon}\dfrac{G_{c}}{8\varepsilon}+\dfrac{2\mathcal{W}_{\texttt{hs}}}{\sqrt{3}\sigma_{\texttt{hs}}}-\dfrac{2\sqrt{3}\mathcal{W}_{\texttt{ts}}}{\sigma_{\texttt{ts}}},

\displaystyle\mathcal{W}_{\texttt{ts}}=\dfrac{\sigma_{\texttt{ts}}^{2}}{2E},\qquad\mathcal{W}_{\texttt{hs}}=\dfrac{\sigma_{\texttt{hs}}^{2}}{2\kappa},

(9)

\left\{\begin{array}[]{l}\mathcal{I}_{1}={\rm tr}\,\mbox{\boldmath$\sigma$}=3\kappa v^{2}{\rm tr}\,{\bf E\vspace{0.2cm}}\\ \mathcal{J}_{2}=\dfrac{1}{2}\,{\rm tr}\,\mbox{\boldmath$\sigma$}^{2}_{D}=2\mu^{2}v^{4}{\rm tr}\,{\bf E}^{2}_{D},\quad{\bf E}_{D}={\bf E}-\dfrac{1}{3}\left({\rm tr}\,{\bf E}\right){\bf I},\end{array}\right.,

and where we have made use of the classical connections $E=\mu(3\lambda+2\mu)/(\lambda+\mu)$ and $\kappa=\lambda+\frac{2}{3}\mu$ between the Young’s modulus $E$ and bulk modulus $\kappa$ and the Lamé constants $\mu$ and $\lambda$ ; the specific forms of the non-negative functions $f$ and $g$ in (8)₂ can be found in Section 3.4 of (Kamarei et al., 2024).

The constitutive prescription (8) for the driving force $c_{\texttt{e}}$ and coefficient $\delta^{\varepsilon}$ leads to a complete macroscopic theory of fracture (6)-(7), one that over the past lustrum has been validated through direct comparisons with experiments on a wide range of nominally elastic brittle solids under a broad range of loading conditions (Kumar et al., 2018a; 2018b; 2020; 2022; 2024; Kumar and Lopez-Pamies, 2021; Kamarei et al., 2024).

The central reason for the apparent status of the PDEs (6)-(7) as a complete theory of fracture nucleation and propagation in a nominally elastic brittle solid subjected to monotonic and quasistatic loading conditions is that, by construction, they encode the competition between the elastic energy (1), the strength surface (2), and the critical energy release rate (3) of the solid in a way that is consistent with the large body of experimental evidence on fracture that has been amassed since the end of the 1800s.

The main objective of this Note is to show that the competition described by the PDEs (6)-(7) can be recast variationally as the minimization of two different functionals. We do so in Section 2. The resulting variational structure happens to be of the same type that has recently allowed to prove existence of solutions in the variational approach to sharp Griffith fracture (Francfort and Marigo, 1998) accounting for boundary loads (Larsen, 2021). It is also of the same type of variational structure that describes the original phase-field theory without material strength when implemented as an alternating minimization (Bourdin et al., 2000; 2006). We devote Section 3 to explaining this overarching connection. There, we also point to some misconceptions in the literature on the interpretation of the original phase-field theory as typically implemented. We conclude in Section 4 by recording a number of final comments.

2 The variational formulation of the phase-field fracture theory (6)-(7)

2.1 The deformation energy functional $\mathcal{E}^{\varepsilon}_{d}({\bf u}_{k};v_{k})$

For any fixed phase field $v_{k}({\bf X})\in[0,1]$ , a standard calculation suffices to show that equations (6) are nothing more than the Euler-Lagrange equations associated with variations in

{\bf u}_{k}\in\mathcal{A}_{u}=\left\{{\bf u}_{k}\in H^{1}(\mathrm{\Omega};\mathbb{R}^{3}):{\bf u}_{k}({\bf X})=\overline{{\bf u}}({\bf X},t_{k}),\;{\bf X}\in\partial\mathrm{\Omega}_{\mathcal{D}}\right\}

of the deformation energy functional

\mathcal{E}^{\varepsilon}_{d}({\bf u}_{k};v_{k}):=\displaystyle\int_{\mathrm{\Omega}}(v_{k}^{2}+\eta_{\varepsilon})W({\bf E}({\bf u}_{k})){\rm d}{\bf X}-\displaystyle\int_{\mathrm{\Omega}}\textbf{b}({\bf X},t_{k})\cdot{\bf u}_{k}\,{\rm d}{\bf X}-\displaystyle\int_{\partial\mathrm{\Omega}_{\mathcal{N}}}\overline{\textbf{s}}({\bf X},t_{k})\cdot{\bf u}_{k}\,{\rm d}{\bf X}.

(10)

Assuming that $\textbf{b}({\bf X},t_{k})\in L^{2}(\mathrm{\Omega};\mathbb{R}^{3})$ , $\overline{\textbf{{\bf u}}}({\bf X},t_{k})\in L^{2}(\partial\mathrm{\Omega}_{\mathcal{D}};\mathbb{R}^{3})$ , and $\overline{\textbf{s}}({\bf X},t_{k})\in L^{2}(\partial\mathrm{\Omega}_{\mathcal{N}};\mathbb{R}^{3})$ , and noting that $\eta_{\varepsilon}>0$ , we can invoke the classical theorems of existence and uniqueness in linear elastostatics (see, e.g., Fichera, 1973) to readily establish that

({\bf u}^{\varepsilon}_{k};v_{k})=\underset{\begin{subarray}{c}{\bf u}_{k}\in\mathcal{A}_{u}\end{subarray}}{\arg\min}\,\mathcal{E}^{\varepsilon}_{d}({\bf u}_{k};v_{k})

(11)

exists, is unique, and is the solution of equations (6).

2.2 The fracture functional $\mathcal{E}^{\varepsilon}_{f}(v_{k};{\bf u}_{k})$

Define the “undamaged” driving force

\widehat{c}_{\texttt{e}}({\bf E}({\bf u}))=\mu\alpha_{2}\sqrt{2\,{\rm tr}\,{\bf E}^{2}_{D}({\bf u})}+3\kappa\alpha_{1}{\rm tr}\,{\bf E}({\bf u})\,\textrm{ so that }\,c_{\texttt{e}}({\bf X},t)=v^{2}\widehat{c}_{\texttt{e}}({\bf E}({\bf u})).

In terms of this driving force, equations (7) can be rewritten as

\left\{\begin{array}[]{lll}\varepsilon\,\delta^{\varepsilon}G_{c}\bigtriangleup v_{k}=\dfrac{8}{3}v_{k}W({\bf E}({\bf u}_{k}))+\dfrac{4}{3}v_{k}^{2}\widehat{c}_{\texttt{e}}({\bf E}({\bf u}_{k}))-\dfrac{\delta^{\varepsilon}G_{c}}{2\varepsilon},&\,\mbox{ if }v_{k}({\bf X})<v_{k-1}({\bf X}),&{\bf X}\in\mathrm{\Omega}\vspace{0.2cm}\\ \varepsilon\,\delta^{\varepsilon}G_{c}\bigtriangleup v_{k}\geq\dfrac{8}{3}v_{k}W({\bf E}({\bf u}_{k}))+\dfrac{4}{3}v_{k}^{2}\widehat{c}_{\texttt{e}}({\bf E}({\bf u}_{k}))-\dfrac{\delta^{\varepsilon}G_{c}}{2\varepsilon},&\,\mbox{ if }v_{k}({\bf X})=1\;\mbox{ or }\;v_{k}({\bf X})=v_{k-1}({\bf X})>0,&{\bf X}\in\mathrm{\Omega}\vspace{0.2cm}\\ v_{k}({\bf X})=0,&\,\mbox{ if }v_{k-1}({\bf X})=0,&{\bf X}\in\mathrm{\Omega}\vspace{0.2cm}\\ \nabla v_{k}\cdot{\bf N}=0,&&\,{\bf X}\in\partial\mathrm{\Omega}.\end{array}\right.

(12)

For any fixed regularization length $\varepsilon>0$ and any fixed displacement field ${\bf u}_{k}({\bf X})\in H^{1}(\mathrm{\Omega};\mathbb{R}^{3})$ , yet another standard calculation suffices to show that equations (12) are the Euler-Lagrange equations associated with variations in

v_{k}\in\mathcal{A}_{v}=\left\{v_{k}\in H^{1}(\mathrm{\Omega};\mathbb{R}):0\leq v_{k}\leq 1,\,v_{k}\leq v_{k-1}\right\}

(13)

of the fracture functional

\mathcal{E}^{\varepsilon}_{f}(v_{k};{\bf u}_{k}):=\displaystyle\int_{\mathrm{\Omega}}v_{k}^{2}W({\bf E}({\bf u}_{k})){\rm d}{\bf X}+\displaystyle\int_{\mathrm{\Omega}}\dfrac{v_{k}^{3}}{3}\,\widehat{c}_{\texttt{e}}({\bf E}({\bf u}_{k})){\rm d}{\bf X}+\dfrac{3\,\delta^{\varepsilon}G_{c}}{8}\int_{\mathrm{\Omega}}\left(\dfrac{1-v_{k}}{\varepsilon}+\varepsilon\nabla v_{k}\cdot\nabla v_{k}\right)\,{\rm d}{\bf X}.

(14)

Since, as outlined next, $\mathcal{E}^{\varepsilon}_{f}$ has minimizers in the admissible set (13), we have that

(v^{\varepsilon}_{k};{\bf u}_{k})\in\underset{\begin{subarray}{c}v_{k}\in\mathcal{A}_{v}\end{subarray}}{\arg\min}\,\mathcal{E}^{\varepsilon}_{f}(v_{k};{\bf u}_{k})

(15)

exists and is a solution of equations (12).

We now briefly describe the existence of minimizers of $\mathcal{E}^{\varepsilon}_{f}$ . We suppose that $v_{k-1}=1$ for simplicity, the general case is similar. The idea is a straightforward application of the Direct Method in the Calculus of Variations. We consider a sequence $v_{n}$ such that $\mathcal{E}^{\varepsilon}_{f}(v_{n};{\bf u}_{k})\rightarrow\inf\mathcal{E}^{\varepsilon}_{f}(\cdot;{\bf u}_{k})$ . Note that

-\frac{1}{3}\displaystyle\int_{\mathrm{\Omega}}\,|\widehat{c}_{\texttt{e}}({\bf E}({\bf u}_{k}))|{\rm d}{\bf X}+\dfrac{3\,\delta^{\varepsilon}G_{c}}{8}\int_{\mathrm{\Omega}}\varepsilon\nabla v_{n}\cdot\nabla v_{n}\,{\rm d}{\bf X}\leq\mathcal{E}^{\varepsilon}_{f}(v_{n};{\bf u}_{k})\leq\displaystyle\int_{\mathrm{\Omega}}(W({\bf E}({\bf u}_{k}))+\,|\widehat{c}_{\texttt{e}}({\bf E}({\bf u}_{k}))|){\rm d}{\bf X}<\infty.

It follows that $\{\nabla v_{n}\}$ is bounded in $L^{2}(\mathrm{\Omega};\mathbb{R}^{3})$ , and since $|v_{n}|\leq 1$ , $\{v_{n}\}$ is bounded in $H^{1}(\mathrm{\Omega};\mathbb{R})$ . Therefore, a subsequence, not relabeled, satisfies $v_{n}\rightharpoonup v_{\infty}$ in $H^{1}(\mathrm{\Omega};\mathbb{R})$ for some $v_{\infty}\in H^{1}(\mathrm{\Omega};\mathbb{R})$ , and so $v_{n}\rightarrow v_{\infty}$ strongly in $L^{2}(\mathrm{\Omega};\mathbb{R})$ . As $\nabla v_{n}\rightharpoonup\nabla v_{\infty}$ in $L^{2}(\mathrm{\Omega};\mathbb{R}^{3})$ , we get

\int_{\mathrm{\Omega}}\left(\varepsilon\nabla v_{\infty}\cdot\nabla v_{\infty}\right)\,{\rm d}{\bf X}\leq\liminf_{n\rightarrow\infty}\int_{\mathrm{\Omega}}\left(\varepsilon\nabla v_{n}\cdot\nabla v_{n}\right)\,{\rm d}{\bf X}

by the convexity of this term.

The first two terms in $\mathcal{E}^{\varepsilon}_{f}$ are of the form $\int_{\mathrm{\Omega}}v_{k}^{p}F\,{\rm d}{\bf X}$ , where $F$ is integrable and $p=2,3$ . Since $v_{n}\rightarrow v_{\infty}$ strongly in $L^{2}(\mathrm{\Omega};\mathbb{R})$ , for a subsequence, not relabeled, $v_{n}\rightarrow v_{\infty}$ a.e. Therefore, $v_{n}^{p}F\rightarrow v_{\infty}^{p}F$ a.e., and since $v_{n}\in[0,1]$ , we have $v_{n}^{p}|F|\leq|F|$ . The Dominated Convergence Theorem then implies

\int_{\mathrm{\Omega}}v_{n}^{p}F\,{\rm d}{\bf X}\rightarrow\int_{\mathrm{\Omega}}v_{\infty}^{p}F\,{\rm d}{\bf X}.

The term in $(1-v_{k})$ converges similarly.

Putting all the terms together, we get

\mathcal{E}^{\varepsilon}_{f}(v_{\infty};{\bf u}_{k})\leq\liminf_{n\rightarrow\infty}\mathcal{E}^{\varepsilon}_{f}(v_{n};{\bf u}_{k})=\inf\mathcal{E}^{\varepsilon}_{f}(\cdot;{\bf u}_{k}),

and so $v_{\infty}$ is a minimizer. Finally, note that for any irreversibility constraint (of the type included in (13)) such that $v_{n}\leq V$ a.e. for some $V$ , then since $v_{n}\rightarrow v_{\infty}$ a.e., we also get $v_{\infty}\leq V$ a.e.

Remark 8

Contrary to the uniqueness of the minimizer (11) of the deformation energy functional (10), the minimizer (15) of the fracture functional (14) need not be unique.

2.3 The variational principle

Having defined the deformation energy and fracture functionals (10) and (14) and having established that they have minimizers, we are now ready to recast the phase-field fracture theory described by the PDEs (6)-(7) as a variational principle.

The idea — as in the standard implementation of the original phase-field fracture theory (see below) — is to alternately minimize the deformation energy and fracture functionals (10) and (14) at every time step $t_{k}$ until reaching a converged solution pair $(\mathfrak{u}_{k}^{\varepsilon},\mathfrak{v}_{k}^{\varepsilon})$ . Doing so leads to sequences ${\bf u}_{i}$ , $v_{i}$ satisfying

{\bf u}_{i}\mbox{ minimizes }\mathcal{E}^{\varepsilon}_{d}(\cdot\;;v_{i-1})

and

v_{i}\mbox{ minimizes }\mathcal{E}^{\varepsilon}_{f}(\cdot\;;{\bf u}_{i}).

The expectation — as, again, in the standard implementation of the original phase-field fracture theory (see below) — is that the limits $\mathfrak{u}_{k}^{\varepsilon}$ and $\mathfrak{v}_{k}^{\varepsilon}$ of ${\bf u}_{i}$ and $v_{i}$ , respectively, satisfy the minimality conditions

\mathfrak{u}_{k}^{\varepsilon}\mbox{ minimizes }\mathcal{E}^{\varepsilon}_{d}(\cdot\;;\mathfrak{v}_{k}^{\varepsilon})

and

\mathfrak{v}_{k}^{\varepsilon}\mbox{ minimizes }\mathcal{E}^{\varepsilon}_{f}(\cdot\;;\mathfrak{u}_{k}^{\varepsilon})

and therefore satisfy the PDEs (6)-(7). At present, there is no mathematical proof. However, as noted in the Introduction, the verity of this expectation has been supported by numerous simulations of experiments on a wide range of nominally elastic brittle solids under a broad range of loading conditions, which have also served as validation results for the theory. In the sequel, we provide some insight into why this is the case.

Remark 9

Larsen (2021) has recently shown that a similar variational principle based on the minimization of two different functionals in the setting of sharp — as opposed to phase-field — fracture is necessary in order to have existence of solutions in the variational approach to fracture accounting for boundary loads.

2.4 Some remarks on the competitions set up by the variational principle

On the one hand, the minimization of the deformation energy functional (10) is nothing more than the classical variational statement that forces in the body satisfy balance of linear and angular momenta for the specific case when the body is made of a linear elastic solid.

On the other hand, the minimization of the fracture functional (14) states that whether cracks nucleate and/or propagate depends on a competition among all three intrinsic material properties (1)-(3) of the solid at states when the body is in elastostatic equilibrium.

Uniform fields.

In particular, for the basic case when the body is subjected to a uniform strain ${\bf E}({\bf u}_{k})=\overline{{\bf E}}_{k}$ , with help of the notation

\displaystyle\overline{{\bf S}}_{k}=\dfrac{\partial W}{\partial{\bf E}}(\overline{{\bf E}}_{k})

for the corresponding uniform “undamaged” stress, noting that

\displaystyle\widehat{c}_{\texttt{e}}(\overline{{\bf E}}_{k})=\dfrac{3\delta^{\varepsilon}G_{c}}{8\varepsilon}\dfrac{(3\sigma_{\texttt{hs}}-\sigma_{\texttt{ts}})}{\sqrt{3}\sigma_{\texttt{hs}}\sigma_{\texttt{ts}}}\left(\mathcal{F}(\overline{{\bf S}}_{k})+\dfrac{\sqrt{3}\sigma_{\texttt{hs}}\sigma_{\texttt{ts}}}{3\sigma_{\texttt{hs}}-\sigma_{\texttt{ts}}}\right)+O(\varepsilon^{0}),

where we recall that the strength function $\mathcal{F}$ is defined in (4), the fracture functional (14) specializes to

	$\displaystyle\mathcal{E}^{\varepsilon}_{f}(v_{k};{\bf u}_{k})=$	$\displaystyle\displaystyle\int_{\mathrm{\Omega}}v_{k}^{2}W(\overline{{\bf E}}_{k}){\rm d}{\bf X}+\displaystyle\int_{\mathrm{\Omega}}\dfrac{v_{k}^{3}}{3}\left[\dfrac{3\delta^{\varepsilon}G_{c}}{8\varepsilon}\dfrac{(3\sigma_{\texttt{hs}}-\sigma_{\texttt{ts}})}{\sqrt{3}\sigma_{\texttt{hs}}\sigma_{\texttt{ts}}}\left(\mathcal{F}(\overline{{\bf S}}_{k})+\dfrac{\sqrt{3}\sigma_{\texttt{hs}}\sigma_{\texttt{ts}}}{3\sigma_{\texttt{hs}}-\sigma_{\texttt{ts}}}\right)+O(\varepsilon^{0})\right]{\rm d}{\bf X}+\vspace{0.2cm}$
		$\displaystyle\dfrac{3\,\delta^{\varepsilon}G_{c}}{8}\int_{\mathrm{\Omega}}\left(\dfrac{1-v_{k}}{\varepsilon}+\varepsilon\nabla v_{k}\cdot\nabla v_{k}\right)\,{\rm d}{\bf X}.$		(16)

For sufficiently small regularization length $\varepsilon$ , in regions where $v_{k}=1$ , the only terms that compete in the minimization of this functional are the second (i.e., the strength) and the third (i.e., the fracture energy) integrals. The first integral (i.e., the elastic energy) is inconsequential in this case.

In view of the definition of the strength function $\mathcal{F}$ in (4), the second integral in (16) starts at a value of 0 when $\overline{{\bf E}}_{k}=\textbf{0}$ . As the strain $\overline{{\bf E}}_{k}$ deviates from 0 along the given loading path, so does the second integral, which may become large enough to compete with the third integral and the minimizing displacement (11) to force the localization of the phase field $v_{k}$ near $v_{k}=0$ and hence the nucleation of fracture. As this localization of $v_{k}$ occurs, the first integral in (16) may become of comparable order to the other two integrals in a way that further enhances localization. Consistent with experimental observations, the numerical experiments referenced in the Introduction indicate that this localization of $v_{k}$ (when the strain is initially uniform in the body) happens when

\displaystyle\mathcal{F}(\overline{{\bf S}}_{k})=0,

that is, when the strength surface (4) of the solid is first violated.

Importantly — in stark contrast to the first integral which is non-negative — the second integral in (16) can be positive or negative. When it is negative, or when it is positive but not sufficiently large, there is no incentive to localize the phase field $v_{k}$ and hence the nucleation of fracture does not occur. Consistent with experimental observations, the numerical experiments referenced in the Introduction indicate that fracture nucleation does not occur so long as

\displaystyle(3\sigma_{\texttt{hs}}-\sigma_{\texttt{ts}})\mathcal{F}(\overline{{\bf S}}_{k})\leq 0.

Non-uniform fields due to the presence of large cracks.

For the opposite basic case when the body contains large²²2“Large” refers to large relative to the characteristic size of the underlying heterogeneities in the solid under investigation. By the same token, “small” refers to sizes that are of the same order or just moderately larger than the sizes of the heterogeneities. In practice, as a rule of thumb, “large” refers to large relative to the material length scale $G_{c}/\mathcal{W}_{\texttt{ts}}$ . cracks, the numerical experiments indicate that the competition set up by the fracture functional (14) describes that such cracks grow according to Griffith’s criticality condition (Griffith, 1921), consistent with experimental observations.

Importantly, all three integrals in (14) enter the competition that describes the growth of large cracks. In particular, because of its different sign for tensile and compressive stresses, the second integral (i.e., the strength) prevents the growth of large cracks under compressive loads, consistent with experimental observations (Liu and Kumar, 2024).

Arbitrary non-uniform fields.

Irrespective of the presence of large cracks, the stress field in a body is typically highly non-uniform. In this general case, consistent yet again with experimental observations, numerical experiments indicate that the competition set up by the fracture functional (14) describes that fracture nucleation is governed neither solely by the strength of the solid nor solely by Griffith’s criticality condition, but by an “interpolation” between the two.

What is more, the results indicate that a necessary condition for the nucleation of fracture at a material point ${\bf X}$ is that the strength surface of the solid be violated at that point. Such a condition is not sufficient, however, as indicated by the observation that stress singularities will always give rise to a stress state that violates the strength surface without necessarily resulting in the nucleation of fracture.

3 A comment on the variational structure of the original formulation of Griffith phase-field fracture without material strength

In the original phase-field fracture formulation (see, e.g., Bourdin et al., 2000; 2006), one seeks a minimizing pair $({\bf u}_{k}^{\varepsilon},v_{k}^{\varepsilon})$ of the energy functional

\mathcal{E}^{\varepsilon}({\bf u}_{k},v_{k}):=\displaystyle\int_{\mathrm{\Omega}}v_{k}^{2}W({\bf E}({\bf u}_{k})){\rm d}{\bf X}-\displaystyle\int_{\mathrm{\Omega}}\textbf{b}({\bf X},t_{k})\cdot{\bf u}_{k}\,{\rm d}{\bf X}-\displaystyle\int_{\partial\mathrm{\Omega}_{\mathcal{N}}}\overline{\textbf{s}}({\bf X},t_{k})\cdot{\bf u}_{k}\,{\rm d}{\bf X}+\dfrac{3\,G_{c}}{8}\int_{\mathrm{\Omega}}\left(\dfrac{1-v_{k}}{\varepsilon}+\varepsilon\nabla v_{k}\cdot\nabla v_{k}\right)\,{\rm d}{\bf X}.

This energy functional is not convex and hence difficult to minimize. It is separately convex, however, and so standard implementations involve alternating minimization. Precisely, sequences of pairs $({\bf u}_{i},v_{i})$ are found, satisfying

{\bf u}_{i}\mbox{ minimizes }\mathcal{E}^{\varepsilon}(\cdot\;,v_{i-1})

and

v_{i}\mbox{ minimizes }\mathcal{E}^{\varepsilon}({\bf u}_{i},\;\cdot).

Looking at the terms that are independent of $v_{k}$ and ${\bf u}_{k}$ respectively, this is equivalent to saying

{\bf u}_{i}\mbox{ minimizes }\mathcal{E}^{\varepsilon}_{d}(\cdot\;;v_{i-1})

and

v_{i}\mbox{ minimizes }\mathcal{E}^{\varepsilon}_{G}(\cdot\;;{\bf u}_{i}),

where we recall that the deformation energy functional $\mathcal{E}^{\varepsilon}_{d}$ is defined by (10), while

\mathcal{E}^{\varepsilon}_{G}(v_{k};{\bf u}_{k}):=\displaystyle\int_{\mathrm{\Omega}}v_{k}^{2}W({\bf E}({\bf u}_{k})){\rm d}{\bf X}+\dfrac{3\,G_{c}}{8}\int_{\mathrm{\Omega}}\left(\dfrac{1-v_{k}}{\varepsilon}+\varepsilon\nabla v_{k}\cdot\nabla v_{k}\right)\,{\rm d}{\bf X}

(17)

defines the Griffith energy functional; note that, as opposed to the fracture functional (14), the Griffith energy functional (17) does not account for the strength of the solid. As in the previous section, the expectation, without proof, is that the limits $\mathfrak{u}_{k}^{\varepsilon}$ and $\mathfrak{v}_{k}^{\varepsilon}$ of ${\bf u}_{i}$ and $v_{i}$ , respectively, satisfy the minimality conditions

\mathfrak{u}^{\varepsilon}_{k}\mbox{ minimizes }\mathcal{E}^{\varepsilon}_{d}(\cdot\;;\mathfrak{v}_{k}^{\varepsilon})

and

\mathfrak{v}_{k}^{\varepsilon}\mbox{ minimizes }\mathcal{E}^{\varepsilon}_{G}(\cdot\;;\mathfrak{u}_{k}^{\varepsilon}).

At this stage, it is plain that the solution pairs $({\bf u}_{k}^{\varepsilon},v_{k}^{\varepsilon})$ described by the phase-field fracture theory of Kumar, Francfort, and Lopez-Pamies (2018) and the original phase-field fracture theory are described by the same type of variational principle. The sole difference between them is that the original phase-field fracture theory is based on the Griffith energy functional (17), which does not account for the strength surface $\mathcal{F}(\mbox{\boldmath$\sigma$})=0$ of the solid, while the phase-field fracture theory of Kumar, Francfort, and Lopez-Pamies (2018) is based on the fracture functional (14), which accounts for all three material properties (1)-(3) required for a complete theory of fracture.

Remark 10

We close this section by emphasizing that for neither setting, the original phase-field fracture theory and that of Kumar, Francfort, and Lopez-Pamies (2018), is there a mathematical proof (yet) that guarantees that large cracks grow according to Griffith’s criticality condition. Whether they do can only be verified (thus far) numerically. For the case of the original phase-field fracture theory, this is because alternating minimization does not necessarily lead to the global minimization of $\mathcal{E}^{\varepsilon}$ . What is more, as discussed by Larsen (2023), it is possible to alter the energy functional $\mathcal{E}^{\varepsilon}$ so that both fracture nucleation and propagation are completely prevented in the alternating minimization described above, while still having $\mathcal{E}^{\varepsilon}$ $\mathrm{\Gamma}$ -converge to the energy functional of sharp Griffith fracture.

4 Final comments

In this Note, we have shown that the phase-field fracture theory of Kumar, Francfort, and Lopez-Pamies (2018) can be formulated as the minimization of two different functionals: ( $i$ ) a deformation energy functional $\mathcal{E}_{d}^{\varepsilon}$ and ( $ii$ ) a fracture functional $\mathcal{E}_{f}^{\varepsilon}$ . While the minimization of the deformation energy functional $\mathcal{E}_{d}^{\varepsilon}$ determines the deformation of the body, the minimization of the fracture functional $\mathcal{E}_{f}^{\varepsilon}$ determines where and when cracks nucleate and propagate.

We have also shown that the original phase-field fracture theory (Bourdin et al., 2000; 2008), as typically implemented in terms of an alternating minimization procedure, is described by the same variational formulation, with the key difference that the fate of cracks is determined not by the fracture functional $\mathcal{E}_{f}^{\varepsilon}$ , but by a Griffith energy functional $\mathcal{E}_{G}^{\varepsilon}$ , which does not account for the strength of solids.

While the focus of this Note has been on the basic case of fracture in nominally linear elastic brittle solids under monotonic and quasistatic loading conditions, it is apparent that the proposed variational approach — in terms of the minimization of two different functionals — has the flexibility to accommodate additional physical phenomena, such as non-monotonic and dynamic loading conditions, as well as possibly the inelasticity of actual solids. From a mathematical point of view, the proposed variational approach may also provide a fruitful path to advance the analysis of fracture theories at large.

Acknowledgements

Support for this work by the National Science Foundation through the Grants DMS–2206114 and DMS–2308169 is gratefully acknowledged. This work began during the 2024 workshop “Fracture as an emergent phenomenon” at the Mathematisches Forschungsinstitut Oberwolfach, who we thank for hosting us.

References

(1) Bourdin B, Francfort GA, Marigo JJ (2000). Numerical experiments in revisited brittle fracture. Journal of the Mechanics and Physics of Solids 48: 797–826.
(2) Bourdin B, Francfort GA, Marigo JJ (2008). The variational approach to fracture. Journal of Elasticity 91: 5–148.
(3) Drucker DC, Prager W (1952). Soil mechanics and plastic analysis for limit design. Quarterly of Applied Mathematics 10: 157–165.
(4) Fichera G (1973). Existence Theorems in Elasticity. In: Truesdell, C. (eds) Linear Theories of Elasticity and Thermoelasticity. Springer.
(5) Francfort GA, Marigo JJ (1998). Revisiting brittle fracture as an energy minimization problem. Journal of the Mechanics and Physics of Solids 46: 1319–1342.
(6) Griffith AA (1921). The phenomena of rupture and flow in solids. Philosophical Transactions of the Royal Society of London. Series A 221: 163–198.
(7) Kamarei F, Kumar A, Lopez-Pamies O (2024). The poker-chip experiments of synthetic elastomers. In preparation.
(8) Kumar A, Francfort GA, Lopez-Pamies O (2018a). Fracture and healing of elastomers: A phase-transition theory and numerical implementation. Journal of the Mechanics and Physics of Solids 112: 523–551.
(9) Kumar A, Ravi-Chandar K, Lopez-Pamies O (2018b). The configurational-forces view of fracture and healing in elastomers as a phase transition. International Journal of Fracture 213: 1–16.
(10) Kumar A, Lopez-Pamies O (2020). The phase-field approach to self-healable fracture of elastomers: A model accounting for fracture nucleation at large, with application to a class of conspicuous experiments. Theoretical and Applied Fracture Mechanics 107: 102550.
(11) Kumar A, Bourdin B., Francfort GA, Lopez-Pamies O (2020). Revisiting nucleation in the phase-field approach to brittle fracture. Journal of the Mechanics and Physics of Solids 142: 104027.
(12) Kumar A, Lopez-Pamies O (2021). The poker-chip experiments of Gent and Lindley (1959) explained. Journal of the Mechanics and Physics of Solids 150: 104359.
(13) Kumar A, Ravi-Chandar K, Lopez-Pamies O (2022). The revisited phase-field approach to brittle fracture: Application to indentation and notch problems. International Journal of Fracture 237: 83–100.
(14) Kumar A, Liu Y, Dolbow JE, Lopez-Pamies O (2024). The strength of the Brazilian fracture test. Journal of the Mechanics and Physics of Solids 182:105473.
(15) Larsen CJ (2021). Variational fracture with boundary loads. Applied Mathematics Letters 121: 107437.
(16) Larsen CJ (2023). Variational phase-field fracture with controlled nucleation. Mechanics Research Communications 128: 104059.
(17) Liu C, Kumar A (2024). Emergence of tension-compression asymmetry from a complete phase-field approach to brittle fracture. Submitted.
(18) Lopez-Pamies, O (2023). Journal Club: Strength revisited: One of three basic ingredients needed for a complete macroscopic theory of fracture. https://imechanica.org/node/26641.

A Variational Formulation of Griffith Phase-Field Fracture with Material Strength