Phase-Field Models for Fatigue Crack Growth

TL;DR

This paper presents a phase-field modeling approach for fatigue crack growth that captures sub-threshold fracture growth under cyclic loads, reproduces the Paris law, and can simulate complex crack geometries.

Contribution

The paper introduces a phenomenological phase-field model incorporating cyclic degradation to simulate fatigue crack growth, including multiple cracks and complex geometries.

Findings

Reproduces the Paris law with high exponents for brittle fatigue growth.

The Paris law exponent decreases with Ginzburg-Landau type phase dynamics.

Model can simulate crack growth in complex geometries in 2D and 3D.

Abstract

We introduce a class of models based on near crack tip degradation of materials that can account for fracture growth under cyclic loads below the Griffith threshold. We incorporate the gradual degradation due to a cyclic load through a flow equation that decreases spatially varying parameters controlling the fracture toughness in the vicinity of the crack tip, with the phase and displacement fields relaxed to an energy minimum at each time step. Though our approach is phenomenological, it naturally reproduces the Paris law with high exponents that are characteristic of brittle fatigue crack growth. We show that the exponent decreases when the phase field dynamics is of the Ginzburg-Landau type with a relaxation time comparable to the cyclic loading period, or when degradation occurs on a scale larger than the process zone. In addition to reproducing the Paris law, our approach can be used to model the growth of multiple cracks in arbitrarily complex geometries under varied loading conditions as illustrated by a few numerical examples in two and three dimensions.

Tables

Table 1. Table 1: The En-passant simulation parameters

$L_y/L_x$	$L_x/L_c$	$L_x/W$	$L_x/\xi$	$\nu$	${\delta }_{\max }/\sqrt{G_c\xi /E}$
$2$	$2.5$	$20$	$200$	$0.1$	$14$

Table 2. Table 2: Simulation parameters for 3-D penny-shaped cracks coalescence

$L/\xi$	$R_c/\xi$	$b/\xi$	$\nu$	${\gamma }_{\min }$	${\delta }_{\max }/\sqrt{G_c\xi /E}$
$50$	$5$	$20$	$1/3$	$0.5$	$6.36$

Table 3. Table 3: Simulation parameters for the 3-D inclined penny-shaped crack

$L/\xi$	$R_c/\xi$	$\nu$	$\theta$	${\gamma }_{\min }$	${\delta }_{\max }/\sqrt{G_c\xi /E}$
$50$	$5$	$1/3$	$\pi /3$	$0.1$	$5.66$

Equations

$$\mathcal{F}(u,\Gamma ):={\int }_{\Omega \setminus \Gamma }W(u)dV-{\int }_{{\partial }_N\Omega }\left(f\cdot u\right)ds+G_c{\mathcal{H}}^{n-1}(\Gamma )$$

$${\mathcal{F}}_{\xi }(u,\varphi ):=\mathcal{B}(u,\varphi )+{\mathcal{S}}_{\xi }(\varphi )$$

$${\mathcal{B}}_{\xi }(u,\varphi ):={\int }_{\Omega }g(\varphi )W(u)dV-{\int }_{{\partial }_N\Omega }\left(f_i{u}_i\right)ds$$

$${\mathcal{S}}_{\xi }(\varphi ):=\frac{G_c}{4C_{\varphi }}{\int }_{\Omega }\left(\frac{w(\varphi )}{\xi }+\xi {\left|\nabla \varphi \right|}^2\right)dV$$

$${\mathcal{S}}_{\xi }(\varphi ):={\int }_{\Omega }\left({\mathcal{E}}_{c0}w(\varphi )+\frac{\kappa }{2}{\left|\nabla \varphi \right|}^2\right)dV$$

$$({u}^{*},\Gamma^{*})=\underset{\scriptstyle{\begin{array}[]{l}{u}\in\mathcal{A}_{u}\\ \Gamma\in\mathcal{A}_{\Gamma}\end{array}}}{\operatorname*{\mathrm{argmin}\,}}\mathcal{F}({u},\Gamma)$$

$$({u}^{*},\phi^{*})=\underset{\scriptstyle{\begin{array}[]{l}{u}\in\mathcal{A}_{u}\\ \phi\in\mathcal{A}_{\phi}\end{array}}}{\operatorname*{\mathrm{argmin}\,}}\mathcal{F}_{\xi}({u},\phi)$$

$$\frac{\partial {\mathcal{F}}_{\xi }}{\partial u}(\tilde{u})=0\forall \tilde{u}\in {\mathcal{A}}_u$$

$$\frac{\partial {\mathcal{F}}_{\xi }}{\partial \varphi }(\tilde{\varphi })=0\forall \tilde{\varphi }\in {\mathcal{A}}_{\varphi }$$

$${\tau }_{\varphi }\frac{\partial \varphi }{\partial t}=-\frac{1}{{\mathcal{E}}_{c0}}\min \left(0,\frac{\delta {\mathcal{F}}_{\xi }}{\delta \varphi }\right)$$

$${\tau }_{\varphi }{\int }_{\Omega }\left(\frac{\partial \varphi }{\partial t}\tilde{\varphi }\right)dV=-\frac{1}{{\mathcal{E}}_{c0}}\frac{\partial {\mathcal{F}}_{\xi }}{\partial \varphi }(\tilde{\varphi })\forall \tilde{\varphi }\in {\mathcal{A}}_{\varphi }$$

$${\mathcal{S}}_{\xi }(\varphi ;{\mathcal{E}}_c)$$

$$\frac{1}{{\chi }_c}\frac{\partial {\mathcal{E}}_c}{\partial t}=\min \left(0,-H({\mathcal{E}}_c-{\mathcal{E}}_{c\min })w(\varphi )\frac{\partial }{\partial t}\left(g(\varphi )W^{+}(u)\right)\right)$$

$$\mathrm{H}(x)=\left\{\begin{array}[]{rcl}1&&x>0\\ 0&&x\leq 0\end{array}\right.$$

$$W^{+}({u})=\left\{\begin{array}[]{rcl}W({u})&&e_{kk}({u})\geq 0\\ 0&&\mathrm{otherwise}\end{array}\right.$$

$$G_c=2\sqrt{2\kappa {\mathcal{E}}_c}{\int }_0^1\sqrt{w(\varphi )}d\varphi =2\sqrt{2\kappa {\mathcal{E}}_c}{C}_{\varphi }$$

$$\frac{d{\mathcal{F}}_{\xi }}{dt}(u,\varphi ,{\mathcal{E}}_c)$$

$$={\int }_{\Omega }\left(w(\varphi )\frac{\partial {\mathcal{E}}_c}{\partial t}\right)dV\le 0$$

$$\frac{d{\mathcal{F}}_{\xi }}{dt}(u,\varphi ,{\mathcal{E}}_c)$$

$$={\int }_{\Omega }\left(-\frac{1}{2{\mathcal{E}}_{c0}{\tau }_{\varphi }}{\left(\frac{\delta {\mathcal{F}}_{\xi }}{\delta \varphi }\right)}^2+w(\varphi )\frac{\partial {\mathcal{E}}_c}{\partial t}\right)dV\le 0$$

$${\mathcal{S}}_{\xi }(\varphi ;\gamma ):=\frac{G_c}{4C_{\varphi }}{\int }_{\Omega }\gamma \left(\frac{w(\varphi )}{\xi }+\xi {\left|\nabla \varphi \right|}^2\right)dV$$

$$\frac{1}{{\chi }_c}\frac{\partial \gamma }{\partial t}=\min \left(0,-\mathrm{H}(\gamma -{\gamma }_{\min })w(\varphi )\frac{\partial }{\partial t}\left(g(\varphi )\frac{W^{+}}{{\mathcal{E}}_{c0}}\right)\right)$$

$$\frac{d{\mathcal{F}}_{\xi }}{dt}(u,\varphi )$$

$$={\int }_{\Omega }\left(\frac{G_c}{4C_{\varphi }}\left(\frac{w(\varphi )}{\xi }+\xi {\left|\nabla \varphi \right|}^2\right)\frac{\partial \gamma }{\partial t}\right)dV\le 0$$

$$\frac{d{\mathcal{F}}_{\xi }}{dt}(u,\varphi )$$

$$={\int }_{\Omega }\left(-\frac{1}{{\mathcal{E}}_{c0}{\tau }_{\varphi }}{\left(\frac{\delta {\mathcal{F}}_{\xi }}{\delta \varphi }\right)}^2+\frac{G_c}{4C_{\varphi }}\left(\frac{w(\varphi )}{\xi }+\xi {\left|\nabla \varphi \right|}^2\right)\frac{\partial \gamma }{\partial t}\right)dV\le 0$$

$${\bar{\mathcal{F}}}_{\xi }(u,\varphi ):=\frac{{\mathcal{F}}_{\xi }(u,\varphi )}{G_c{\xi }^{n-1}}=$$

$$-$$

$${{\mathcal{E}}_c}_{n+1}=\max \left({\mathcal{E}}_{c\min },\min \left({{\mathcal{E}}_c}_n,{{\mathcal{E}}_c}_n-{\chi }_c\left(g({\varphi }_{n+1})W^{+}(u_{n+1})-g({\varphi }_n)W^{+}(u_n)\right)\right)\right)$$

$${\gamma }_{n+1}=\max \left({\gamma }_{\min },\min \left({\gamma }_n,{\gamma }_n-{\chi }_c\left(g({\varphi }_{n+1})W^{+}(u_{n+1})-g({\varphi }_n)W^{+}(u_n)\right)\right)\right)$$

$${\int }_{\Omega }\left(\frac{{\varphi }_t-{\varphi }_{t-1}}{\delta t}\right)\tilde{\varphi }dV=-\frac{\partial {\mathcal{F}}_{\xi }}{\partial \varphi }(\tilde{\varphi })\forall \tilde{\varphi }\in {\mathcal{A}}_{\varphi }$$

$$\frac{1}{{\chi }_c}\frac{\partial \gamma }{\partial t}=\min \left(0,-\mathrm{H}(\gamma -{\gamma }_{\min })\mathrm{H}(W^{+}-{\sigma }_c)\frac{\partial }{\partial t}\left(g(\varphi )\frac{W^{+}}{{\mathcal{E}}_{c0}}\right)\right)$$

$$r_c=\frac{1}{2\pi }{\left(\frac{K_{IC}}{{\sigma }_c}\right)}^2$$

$$\gamma =\max ({\gamma }_{\min },1-{\chi }_c{e}_{eq}^p)$$

Full text

1 Introduction

More than one and a half centuries after the first systematic investigations by August Wöhler [60], fatigue fracture remains the major mode of failure of mechanical components. The pioneering experimental observations of Paris [47] have offered a concrete connection between linear elastic fracture mechanics (LEFM) and fatigue crack propagation for long cracks under constant Mode-I loading. This is further evidenced by the large number of articles offering ad hoc modifications of the power law [46, 24, 58, 45]. In this article, we aim to present a phenomenological continuum model for fatigue crack growth within the framework of the phase-field approach to fracture. This model is shown to naturally reproduce the Paris law $da/dN\sim\Delta K^{m}$, which relates the crack growth per cycle to the range of variation $\Delta K$ of the mode-I stress intensity factor during cyclic loading. In addition, this model allows us to perform component level simulations of fatigue crack growth.

Paris-law type scaling has been widely observed experimentally both in ductile materials [52, 56] where the exponent $m\sim 2\textup{--}4$ and in brittle materials [13, 54] where $m\sim 10\textup{--}40$ is typically much higher. Since its inception, many attempts have been made to explain the different observed exponents of the Paris power-law. Rice [51] shows that simple continuum level analysis of fatigue crack growth based on perfect plasticity also predicts $m=2$. Simple dimensional analysis with complete self-similarity assumption also predicts a second order Paris type power-law (*i.e., *$m=2$), whereas higher $m>2$ exponents can also be justified using incomplete self-similarity arguments [53, 11]. Classical continuum and meso-scale simulations have also been used to obtain $m$. Bazant and Hubler [12] obtain $m=4$ through closed-form analysis of cyclic creep in quasi-brittle materials. Deshpande et al. [22] obtain $m\sim 4-8$ from meso-scale simulations of discrete dislocation dynamics and a simple straight cohesive crack under pure mode-I loading. The interested reader can find a comprehensive review of different theoretical models of fatigue fracture in [30].

What distinguishes our class of models is not only that the Paris law with a variable exponent emerges naturally under simple mode-I loading, but also that they enable direct numerical simulation of fatigue crack growth without a need for micro/meso-scale simulations and a priori assumptions. Within the phase-field framework, Amor et al. [8] developed a fatigue model by introducing a dissipation energy that depends on crack length. This model has the advantage of directly modeling the Paris law with a prescribed exponent, but is limited to describing the growth of a single crack. Nguyen et al. [44] developed a cohesive model of fatigue where the cohesive stiffness at the crack tip was reduced as a function of load rate. Abdelmoula et al. [1] introduced a similar cohesive law where the accumulation of the opening displacement was tracked by a counter variable, thereby allowing for the fatigue crack to propagate under any load. They further showed that a Paris power-law $m=4$ emerges naturally from their formulation for small loads. More recently, Alessi et al. extended this type of counter variable approach in a one-dimensional phenomenological phase-field model of fatigue [3]. Here we develop a class of phase field models that phenomenologically describe fatigue crack growth by degradation of the local fracture toughness as a result of cyclic loading following a qualitatively similar approach as Abdelmoula et al. [1] in a cohesive zone framework. The model has the advantage of describing the growth of a collection of interacting cracks in complex geometries under arbitrary loads, including fatigue and continuous growth below and above the Griffith threshold, respectively.

We can justify our phenomenological assumptions by using the experimental observations of fatigue in brittle materials. Even though dislocation processes that underlie ductile fatigue in crystalline solids are minimal or absent in brittle materials, failure under cyclic sub-critical loads has been widely reported for nominally brittle materials such as alumina, TZP, silicon nitride, and silica glasses as well as in ceramics and polymers (*e.g., *see [55] for a complete survey of fatigue in ceramics). This crucial difference is also highlighted by examining postmortem crack surfaces. Unlike in ductile fatigue of metals where there are distinct fractographical differences between the fracture surfaces of cracks in quasi-static monotonic failure and the surface of fatigue cracks, in most brittle solids there are few discernible features due to fatigue crack propagation [21]. Nevertheless, fatigue fracture propagation in brittle materials can be explained by near-crack-tip phenomena which facilitate the sub-critical crack propagation. For example, while bulk silicon is not susceptible to fatigue failure, cyclic failure was observed in micron scale single crystal silicons [20, 31, 4] due to the creation of a silica layer which thickens at high stress sites and is prone to environmentally assisted corrosive cracking [42, 5].

Our work builds on a growing body of research aimed at expanding the phase-field approach to fracture. Since their inception [25, 16, 32], phase-field models of fracture have proven to be powerful tools in modeling a wide variety of fracture phenomena. These models have been validated by theoretical analyses [29] and comparisons of predicted and observed crack paths in non-trivial geometries [40]. They have been used to reproduce complex experimental observations in brittle fracture including thin-film fracture [39], thermal fracture [17], mixed mode fracture [19], and dynamic fracture [19, 18, 37], and fracture in colloidal systems [49] to give only a few examples. The approach has also been extended to ductile fracture [41, 6, 15].

This article is organized as follows: we first review the classical phase-field models of brittle fracture in section 2. In section 3, we introduce two models of brittle fatigue fracture. We present our numerical simulation algorithms and strategy in section 4. We show in section 5.1 how the Paris law naturally emerges from our models. We then show two sets of simulations for complex crack paths in 2-D in section 5.2 and 3-D in section 5.4. Finally in section 6 we show how extending the degradation zone beyond the fracture process zone results in smaller Paris law exponents.

2 Phase-field models of brittle fracture

Brittle fracture is signified by the Griffith criterion that states that the crack propagates if and only if the energy release rate $G$ (*i.e., *the derivative of the stored elastic energy with respect to crack length under constant loading) from the tip of the crack is larger than the fracture toughness of the material $G_{c}$.

We define the free energy $\mathcal{F}$ by combining the elastic (bulk) and fracture (surface) energies for an elastic domain $\Omega\subset\mathbb{R}^{n}$ containing fracture set $\mathbb{R}^{n-1}\supset\Gamma\subset\Omega$ and due to deformation ${u}=\{u_{i}\}\,,i=1\dots N$ as

[TABLE]

where $W({u})=({\mathcal{C}}_{ijkl}\,e_{kl}({u})e_{ij}({u}))/2$ is the elastic energy density, $e_{ij}({u})=(u_{i,j}+u_{j,i})/2$ is the linear strain (*i.e., *symmetric gradient) of the displacement vector ${u}$, where for Lame parameters $\lambda,\mu$ we write the fourth order linear elastic constitutive tensor (*i.e., *Hooke’s law) as ${\mathcal{C}}_{ijkl}=\lambda\delta_{ij}\delta_{kl}+\mu(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{kj})$. Furthermore, $\mathcal{H}^{n-1}(\Gamma)$ is the Hausdorff measure of the crack set *i.e., *the aggregate length of cracks in two dimensions or aggregate surface area of cracks in three dimensions. The classical elasticity boundary conditions are given as external forces $f$ acting on part of its boundary $\partial_{N}\Omega\subset\partial\Omega$ and Dirichlet condition ${u}\equiv\hat{{u}}$ on $\partial_{N}\Omega\subset\partial\Omega\setminus\Gamma$.

The phase-field method approximates the above energy by replacing the crack set $\Gamma$ with a phase field $\phi\in[0,1]$, where $\phi=1$ almost everywhere and converges to $\phi\approx 0$ at $\Gamma$. We can write the approximate free energy as:

[TABLE]

where

[TABLE]

and where we introduce the phase-field length scale $\xi$. $g(\phi)$, $w(\phi)$ are monotonically increasing and decreasing functions of $\phi$ respectively where $g(1)=w(0)=1-g(0)=1-w(1)=1$, and $C_{\phi}=\int_{0}^{1}\sqrt{w(\phi)}\,d\phi$ is a scaling constant. Throughout this manuscript we will use the formulation developed by Karma, Kessler and Levine (KKL) [32] and set $g(\phi)=4\phi^{3}-3\phi^{4}$ and $w(\phi)=1-g(\phi)$. This choice allows for the propagation of a single crack and prohibits the creation of any new cracks where $\phi=1$. Rigorous mathematical results [16] show that the minimizers of (7) converge to those of (6) when $\xi\rightarrow 0$, and crack propagation laws traditionally used in the LEFM framework (Griffith criterion $G=G_{c}$ and the principle of local symmetry) have also been rigorously derived for this model [29]. However, there is an ever-growing scientific consensus that $\xi$ represents an important material property that is needed for faithfully reproducing the experimental observations [17, 40, 57]. Indeed, for this class of functionals defined by (2)-(4) one can obtain a physical process zone size using the maximum tensile strength and fracture toughness of the material [50].

We should also note that we can write (4) equivalently as [32]:

[TABLE]

where $\mathcal{E}_{c0}=G_{c}/4C_{\phi}\xi$ is interpreted in the KKL model as a critical energy density (stress) above which the $\phi=0$ state becomes energetically favored and $\kappa=G_{c}\xi/2C_{\phi}$ is the coefficient of the gradient square term that together with $\mathcal{E}_{c}$ determines the size of the process zone and fracture energy.

2.1 Quasi-static evolution

In a quasi-static (*i.e., *rate independent) setting where we assume that fracture and elasticity have much smaller time-scales than loading (and other variables), we can pose the problem of fracture as a minimization problem. We define the problem of fracture as:

[TABLE]

where $\mathcal{K}_{A}=\{{u}|\,{u}(x)=\hat{{u}}(x),\,\forall x\in\partial_{D}\Omega\}$ is the set of displacement that satisfies the boundary conditions (*i.e., *the admissible displacements) and $\mathcal{A}_{\Gamma}=\{\Gamma|\,\Gamma(t^{+})\supset\Gamma(t^{-}),\,\forall t^{+}>t^{-}\}$ is the admissible fracture set. The energy minimization problem assumes that the system evolves in such a manner as to minimize the free energy (1). It is easy to show that for a single connected crack (7) is equal to the Griffith fracture theory which states that quasi-static cracks only propagate at $G=G_{c}$.

This can be translated to its phase-field equivalent problem (6) as:

[TABLE]

where $\mathcal{A}_{\phi}=\{\phi|\,0\leq\phi(t^{+})\leq\phi(t^{-})\,\forall t^{+}>t^{-}\}$. We can then write the governing equations for ${u}$ and $\phi$ as Euler-Lagrange equations of (2):

[TABLE]

where $\dfrac{\partial\mathcal{F}}{\partial x}(\tilde{x})$ is the Gâteaux derivative of free energy $\mathcal{F}$ with respect to $x$ in $\tilde{x}$ direction.

2.2 Ginzburg-Landau dynamics

It is also possible to pose the problem of fracture as a gradient flow toward the local minimum energy basins while assuming a rate independent elasticity as (8) (*i.e., *ignoring the elastic kinetic effects). In this setting we write the evolution equation of the phase field $\phi$ in the standard Ginzburg-Landau form [29]:

[TABLE]

where $\dfrac{\delta\mathcal{F}_{\xi}}{\delta\phi}$ is the Fréchet derivative of free energy $\mathcal{F}_{\xi}$ with respect to $\phi$. Integrating both sides, (11) can be equivalently written as

[TABLE]

where $\tau_{\phi}$ is the time scale associated with the evolution of the phase field $\phi$. Not surprisingly, we can recover the quasi-static Euler-Lagrange equation (9) when $\tau_{\phi}\rightarrow 0$.

3 Phase-field models of fatigue fracture

Since the classical phase-field method essentially approximates the Griffith problem, it does not allow for fracture propagation at sub-critical loads *i.e., *$G<G_{c}$ ($K_{I}<K_{IC}$). Since fatigue fracture, by definition, happens at sub-critical loads, it cannot be modeled with classical phase-field formulations. Although there is an abundance of literature on the physical mechanisms of fatigue fracture initiation and propagation, here we aim to model fatigue on a phenomenological level as a dissipative process that lowers the effective near-tip fracture toughness. This allows sub-critical cracks to propagate under cyclic loading without changing the monotonic fracture toughness. The two models introduced subsequently are based on this observation.

3.1 Model-I: the variable critical energy model

In the first model we achieve sub-critical propagation by lowering near-tip fracture toughness by reducing the critical energy density $\mathcal{E}_{c}$ around the crack tip. By lowering $\mathcal{E}_{c}$ near the crack tip, in essence, we lower the energy barrier for propagation of the fracture phase-field. We rewrite the surface energy (5) with spatially varying $\mathcal{E}_{c}$:

[TABLE]

where $\mathcal{E}_{c}(t=0)=\mathcal{E}_{c0}$ and $\mathcal{E}_{c}(x)$ evolves according to

[TABLE]

where $\chi_{c}$ is a scaling factor which controls the evolution of the critical energy, and $\mathcal{E}_{c\min}\geq 0$ signifies the lowest energetic barrier against fracture propagation. We have set $\mathcal{E}_{c\min}=0$ for all the simulations presented in this article. In addition, we define the Heaviside step function:

[TABLE]

Using the above formulation, one can immediately see that the critical energy density remains constant under a constant load but decreases when the load is increasing in time. The stress divergence at the crack tip, which is cut off on a scale $\xi$ in the phase-field model, naturally decreases $\mathcal{E}_{c}$ near the tip thereby allowing for fatigue propagation. However, as a first step in describing fatigue crack growth of brittle materials, we keep the degradation process localized to the process zone by introducing the multiplicative factor $w(\phi)$ in (13). Furthermore, we define the positive energy potential as

[TABLE]

This allows the critical energy to only decrease in tensile stresses, thereby suppressing the propagation of compressive fatigue cracks.

Although the above formulation limits the propagation of compressive fatigue cracks, it does not limit inter-penetration of the created fatigue fracture faces under compressive loads. It is straightforward to extend this formulation to account for non-interpenetration by additively splitting the elastic energy potential into volumetric and deviatoric parts [7]. Furthermore, it is possible to define the positive energy potential (15) to allow for shear degradation both in tension and compression using the same additive split introduced in [7] .

To elucidate how this model allows for sub-critical propagation, we should note that in the elliptic approximations of brittle fracture, the fracture toughness can be calculated as

[TABLE]

As a result, by decreasing the critical energy density in front of the crack, the near-tip fracture toughness (non-uniformly) decreases allowing for propagation of the fracture. We should note, however, that the energy release rate for the crack to propagate (at process zone scale $\xi$) in the resulting non-uniform $\mathcal{E}_{c}$ field is non-trivial. Therefore, obtaining a closed-form solution for the equivalent near-tip fracture toughness may not be possible.

Reducing the critical energy density $\mathcal{E}_{c}$ in our formulation is very similar to the lowering of cohesive stiffness in [44]. Indeed, as pointed out by Pham and Marigo [50], the process zone size $\xi$ and the critical energy density $\mathcal{E}_{c}$ in the phase-field formulation are analogous to the opening displacement and maximum stress in the cohesive zone formulation, respectively.

3.1.1 Thermodynamical considerations

It is straightforward to calculate the energy dissipation with respect to time for (12). Taking the time derivative of the approximate free energy $\mathcal{F}_{\xi}$ and using the governing equations (8)-(9) for the quasi-static evolution and further using the fact that the local phenomenological degradation law (13) implies $\partial\mathcal{E}_{c}/\partial t\leq 0$, we obtain:

[TABLE]

which ensures the dissipative nature of the quasi-static algorithm taking into account that (9) is equivalent to $\delta F_{\xi}/\delta\phi=0$.

Similarly, if the evolution of the phase field $\phi$ is governed by gradient dynamics as defined in (10), we can write:

[TABLE]

Therefore, the proposed algorithm is dissipative.

3.2 Model-II: $\gamma$-model

It is also possible to reduce the spatial fracture toughness defining a variable fracture toughness ratio $\gamma$ in (4).

[TABLE]

where $\gamma(x)$ evolves as

[TABLE]

The minimum relative fracture toughness $\gamma_{\min}>0$ is introduced (in the same spirit as $\mathcal{E}_{c\min}$ (13)) as the lowest relative fracture toughness achievable. It should be noted that (if $\gamma_{\min}=0$) the limit $\gamma\rightarrow 0$ of (19) can result in $\phi=0$ at a dense subdomain of $\Omega$ which results in an ill-defined elasticity problem.

Similar to the first model, it is worth emphasizing that unlike the classical phase-field formulation with uniform fracture toughness, the resulting fracture toughness near the crack tip defined by (19) is non-trivially related to $G_{c}$ because $\gamma$ is spatially varying in the crack tip region.

3.2.1 Thermodynamical considerations

Similar to section 3.1.1 taking the derivative of the approximate free energy $\mathcal{F}_{\xi}$ with respect to time and taking in account the quasi-static governing equations (8)–(9) we can write:

[TABLE]

once again noticing that ${\partial\gamma}/{\partial t}\leq 0$ as imposed by the phenomenological degradation law (20).

The same procedure can be duplicated for the gradient dynamics:

[TABLE]

4 Numerical implementation

4.1 Dimensionless parameters

In all sections that follow, we define dimensionless spatial coordinates using one of the sample dimensions $\bar{{x}}={x}/\xi$ and displacements by $\bar{{u}}={u}/h_{0}$ where we define $h_{0}=\sqrt{G_{c}\,\xi/E}=K_{IC}\sqrt{\xi}/E$. We can write the dimensionless free energy form of (2) as:

[TABLE]

where $\bar{{\mathcal{C}}}_{ijkl}={\mathcal{C}}_{ijkl}/E$, and $\bar{f}_{i}=f_{i}\sqrt{\xi/G_{c0}E}$. Furthermore, for the dynamic phase-field evolution (11), we use $\tau_{\phi}$ to adimensionalize time *i.e., *$\bar{t}=t/\tau_{\phi}$.

4.2 Numerical algorithms

In this article, the spatial discretization is done using the Galerkin finite element method. The critical energy field $\mathcal{E}_{c}$ and the fracture toughness ratio parameter $\gamma$ are discretized as discontinuous constants across each element (*i.e., *DG-0), and components of ${u}$ and $\phi$ are discretized using Lagrange basis functions.

4.2.1 Alternate minimum algorithm for quasi-static evolution

Our minimization strategy for (2) is now classical [16]. At each time step, it is achieved by alternating minimizations with respect to $u$ and $\phi$ until convergence, leveraging the separate convexity of the regularized energy with respect to each field. As stated before, the solution for displacement ${u}$ and phase field $\phi$ is calculated by solving their respective Euler-Lagrange equations (8)-(9). The first step to obtain ${u}$ using the alternate minimization algorithm is a simple convex problem implemented by solving the elastic equilibrium. We use a bounded reduced space Newton minimization scheme [14] for the second step to solve for $\phi$ . The phase-field method requires the spatial resolution of discretization to be smaller than the process zone length $\xi$.

As presented in algorithm 1, the variables of the fatigue models $\mathcal{E}_{c}$ or $\gamma$ are only updated outside the alternate minimization loop after the alternate minimization algorithm converges at each time step i.e.,

[TABLE]

or

[TABLE]

The resulting problems are often very large and require the use of a parallel programming paradigm and the complex numerical tools therein. Our implementation [40] relies on the distributed data structures provided by libMesh [33] and linear algebra provided by PETSc [10, 9].

4.2.2 Dynamic algorithm

The algorithm 2 for dynamic evolution is much simpler and only ensures elastic equilibrium at each time step. We have implemented an implicit forward Euler time stepping algorithm for $\phi$ where at each time step we solve:

[TABLE]

using the bounded reduced space Newton minimization scheme implemented in PETSc [9].

5 Numerical simulations

5.1 Paris law

The first validation on any continuum fatigue model is to assess its ability to simulate the Paris law. This law states that the increase in the length of a fatigue crack per unit cycle relates to the range of applied stress intensity factors as a power-law *i.e., *$da/dN\sim(\Delta K)^{n}$ where $a$ is the crack length and $N$ is the number of cycles. For this purpose, we set up a simple problem illustrated in Fig. 1 for a rectangular domain $\Omega=(-L_{x}/2,L_{x}/2)\times(-L_{y}/2,L_{y}/2)$ containing an initial crack $\Gamma_{0}=(-L_{x}/2,0)\times{0}$. We impose an oscillatory vertical displacement on the top and bottom $u_{y}(\pm L_{y})=\delta(t)$ where $\delta(t)$ is oscillating in a sawtooth fashion by linearly increasing between [math] and $\delta_{\max}$ when $t$ increases from $nT$ to $nT+T/2$ and linearly decreasing from $\delta_{\max}$ to [math] when $t$ increases from $nT+T/2$ to $(n+1)T$, where the integer $n$ is the cycle number and $T$ is the period. Additionaly, we set the horizontal displacement $u_{x}(\pm L_{y})=0$. All the simulations in this section are done for $L_{y}=50\xi$ and $L_{x}=200\xi$.

The energy release rate from the tip of a straight crack in the geometry shown in Fig. 1 can be approximated by $G\simeq 2E\delta^{2}/(1-\nu^{2})L_{y}$. Using this formula, we performed the numerical simulation at different $\delta_{\max}$ corresponding to different $G/G_{c}$ ($K_{I}/K_{IC}$) ratios where we set $\nu=1/3$. We note that here $\Delta K/K_{IC}=K_{I}/K_{IC}$ since the material is completely unloaded at each cycle.

Fig. 2 shows the evolution of the crack ($\phi=0.5$ contour line) over 40 cycles (after the initial transient regime) as well as the spatial distribution of $\gamma$ at the last step. The translational invariance of the geometry and the phase-field formulation result in a steady state propagation of the initial crack with $da/dN$ approaching a constant after a transient regime lasting a few cycles.

Figs. 3–4 illustrates the emergent Paris law for different values of $\chi_{c}$ for the two models. We can observe that the exponent of the power-law decreases with increasing $\chi_{c}$. This non-trivial behavior can be interpreted as resulting from the fact that, for sufficiently large $\chi_{c}$, the degradation variable $\gamma$ reaches its minimum value $\gamma_{\min}$ in a finite region around the crack tip, as found by inspection of the numerical results. The size of this region increases with the amplitude $K_{I}/K_{IC}$ of cyclic loading. Therefore, as $\chi_{c}$ increases, the size of the region around the crack tip where $\gamma>\gamma_{\min}$ decreases, thereby decreasing the sensitivity of degradation and crack growth to load $K_{I}/K_{IC}$.

Figs. 3–4 show convergence of the rate of crack growth $da/dN$ with increasing number of time steps per cycle ($T/\delta t$ varying from 2 to 200). Convergence is faster at lower $\chi_{c}$ for a given load $\Delta K/K_{IC}$ since a lower degradation rate can be resolved with a smaller number time steps per cycle. As a result, the rate of crack growth is converged for all loads at $\chi_{c}=1$ but only a smaller range of loads for larger $\chi_{c}$. We note that, even though resolving the continuum limit of the time-dependent evolution equations (13) and (20) controlling the degradation process requires the time step to be chosen small enough for a given $\chi_{c}$ and load, the corresponding discrete evolution equations (24) and (25) still produce Paris law type behavior over a finite range of loads. Hence, similarly as iterated maps of dynamical systems, they provide discrete time models of fatigue crack growth with exponents that decrease with decreasing $\delta t/T$ as seen in Figs. 3–4.

Using the dynamic formulation (11) allows for further probing of the models where we set $\delta t=\tau_{\phi}/5$. As expected, when the relaxation time scale of the phase-field is much smaller than the cyclic loading period ($\tau_{\phi}\ll T$), the Paris-law exponent converges to its quasi-static value as depicted in Fig. 5–6. Interestingly, away from this limit (i.e. for $\tau_{\phi}\sim T$), the Paris law exponent is significantly lower than its quasi-static value. This result suggests that dissipation processes that slow down crack growth on the time scale of cyclic loading contribute to decreasing the sensitivity of the crack growth rate to the local fracture toughness degradation rate and hence lowering the Paris law exponent.

5.2 Crack interaction

In the next set of simulations, we examine the problem of crack interactions for two offset cracks approaching each other under tensile loading. The naturally occurring En-passant configuration [2, 48, 59], has been the focus of a series of theoretical [35, 38, 26] and experimental studies [23, 34]. The geometry of the problem is illustrated in Fig. 7 where we apply a vertical displacement ${u}(\pm L_{y})=(0,\delta(t))$. Table 1 summarizes the parameters used for the En-passant numerical simulations. All simulations were performed using a mesh consisting of $\sim 520k$ delaunay triangular elements with an average size of $\xi/5$.

It is essential to first perform a quasi-static simulation without degradation to calculate the initial energy release rate $G(\delta)$ (or $K_{I}$) needed for the fatigue simulations. Fig. 8 depicts the results of this En-passant quasi-static simulation where the load $\delta$ was increased monotonically ($\delta(t)=kt$). The setup of the problem is such that each crack first deviates from its original plane and then rotates toward the opposing crack. The energy release rate from the crack tip increases at first, destabilizing the crack at $\delta\simeq 20\sqrt{G_{c}\xi/E}$. As the cracks turn, they lose energy (the energy release rate decreases) so that the complete fracture happens at $\delta\simeq 42.4\sqrt{G_{c}\xi/E}$.

We performed two fatigue simulations corresponding to each model, where $\delta(t)$ oscillates as before in a sawtooth fashion between [math] and $\delta_{\max}$ in period $T$. Here we set $\delta_{\max}=14\sqrt{G_{c}\xi/E}$ which corresponds to $K_{I}/K_{IC}\simeq 0.7$. We also set the proportionality factors for both simulations to $\chi_{c}=1$. Since our models reduce to their quasi-static counterparts at $\mathcal{E}_{c}=\mathcal{E}_{c0}$ or $\gamma=1$, they capture the transition between the fatigue fracture and the quasi-static regime naturally. Both simulations show that the crack initially propagates at the sub critical load level but then destabilizes once $G(\delta_{\max})\geq G_{c}$ and propagates abruptly (depicted in Fig. 9 frames 3-4 and in Fig. 10 frames 4-5). Finally, as the energy release rate from the tip of the turning cracks decreases, the fatigue propagation slows down significantly.

Fig. 11 compares the crack path of the quasi-static evolution under monotonic loading with the crack paths obtained from the critical energy and the $\gamma$-model under cyclic loading. The three macroscopic crack paths are in good agreement up to a spatial resolution comparable to the process zone size. Also, the abrupt propagation segment of both fatigue crack paths remains unchanged compared to the quasi-static fracture path.

5.3 3-D penny-shaped cracks coalescence

In the first three-dimensional set of simulations, we model the coalescence of two 3-D penny-shaped cracks. Two penny-shaped cracks with radius $R_{c}$ are initialized at distance $b$ in an $L^{3}$ cube where we apply displacement boundary conditions ${u}=(0,0,\delta(t))$ at $(x,y,\pm L)$. A similar problem was originally investigated by Legrand et al. [36] using a semi-analytic method where the authors used a perturbation method along with the principle of local symmetry (i.e. vanishing mode II stress intensity factor at the crack tip) to predict the propagation and coalescence of the cracks. We performed this simulation using a mesh containing $16.5\,M$ tetrahedral elements and the $\gamma$-model 3.2. Since the computational cost of performing a quasi-static simulation (over many cycles) is prohibitive, we opted for the dynamic phase-field formulation 2.2 with the loading period $T/\tau_{\phi}=10$ and $\chi_{c}=10$. The simulation parameters are summarized in table 3. In this geometry, the energy release rate increases as the cracks grow and speeds up the fatigue crack propagation. Furthermore, the energy release rate increase is more dramatic at the coalescence site resulting in the cracks merging over 5 cycles between $N=50-55$. This is similar to the behavior predicted by Legrand et al. (see for example Fig. 5 in [36]).

5.4 3-D inclined penny-shaped crack

In the last set of numerical simulations, we model the propagation of a 3-D penny shaped crack. A tilted penny-shaped crack with radius $R_{c}$ is initialized in the middle of an $L^{3}$ cube with its normal $\boldsymbol{n}_{c}=\left(\cos(\theta_{c}),\cos(\theta_{c}),\sin(\theta_{c})\right)$. The loading ${u}=(0,0,\delta(t))$ is applied at $(x,y,\pm L)$. A similar problem was originally simulated by Gravouil et al. [27] using a level set method in a quasi-static setting for $\theta=\pi/4$. We performed this simulation using a mesh containing $16\,M$ tetrahedral elements. The simulation was performed using the $\gamma$-model 3.2 and the dynamic phase-field formulation 2.2 (similar to the previous section to limit the computational cost). We chose the loading period $T=10\,\tau_{\phi}$ and $\chi_{c}=10$. The simulation parameters are summarized in table 3.

Fig. 15 shows the crack front obtained where, similar to brittle fracture, the crack front twists and turns towards the symmetry line $z=0$ where its $K_{II}=0$. As shown in this figure, the propagation of the fatigue crack accelerates over time. This is due to the increase of the energy release rate from the crack tip, which can be easily understood since the area under tension decreases as the crack propagates. Fig. 16 shows how the relative fracture toughness ratio $\gamma$ changes spatially around the crack front. Like the En-passant simulations in section 1, the change in $\gamma$ is limited to the immediate vicinity of the crack front.

Finally, Fig. 17 shows a comparison between the quasi-static and fatigue crack paths. As in the case of the En-passant simulations 5.2, the crack path remains unaltered with respect to the quasi-static path.

6 Towards a ductile phase-field fatigue formulation

In all previous sections, the fatigue degradation was limited to the process zone around the crack tip. However, in K-dominant fatigue of ductile materials, the crack-tip plastic zone has a scale larger than the fracture process zone but much smaller than the specimen dimensions. Furthermore, ductile fracture is driven by void nucleation, growth, and coalescence mediated by the plasticity. Locally, voids reduce the load-bearing portion of the crack plane similar to a damage variable in the continuum description. Indeed, such modifications have been suggested and implemented in the celebrated Gurson-Tvergaard-Needleman (GTN) [28, 43] model to account for the effect of plasticity on fracture [61]. Thus, one can posit that in ductile fracture, the fracture toughness locally reduces due to plasticity. Such a ductile phase-field fracture formulation has been proposed in [6] where $g(\phi)$ was modified in (3) to include the equivalent plastic strain in a way that lowers the fracture toughness locally in plastically deforming areas and retards the creation of cracks in the elastic regions.

Here we suggest a simple extension of the previously developed fracture toughness model to account for degradation on a scale larger than the fracture process zone. We modify the evolution equation for $\gamma$ (20) to allow for fatigue degradation of fracture toughness to take place at any point where the tensile elastic energy potential is greater than a critical stress $W^{+}>\sigma_{c}$ i.e.,

[TABLE]

Fig. 19 shows the results of the simulations for a simple mode-I geometry depicted earlier in figure 1. For all simulations shown, we set $\chi_{c}=1$, $\nu=1/3$ and $L_{y}/L_{x}=1/4$. As one can see, the Paris law is reproduced with a smaller exponent for smaller $\sigma_{c}$. Comparing $\sigma_{c}$ to the divergent stresses at the crack tip, one can derive the fatigue degradation zone’s radius as

[TABLE]

Fig. 19 shows that the modified degradation law (27) reproduces the Paris law with an exponent that decreases with increasing radius $r_{c}/\xi$ of the degradation zone. Importantly, when $r_{c}\gg\xi$, this exponent ($m\approx 4.64$ when $r_{c}/\xi=32$) is substantially smaller than the exponent ($m\approx 8.22$) obtained when degradation is localized to the process zone (Fig. 4). As expected, when $r_{c}$ is comparable to $\xi$, the exponent is similar to the one previously obtained.

By extending this approach, we have developed a ductile fatigue model where the degradation of the fracture toughness is controlled by the accumulated plastic strain $e^{p}_{eq}$, which is achieved by replacing (27) by

[TABLE]

This model reproduces non-trivial features of ductile fatigue such as crack closure. Both the model and results will be presented elsewhere.

7 Conclusions

We introduced two phase-field models of fatigue crack growth. In our models, the material degrades near the crack tip as a result of cyclic loading which allows for sub-critical (*i.e., *$G<G_{c}$) propagation of cracks. Crucially, this class of models preserves the fracture toughness as a material property globally and only affects the near-tip properties on the process zone scale. In other words, if the crack is loaded beyond $G_{c}$($K_{IC0}$), it propagates at a constant (non-cyclical) load. This allows us to perform uniform simulations which seamlessly switch between monotonic $G>G_{c}$ and fatigue $G<G_{c}$ crack propagations (see, for example, figure 11 in section 5.2).

The most intriguing finding of this article is the emergence of the Paris law. We showed (using numerical simulations) that the Paris power-law emerges naturally from our formulation under simple mode-I loading section 5.1 in both quasi-static simulations and in simulations where the phase-field $\phi$ follows Ginzburg-Landau dynamics. In the quasi-static simulations, we recovered the Paris law with high exponents typical of fatigue crack growth in brittle materials [54]. Additionally, we were able to capture the Paris law over a wider range of exponents in Ginzburg-Landau gradient dynamics by changing the ratio of the cyclic loading period over the charachteristic relaxation time of the phase-field $T/\tau_{\phi}$.

Here, we should also emphasize that our two models should be thought of as two specific cases of a more general class of fatigue phase-field models where the fracture toughness is degraded locally around the crack tip as a function of the elastic power (the time derivative of elastic energy density) stored at each point in the system. For example, our choice of linear degradation of $\mathcal{E}_{c}$ and $\gamma$ are the simplest cases, and replacing $H(\mathcal{E}_{c}-\mathcal{E}_{c\min})$ with $(\mathcal{E}_{c}-\mathcal{E}_{c\min})$ in (13) (or, similarly, replacing $H(\gamma_{c}-\gamma_{\min})$ with $(\gamma-\gamma_{\min})$ in (20)) produces similar results to those presented in this article.

The En-passant simulations 5.2 and the inclined penny-shaped 3-D crack simulations section 5.4 suggest that, in an isotropic elastic solid, the fatigue crack propagates on a path similar to quasi-static fracture *i.e., *at $K_{II}=0$ [29]. Even though this observation agrees with current wisdom, it is not formally derived here from an analysis of the phase-field models, as previously done for quasi-static propagation under constant loading [29].

Finally, in section 6, using a simple extension of our brittle model, we showed that in a ductile material where the degradation occurs at the crack-tip plastic zone (on a scale larger than the process zone $\xi$, but smaller than the system size) the Paris law exponent is reduced ($m\simeq 4.64$ in 19). This key finding paves the way for the introduction of ductile fatigue phase-field models where degradation is driven by plastic deformation.

8 Acknowledgments

The authors’ work was supported by the U.S. Navy’s SBIR/STTR office, under the auspices of contract N68335-16-C-0206. The majority of the numerical simulations were performed using resources of Northeastern University’s Discovery cluster located in Massachusetts Green High Performance Computing Center (MGHPCC) in Holyoke, MA.