High-accuracy phase-field models for brittle fracture based on a new family of degradation functions

TL;DR

This paper introduces a new family of degradation functions for phase-field models of brittle fracture, significantly improving accuracy in predicting crack initiation and growth while maintaining elastic response.

Contribution

A novel parametric family of degradation functions is proposed, enhancing phase-field fracture models' predictive accuracy and preserving elastic behavior before fracture occurs.

Findings

Superiority over classical quadratic degradation functions in numerical tests

Improved prediction of critical loads for crack nucleation

Enhanced modeling of crack propagation dynamics

Abstract

Phase-field approaches to fracture based on energy minimization principles have been rapidly gaining popularity in recent years, and are particularly well-suited for simulating crack initiation and growth in complex fracture networks. In the phase-field framework, the surface energy associated with crack formation is calculated by evaluating a functional defined in terms of a scalar order parameter and its gradients, which in turn describe the fractures in a diffuse sense following a prescribed regularization length scale. Imposing stationarity of the total energy leads to a coupled system of partial differential equations, one enforcing stress equilibrium and another governing phase-field evolution. The two equations are coupled through an energy degradation function that models the loss of stiffness in the bulk material as it undergoes damage. In the present work, we introduce a new parametric family of degradation functions aimed at increasing the accuracy of phase-field models in predicting critical loads associated with crack nucleation as well as the propagation of existing fractures. An additional goal is the preservation of linear elastic response in the bulk material prior to fracture. Through the analysis of several numerical examples, we demonstrate the superiority of the proposed family of functions to the classical quadratic degradation function that is used most often in the literature.

Tables

Table 1. Table 1 : Simulated critical internal bending moment for the four-point bending specimen. R ​ E 𝑅 𝐸 RE denotes relative error.

Simulation	Description	$\ell$	$\ell /h$	$R{E}_{P_c^{\text{calib}}}$	$U_c$	$M_c$	$R{E}_{M_c}$
run		(mm)		(%)	(mm)	(N-mm)	(%)
1	quadratic	0.94	2.0	$0.09$	$-0.0321$	187.80	$3.99$
2	$n=5.26$	0.94	2.0	$0.21$	$-0.0298$	192.75	$6.72$
3	$n=5.314$	0.5	2.0	$0.18$	$-0.0285$	183.75	$1.74$
4	$n=5.325$	0.3	2.0	$0.16$	$-0.0280$	180.25	$0.194$

Table 2. Table 2 : Material properties for the bi-material specimen.

Region	$A$	$B$
$E$	100 GPa	200 GPa
$\nu$	0.2	0.2
${𝒢}_c$	0.1 N/mm	$-$
${\sigma }_c$	70 MPa	$-$

Table 3. Table 3 : Details of simulation results pertaining to the bi-material problem: magnitude of applied displacement at crack initiation ( U c subscript 𝑈 𝑐 U_{c} ), total vertical force at top boundary ( F c subscript 𝐹 𝑐 F_{c} ), maximum tensile stress ( σ max subscript 𝜎 max \sigma_{\text{max}} ), and phase-field value in the critical element ( ϕ c subscript italic-ϕ 𝑐 \phi_{c} ).

Simulation	Description	$\ell$ (mm)	$U_c$ (mm)	$F_c$ (N)	${\sigma }_{\text{max}}$ (MPa)	${\varphi }_c$
1	quadratic	$0.31$	0.034	5142	71.02	0.4430
2	quadratic	$1.25$	0.023	3415	48.06	0.4611
3	quadratic	$5.0$	0.019	2682	39.71	0.4135
4	$n=4.4$	$1.25$	0.023	3562	69.71	0.1103
5	$n=2.95$	$5.0$	0.023	3553	70.17	0.1235

Equations

$$G-{\mathcal{G}}_c$$

$$\dot{a}$$

$$\left(G-{\mathcal{G}}_c\right)\dot{a}$$

$$\Psi \left(\mathbf{u},\Gamma \right)={\int }_{\Omega \backslash \Gamma }\frac{1}{2}\mathbf{\epsilon }\left(\mathbf{u}\right):{\mathbb{C}}^e:\mathbf{\epsilon }\left(\mathbf{u}\right)\mathrm{d}\Omega +{\mathcal{G}}_c{\mathcal{H}}^{n-1}\left(\Gamma \right)$$

$${\Gamma }_{t+\Delta t}\supseteq {\Gamma }_t$$

$$\Psi \left(\mathbf{u},\varphi \right)={\int }_{\Omega }\frac{1}{2}\left[{\left(1-\varphi \right)}^2+\kappa \right]\mathbf{\epsilon }\left(\mathbf{u}\right):{\mathbb{C}}^e:\mathbf{\epsilon }\left(\mathbf{u}\right)\mathrm{d}\Omega +{\mathcal{G}}_c{\int }_{\Omega }\left(\frac{1}{2\ell }{\varphi }^2+\frac{\ell }{2}\Vert \nabla \varphi {\Vert }^2\right)\mathrm{d}\Omega$$

$$\varphi \left(x\right)=\{\begin{array}{ll}1,& x=a\\ 0,& \text{otherwise}.\end{array}$$

$$\Phi=\Set{\phi\in H^{1}\left(\mathbb{R};\left[0,1\right]\right)}{\begin{array}[]{l}\phi\left(a\right)=1\\ \phi^{\prime}\left(x\right)>0\text{ if }x<a\\ \phi^{\prime}\left(x\right)<0\text{ if }x>a\\ \phi\left(x\right)\rightarrow 0\text{ as }x\rightarrow\pm\infty\end{array}}.$$

$$\varphi \left(x\right)$$

$$W_{\ell }\left(\mathbf{\epsilon }\left(\mathbf{u}\right),\varphi ,\nabla \varphi \right)=\frac{1}{2}\mathbf{\epsilon }:\mathbb{C}\left(\varphi \right):\mathbf{\epsilon }+w\left(\varphi \right)+\frac{1}{2}{w}_1{\ell }^2\nabla \varphi \cdot \nabla \varphi$$

$$\mathbb{C}\left(\varphi \right)=g\left(\varphi \right){\mathbb{C}}^e,$$

$$\Pi =\Psi -W={\int }_{\Omega }\left[g\left(\varphi \right)\psi \left(\mathbf{\epsilon }\right)+{\mathcal{G}}_c\left(\frac{1}{2\ell }{\varphi }^2+\frac{\ell }{2}\nabla \varphi \cdot \nabla \varphi \right)\right]\mathrm{d}\Omega -{\int }_{\Omega }\mathbf{b}\cdot \mathbf{u}\mathrm{d}\Omega -{\int }_{\partial {\Omega }^t}\mathbf{t}\cdot \mathbf{u}\mathrm{d}S$$

$$\begin{gathered} \delta \Pi ={\int }_{\Omega }g\left(\varphi \right)\frac{\partial \psi }{\partial \mathbf{\epsilon }}:\delta \mathbf{\epsilon }\mathrm{d}\Omega -{\int }_{\Omega }\mathbf{b}\cdot \delta \mathbf{u}\mathrm{d}\Omega -{\int }_{\partial {\Omega }^t}\mathbf{t}\cdot \delta \mathbf{u}\mathrm{d}S \\ +{\int }_{\Omega }{g}^{\prime }\left(\varphi \right)\psi \left(\mathbf{\epsilon }\right)\delta \varphi \mathrm{d}\Omega +{\mathcal{G}}_c{\int }_{\Omega }\left(\frac{1}{\ell }\varphi \delta \varphi +\ell \nabla \varphi \cdot \delta \nabla \varphi \right)\mathrm{d}\Omega =0. \end{gathered}$$

$${\int }_{\Omega }g\left(\varphi \right)\frac{\partial \psi }{\partial \mathbf{\epsilon }}:\delta \mathbf{\epsilon }\mathrm{d}\Omega ={\int }_{\Omega }\mathbf{b}\cdot \delta \mathbf{u}\mathrm{d}\Omega +{\int }_{\partial {\Omega }^t}\mathbf{t}\cdot \delta \mathbf{u}\mathrm{d}S$$

$${\int }_{\Omega }{g}^{\prime }\left(\varphi \right)\psi \left(\mathbf{\epsilon }\right)\delta \varphi \mathrm{d}\Omega +{\mathcal{G}}_c{\int }_{\Omega }\left(\frac{1}{\ell }\varphi \delta \varphi +\ell \nabla \varphi \cdot \nabla \delta \varphi \right)\mathrm{d}\Omega =0.$$

$$\nabla \cdot \left[g\left(\varphi \right)\mathbf{\sigma }\right]+\mathbf{b}$$

$$g\left(\varphi \right)\mathbf{\sigma }\cdot \mathbf{n}$$

$$\mathbf{u}$$

$${\mathcal{G}}_c{\ell }_0{\nabla }^2\varphi -\frac{{\mathcal{G}}_c}{{\ell }_0}\varphi$$

$$\nabla \varphi \cdot \mathbf{n}$$

$$\varphi {\left(\mathbf{x}\right)}_{t+\Delta t}\ge \varphi {\left(\mathbf{x}\right)}_t\forall \mathbf{x}\in \Omega .$$

$$\mathcal{H}\left(\mathbf{x},t\right)=\underset{s\in \left[0,t\right]}{\max }\psi \left(\mathbf{\epsilon }\left(\mathbf{x},s\right)\right)$$

$$\mathcal{H}\left(\mathbf{x},t\right)=\{\begin{array}{ll}\underset{s\in \left[0,t\right]}{\max }\psi \left(\mathbf{\epsilon }\left(\mathbf{x},s\right)\right)& \text{if }\varphi >{\varphi }_c\\ \psi \left(\mathbf{\epsilon }\left(\mathbf{x},t\right)\right)& \text{otherwise.}\end{array}$$

$$\mathbf{u}=\sum _{I=1}^m{\mathbf{N}}_I^{\mathbf{u}}{\mathbf{u}}_I\text{and}\varphi =\sum _{I=1}^m{N}_I{\varphi }_I$$

$${\mathbf{N}}_I=\left[\begin{array}{cc}{N}_I& 0\\ 0& N_I\end{array}\right]$$

$$\mathbf{\epsilon }=\sum _{I=1}^m{\mathbf{B}}_I^{\mathbf{u}}{\mathbf{u}}_I\text{and}\nabla \varphi =\sum _{I=1}^m{\mathbf{B}}_I^{\varphi }{\varphi }_I$$

$${\mathbf{B}}_I^{\mathbf{u}}=\left[\begin{array}{cc}{N}_{I,x}& 0\\ 0& N_{I,y}\\ N_{I,y}& N_{I,x}\end{array}\right]\text{and}{\mathbf{B}}_I^{\varphi }=\left[\begin{array}{c}{N}_{I,x}\\ N_{I,y}\end{array}\right].$$

$${\mathbf{r}}_I^{\mathbf{u}}$$

$$r_I^{\varphi }$$

$${\left\{\varphi \right\}}_i^{m+1}={\left\{\varphi \right\}}_i^m-{\left[{\mathbf{K}}^{\varphi \varphi }\left({\varphi }_i^m\right)\right]}^{-1}{\left\{r^{\varphi }\right\}}_i^m$$

$${\mathbf{K}}_{IJ}^{\varphi \varphi }\left({\varphi }_i^m\right)={\int }_{\Omega }\left[{\mathcal{G}}_c\ell {\mathbf{B}}_I^{\varphi T}{\mathbf{B}}_J^{\varphi }+\left(\frac{{\mathcal{G}}_c}{\ell }+\mathcal{H}{g}^{\prime \prime }\left({\varphi }_i^m\right)\right)N_I{N}_J\right]\mathrm{d}\Omega$$

$$g^{\prime }\left({\varphi }_1\right)=g^{\prime }\left({\varphi }_2\right)+g^{\prime \prime }\left({\varphi }_2\right)\left({\varphi }_1-{\varphi }_2\right)+\frac{g^{\prime \prime \prime }\left({\varphi }_2\right)}{2}{\left({\varphi }_1-{\varphi }_2\right)}^2+\dots$$

$$0\approx g^{\prime }\left(\varphi \right)+g^{\prime \prime }\left(\varphi \right)\left(1-\varphi \right)$$

Full text

1 Introduction

The accurate simulation of fracture evolution in solids is a major challenge for computational algorithms, in large part due to crack paths that are generally unknown a priori. In this regard, phase-field approaches have shown great potential with their ability to automatically determine the direction of crack propagation through minimization of an energy functional. The phase-field framework naturally handles the emergence of phenomena such as crack nucleation and branching without the need to introduce additional criteria. In particular, formulations derived from the variational theory of Francfort and Marigo (1998) have received a lot of attention from the applied mechanics community due to its strong ties to Griffith’s theory for brittle fracture. Phase-field models belong to the category of continuum approaches for fracture propagation, utilizing a diffuse representation of cracks in place of actual discontinuities. The amount of crack regularization is controlled via a prescribed length scale $\ell$, which constitutes an additional parameter of the model.

The aim of the present work is to address two long standing issues associated with the phase-field formulation that arise in conjunction with use of the now-classical quadratic degradation function. The first has to do with premature stiffness degradation stemming from the evolution of damage around regions of stress concentration. The second and more serious issue deals with the observed dependence of simulated failure loads on the phase-field regularization parameter, a phenomenon that has largely gone unexplored in the literature until very recently. The problem is most noticeable in cases of crack growth emanating from a notch and undermines the usefulness of phase-field approaches in solving problems that involve crack initiation at a priori unknown locations (Klinsmann et al., 2015). The parameter $\ell$ was initially introduced by Bourdin et al. (2000) as a purely mathematical construct allows for the Griffith energy corresponding to a discrete crack to be recovered in the limit as $\ell$ goes to zero, in the sense of $\Gamma$-convergence (Braides, 2006). Of the two aforementioned issues, the first may be remedied by making use of alternative formulations as discussed in Pham et al. (2011). On the other hand, the dependency of mechanical response on $\ell$ is not yet fully understood. It has been suggested recently (e.g. Bourdin et al., 2014) that the latter should be viewed as a material parameter that is closely connected to the crack nucleation stress. While we consider the two points raised above as distinct issues, we nonetheless recognize that they are also closely inter-related, in particular because the first often exacerbates the second.

Our contribution in the following work is twofold. First, we provide a conceptual explanation of how the choice of length scale can result to either delay or acceleration of failure under quasi-static conditions. Secondly, we introduce a new family of degradation functions that allows for correctly reproducing the onset of failure for reasonably chosen arbitrary values of the regularization parameter. The latter point is important since the problem of regularization-dependent material response is not solved in the alternative formulations previously mentioned, and furthermore is not only confined to brittle fracture as demonstrated in the numerical results of Areias et al. (2016) on cracking in elastoplastic materials. Enthusiasm in the relatively new phase-field paradigm has led to a number of multi-physics applications, which include cracking in piezoelectric solids (Miehe et al., 2010b), fluid-driven fracture propagation (Mikelić et al., 2015; Miehe et al., 2015b), thermal shock-induced cracks (Bourdin et al., 2014) and fragmentation of battery electrode particles (Miehe et al., 2015a). This underscores the need for quantitative accuracy with regard to the fracture model, particularly in the case of crack nucleation which is often the critical failure mechanism for many such applications.

The remainder of this paper is structured as follows: We begin in Section 2 with a discussion of important fundamental concepts underlying the phase-field approach as well as our motivation for pursuing the current research direction. Specifics regarding the formulation used in the present work are given in Section 3 which also includes important details with regard to numerics; in particular for the case where the phase-field evolution equation is nonlinear, we outline a linearization scheme based on a truncated Taylor series approximation that avoids the implementation of nested loops in the solution scheme. The following two sections contain the main novelties of the current work: section 4 begins with a numerical example demonstrating the apparent contradiction that is often seen between simulation results and what is expected from the $\Gamma$-convergence property of the fracture phase-field model. This is followed by a discussion which aims to explain why the latter alone is not sufficient to ensure the proper behavior of the model. In Section 5, we propose a new parametric family of degradation functions that aims to increase the accuracy of phase-field simulations by addressing key issues discussed in the previous section. Superiority of the resulting formulation over the classical model employing quadratic degradation is demonstrated via several numerical examples in Section 6. Finally, concluding remarks and outlook are given in Section 7.

2 Theoretical aspects of phase-field modelling

The phase-field framework was first introduced by Fix (1983) and Langer (1986) for modeling phase transitions in materials, and later extended to free discontinuity problems by Ambrosio and Tortorelli (1990) who worked on image segmentation. Its specific application to crack propagation in solids is much more recent, and is the result of independent work by researchers coming from the fields of physics (Aranson et al., 2000; Karma et al., 2001) and applied mechanics (Bourdin et al., 2000). We adapt the latter viewpoint in this study, and furthermore note that while the original formulation introduced by Bourdin et al. has remained virtually unchanged in current usage, the argument on what constitutes proper solutions to the resulting equations as well as the meaning of key quantities is far from resolved. In view of this, we begin with a short review of theory pertaining to the phase-field formulation for brittle fracture along with a discussion of significant developments in the field in order to provide context for the present work.

2.1 Brittle fracture: from Griffith to Francfort-Marigo

Griffith (1921) can be credited as being the first to formally state the thermodynamic principles governing the propagation of fractures in brittle materials that has become the foundation of modern linear elastic fracture mechanics. According to Griffith’s theory, an existing crack will propagate when the rate of energy release $G$ associated with crack extension exceeds a critical value equal to the material fracture toughness, $\mathcal{G}_{c}$. This can be expressed via the following set of Kuhn-Tucker conditions (Negri and Ortner, 2008):

[TABLE]

where $\dot{a}$ denotes the rate of crack length increase. The first inequality precludes the case of unstable cracking where $G>\mathcal{G}_{c}$, while the second is an irreversibility constraint that prevents unphysical healing of fractures. Finally, condition (1c) implies that $G$ must be equal to $\mathcal{G}_{c}$ when the crack is growing, and conversely that a crack cannot extend when $G<\mathcal{G}_{c}$. An important shortcoming of Griffith’s theory is its inability to accommodate crack nucleation or predict the branching of fractures. An extension of the framework was developed by Francfort and Marigo (1998) in the form of a variational theory of fracture, which is aimed at overcoming the earlier drawbacks through adoption of an energy minimization paradigm. It stipulates that the total potential energy corresponding to a linear elastic body $\Omega$ containing a set of crack points $\Gamma$ can be written as a sum of bulk and surface terms. That is,

[TABLE]

where $\bm{\varepsilon}\left(\bm{u}\right)=\frac{1}{2}\left[\nabla\otimes\bm{u}+\bm{u}\otimes\nabla\right]$ is the symmetric small-strain tensor, $\mathbb{C}^{e}$ is the standard linear isotropic elasticity tensor, and $\mathcal{H}^{n-1}$ is the $\left(n-1\right)$-dimensional Hausdorff measure giving the surface area associated with $\Gamma$. Equation (2) is referred to as the Griffith functional, and it is assumed that for some given boundary conditions on $\Omega$, the unknown displacements $\bm{u}$ as well as the crack set $\Gamma$ can be obtained via a global minimization of said functional subject to the irreversibility condition

[TABLE]

that is comparable to (1b). In contrast, Griffith requires only stationarity of (2). Furthermore in the variational theory, $\Gamma$ is not restricted to consist of a single crack. Thus a body that is initially without flaw ($\Gamma=\emptyset$) may nucleate a crack if the resulting configuration has lower total energy compared to one where no crack forms. Similarly, a crack is allowed to branch if this leads to a lower potential energy than simple extension. The directions of advance are naturally obtained as those leading to minimum increase in (2), so that in theory the energy minimization framework is able to handle crack initiation and branching without the need to introduce additional criteria.

2.2 Phase-field and gradient damage models

The main difficulty in performing a direct minimization of the Griffith functional is that $\bm{u}$ is generally discontinuous across $\Gamma$, so that (2) contains a locus of jump sets whose locations are a priori unknown and which generally do not align with the predefined domain discretization that is utilized in numerical simulations. To render the problem tractable, Bourdin et al. (2000) adopted a strategy wherein the minimization is instead performed on an approximation of the Griffith functional having regularized jump sets so that $\bm{u}$ is continuous over the entire domain. This was inspired by the earlier work of Ambrosio and Tortorelli (1990, 1992) who solved a similar problem in image segmentation involving the functional of Mumford and Shah (1989). A scalar order parameter known as the crack phase-field is introduced to interpolate between fractured and intact regions, with the amount of regularization controlled by a characteristic length parameter denoted by $\ell$. The validity of such as strategy rests on whether the regularized approximation tends towards the original functional as $\ell$ goes to zero, in the sense of $\Gamma$-convergence (Braides, 2006). Although an additional equation governing the phase-field evolution must now be solved along with the linear momentum equation, the main advantage of this approach is that numerical solutions may be obtained via classical finite element algorithms as both $\bm{u}$ and the phase-field are continuous. Bourdin et al. (2000)’s regularization of (2) took the form

[TABLE]

where $\phi$ is the phase-field that takes on values between 0 and 1, corresponding respectively to fully intact and broken states. On the other hand, $\kappa$ is a small positive constant meant to ensure positivity of the bulk energy as $\phi\rightarrow 1$. The two important features of the above expression are (a) the replacement of $\mathcal{H}^{n-1}\left(\Gamma\right)$ by an elliptic functional that calculates the combined length of all the cracks, and (b) the coefficient $\left(1-\phi\right)^{2}+\kappa$ known as the energy degradation function that penalizes the material stiffness according to the value of $\phi$. Equation 4 is essentially a direct adaptation of the earlier functional of Ambrosio and Tortorelli (1992), for which a proof of $\Gamma$-convergence was subsequently given by Chambolle (2004).

2.2.1 Alternative variational problems

It was later suggested by Miehe et al. (2010c) that development of elliptic approximations to $\mathcal{H}^{n-1}\left(\Gamma\right)$ could also be motivated from a more physical standpoint by considering the 1-dimensional case of an infinitely long bar of uniform cross section. Assuming that the bar is aligned with the $x$-axis and that a single crack fully cuts the bar at $x=a$, the phase-field profile corresponding to this sharp crack is none other than the discontinuous scalar function

[TABLE]

A regularized approximation of the above can then be made via a function $\phi\in\Phi$, where

[TABLE]

Some candidate functions are

[TABLE]

We note that the function (7) utilized by Miehe et al. (2010c) is in fact the solution obtained by minimizing the functional of Bourdin et al. (2000) in one dimension. On the other hand, (8) is a compactly supported function that leads to a variational inequality problem (Pham et al., 2011), while (9) is obtained by solving a 4th order governing equation and was introduced by Borden et al. (2012) with the aim of loosening the mesh size requirements with respect to $\ell$ at the same time taking advantage of numerical methods that can adequately model $C^{1}$-continuous solutions. Figure 1 shows a comparison of the three functions mentioned above. It can be observed that for a given value of $\ell$, the amount of crack diffusion is additionally dependent on the specific form of $\phi\left(x\right)$.

Higher order phase-field formulations for fracture have so far not achieved the same popularity as their lower order counterparts. This is mainly due to the higher order of continuity that they require in conjunction with numerical solutions, which is expensive to obtain with traditional frameworks such as finite elements. In addition, $\Gamma$-convergence is yet to be proven for these formulations. On the other hand, Li et al. (2015) point out that the incorporation of general anisotropic effects relating to the surface energy requires a formulation that is at least 4th order.

2.2.2 Damage

The connection between phase-field approaches and nonlocal damage models was explored by Pham et al. (2011), who noted that elliptic functionals approximating (2) can be seen as specific cases of the integral of a general state function pertaining to a gradient damage model:

[TABLE]

where $w\left(\phi\right)$ is a monotonically increasing function in the interval $\left[0,1\right]$ with $w\left(0\right)=0$ and $w\left(1\right)=w_{1}$. They suggest using the linear form $w\left(\phi\right)=w_{1}\phi$ which leads to models having a real elastic phase with no premature decrease in material stiffness, at the cost of solving a variational inequality problem for the damage evolution. In contrast, the quadratic form of $w\left(\phi\right)$ employed in (4) leads to material behavior having no real elastic phase, with damage already occurring at the onset of loading (Amor et al., 2009). Bulk energy release resulting from evolution of the phase-field is facilitated through a damage-dependent elasticity tensor $\mathbb{C}\left(\phi\right)$ as seen in (10). The simplest form which leads to isotropic behavior consists of the multiplicative ansatz

[TABLE]

in which $g\left(\phi\right)$ is the energy degradation function mentioned previously that is non-negative in the interval $\left[0,1\right]$ with $g\left(0\right)=1$ and $g\left(1\right)=g^{\prime}\left(1\right)=0$. Such form however has limited applicability, since it allows for unphysical compressive cracking based on $\mathcal{G}_{c}$. Lancioni and Royer-Carfagni (2009) adopted the above model to shear cracking by having $g\left(\phi\right)$ act only on the deviatoric portion of the strain. This was subsequently improved upon by Amor et al. (2009) who proposed that crack growth be driven also by the spherical part of the energy when the volumetric strain is positive. This results in more realistic anisotropic behavior where the material is allowed to crack in tension and shear, but not in compression. An alternative formulation was introduced by Miehe et al. (2010a, c) in which degradation occurs only on tensile components of the principal strain tensor, leading to pure mode-I cracks.

2.3 Going back to Griffith

The functional $\Psi\left(\bm{u},\phi\right)$ in (4) is neither linear nor convex, which makes the task of finding global minimizers non-trivial. Bourdin et al. (2000) proposed an alternate minimization algorithm that takes advantage of the convexity of $\Psi$ with respect to either $\bm{u}$ or $\phi$ when the other is held constant. Nonetheless, solution schemes based on descent algorithms only converge to local minimizers or saddle points and are by themselves inadequate for obtaining global minimizers. It was found that naive application of the alternate minimization often resulted in solutions that exhibited unphysical dips in the total energy, particularly when crack growth is brutal. In an attempt to remedy this behavior, a heuristic backtracking scheme was developed by Bourdin et al. (2008) (with subsequent improvements by Mesgarnejad et al. (2015)) based on an additional optimality condition that enforces monotonic evolution of the total energy when the load is also monotonically increasing. From a phenomenological standpoint however, the main objection to using global energy minimization is that it allows for evolutions where the current configuration jumps over arbitrarily large energy barriers in order to reach the new configuration corresponding to the global minimizer (Negri and Ortner, 2008). While this enables the strict preservation of energy conservation, it may also result in unphysical response where cracks propagate at lower energy release rates than $\mathcal{G}_{c}$ which violates Griffith’s criterion. Recently, Larsen (2010) introduced the notion of $\varepsilon$-stability as a stepping stone towards formulations that can predict crack paths based on local minimality, which is in turn closer to Griffith’s original idea. As with Griffith’s model, such solutions will also exhibit dissipation in the total energy in cases where crack propogation occurs in a brutal manner. However as pointed out by Negri and Ortner (2008), brutal cracking is primarily a dynamic phenomenon which explains why the total energy cannot be completely accounted for in a quasi-static framework.

2.4 On the treatment of $\ell$ as a material parameter

If we settle for invoking local versus global minimality, then classical solution schemes that were previously deemed inadequate are now robust without the need to perform backtracking. We are left with the non-convexity of Griffith’s functional, but this is easily dealt with via the alternate minimization algorithm as earlier mentioned. On the other hand, we end up as well with Griffith’s original conundrum concerning crack initiation. It turns out that the saving grace is none other than the regularization of the functional, which as will be discussed in the later sections actually allows for crack initiation in the absence of stress singularities provided that the characteristic length $\ell$ associated with the regularization is finite. This also brings us more in line with nonlocal damage theory where the thickness of the localization zone is usually a constant parameter associated with some physical internal length. In this context, the treatment of $\ell$ as a material parameter as suggested in Mesgarnejad et al. (2015) and Nguyen et al. (2016) makes a lot of sense since it can be shown that under the assumption of uniform damage, the peak stress $\sigma_{c}$ (which can be interpreted as the critical stress for crack nucleation) is dependent on $\ell$. However it has also been observed (e.g. Klinsmann et al., 2015) that propagation cracks explicitly modeled in the mesh exhibits a dependence on $\ell$ as well. Thus if $\sigma_{c}$ and $\mathcal{G}_{c}$ are assumed to be two independent material parameters, the interdependence of each with $\ell$ means that we can adjust the regularization based on one or the other but generally not both at the same time. Hence the idea of only relying on $\ell$ to calibrate the model is not entirely adequate. Furthermore in heterogeneous media, $\ell$ may take on a different value for each medium so that the diffuse crack becomes thicker or thinner as it passes from one material to the next. This is especially undesirable in a multiphysics setting in which the model for some overlying physical process depends directly on the phase-field. From (10) and (11) we can see that the degradation function is the remaining component through which we can rectify the model, and it is in fact this realization that has motivated the present work.

3 Governing equations and numerical implementation

For the remainder of this study, we have chosen to adopt the functional of Bourdin et al. (2000) by reason of its simplicity. Incorporating the work done by external forces, the regularized total potential energy for a given body $\Omega$ subject to boundary conditions is given by

[TABLE]

where $\psi\left(\bm{\varepsilon}\right)=\frac{1}{2}\bm{\varepsilon}\!:\!\mathbb{C}^{e}\!:\!\bm{\varepsilon}$ is the Helmholtz free energy density and $\partial\Omega^{t}$ denotes the part of the boundary for which Neumann (i.e. traction) conditions are prescribed. The quantity $\bm{b}$ represents the body force, while $\bm{t}$ is the vector of prescribed tractions acting on the Neumann boundary. By imposing the stationarity of $\Pi$ we obtain the variational equation

[TABLE]

The above equality must hold for arbitrary values of $\delta\bm{u}$ and $\delta\phi$, implying that

[TABLE]

These constitute the weak form of the governing equations. Noting that $\bm{\sigma}=\partial\psi/\partial\bm{\varepsilon}$, the equivalent strong formulation may be obtained by applying Gauss’ divergence theorem yielding the following coupled system:

[TABLE]

Equations 15a, 15b and 15c comprise the linear momentum equation and its corresponding boundary conditions, while (15d) is the phase-field evolution equation with the associated boundary condition given by (15e). As mentioned earlier, $\phi$ should go to zero away from the crack. Thus it is tacitly assumed that the domain is of sufficient size to provide adequate separation between the regularized crack and the boundary, allowing the phase-field to decay to values that are small enough to approximate this condition.

Extension of the irreversibility condition (3) to the regularized case is not immediately obvious, since the intermediate states $0<\phi<1$ do not have a straightforward physical interpretation. The natural course from the perspective of damage mechanics is to enforce the condition

[TABLE]

This can be imposed via an additional penalty term in the phase-field evolution equation (Miehe et al., 2010c), or alternatively through a history variable $\mathcal{H}$ that replaces the quantity $\psi\left(\bm{\varepsilon}\right)$ in (15d), defined as (Miehe et al., 2010a)

[TABLE]

A closer look at the phase-field localization process however reveals that (16) may not be the best extension of (3). For example in the 1-dimensional case Kuhn et al. (2015) demonstrate that localization of the phase-field at the center of a diffuse crack involves both the growth of $\phi$ near the crack tip, as well as a decrease of the same in adjacent regions, which enables the phase-field to correctly settle to the exponential profile given in (7). Based on this observed behavior, a strict imposition of (16) may lead to an overestimate of the crack length. Consequently, we introduce a modified version of (17) in which imposes irreversibility only when $\phi$ exceeds a certain threshold, i.e.

[TABLE]

The parameter $\phi_{c}$ represents the maximum value for damage that is allowed to heal during unloading. For a material point undergoing damage $\phi\leq\phi_{c}$, the resulting stress-strain curves will be nonlinear, with the amount of departure from linearity dependent on the specific form of the degradation function. Nonetheless having $\phi_{c}>0$ allows the stress paths for subsequent unloading and reloading to coincide with the initial loading curve, which cannot be achieved otherwise. For the present work, we have chosen to set $\phi_{c}=0.5$ in all of our simulations.

The coupled system described in (15) is implemented in a classical finite element framework, with the primary unknowns being the displacement $\bm{u}$ and phase-field $\phi$. In a 2-dimensional setting, these are expressed in terms of the corresponding nodal degrees of freedom as

[TABLE]

wherein

[TABLE]

with $N_{I}=N_{I}\left(\bm{x}\right)$ denoting the shape function associated with node $I$, and $\bm{u}_{I}$ and $\phi_{I}$ the respective displacement and phase-field degrees of freedom at node $I$. The strain and phase-field gradient are given by

[TABLE]

in which

[TABLE]

The former is the symmetrized gradient matrix associated with the Voigt form of the strain tensor. The test functions and their corresponding derivatives can be obtained from the above expressions by replacing $\bm{u}_{I}$ and $\phi_{I}$ with $\delta\bm{u}_{I}$ and $\delta\phi_{I}$ respectively. Due to the arbitrariness of the latter two quantities, numerical approximation of the weak form in (14) yields the following nonlinear system of equations at each node $I$:

[TABLE]

where $\mathcal{H}$ is the threshold-based history variable defined in (18). Due to the non-convexity of (12), the above coupled system is solved using the alternate minimization algorithm outlined in Bourdin et al. (2000). This involves cycling between (23a) and (23b): the linear momentum equation is solved first using values of $\phi$ from the previous iteration. Next, the phase-field evolution equation is solved using the newly obtained values of $\bm{u}$. The iterations are carried out repeatedly until the prescribed criteria on the size of residuals and inter-iteration corrections on the unknowns are met.

For degradation functions in which the derivative $g^{\prime}\left(\phi\right)$ is nonlinear, the subsystem represented by (23b) also becomes nonlinear, resulting in the need to perform nested iterations. In our experience, a naive implementation of the alternate minimization algorithm wherein the linearized phase-field equation is solved only once before going back to linear momentum is generally unstable and leads to incorrect results. That is, in the subsystem

[TABLE]

corresponding to the $\left(m+1\right)$th iteration within the $i$th time step, use of the exact Jacobian given by

[TABLE]

often produces an incorrect evolution of the phase-field and eventual blow-up. However, we have found that such an approach can be made to work by replacing $g^{\prime\prime}\left(\phi_{i}^{m}\right)$ with an approximate expression derived from low order terms in a Taylor expansion. Assuming that the degradation function is smooth, we can obtain $g^{\prime}\left(\phi_{1}\right)$ as an infinite sum of terms involving higher order derivatives of $g$ evaluated at $\phi_{2}\in\left[0,1\right]$:

[TABLE]

Now let $\phi_{1}=1$ so that $g^{\prime}\left(\phi_{1}\right)=0$. Dropping higher order terms as well as the subscript on $\phi_{2}$, we obtain

[TABLE]

which then gives us our approximation for the 2nd derivative of $g$:

[TABLE]

4 A curious case of crack nucleation

Pre-existing cracks maybe accounted for in the phase-field model in two ways. The first is through initialization of the phase-field profile, the second by modeling the crack faces directly as internal boundaries in the discretized geometry. Sicsic and Marigo (2013) provide analytical results showing that the growth of fully developed fractures described via the phase-field obeys Griffith’s law as the regularization parameter goes to zero. On the other hand it has been shown in numerical experiments that the same is not generally true for the extension of cracks built into the mesh. Recent studies (e.g Klinsmann et al., 2015; Nguyen et al., 2016) have observed that the simulated critical energy release rate (or analogously, the peak load) overshoots the correct value for sufficiently small $\ell$. This inconsistency becomes more understandable upon the realization that propagation of mesh-described cracks is actually a manifestation of nucleation rather than extension in the context of phase-field approaches, i.e. crack formation at a region where the phase-field is uniformly zero. Similar behavior can be observed in the case of crack nucleation at a notch or reentrant corner; the extension of mesh-modeled cracks is in fact a limiting case of the former where the notch angle is zero. Thus one can observe that the Francfort-Marigo-Bourdin phase-field model gives rise to three distinct types of simulated material response in connection with fracture: (a) propagation of phase-field-described fractures which is relatively well-understood, (b) crack nucleation in the absence of stress singularities which will be discussed in Section 5.4, and (c) quasi-nucleation behavior associated with the extension of mesh-modeled cracks, which is our immediate concern in this section.

4.1 Preliminary numerical example

To illustrate the dependence of the material response on the phase-field length scale, we simulate fracture propagation in a homogeneous specimen containing a center crack and subjected to tensile loading as shown in Fig. 2(a).

Assuming that $H$ is taken large enough such that the tensile stresses at the boundary are acceptably uniform, the mode-I stress intensity factor can be computed for finite values of the ratio $a/b$ as

[TABLE]

where $F\left(a/b\right)$ is a shape factor given by

[TABLE]

The above formula has a reported accuracy of $0.1\%$ or better for any $a/b$ (Tada et al., 2000). We note that $a/b=0$ and $H=\infty$ corresponds to the original fracture problem of Griffith (1921), for which $F\left(0\right)=1$. For the plane strain case, the critical stress intensity factor and strain energy release rate are related by

[TABLE]

Combining the above with (29), we obtain the following expression for the critical failure load:

[TABLE]

The actual computational domain is shown in Fig. 2(b) along with the relevant boundary conditions; due to symmetry, only half the geometry needs to be considered. Note that the initial crack $\Gamma$ of length $a$ in the computational domain is modeled as part of the geometry, and no initialization of the phase-field is performed to account for its presence. Actual dimensions used are $a=10$mm, $b=2a$ and $H=10a$. Thus $a/b=0.5$ and we obtain $F\left(a/b\right)=1.1862$ from (30). The material constants are $E=70,000$ MPa, $\nu=0.22$ and $\mathcal{G}_{c}=0.007$ N/mm. From the preceding equation, we obtain the critical failure load as $P_{c}=68.26$ N. The simulation is carried out using monotonic displacement control with the specimen gradually stretched in increments of $\Delta U=2.5\times 10^{-4}$ mm until failure occurs in the form of brutal cracking. In order to obtain a precise determination of the failure load, the increase in boundary displacements is carried out using smaller increments of $\Delta U=2.5\times 10^{-5}$ mm as failure is approached. Figure 2(c) shows the dependence of simulation results on the phase-field characteristic length as well as the relative mesh refinement, $\ell/h^{e}$. It can be seen that larger values of $\ell$ lead to an increase in deviation from linear behavior, whereas a smaller $\ell$ drives the peak load upwards. Furthermore there is an apparent lack of convergence with respect to the regularization parameter, since the load-displacement curve overshoots the true failure load when values of $\ell$ smaller than some threshold are used. As can be observed, the severity of this phenomenon is also influenced by the mesh refinement, and in particlar is greater for coarser meshes relative to $\ell$. We have found that different sets of material parameters give qualitatively the same behavior as what we have shown.

4.2 Exploring the overshoot phenomenon

The apparent lack of convergence in the material response with respect to $\ell$ in the above numerical example seems to contradict the $\Gamma$-convergence property of (4), however this can be explained by the fact that Griffith’s criterion does not actually involve the energy functional directly but rather its gradients (Fréchet derivatives). This nuance has no corresponding counterpart in image segmentation, and makes phase-field simulation of brittle fracture a fundamentally different problem from the former, despite the similarity of the Griffith energy to the Mumford-Shah functional. Thus while $\Gamma$-convergence of the regularized approximation to the sharp-boundary functional is by itself sufficient to produce physically meaningful results in an imaging context, this is no longer the case for brittle fracture.

At present, our understanding of the above phenomenon relies on numerical evidence obtained from analyzing problems such as the one presented in Section 4.1. To elucidate further, recall that for a material which fractures according to Griffith’s theory as summarized in (1), the following inequality applies with regard to energy increments:

[TABLE]

for some arbitrary small crack extension $\delta\Gamma>0$. Now $-\delta\Psi_{b}^{e}=G\delta\Gamma$ where $G$ is the energy release rate at the crack tip, so the strict inequality $-\delta\Psi_{b}^{e}<\mathcal{G}_{c}\delta\Gamma$ means that the crack must be stationary due to (1c). If $-\delta\Psi_{b}^{e}=\mathcal{G}_{c}\delta\Gamma$, then a positive $\delta\Gamma$ is admissible and the crack can propagate stably. On the other hand, the reverse inequality $-\delta\Psi_{b}^{e}>\mathcal{G}_{c}\delta\Gamma$ is generally understood as corresponding to brutal cracking. In a quasi-static framework where dynamic effects are disregarded, fracture propagation simply continues until a state is reached wherein the condition $\delta\Psi_{b}^{e}<\mathcal{G}_{c}\delta\Gamma$ is once again satisfied, resulting in arrest of the crack. It can be shown that in many cases, such a condition cannot be satisfied for any length of crack advance which results in the fracture cutting through the entire width of the domain.

Similar behavior is manifested by the evolution of $\phi$ in the diffuse-crack model during brutal crack propagation, and can observed by scrutinizing successive iterations within the relevant time step. In contrast to the original theory however, the phase-field model contains only the equality part of Griffith’s criterion in the phase-field evolution equation. That is, (14b) can be written as

[TABLE]

wherein

[TABLE]

for some positive $\delta\Gamma$ that arises from an arbitrary incremental evolution of the phase-field, denoted by $\delta\phi$. The implications of this are immediately obvious when on looks at the strong form of the phase-field equation in (15d): assuming that $g^{\prime}\left(\phi\right)<0$ everywhere except at $\phi=1$ (and this is in fact necessary for damage to evolve at all), then it is clear that $\phi$ must begin moving away from its initial value of 0 from the moment that nonzero stress is induced in the material. Furthermore, let $\eta_{b}$ denote the error arising from using $\Psi_{b}^{\text{app}}$ in place of $\Psi_{b}^{e}$, i.e.

[TABLE]

Plugging the above into (34) and writing $\delta\Psi_{b}^{e}$ in terms of $G$, we obtain the following relation:

[TABLE]

From the previous numerical example, we can infer that at some critical loading $U_{s}$ brutal propagation of the crack will occur, presumably because now $-\Psi_{b}^{\text{app}}>\mathcal{G}_{c}\delta\Gamma$ for any $\delta\phi$. The equivalent condition in terms of $\eta_{b}$ and $G$ is given by

[TABLE]

The different curves in Fig. 2(c) demonstrate how the actual value of $U_{s}$ depends on $\ell$. The key idea here is that both $\delta\Gamma$ and $\delta\eta_{b}$ are influenced by $\ell$, but in varying degrees from one another. In particular, it is no longer just the quantity $G-\mathcal{G}_{c}$ that determines the onset of brutal cracking; as can be observed from Fig. 2(c), both undershoot and overshoot of the correct failure load are possible. The challenge is to have (39) occur at the precise moment that $G$ exceeds $\mathcal{G}_{c}$, so that brutal fracture occurs at the correct magnitude of loading.

The current prevailing thought is that one can achieve the above scenario by some “correct” choice of the regularization parameter. However, the need to specify $\ell$ (and obviously $\ell>0$) brings into question the benefit of having regularized approximations $\Gamma$-converge to the Griffith energy at all as $\ell$ goes to zero. One can argue that the removal of such a requirement is not a disadvantage since it lends more flexibility to the phase-field framework and likewise opens the door to other interesting and more exotic approximations of (2), such as the higher order formulations by Borden et al. (2014) and Li et al. (2015) that have so far not been proven to be $\Gamma$-convergent to Griffith’s energy. More importantly, the main problem with relying on calibrating $\ell$ in order to obtain the correct instance of failure is that such a strategy is not guaranteed to succeed in all possible cases, in particular when the setup is very different from the one analyzed above. This is demonstrated in Section 6, where we study a problem for which the aforementioned technique does not work at all, at least within practical limitations.

4.3 Preserving linearity in the material response

An important consequence of (38) is that the material response of the regularized model inevitably drifts from linear elastic behavior prior to fracture, since growth of $G$ as a result of increasing $U$ must be matched by a corresponding increase in the incremental error term $\delta\eta_{b}$. Since $\eta_{b}$ represents the discrepancy between approximate and the exact bulk energies, an ever-increasing increment in the error term means that the simulated material behavior deviates further and further from linear elasticity with increasing $U$ as evident in Fig. 2(c). Some control on $\delta\eta_{b}$ can be exercised through the factor $\delta\Gamma$, i.e. we keep $\delta\eta_{b}$ small by keeping $\delta\Gamma$ small as well. However since $\delta\phi$ is arbitrary, we can accomplish this only by careful construction of either the degradation function (which affects the bulk energy), the crack length functional, or possibly both.

5 A new family of degradation functions

The ideas presented in Sections 4.2 and 4.3 can be combined together to give us a set of properties for what we would consider an accurate phase-field model with regard to the extension of mesh-described cracks:

The simulated critical displacement should preferably be close to the correct value, and 2. 2.

the accumulated error $\eta_{b}$ should be kept small prior to the occurrence of brutal fracture.

Item (b) is quite straightforward, and is achieved by having brutal fracture occur at low values of the phase-field. Such behavior is readily observed with the alternatives to quadratic degradation that have appeared in the literature, for instance the quartic function

[TABLE]

utilized by Karma et al. (2001) in conjunction with their own phase-field theory, and the cubic function

[TABLE]

analyzed by Borden (2012), in which the quantity $s$ controls the slope of the degradation function at the unbroken state. Kuhn et al. (2015) have shown that all three functions have similar post-failure behavior in stable crack growth, i.e. their differences lie primarily in the prediction of the level of strain or stress at which crack propagation occurs, and also in the amount of stiffness reduction observed prior to the onset of cracking.

Item (a) is more difficult to satisfy, in particular since quantities pertaining to the bulk energy are also dependent on material properties. The degradation function must then be parametric, in order to have the means of compensating for different values of these properties. We can see that none of the different functions mentioned above possess the latter property, so that one is instead forced to rely on tweaking $\ell$ as is done with the quadratic degradation function. From a conceptual standpoint this is not entirely satisfactory, since $\ell$ as a parameter belongs to the crack functional term and not the bulk energy. Furthermore a change in the regularization parameter leads to corresponding changes in both the bulk and surface terms. It is our view that it is better to introduce parameters directly into the degradation function. In doing so one is able to alter the behavior of $\delta\eta_{b}$ independently of $\delta\Gamma$. Furthermore the resulting formulation does not force the interpretation of $\ell$ as a de facto material parameter, but is rather nearer to the original concept of Bourdin et al. (2000) where $\ell$ is purely a mathematical construct that arises in connection with the regularization of discontinuities.

5.1 Exponential-type degradation

Consider now the family of degradation functions defined by the 3-parameter function

[TABLE]

where $k$, $n$ and $w$ are real numbers such that $k>0$, $n\geq 2$ and $w\in\left[0,1\right]$. The function $f_{c}$ is a corrector term whose role shall be explored in the later discussions. For now let us assume that $w=0$ so that (42) has effectively only 2 free parameters. The resulting expression has the following properties:

$g_{e}\left(\phi\right)$ is monotonically decreasing , 2. 2.

$g_{e}\left(0\right)=1$, $g_{e}\left(1\right)=0$, 3. 3.

$g_{e}^{\prime}\left(0\right)<0$, $g_{e}^{\prime}\left(1\right)=0$.

In choosing the form of (42) we have aimed for a minimal but sufficient number of parameters that allows us to have some control in the overall shape of the function in order to restore proper balance between bulk and surface energy increments, as well as suppress unphysical stiffness reduction prior to fracture. Note that one obtains the function $\left(1-\phi\right)^{n}$ in the limit as $k$ approaches 0, as shown in Fig. 3.

On the other hand, increasing $n$ has the effect of flattening $g_{e}\left(\phi\right)$ as $\phi$ goes to 1, shown in Fig. 4. The parameters $k$ and $n$ must be chosen such that crack propagation occurs at the right energy release rate for some given $E$, $\mathcal{G}_{c}$ and $\ell$. Of prime importance here is the shape of $g_{e}\left(\phi\right)$ at the vicinity of $\phi=1$ which controls the amount of elastic bulk energy that is dissipated at the diffuse crack tip. On the other hand, spurious stiffness reduction prior to fracture is connected to the behavior of $g^{\prime}\left(\phi\right)$ at $\phi=0$; we want to keep $g^{\prime}\left(0\right)$ small which is equivalent to setting $k$ to be large. However, choosing an excessively large value for $k$ also results in undesirable stress-strain behavior. In the following analysis, we show that it is possible to eliminate one parameter in (42) by selecting the largest values of $k$ (given some $n$) for which the resulting stress-strain relationships is considered acceptable.

5.2 Analytic model behavior in 1D

In order to study the effect the parameters $k$ and $n$ in our proposed family of functions, we take a look at the 1-dimensional case of a materially homogeneous bar with uniform cross section and length equal to $2L$. The bar is subjected to the boundary conditions $u\left(\pm L\right)=\pm u_{0}$ and $\phi^{\prime}\left(\pm L\right)=0$ as shown in Fig. 5.

Assuming zero body forces, the governing equations in (15) reduce to

[TABLE]

in which $\varepsilon=du/dx$, $\sigma=E\varepsilon$ and $\psi=\frac{1}{2}\sigma\varepsilon$. We focus on spatially homogeneous solutions for the phase-field, $\phi\left(x\right)\equiv\phi_{0}$ which implies that the stress is also spatially uniform, i.e. $\sigma\equiv\sigma_{0}=E\varepsilon_{0}$. As $\phi$ is no longer a function of $x$, (43b) simplifies to

[TABLE]

While it is physically more correct to express $\phi$ as a function of $\varepsilon$ (since crack formation is driven by the mechanical response), for complicated forms of $g\left(\phi\right)$ it becomes more convenient to adopt the opposite order of dependence. Hence we obtain

[TABLE]

with the corresponding derivative given by

[TABLE]

Consequently the derivative of the damaged-reduced stress can be obtained with respect to the phase-field as

[TABLE]

The effective stress-strain curve accounting for damage due to the phase-field can then be defined as

[TABLE]

Combining the last equation above with (45) and (46), we obtain after further manipulation the expression

[TABLE]

Now $g^{\prime}\left(0\right)<0$ by construction for $\phi<1$, and if $g\left(\phi\right)$ has monotonically increasing slope (i.e. $g^{\prime\prime}\left(\phi\right)\geq 0$) then the above expression is well defined for $\phi\in\left[0,1\right]$. However for degradation functions of the form given by (42), the existence of an inflection point means that the denominator in (49) may become zero at some point, implying the existence of a vertical tangent in the $\sigma$–$\varepsilon$ curve and possibly also snap-back behavior. This phenomenon is more pronounced for larger values of $k$, as illustrated in Fig. 6.

However it can be seen that a high value of $k$ also acts to suppress the undesired deviation from linear elastic behavior. Hence we want to choose this parameter as large as possible in order to minimize the said effect, but still small enough so as not to generate snap-back. This implies that $k=k\left(n\right)$, and the specific relationship is found by considering the limiting case where the denominator in (49) goes to zero. This yields the expression

[TABLE]

where

[TABLE]

with the relevant calculations given in A. Plugging the above results into (42) gives the reduced-parameter degradation function

[TABLE]

where again for the meantime we take $w=0$. The profile of $g_{s}\left(\phi\right)$ and its derivative are shown in Fig. 7 for several $n$.

Due to the fact that $g^{\prime}_{s}\left(0\right)<0$, growth of the phase-field takes place naturally in the presence of local stress gradients. This is in contrast to degradation functions where $g^{\prime}\left(0\right)=0$, for which special solution procedures are required to trigger the evolution of $\phi$ away from an undamaged state. Furthermore it has been shown (e.g. Kuhn et al., 2015) that for certain configurations of polynomial degradation functions, Eq. (45) may predict inadmissible values of the phase-field (e.g. $\phi\notin\left[0,1\right]$) at low strains implying a bifurcation-type behavior with $\phi$ remaining at the undamaged state until the point of bifurcation. However said point does not generally coincide with the onset of fracture, so that some stiffness reduction still occurs prior to the realization of peak loads. On the other hand, the family of degradation functions represented by Eq. (52) give rise to smooth $\varepsilon$-$\phi$ and $\sigma$-$\phi$ relationships as shown in Fig. 8.

It follows that for these type of functions, the material stress-strain behavior will exhibit elastic stiffness reduction, albeit in much more reduced magnitudes compared to the quadratic degradation function (see Fig. 9).

Likewise an important result is that for some given (fully determined) degradation function, the resulting normalized $\sigma$-$\varepsilon$ curve is unique so that the actual failure stress and strain are dependent on $\ell$ as well as the material parameters. This implies that there is no single degradation function that works for the entire range of values of $E$, $\mathcal{G}$ and $\ell$. Otherwise, the regularization parameter $\ell$ cannot be freely chosen but rather must be determined from the other material parameters in order to give the correct failure stress. The latter condition imposes a severe limitation on the phase-field method, particularly when viewed in the context of multi-physics simulations where one might desire to have control over the amount of crack regularization in order to satisfy requirements stemming from physics external to the mechanics and fracture propagation.

5.3 Role of $w$ and $f_{c}\left(\phi\right)$

Unfortunately, the simplified form of (52) with $w=0$ is not entirely adequate due to the fact that for higher values of $n$, the flattened shape of $g_{s}\left(\phi\right)$ means that near-total annihilation of the material stiffness already occurs at values of $\phi$ significantly less than 1. As a consequence, the phase-field stagnates below unity even though the material is fully damaged. As a result, calculation of crack lengths via evaluation of $\Gamma\left(\phi\right)$ will yield incorrect results. The above shortcoming can be remedied through $f_{c}\left(\phi\right)$, which acts as a correction term influencing how $g\left(\phi\right)$ goes to zero as $\phi\rightarrow 1$. Its purpose is to impart a residual gradient to $g_{s}\left(\phi\right)$ that is independent of $n$, so that $g_{s}^{\prime}\left(\phi\right)$ is always sufficiently below zero for $\phi<1$. A suitable expression satisfying the properties enumerated in Section 5.1 is

[TABLE]

In order to fully determine the constants $a_{2}$ and $a_{3}$, we impose two conditions. The first is that $f_{c}^{\prime}\left(\phi^{\star}\right)-\phi^{\star}f_{c}^{\prime\prime}\left(\phi^{\star}\right)=0$ in order to retain validity of expressions obtained based on $\phi^{\star}$ in A. The second is that $f_{c}\left(0\right)=1$. This yields the following expressions for the constants:

[TABLE]

We note that $f_{c}\left(\phi\right)$ itself is in general not a degradation function since for sufficiently large $\phi^{\star}$ it may be that $f_{c}\left(\phi\right)>1$ at certain values of $\phi$. Thus $w$ should be kept small, otherwise the correction term dominates. We have found that setting $w=0.1$ imparts a sufficient residual in the gradient of $g_{s}\left(\phi\right)$ while still satisfying that requirements given in Section 5.1. The resulting plots for $g\left(\phi\right)$ and $g^{\prime}\left(\phi\right)$ are shown in Fig. 10.

An additional benefit of having the correction term in the form of (53) is that for small values of $w$ the material response prior to fracture is closer to linear.

5.4 Fracture initiation based on tensile strength

The ability to initiate cracks in the absence of stress singularities requires the notion of strength in the form of a critical tensile stress $\sigma_{c}$ that is absent in the original theory of Griffith. Following the approach of Pham et al. (2011) and Bourdin et al. (2014), let us assume that this coincides with the peak stress in the stress-strain curve associated with the 1-d homogeneous-stress model described above. For the case of the quadratic degradation function, the peak stress is reached at a phase-field value of 0.25, leading to the relation (Nguyen et al., 2016)

[TABLE]

For the family of degradation functions defined by (52), the above expression is further dependent on $n$ as evident from Fig. 8(b). An explicit expression relating $n$ to the material parameters including $\ell$ is not easily obtained owing to the complicated form of (52). Instead we can utilize an approximate expression made by fitting a function to numerical evaluations of the peak stress for different values of $n$ as shown in Fig. 11.

This function is expressed in terms of the dimensionless quantity $\sigma_{nd}=\sigma_{c}\sqrt{\ell/E\mathcal{G}_{c}}$ and is of the form

[TABLE]

With the weighting factor $w$ set to 0.1, the resulting values of the coefficients $c_{0}$ to $c_{3}$ are as follows:

[TABLE]

It should be emphasized however that (56) much like (55) is valid only for the case where there are no stress gradients, and therefore has very limited applicability to cases where fracture nucleates from a stress concentration. Furthermore, these two equations do not account for the dependence that $n$ or $\ell$ must have on the mesh refinement when stresses are no longer uniform. On the other hand when stress concentrations are finite, it is straightforward to check via inspection of numerical results whether critical stresses have been exceeded, and thus model calibration in such a case is much easier compared to one where the fracture emanates from a stress singularity.

6 Numerical Examples

In this section we examine the performance of the proposed single-parameter degradation function relative to the conventional quadratic model through several examples. Our particular interest is in examining its ability to accurately capture the onset of fracture in the case of (a) a phase-field crack initiating at a location of stress singularity, and (b) one where a crack nucleates due to a nonsingular stress concentration reaching the prescribed material strength. The first numerical example is a recalculation of the problem presented in Section 4.1 using the new degradation function. It demonstrates how to determine the proper value of the parameter $n$ and also explores the effect of mesh refinement. The second example provides numerical evidence that the parameter tuning for $n$ becomes increasingly robust as the phase field parameter $\ell$ is reduced. The third example deals with a problem featuring strength-based crack initiation and also subsequent branching in a bi-material specimen. It highlights the need to carefully scrutinize numerical results and also the danger in blindly utilizing ready-made formulas for determining $\ell$ or $n$ which do not account for the specific local stress distributions in the problem at hand. In the final example, we investigate the new degradation function’s potential to accurately model the stable propagation of an initial crack that is explicitly modeled in the geometry.

Numerical computations were carried out within a finite element framework implemented in our in-house C++ code, which utilizes OpenMP to achieve shared-memory parallelization on a desktop machine having a multi-core processor. For all problems, the relevant domains are discretized using 3-node triangles having linear shape functions and assume plane strain conditions. Our code allows the combination of elements having a different number of primary unknowns, and this feature is utilized in some of the examples below. In such cases, additional boundary conditions have to be implemented at element interfaces in order to have proper closure of the governing equations. In using (52), we have set $w=0.1$ leaving $n$ as the sole free parameter subject to calibration/tuning. The coupled system of equations is solved using the alternate minimization algorithm, where we apply the linear approximation described at the end of Section 3 for the portion of the Jacobian matrix pertaining to the phase-field equations. With the aforementioned technique, very little difference is observed in computation times (e.g. number of iterations per step) between the simulations which utilize the quadratic degradation function and those which make use of our proposed alternative that is significantly more nonlinear.

6.1 Brutal crack propagation in center-cracked specimen

We revisit the brutal cracking problem of Section 4.1 involving a center-cracked specimen loaded in tension. As the analytical failure load is known for such a setup, it is useful not only for comparing the effect of our proposed single-parameter degradation function on the model behavior versus the original quadratic, but also as a means of calibrating the former by determining the proper value of $n$. Using the same specimen dimensions and material properties as before, along with a characteristic length of $\ell=0.5$ mm and critical mesh refinement ratio of $\ell/h^{e}=2$ (see Fig. 12 for detail of meshing in the crack vicinity), we resolve the problem utilizing our newly proposed degradation function given by (52) with $w=0.1$ as earlier recommended.

Initially, the prescribed upward displacement at the top boundary (see Fig. 2(b)) is increased using constant increments of $\Delta u_{\text{coarse}}=2.5\times 10^{-4}$ mm to determine the approximate displacement $u_{\text{crit}}$ at which failure occurs, after which the simulation is rerun with displacement increments refined to $2.5\times 10^{-5}$ mm between $u_{\text{crit}}+\Delta u_{\text{coarse}}$ in order to achieve higher precision in the simulated failure load. Fig. 13 shows the results obtained from using different values of $n$ between 4.5 and 6. The proper value of the degradation function parameter corresponding to the desired critical load of $P_{c}=68.26$ N is obtained via polynomial curve fitting, which yields a value of $n=5.314$. Incorporating this into the simulation produces a failure load of 68.38 N, representing an error of 0.18% with respect to the benchmark solution. While accuracy of the calculated load may be further improved by employing smaller $\Delta u_{\text{fine}}$ in addition to adjusting the value of $n$, the curve fitting procedure employed above nonetheless serves as a simple and straightforward means of achieving a reasonably accurate calibration of our proposed degradation function. An important property of (52) evident from Fig. 13 is that a higher value of $n$ always leads to lower simulated failure load. Plots of the load-displacement curves for different $n$ are shown in Fig. 14. We observe that the results are reasonably robust in terms of the exponent $n$ in that all choices of $n$ lead to an accurate representation of the linear regime prior to onset of fracture, in contrast to the classical quadratic degradation function. Furthermore, for the specific value of $\ell$ employed in the simulations, even an inaccurate calibration of $n$ having around 10% deviation from the optimal value still leads to a more accurate failure load than predicted by the quadratic model.

We also investigate the influence of the mesh refinement on the numerical results, as it is well known that the discretization of the domain close to the cracks must satisfy certain requirements on element sizes with respect to $\ell$ in order to properly resolve the exponential character of the phase-field. Specifically, $h<\beta\ell$ where $h$ is the length of element edges at the fracture vicinity and $\beta$ is a factor typically set to 1/2 in the literature based on results from Miehe et al. (2010c). However this estimate was based on a setup where the crack is aligned with element edges, allowing for the natural reproduction of the gradient discontinuity that occurs at $\phi=1$. In practice, the peak of the phase-field profile must occur at element Gauss points in order to effect a full degradation of the material stiffness. This implies that for constant gradient elements, this peak actually exists as a plateau of width $h$, which is an additional source of error when calculating the functional $\Gamma\left(\phi\right)$. Hence it may be necessary to choose a smaller value of $\beta$. Keeping the value of $\ell=0.5$ mm constant for the above problem, we determine $n$ for different values of the effective element size at the critical zone. The resulting plot is shown in Fig. 15. It can be observed that the change in $n$ becomes significantly smaller for $h\leq\ell/10$, indicating that we see numerical convergence with respect to the ratio $\ell/h$. Unfortunately, a full convergence study of $n$ with respect to mesh refinement is limited by the accuracy of the approximate analytical solution in equation (30).

6.2 Four-point bending test

For the second example, we simulate fracture propagation in a beam having an initial crack of length $a$ and subjected to four-point bending as shown in Fig. 16.

Our aim is to investigate the robustness of the degradation function parameter $n$ obtained in the previous section, by solving auxiliary problems that involve loading configurations fundamentally different to those in the main problem. To this end, we use the same values for the material parameters as given in Section 4.1. Likewise, we treat Example 6.1 as a prior calibration step.

The particular loading configuration investigated in this section produces a uniform internal moment between the inner applied loads, and by setting $a=10$ mm, $b=2a$ and $L_{1}=10a$ for the specimen dimensions we end with what is essentially the same computational domain as the previous example, albeit subjected to pure bending in the central beam portion of length $2L_{1}$. The internal bending moment at this region has a magnitude of $\left(L_{2}-L_{1}\right)P$, and for the current example we have set the moment arm $L_{2}-L_{1}$ equal to 50 mm. However since the loading consists of concentrated forces and support reactions, we have found it necessary to model as non-fracturing the beam portions where the forces are applied in order to avoid spurious damage evolution at these locations. The regions colored white in Fig. 16 indicate portions of the domain that are modeled as linear elastic with only the displacement field $\bm{u}$ as the primary unknown, whereas the gray region has both $\bm{u}$ and $\phi$. Thus a boundary condition for the phase-field must be specified at the interface between fracturing and non-fracturing regions. For the current example, this is the Neumann condition $\nabla\phi\cdot\bm{n}=0$, with $\bm{n}$ denoting the unit normal vector to the interface.

A semi-analytical solution for the critical moment corresponding to an energy release rate of $\mathcal{G}_{c}$ at the crack tip can be computed as

[TABLE]

where for pure bending the shape factor $F\left(a/b\right)$ has the form

[TABLE]

with a reported accuracy of 0.2% in the stress intensity factor $K_{I}$ for $a/b\leq 0.6$ (Tada et al., 2000). For the current specimen, $a/b=0.5$ and the above formula gives $F\left(a/b\right)=1.4945$. Plugging this into (58) yields a critical bending moment of 180.60 N-mm, which we designate as the benchmark solution for the problem.

Four simulation runs were carried out for comparison with the benchmark solution given above. For the first, we employ the standard quadratic degradation function with $\ell$ calibrated to a value of 0.94 by matching the simulated critical load to the benchmark solution for the center-cracked specimen (see Section 4.1). The second simulation run makes use of our proposed degradation function, where we have set $\ell=0.94$ in order to compare results of different degradation functions given the same regularization length scale. The corresponding value of $n$ for this case is found to be 5.26 based on calibration runs using the CC-specimen setup. In the third run, we set $\ell=0.5$ which allows us to directly use the result of Section 6.1. The fourth simulation uses $\ell=0.3$, with the obligatory calibration step yielding a value of 5.325 for the parameter $n$. In all four cases, the loading was applied in the form of prescribed downward displacements, first at increments of $\Delta U=-0.0025$ mm per step and then later refined to $-0.0001$ mm per step near the onset of crack propagation.

A summary of results for the four simulations is given in Table 1, where the relative error of a quantity $Q$ with respect to its benchmark value is computed as

[TABLE]

The corresponding load-displacement curves are shown in Fig. 17.

We can see that for the two runs with a coarse length scale of $\ell=0.94$, the relative errors for $M_{c}$ are significantly different from those for $P_{c}$ in the auxiliary problem used for calibration. These discrepancies show that the stress distributions around crack tips have a non-negligible influence on the model behavior, regardless of the form used for the degradation function. This is an unavoidable consequence of the diffuse approaches since the energy release rate at a crack tip is obtained via a nonlocal calculation. One should note that in this case the relative errors themselves are not definitive of a particular model’s accuracy since they are influenced by the size of load increments (i.e., time steps), and also because the model parameter can often simply be re-tuned to give better results although this latter step was not done in the current example. It is however evident from comparing the relative errors obtained during calibration and those for the main setup that parameter values are not automatically transferable from one problem to another, and that for the degradation function proposed in the present work such transferability is affected by the value of $\ell$ used (presumably in relation to the material parameters $E$ and $\mathcal{G}_{c}$).

In contrast, the results of the last three simulations seem to indicate that transferability of values for $n$ improves as $\ell$ is decreased. This is indeed an interesting outcome, and is consistent with the understanding that it is the diffuse representation of the crack tip which influences the energy release rates. Without drawing too broad conclusions from a single example, it appears that the new degradation function proposed here is less sensitive to calibration than the traditional degradation function. Nonetheless it is clear from this exercise that one must be careful in designing calibration procedures for any kind of degradation function, particularly when $\ell$ is large.

6.3 Crack initiation and branching

Our third example involves a bi-material specimen that is loaded in tension as shown in Fig. 18(a).

The material properties corresponding to the regions designated as $A$ and $B$ in the figure are given in Table 2, where it can be seen that material $B$ is stiffer than the other and is also non-fracturing.

We thus adopt the approach employed in the previous example: region $B$ is modeled as a linear elastic material with only displacement degrees of freedom, while in region $A$ we incorporate additional unknowns pertaining to the phase-field. In contrast to Section 6.2 however, here we impose the homogeneous Dirichlet condition $\phi=0$ on the interface separating between the two regions. This is done to ensure that the resulting phase-field profile is meaningful with respect to crack length calculations. Prescribed uniform vertical displacements of magnitude $U=0.05$ mm are applied at the top and bottom boundaries in increments of $\Delta U=0.001$ mm. We compare simulation results obtained from using our proposed new degradation function to that of the classical model employing quadratic energy degradation for two values of the phase-field regularization parameter, namely $\ell=1.25$ mm and $\ell=5$ mm. All four cases utilize the same discretization of the problem domain shown in Fig. 18(b), where the effective size of element edges along the anticipated path of crack propagation have been set to $h=0.4$mm. In addition a fifth simulation run was carried out with $\ell$ set to 0.31 mm on a finer discretization having $h=0.15$ mm; this corresponds to the case where failure occurs at the specified value fo $\sigma_{c}$ in connection with a quadratic degradation model. Force-displacement curves for the five cases are displayed in Fig. 19, while values of specific quantities of interest at crack initiation are listed in Table 3.

Due to the fact that boundary displacements are applied in constant increments without refinement near the instance of fracture initiation as done in the previous example, it is not possible to reproduce exactly the specified critical stress of 70 MPa during crack nucleation. Simulations 1 to 3 were carried out using the classical quadratic degradation function, while 4 and 5 utilize the new single-parameter degradation function given in (52). For the former, an estimate for the required magnitude of $\ell$ corresponding to $\sigma_{c}=70$ MPa may be obtained from (55). This yields $\ell=0.215$, however as observed from Fig. 19, the correct value of the regularization length for the quadratic case is nearer to 0.31. We also note that simulations 4 and 5 produce virtually identical results with respect to the peak load, demonstrating the ability of the proposed new degradation function to properly compensate for different magnitudes of crack regularization. It can be observed that past the initial crack formation which is represented by the sudden drop in the force-displacement curve, all simulations display essentially the same behavior. This is not surprising, since the for all cases, the initial crack traverses the entire width of region $A$, and so the subsequent residual force comes mainly from the resultant stresses in the non-fracturing part of the specimen as shown in Fig. 20.

Meanwhile, the final crack trajectories corresponding to $U=0.05$ mm obtained from simulations 2 to 4 are shown in Fig. 21.

We point the reader to a particular nuance of the current numerical example, namely that it is not immediately obvious simply from looking at the combined force-displacement plots in Fig. 19 which curve represents the correct specimen behavior under the given loading conditions.

An important insight can be found by examining the value of $\phi$ in the critical element at which the stress is maximum. From Table 3 we see that for simulation 1 this equal to 0.443 which corresponds to a degradation factor of $g_{2}\left(0.443\right)=0.310$. This means that just prior to fracture, the critical element has a stiffness of only slightly more than a third of its original value. This leads to a severe under-calculation of the critical stress, with the simulation reporting a value of $\sigma=71.02$ MPa at the critical element whereas a separate simulation assuming linear elastic behavior of the whole domain produces a stress of 106.1 MPa at the same location. This amounts to an overshoot of more than 50% of the true tensile strength. However since the damaged region comprises only a small fraction of the specimen’s total area (see Fig. 22), the linear elastic behavior exhibited by the undamaged region dominates the specimen response leading to deceptively small deviations in the force-displacement curves.

In contrast, simulation 4 has $\phi=0.1103$ at the critical element prior to failure. Coupled with the form of (52) that minimizes stress degradation for small values of the phase-field, we obtain $g_{s}\left(0.1103\right)=0.938$ for $n=4.4$ and $w=0.1$ which gives rise to much less distortion of the stress compared to the classical quadratic degradation function. Indeed at a displacement magnitude of $U_{c}=0.023$ mm, the phase-field model predicts a tensile stress of 69.71 MPa which is much closer to the value of 71.75 MPa obtained from assuming purely elastic material behavior. Additionally, we note that the values of $\ell$ and $n$ which lead to what may be considered as the “correct” model response in the case of using respectively the quadratic and exponential degradation functions are significantly different from the estimates obtained by using the formulas given in Section 5.4. This is due to the incompatibility between actual stress states at the crack initiation region for the current example (which are already localized prior to crack initiation) and the assumption of homogeneous stress pertaining to the 1-d case that was used in deriving the expressions in the aforementioned section.

6.4 Stable crack growth in a homogeneous medium

The problem of a rectangular specimen subjected to so-called surfing boundary conditions was initially used by Hossain et al. (2014) for studying the effective toughness of heterogeneous media and later adopted by Kuhn and Müller (2016) in the context of configurational forces. Details of the specimen geometry together with the initial crack are given in Fig. 23.

In both of the works mentioned, the phase-field profile is initialized such that $\phi=1$ at all points in the crack locus, decaying with the proper gradients towards zero away from the crack. In contrast for the current example, no such initialization is carried out in order to simulate the transition from crack initiation at a location of stress singularity towards propagation of a fracture that is fully described by the phase-field. The Dirichlet boundary conditions are derived from a $K_{I}$-controlled displacement field corresponding to a crack under mode-I loading, given by

[TABLE]

in terms of polar coordinates $r$ and $\theta$, with the crack extending infinitely along the line $\theta=\pi$ from a tip located at $r=0$. The quantity $\kappa$ is Kolosov’s constant which is equal to $3-4\nu$ for the case of plane strain. Crack propagation is achieved by translation of the above coordinate system with respect to the original configuration of the specimen resulting in the horizontal motion of the crack tip. Letting $x_{K_{I}}\left(t\right)=vt$ and $y_{K_{I}}\left(t\right)=0$ be the Cartesian coordinates of the crack tip for some fictitious time $t$ and positive constant $v$, we obtain

[TABLE]

For simplicity, we have chosen to let $v=1$. The material properties used for the specimen are $E=210$ GPa, $\nu=0.3$ and $\mathcal{G}_{c}=2.7$ N/mm. Finally, we set $K_{I}$ to a constant value of $\sqrt{E\mathcal{G}_{c}}$ and run the simulation from $t=5$ up to $t=30$ in increments of $\Delta t=0.5$. The analytical response of the specimen can be understood as follows: for $t\in\left[5,10\right)$, the energy release rate at the crack tip is smaller than $\mathcal{G}_{c}$, so that the crack does not propagate. At $t=10$, this quantity is exactly equal to $\mathcal{G}_{c}$, allowing the crack to growth. Henceforth for $t>10$, the crack tip moves to the $K$-field center denoted by $x_{K_{I}}\left(t\right)$.

To gain insight on the numerical behavior of the fracturing specimen, we conduct a preliminary simulation assuming plain linear elastic response without fracturing, the results for which are shown in Fig. 24(a).

The energy release rate at the crack tip is obtained by calculating the $J$-integral over the contour defined by the specimen boundary, $\partial\Omega$. One can see that this is underestimated in the numerical solution, i.e. the $J$-integral is less than $\mathcal{G}_{c}$ at $t=10$. Consequently, location of the crack tip predicted by the numerical simulation lags behind the analytical location as illustrated by the black and gray dashed lines in Fig. 24(a).

We now examine phase-field model behavior in connection with the quadratic degradation function for three different values of the regularization parameter, namely $\ell=1$ mm, 3 mm and 5 mm. Again this is done by looking at two quantities of interest: the energy release rate at the crack tip represented by the $J$-integral, and the length of crack extension described by the phase-field that is obtained by evaluating the functional $\Gamma\left(\phi\right)=\int_{\Omega}\left(\frac{1}{2\ell}\phi^{2}+\frac{\ell}{2}\nabla\phi\cdot\nabla\phi\right)\,\mathrm{d}\Omega$. The results for different values of $\ell$ are summarized in Fig. 24(b) and exhibit similar behavior. We first observe a zone of premature crack growth where the crack length is seen to increase at rates less than $v$. This preliminary growth is not related to any physical extension of the crack but is in fact due to the evolution of the phase-field profile representing the diffuse crack tip, illustrated in Fig. 25.

This is followed by brutal cracking represented by a sudden increase in the crack length (see Fig. 26), after which the fracture grows stably at a rate more or less equal to $v$.

As expected, the numerical location of the crack tip trails the $K$-field center associated with the applied boundary condition. In addition, the critical energy release rate is overestimated in the region of stable crack growth, i.e. $\mathcal{G}_{c}^{\text{num}}>\mathcal{G}_{c}$. This phenomenon has been previously reported in the literature, and we refer the reader to Hossain et al. (2014) and Kuhn and Müller (2016) for more detailed discussions on the matter. Our main focus at the moment is the overshoot that occurs in the $J$-integral prior to the onset of fracture, which results in even further delay of the actual crack extension. Such behavior is obviously unphysical, and more so does not occur when the initial crack is described by the phase-field as demonstrated in various numerical examples from the aforementioned literature. It is our belief that this artifact is heavily dependent on the specific form of the degradation function. More importantly in this case, the overshoot does not decrease with smaller $\ell$. The results shown in Fig. 24(b) provide evidence that one cannot in general rely on the strategy of calibrating $\ell$ in order to obtain correct model behavior, and furthermore casts doubt on the notion that $\ell$ should be viewed as a material parameter, particularly in connection with the reproduction of $\mathcal{G}_{c}$. As previously mentioned, the overshoot of the critical energy release rate results in a delay of actual crack extension such that by the time it occurs, there is an excess in bulk energy that must be dissipated. Upon the onset of fracture, instantaneous catch-up growth occurs resulting in a finite increase of the crack length as shown in Fig. 26.

On the other hand, simulations utilizing the proposed single-parameter degradation show behavior that more closely reflects the physics as shown in Fig. 27.

As with the previous results pertaining to quadratic degradation, we can observe that for the same value of the degradation parameter $n$, varying the magnitude of $\ell$ has little effect on the amount of spurious overshoot in the energy release rate prior to crack extension. Rather, the parameter $n$ itself is effective in controlling this feature, and thus can be calibrated such that crack extension occurs when the $J$-integral reaches a value of $\mathcal{G}_{c}^{\text{num}}$. Additionally the phase-field remains at very low values before the onset of fracture, resulting in negligible increase of the crack length prior to the actual onset of crack growth as seen in Fig. 28.

Nonetheless, we can observe from Fig. 27 that setting $n$ too low results in an overshoot behavior similar to the case of quadratic degradation, while for sufficiently high values of $n$ a dip occurs in the $J$-integral following onset of crack extension. The latter is also a numerical artifact, with an underlying mechanism that is converse to what occurs for overshooting. That is, $\phi$ experiences a jump in value at the crack tip around which an exponentially decaying profile is enforced by the evolution equation for the phase-field. This leads to the crack extension being too big, so that now the numerical crack tip may be understood to have jumped ahead of $x_{K_{I}}$. The result is a virtual unloading at the crack tip vicinity evidenced by the decrease in maximum tensile stress shown in Fig. 29.

However since the crack tip diffusion is controlled by the phase-field length scale, the aforementioned dip may be reduced by simply decreasing the magnitude of $\ell$.

7 Concluding remarks

In this paper, we have introduced a novel family of energy degradation functions aimed at overcoming major drawbacks of the standard phase-field model in simulating fracture nucleation in brittle materials. A key feature of these functions is their dependence on a set of parameters, permitting us to effect minute changes to their shape. This allows for a more detailed study on how the form of the degradation function influences the phase-field model response independent of the regularization parameter $\ell$. Of particular interest is the discovery that use of the standard quadratic degradation function leads to a delay in the onset of crack propagation, leading to an overshoot in predicted critical loads which cannot be ameliorated by adjusting the value of $\ell$. This finding is remarkable, since such a strategy was previously thought to be adequate for recovering correct failure loads based on prior numerical examples found in the literature. Furthermore, it casts doubt on the idea that $\ell$ should be viewed as a material parameter, at least for the case involving brittle fracture of linear elastic materials. On the other hand with the proposed family of degradation functions, it is possible to obtain significantly more accurate simulations provided that proper calibration of the function parameters is carried out. The computational overhead resulting from the consequent nonlinearity of the phase-field evolution equation can be virtually eliminated by employing suitable linear approximations for the tangent matrices which then allows for straightforward application of the alternate minimization algorithm.

An important consideration for the proposed family of degradation functions is the actual number of independent parameters that must be specified, since this directly affects the difficulty or ease of calibration. In this paper we have chosen to work with a function that has only one parameter to be calibrated out of an initial four, in the belief that more would render the model unappealing for use in an industry setting. Consequently, we do not take full advantage of the flexibility of our model. Furthermore the elimination of extra, unwanted parameters was done based on a rationale that prioritized the preservation of linear elastic response prior to fracture. Looking at results of the numerical examples we can see that this objective has been sufficiently accomplished, however the price to pay is a spurious dip in the bulk energy after the initial crack nucleation which occurs even with proper calibration as seen in the surfing problem. In the current model, this can only be alleviated by reducing the phase-field regularization which in turn increases computational expense due to meshing requirements along the crack trajectories. As an alternative, one can allow damage to occur gradually in the vicinity of the crack nucleation point prior to failure, however this requires a degradation of the bulk energy to preserve energy balance and runs counter to the rationale mentioned above. It is thus outside the scope of the present paper, and will be explored in a future work.

Acknowledgements

This work was funded by the Research Council of Norway through grant no. 228832/E20 and Statoil ASA through the Akademia agreement. The authors likewise express their gratitude to Prof. Blaise Bourdin for suggesting the use of the surfing boundary problem as a test case for model comparison.

Appendix A Derivation of $k\left(n\right)$

Let $g\left(\phi\right)$ be defined according to Eq. (42). The corresponding first and second derivatives are given by

[TABLE]

We are interested in the limiting scenario where Eq. (49) becomes infinite, i.e.

[TABLE]

Additionally, we will consider only the case where $\phi<1$, since setting $\phi=1$ results in a zero numerator in Eq. (49). Plugging Eqs. (63) and (64) into the above expression gives

[TABLE]

which simplifies to

[TABLE]

Solving for $k$ in the above equation, we obtain

[TABLE]

The behavior of $k\left(\phi,n\right)$ is shown in Figure 30 for several values of $n$.

We are interested in the minimum possible value of $k$ for each given $n$, hence the relevant condition to consider is

[TABLE]

Assuming that the denominator does not equal zero, the above equation reduces to

[TABLE]

whereupon we obtain the positive root

[TABLE]

via the quadratic formula. The above result is applicable for $n\neq 2$. For the case where $n=2$, Eq. (67) becomes

[TABLE]

Proceeding similarly to the previous case, we have

[TABLE]

so that the solution is

[TABLE]