1 Introduction

As it is in the nature of dynamics, components of machines are subjected to oscillating loads. The magnitude of these loads may be very low in comparison with the materials ultimate strength, but nevertheless, it might be the very reason for failure of machine parts and finally the cause of catastrophic accidents. The underling physical phenomenon is denoted as mechanical fatigue of materials. This process involves the accumulation of damage by ongoing repetition of loading and unloading, which leads to the formation of microscopic cracks and with continuing cycling even to macroscopic cracks and fatal failure of the structure. Catastrophic accidents from the past decades of e.g., airplanes, trains, cars, motorcycles, turbine rotors, pressure vessels and others make clear that this field must be the focus of investigations and studies. Fatigue of materials has been the focus of numerous textbooks like e.g., those from Suresh, Schijve or Dowling [11, 33, 39]. The progression of the damage caused by fatigue is divided into the phases: (i) nucleation of microstructural damage and growth of microcracks to form a macrocrack and (ii) growth of this macrocrack until final fracture of the structure. A specific approach for the description of fatigue was developed for each of these two stages of the total fatigue life. The total life approach deals with the first phase by means of S–N curves, which reveal the number of load cycles for specific amplitudes to be applied until a macrocrack is initiated and extended to complete fracture of a specimen. It must be noticed that specimens to record S–N lines are usually designed such that the amount of cycles until initiation of a crack is significantly higher (about \(90 \%\) of the total life) than the period of crack growth. The second approach has its focus on the period of the extension of an existing macrocrack to a critical size. The small steps of crack extension of a fatigue crack are characterized by the crack growth rate, which can be evaluated using the so-called Paris law [30]. Within this law, the cyclic load of a constant amplitude sequence is commonly represented by the range of the stress intensity factor \(\Delta K\) (see e.g., [32]). The crack growth rates associated to a certain load range are fit using a power law, which reveals a straight line within double logarithmic scales. Practically, it is crucial which period of the fatigue life is of interest in order to decide which concept to life-time estimation is applied.

From a numerical modeling perspective, the problem is on the one hand the high number of load cycles leading to high computational effort, and on the other hand, the lag within the mechanistic basis of the fatigue phenomenon. Fatigue is affected by a large number of complex processes occurring on different scales. Accordingly, a large number of numerical approaches to fatigue incorporate the phenomenological laws from above. For instance, in [20, 29], Paris’ law is used to simulate fatigue crack growth with finite elements and generalized finite elements, respectively. The law is also incorporated in a widely used tool for fatigue crack growth simulations, namely NASGRO, [15]. However, in some studies (see e.g., [9, 25]), cyclic fatigue life until crack initiation was modeled utilizing the framework of Continuum Damage Mechanics (CDM) from Lemaitre [24]. Fish and Oskay introduced a numerical model in [13], where the damage evolution is integrated by a growing void volume fraction governed by Gurson’s model. In this work, cracks are modeled by the deletion of finite elements after a the void volume fraction exceeds a critical damage value. Generally, the approximation of cracks or crack extension within a conventional simulation setting is problematic, as element deletion or remeshing algorithms must be integrated in the simulation scheme. A very beneficial technique to overcome these issues is provided by the phase field framework. Within this method, cracks are modeled via an additional scalar field indicating broken or intact material. This scalar field is coupled with the material’s stiffness in order to account for the significant stiffness reduction in broken areas. The method has gained much attention in the last decade and a high number of studies in different fields of fracture mechanics were published so far. The topics of the presented models range from dynamic fracture by e.g., Schlüter et al. [34] over brittle fracture by e.g., Borden et al. [5], Miehe et al. [27] and Kuhn and Müller [21], and also ductile fracture, proposed by e.g., Borden et al. [4], Kuhn et al. [22] and Aldakheel et al. [1] to models including effects of anisotropic materials like e.g., the models from Teichtmeister [40] or Schreiber et al. [35]. The field of fatigue crack growth is relatively new in the context of phase field modeling, and to the best of the authors knowledge, it was first dealt with this issue by Alessi et al. [2] for a one-dimensional setting. Within this work, the critical energy release rate of the material becomes degraded according to an accumulated strain value. Following up on this publication, several phase field approaches for fatigue fracture in higher dimensions as e.g., [7, 26, 38] were proposed. Furthermore, a phase field fatigue model for cracking in cyclically loaded concrete was proposed by Wriggers et al. [42]. In [37], we presented an alternative basic phase field approach to handle the issue of cyclic fatigue, which is generally an extension of the work from Kuhn and Müller for brittle fracture [21]. Fatigue is incorporated by an additional energy contribution controlled by accumulated damage, and we showed that crack growth rates evaluated by this model follow Paris’ law. However, so far, we only considered constant amplitude loading for a pulsating load sequence in unidirectional cracks. As in real applications, generally various loading amplitudes occur with different stress ratios, which have significant effect on the crack growth curve, we apply appropriate modifications to the model in order to meet this requirement.

2 Model description

In order to provide a comprehensive description and to point out the motivation for the underlying modification regarding fatigue damage, we start with a brief illustration of the phase field method for fracture and give a characterization of the basic model and its equations. Within the subsequent section, we then continue to introduce the new model for fatigue crack growth, its properties, and the underlying simulation scheme.

2.1 Regularization of brittle fracture using the phase field method

The phase field method very elegantly overcomes the issue of modeling interfaces by coupling fields of physical quantities to an order parameter that indicates the phase at a certain point. Thus, a defuse approximation of sharp interfaces is obtained. In terms of fracture mechanics the underlying pure phases are either intact material or broken material. The range of the order parameter is used to smear cracks over a certain area, which is then referred as transition zone. In order to yield an approximation of the sharp crack, this zone shall be as small as possible.

Within the phase field model for brittle fracture [21], the approximation of cracks is accomplished by means of the continuous transition of the variable \(s({{\varvec{\mathrm {x}}}},t)\) in the range [0, 1], where the value \(s=1\) refers to the initial intact condition of the material, and \(s=0\) refers to broken material. In the context of phase field fracture modeling, the theory by Griffith [16] is a fundamental contribution. According to this theory, the condition for crack extension is the compensation of strain energy release by additional crack energy, related to crack surface area of infinitesimal crack extension. Following up on this consideration, the phase field model uses the regularized formulation of brittle fracture proposed by Bourdin et al. [6], where it is postulated that displacements and crack patterns represented by \(s({{\varvec{\mathrm {x}}}},t)\) minimize the total internal energy

$$\begin{aligned} \begin{aligned} E= \int \limits _{\mathrm {\Omega }} \psi ( {{\varvec{\varepsilon }}},s,\nabla s) \, \mathrm {d}V = \int \limits _{\mathrm {\mathrm {\Omega }}} (g(s)+\eta )W( {{\varvec{\varepsilon }}}) + \mathrm {\Gamma }(s,\nabla s) \; \mathrm {d}V \end{aligned} \end{aligned}$$
(1)

of a loaded and potentially fractured body. The first energy contribution in Eq. (1) accounts for strain energy with the strain energy density of a linear elastic material

$$\begin{aligned} W=\frac{1}{2} {{\varvec{\varepsilon }}} : \left( \mathbb {C} {{\varvec{\varepsilon }}} \right) , \end{aligned}$$
(2)

with the material stiffness tensor \(\mathbb {C}\) and the infinitesimal strain tensor \({{\varvec{\varepsilon }}}\). The strains are computed by the spatial gradient of the displacements \({{\varvec{\mathrm {u}}}}\) via

$$\begin{aligned} {{\varvec{\varepsilon }}}({{\varvec{\mathrm {u}}}})=\frac{1}{2}\left( \nabla {{\varvec{\mathrm {u}}}}+\left( \nabla {{\varvec{\mathrm {u}}}}\right) ^T \right) . \end{aligned}$$
(3)

The connection to Griffith’s theory becomes evident considering the second energy contribution of Eq. (1), which accounts for crack energy characterized by the crack energy density

$$\begin{aligned} \mathrm {\Gamma }=\mathcal {G}_c\left( \frac{(1-s)^2}{4\epsilon }+\epsilon |\nabla s|^2 \right) . \end{aligned}$$
(4)

The critical energy release rate \(\mathcal {G}_c\) is a material parameter with unit energy per surface area. This quantity can be referred to as a material specific measure of the energy release per crack surface increase and within this regularization it applies to a local and also to a nonlocal contribution governed by s and its spatial gradient. According to Eq. (4), crack energy will increase once the phase field order parameter decreases or corresponding gradients evolve. A length scale \(\epsilon \) is incorporated in order to control the width of the transition zone form intact to broken material. Accordingly, this parameter is also referred as regularization parameter. In [21], it is shown that crack energy is correctly estimated as \(\epsilon \) goes to zero. However, the feasible minimum of \(\epsilon \) is certainly a matter of discretization.

The stress \({{\varvec{\sigma }}}\) is obtained by applying the derivative of the potential \(\psi \) with respect to strain. This yields

$$\begin{aligned} {{\varvec{\sigma }}}=\frac{\partial \psi }{\partial {{\varvec{\varepsilon }}}}= (g(s)+\eta )\mathbb {C}{{\varvec{\varepsilon }}} \end{aligned}$$
(5)

and as a consequence of the coupling to the scalar phase field the stress is degraded, governed by the function g(s) with \(s<1\). Furthermore, the application of the parameter \(\eta \) is important. This coefficient is chosen such that \(0<\eta \ll ~1\) in order ensure a marginal residual stiffness and therefore avoid numerical problems concerning the static solution scheme. The most widely used form of a degradation functions is a simple quadratic approach

$$\begin{aligned} g(s)=s^2. \end{aligned}$$
(6)

It may be noticed that alternative approaches for degradation functions were proposed (see e.g., [4, 23]) affecting the material behavior. Thus, applying a certain cubic approach as degradation function, it can be shown that the material does not experience a significant degradation of stiffness just before cracking, which is not the case for Eq. (6). However, this degradation function was used in many studies as it has been proven to provide a robust implementation and reliable results.

As the potential does not only dependent on the order parameter but also on its gradient, within [21], a generalized Ginzburg-Landau equation (see [18]) is applied for the time evolution of the phase field s with

$$\begin{aligned} \begin{aligned} \dot{s}=-M \frac{\delta \psi }{\delta s}=- M \left\{ \frac{\partial \psi }{\partial s}- \nabla \left( \frac{\partial \psi }{\partial \nabla s} \right) \right\} =-\frac{M}{2}\left\{ g'(s) {{\varvec{\varepsilon }}} : \left( \mathbb {C} {{\varvec{\varepsilon }}} \right) - \mathcal {G}_c \left( 4\epsilon \mathrm {\Delta } s+ \frac{1-s}{\epsilon }\right) \right\} \end{aligned} \end{aligned}$$
(7)

utilizing the Laplace operator \(\mathrm {\Delta }\). The scalar mobility parameter M in Eq. (7) may be considered as viscous regularization parameter regarding static conditions with the limit case \(M \rightarrow \infty \) to approximate quasi static conditions.

2.2 Model formulation for cyclic fatigue cracks

2.2.1 Interpretation of the order parameter

A basic property of the phase field model for cyclic fatigue proposed in this work is the interpretation of the order parameter s, since it can be interpreted as either physical damage variable according to CDM or as indicator quantity for the regularization. A formulation in terms a of a damage variable

$$\begin{aligned} D=\frac{\mathrm {\Delta } S_D}{\mathrm {\Delta } S}. \end{aligned}$$
(8)

with the area of the intersection of all microvoids with a certain plane \(\mathrm {\Delta } S_D\) and undamaged area \(\mathrm {\Delta } S\) as introduced by Lemaitre’s framework of CDM ([24]) may in a first place seam likely. This definition restricts the range of the damage variable D to \(0,\dots ,1\). A damage state \(D=0\) indicates undamaged material, while \(D=1\) indicates fully damaged (broken) material. Suppose now a material point is exposed to a constant stress \({{\varvec{\sigma }}}\) . The effective stress tensor \(\hat{{{\varvec{\sigma }}}}\) has to increase as the effective surface area carrying the load is reduced and therefore

$$\begin{aligned} \hat{{{\varvec{\sigma }}}}=\frac{1}{1-D}{{\varvec{\sigma }}}. \end{aligned}$$
(9)

By means of the strain equivalence principle and considering a linear elastic constitutive law, the strain energy density for the damaged material point can be formulated by means of effective stresses

$$\begin{aligned} W_D=\frac{1}{2}\hat{{{\varvec{\sigma }}}} : \left( \mathbb {C}^{-1}\hat{{{\varvec{\sigma }}}}\right) =\frac{1}{2(1-D)^2}{{\varvec{\sigma }}} : \left( \mathbb {C}^{-1}{{\varvec{\sigma }}}\right) . \end{aligned}$$
(10)

The strains can be obtained by differentiation

$$\begin{aligned} {{\varvec{\varepsilon }}}=\frac{\partial W_D}{\partial {{\varvec{\sigma }}}}= \frac{1}{(1-D)^2}\mathbb {C}^{-1}{{\varvec{\sigma }}} \end{aligned}$$
(11)

and an effective compliance tensor can be defined as

$$\begin{aligned} \tilde{\mathbb {C}}^{-1}=\frac{1}{(1-D)^2}\mathbb {C}^{-1}. \end{aligned}$$
(12)

If Eq. (5) is solved for \({{\varvec{\varepsilon }}}\), and Eq. (6) is used as degradation function one obtains

$$\begin{aligned} {{\varvec{\varepsilon }}}=\frac{{{\varvec{\sigma }}}}{s^2}\mathbb {C}^{-1} \qquad \text {with:} \quad \tilde{\mathbb {C}}^{-1}=\frac{1}{s^2}\mathbb {C}^{-1}, \end{aligned}$$
(13)

where the residual stiffness \(\eta \) is neglected. The state of a material in terms of differentiation between intact and broken material can be described by the damage variable D from CDM, and it can also be described by the phase field order parameter, where \(s=1\) for intact material and \(s=0\) for broken material. Further, \(D \in {[0,1]}\), \(s\in {[0,1]}\) and with the definitions for undamaged and fully damaged state the relation \(s=1-D\) can be formulated. Incorporating this relation in Eqs. (12) and (13), one finds the constitutive law of both formulations to coincide. An analog comparison is used in [41] to formulate a damage model in terms the phase field method. However, according to Bourdin [5], the phase field variable is in general not meant to describe damage in terms of a physical criterion like the evolution of microvoids. It results from the basic approach of approximating sharp interfaces by a smooth transition for numerical convenience. Accordingly, even if s may also be interpreted as a kind of nonlocal damage variable, in the following damage refers to fatigue and the connection to s is not explicit. Therefore, a proper estimate to represent the fatigue damage occurring in the material must be adopted.

2.2.2 Phase field approach for cyclic fatigue cracks

The convenient properties of the phase field modeling philosophy result from the formulation of the evolution equations for certain phases, which are commonly derived from potentials. Within the solution scheme a sequence of states is generated, which approximate energetic minima. No additional criteria are required. In terms of quasi static fracture, this means a crack will be predicted once an increase of crack energy is beneficial regarding the local minimum of the total potential energy. Also, knowledge about the crack extension direction is not required a priori.

Naturally, the underlying processes caused by fatigue of materials like e.g., movement of dislocations, formation of microcracks and others do also consume energy. Accordingly, this has to be taken into account by a respective potential. As a matter of fact, the evolution of s described by Eq. (7) will not lead to an evolution of a fatigue crack, no matter how many load cycles are explicitly simulated. This is simply because the driving processes associated with fatigue already occur far below the global ultimate strength of the material, but they accumulate. This accumulation of disruption is not considered by the total internal energy form Eq. (1). Therefore, the basic modification we apply to this functional is to incorporate an additional density contribution accounting for accumulated driving forces caused by the phenomenon of cyclic fatigue. The resulting regularized total internal energy is

$$\begin{aligned} \begin{aligned} E^{\text {f}}=\int \limits _{\mathrm {\mathrm {\Omega }}} (g(s)+\eta )W+ \mathrm {\Gamma } \; \mathrm {d}V + \int \limits _{\mathrm {\Omega }} h(s) P(D_\text {f}) \; \mathrm {d}V, \end{aligned} \end{aligned}$$
(14)

where \(P(D_\mathrm {f})\) represents the additional energies as piecewise defined function

$$\begin{aligned} P(D_\mathrm {f})=q\langle D_\mathrm {f}-D_{_\mathrm {c}} \rangle ^b. \end{aligned}$$
(15)

The \(\langle \cdot \rangle \) are referred to as Macaulay brackets with the definition

$$\begin{aligned} \langle \cdot \rangle ^n= \left\{ \begin{array}{l l l} 0 &{} \,\,\text {for} &{} \,\, (\cdot ) \le 0\; \\[2mm] (\cdot )^n &{}\,\, \text {for} &{} \,\, (\cdot ) > 0.\\ \end{array} \right. \end{aligned}$$
(16)

This additional energy will, appropriate regularization parameters q and b presumed, rapidly increase once the threshold value \(D_\mathrm {c}\) is exceeded. These parameters may also be set in order to obtain a correct behavior in terms of crack growth rates. The current state of fatigue damage is generally determined as sum of previous damage \(D_0\) and the damage increment associated with a certain period of time \(\mathrm {d}D_\text {f}\), which is in general dependent on the load amplitude, the midpoint load of the particular cycle, the load ratio, the stress rate and also on the damage history. Considering the complexity of fatigue, an accurate approach has to be incorporated within the phase field model as estimate for \(\mathrm {d}D_\mathrm {f}\), which is also reasonable in terms of computation. Different approaches have been proposed. These approaches can be characterized between those being of rather phenomenology kind as e.g., in [8, 28], and contributions based on continuum mechanics consistent with the second law of thermodynamics (e.g., [9, 14]). Even if probably any approach could be integrated in the proposed model, it is less convenient to apply a micromechanical approach in this case as one would end up at a CDM-related model. Two basic approaches within the field of fatigue lifetime estimation are Miner’s rule [28] and the fatigue damage model proposed in [8]. Both models can be formulated in terms of the damage increment, which is

$$\begin{aligned} \mathrm {d}D_{\text {fMiner}}=\frac{1}{N_{Fi}} \mathrm {d}N \end{aligned}$$
(17)

according to the Miner’s rule. \(N_{Fi}\) is the number of bearable cycles at a certain stress amplitude and stress ratio, which is given by a S–N curve (see Fig. 1) to indicate the bearable stress magnitude as function of cycle numbers \(N_{Fi}\). The data for these curves are obtained by fatigue experiments.

Fig. 1
figure 1

Schematic illustration of S–N line with explanation of several important quantities

A straight line is obtained using a double logarithmic scale. In contrast, within the fatigue damage model from Chaboche the increment is formulated as

$$\begin{aligned} \mathrm {d}D_{\text {fChaboche}}=D^{\alpha _1} \left( \frac{\sigma _{\text {max}}-\sigma _{\text {mean}}}{f(\sigma _{\text {mean}})},\right) ^{\alpha _2} \mathrm {d}N \end{aligned}$$
(18)

where the parameters \(\alpha _1\), \(\alpha _2\) and the effect of the mean stress \(f(\sigma _{\text {mean}})\) have to be determined in experiments. Separation of variables and integration yields the accumulated damage states

$$\begin{aligned} D_{\text {fMiner}}=\frac{N}{N_{Fi}} \;\; \text {and} \;\; D_{\text {fChaboche}}=\left( \frac{N}{N_{Fi}}\right) ^{\frac{1}{1-\alpha _1}}, \end{aligned}$$
(19)

for single-level loading. Considering Eqs. (17)–(19), one can simply identify the main difference between the two approaches. Within Miner’s rule, the effect of the initial state of damage is neglected, and therefore, the damage accumulation is linear. The approach of Chaboche on the other hand reveals a nonlinear curve if the damage is plotted as function of the number of cycles. The consequences are illustrated by the plot in Fig. 2. For the Miner rule, no sequence effect occurs, whereas the final damage state within Chaboche’s model clearly differs for different load sequences. In the following, the extension of the introduced phase field model will be discussed exemplary in terms of the Miner’s rule due to its convenient way of implementation. However, in [36], we showed that even by adapting this linear accumulation method, the global response in terms of fatigue crack growth predicted by our model shows a significant sequence effect. Furthermore, a direct connection between Miner’s rule of damage accumulation and Paris’ law for the description of crack growth rates can be shown (see [10]). The model proposed in [36, 37] was capable of dealing with constant amplitude loading only, which means cyclic loading constant in amplitude and mean load. However, in reality, components are subjected to load time histories, which vary in amplitude and mean load and are therefore rather considered as complex. An example for a sequence of various amplitude loading is illustrated in Fig. 3. As the amplitude of the particular cycles and also the mean load has a crucial effect on the fatigue crack growth behavior, these quantities have to be taken into account within the damage increment. Thus, we define the ratio of the mean external load with respect to the maximum load following the envelope of the load time history of the external load as the mean load ratio L. According to linear elasticity, the damage increment due to cyclic fatigue is formulated utilizing this ratio L as

$$\begin{aligned} \mathrm {d}D=\frac{\mathrm {d}N}{n^D}\left( \frac{\tilde{\sigma }(1-L)}{f(L)}\right) ^k \end{aligned}$$
(20)

with knee point cycle number \(n^D\) and slope factor k of the S–N curve. Increments from amplitudes below the threshold \(f(L)= A_D(1-L)^\beta \) with parameters \(A_D\) and \(\beta \) are not taken into account for damage accumulation. Integrating this increment in Eq. (14), the total internal Energy \(E^\text {f}\) of a component that may contain cracks caused by quasi static or under critical cyclic loading is

$$\begin{aligned} \begin{aligned} E^{\text {f}}&=\int \limits _{\mathrm {\Omega }} \left[ (g(s)+\eta )\frac{1}{2}{{\varvec{\varepsilon }}} : \left( \mathbb {C} {{\varvec{\varepsilon }}}\right) + \mathcal {G}_c\left( \frac{(1-s)^2}{4\epsilon }+ \epsilon |\nabla s|^2 \right) \right] \mathrm {d}V\\&\quad +\int \limits _{\mathrm {\Omega }} h(s) q\left\langle D_0+\frac{\mathrm {d}N}{n^D}\left( \frac{\tilde{\sigma }(1-L)}{f(L)}\right) ^k-D_\mathrm {c} \right\rangle ^b \mathrm {d}V , \end{aligned} \end{aligned}$$
(21)
Fig. 2
figure 2

Sequence effect of Chaboche’s model and comparison with the Miner’s rule

Fig. 3
figure 3

Schematic illustration of various amplitude load sequence consisting of five different blocks

with \(D_0\) being the previous state of fatigue damage, which within the following transient simulation framework is treated as history variable. The parameters \(A_D\), k, and \(n^D\) can be obtained from respective S–N curves. The aim in this concern must be to get a corresponding Paris law from a number of simulations with different \(\Delta K\) range. This can be achieved as shown in [37]. Considering the energy in Eq. (21), the stress reveals an additional component besides the usual static stresses with

$$\begin{aligned} {{\varvec{\sigma }}}= (g(s)+\eta )\mathbb {C}{{\varvec{\varepsilon }}}+h(s)q b\langle D_\text {f}-D_\text {c} \rangle ^{(b-1)}\frac{\mathrm {d} D_\text {f}}{\mathrm {d} {{\varvec{\varepsilon }}}}. \end{aligned}$$
(22)

The second term in Eq. (22) is considered as component accounting for micro-stresses caused by the fatigue mechanism. Note that, a small deformation framework with a linear elastic material law is used within this work. Accordingly, the model may be restricted to the range of high cycle fatigue, where \(N_F\ge 1000\). Within fatigue, the number of cycles to reach a certain state of damage or crack length is the essential. Therefore, the phase field evolution Eq. (7) is transferred into the cycle domain by chain rule differentiation. The evolution can then be formulated in terms of number of cycles as

$$\begin{aligned} \begin{aligned} \frac{\mathrm {d}s}{\mathrm {d}N}&=-\hat{M} \left\{ g'(s) \frac{1}{2} {{\varvec{\varepsilon }}} : \left( \mathbb {C}{{\varvec{\varepsilon }}}\right) +h'(s)P(D_\text {f}) -\mathcal {G}_c \left( 2 \epsilon \mathrm {\Delta } s -\frac{s-1}{ 2 \epsilon } \right) \right\} , \end{aligned} \end{aligned}$$
(23)

Generally, the large number of cycles, which are to be applied within the fatigue range before a noticeable damage occurs, are an issue for simulations within this field. In order to keep simulation times within an acceptable range, cycles of similar quantities are collected into blocks and the additional fatigue damage handled within one simulation step refers to these blocks. The block size is chosen adaptively to control the rate of damage. The idea of this scheme is referred to as “cycle jump” and was proposed in [14]. In detail, the cycle scheme we applied within our simulations deals with a fixed block size (cycles per simulation step). In order to control this block size and also to ensure a proper convergence behavior even for situations of a rapid increase of the energy contribution P adaptive step size adjustment is applied. The state of fatigue damage of the previous simulation step is treated as history variable and accordingly the damage at a certain state or simulation step i, respectively, is

$$\begin{aligned} D_i=D_{i-1}+\mathrm {d}D_{\text {f}}. \end{aligned}$$
(24)

Also worth to note is that an arbitrary load time sequence is handled within this simulation scheme using a proportional load command, which follows the envelope curve illustrated in Fig. 3. In case of sequences consisting of cycles, which are not only pulsating the ratio L is varied with the particular cycle numbers.

3 Numerical implementation

The phase field model for fatigue failure proposed in the previous sections is discretized and implemented into a finite element framework. With the evolution equation for the phase field Eq.  (23), the governing set of partial differential equations for the unknown displacements \({{\varvec{\mathrm {u}}}}\) and the order parameter s is

$$\begin{aligned} \begin{aligned} {{\varvec{0}}}&= \text {div}\ {{\varvec{\sigma }}}\\ \frac{\mathcal {G}_c}{2 \epsilon }&=\frac{1}{\hat{M}} \frac{\mathrm {d}s}{\mathrm {d}N}+2s(W+P) +s \frac{\mathcal {G}_c}{2\epsilon } - \mathcal {G}_c 2 \epsilon \mathrm {\Delta }s , \end{aligned} \end{aligned}$$
(25)

where volume forces are neglected within the mechanical equilibrium and the quadratic degradation function Eq. (6) is incorporated. The respective weak forms are obtained if both equations are multiplied with virtual quantities \(\delta {{\varvec{\mathrm {u}}}}\) and \(\delta s\), respectively. This yields

$$\begin{aligned} 0=-\int \limits _{\mathrm {\Omega }} {{\varvec{\sigma }}} : \delta {{\varvec{\varepsilon }}} \; \mathrm {d}V+\int \limits _{\partial \mathrm {\Omega }_t} \delta {{\varvec{\mathrm {u}}}} \cdot {{\varvec{\mathrm {t}^*}}} \; \mathrm {d}S \end{aligned}$$
(26)

with vector \({{\varvec{\mathrm {t}}}}^*\) of prescribed tractions that act on the boundary \(\partial \mathrm {\Omega }_t\) and

$$\begin{aligned} \begin{aligned} 0=-\int \limits _{\mathrm {\Omega }} \delta s \frac{\mathrm {d}s}{\hat{M}\mathrm {d}N} \;\; \mathrm {d}V -\int \limits _{\mathrm {\Omega }} \delta s \left( s {{\varvec{\varepsilon }}} : \left( \mathbb {C} {{\varvec{\varepsilon }}}\right) +2P +(s-1) \frac{\mathcal {G}_c}{2\epsilon } \right) \;\; \mathrm {d}V -\int \limits _{\mathrm {\Omega }} 2 \mathcal {G}_c \epsilon \nabla s \cdot \nabla \delta s \;\; \mathrm {d}V, \end{aligned} \end{aligned}$$
(27)

where integration by parts was applied. The following derivations are made in terms of a 2d setting and under using Voigt’s notation for symmetric tensors. The field quantities \({{\varvec{\mathrm {u}}}}\) and s are discretized by means of shape functions \(H_I\) with \(I=1,\dots ,n\) and n the number of nodes per element. Accordingly,

$$\begin{aligned} {{\varvec{\mathrm {u}}}}({{\varvec{\mathrm {x}}}})=\sum _{I=1}^n H_I({{\varvec{\mathrm {x}}}}) {{\varvec{\mathrm {u}}}}_I \quad \mathrm{and} \quad s({{\varvec{\mathrm {x}}}})=\sum _{I=1}^n H_I({{\varvec{\mathrm {x}}}}) s_I, \end{aligned}$$
(28)

with nodal quantities \({{\varvec{\mathrm {u}}}}_I\) and \(s_I\). The strain matrix \({{\varvec{\varepsilon }}}\) and the spatial gradient \(\nabla s\) are discretized by

$$\begin{aligned} {{\varvec{\varepsilon }}}({{\varvec{\mathrm {x}}}})=\sum _{I=1}^n {{\varvec{B}}}_I^{{{\varvec{{{\varvec{\mathrm {u}}}}}}}}({{\varvec{\mathrm {x}}}}){{{\varvec{\mathrm {u}}}}}_I \quad \mathrm{and}\quad \nabla s ({{\varvec{\mathrm {x}}}})=\sum _{I=1}^n {{\varvec{B}}}_I^s({{\varvec{\mathrm {x}}}}){s}_I, \end{aligned}$$
(29)

using the respective operator matrices

$$\begin{aligned} {{\varvec{B}}}_I^{{{\varvec{{{\varvec{\mathrm {u}}}}}}}}({{\varvec{\mathrm {x}}}})=\begin{bmatrix} H({{\varvec{\mathrm {x}}}})_{I,x} &{} 0\\ 0 &{} H({{\varvec{\mathrm {x}}}})_{I,y}\\ H({{\varvec{\mathrm {x}}}})_{I,y} &{}H({{\varvec{\mathrm {x}}}})_{I,x} \end{bmatrix} \quad \mathrm{and} \quad {{\varvec{B}}}^s_I =\begin{bmatrix} H({{\varvec{\mathrm {x}}}})_{I,x}\\ H({{\varvec{\mathrm {x}}}})_{I,y} \end{bmatrix}. \end{aligned}$$
(30)

The internal forces corresponding to the displacements \({{\varvec{\mathrm {u}}}}\) and s at node I are derived from Eqs. (26) and (27) as

$$\begin{aligned} {{\varvec{F}}}^{{\varvec{{{\varvec{\mathrm {u}}}}}}}_I=\int \limits _{\mathrm {\Omega }} \left( {{\varvec{B}}}_I^{{\varvec{{{\varvec{\mathrm {u}}}}}}}\right) ^T {{\varvec{\sigma }}} \;\mathrm {d}V \end{aligned}$$
(31)

and

$$\begin{aligned} \begin{aligned} F^s_I =\int \limits _{\mathrm {\Omega }} H_I \hat{M}^{-1} \frac{\mathrm {d}s}{\mathrm {d}N} \; \mathrm {d}V + \int \limits _{\mathrm {\Omega }} H_I\left( s {{\varvec{\varepsilon }}}^T \mathbb {C} {{\varvec{\varepsilon }}} + 2sP(D_\text {f})+\mathcal {G}_c \frac{s-1}{2\epsilon }\right) \; \mathrm {d}V + \int \limits _{\mathrm {\Omega }} 2\mathcal {G}_c\epsilon \left( {{\varvec{B}}}_I^s\right) ^T \nabla s \; \mathrm {d}V, \end{aligned} \end{aligned}$$
(32)

where the same shape functions \(H_I\) are used for virtual quantities. In order to find a solution of the nonlinear system, the Newton-Rapson scheme is applied. The associated components of the current tangential stiffness matrix are obtained by derivation of the residual with respect to the unknown field quantities as

$$\begin{aligned} {{\varvec{K}}}_{IJ}^{{{\varvec{{{\varvec{\mathrm {uu}}}}}}}}&=\frac{\partial {{\varvec{F}}}_I^{{\varvec{{{\varvec{\mathrm {u}}}}}}}}{\partial {{\varvec{\mathrm {u}}}}_{J}}=\int \limits _{\mathrm {\Omega }} \biggl \lbrace {{\varvec{B}}}_I^{{\varvec{{{\varvec{\mathrm {u}}}}}}}(s^2+ \eta )\mathbb {C} {{\varvec{B}}}_J^{{{\varvec{{{\varvec{\mathrm {u}}}}}}}} + s^2 {{\varvec{B}}}_I^{{\varvec{{{\varvec{\mathrm {u}}}}}}} \frac{\mathrm {d} ^2 P}{\mathrm {d {{\varvec{\varepsilon }}}}^2} {{\varvec{B}}}_J^{{{\varvec{{{\varvec{\mathrm {u}}}}}}}} \biggl \rbrace \; \mathrm {d}V, \end{aligned}$$
(33)
$$\begin{aligned} {{\varvec{K}}}_{IJ}^{s{{\varvec{{{\varvec{\mathrm {u}}}}}}}}&=\frac{\partial {{\varvec{F}}}_I^s}{\partial {{\varvec{\mathrm {u}}}}_{J}}= \int \limits _{\mathrm {\Omega }} H_I 2 s \left( \mathbb {C} {{\varvec{\varepsilon }}} + \frac{\mathrm {d}P}{\mathrm {d {{\varvec{\varepsilon }}}}} \right) ^T {{\varvec{B}}}_J^{{\varvec{{{\varvec{\mathrm {u}}}}}}} \; \mathrm {d}V, \end{aligned}$$
(34)
$$\begin{aligned} {{\varvec{K}}}_{IJ}^{{{\varvec{{{\varvec{\mathrm {u}}}}}}}s}&=\frac{\partial {{\varvec{F}}}_I^{{\varvec{{{\varvec{\mathrm {u}}}}}}}}{\partial s_{J}}= \int \limits _{\mathrm {\Omega }} \left( {{\varvec{B}}}_I^{{\varvec{{{\varvec{\mathrm {u}}}}}}} \right) ^T 2 s H_J \left( \mathbb {C} {{\varvec{\varepsilon }}} + \frac{\mathrm {d}P }{\mathrm {d {{\varvec{\varepsilon }}}}} \right) \mathrm {d}V, \end{aligned}$$
(35)
$$\begin{aligned} {{\varvec{K}}}_{IJ}^{ss}&=\frac{\partial {{\varvec{F}}}_I^s}{\partial s_{J}}= \int \limits _{\mathrm {\Omega }} \biggl \lbrace H_I H_J \left( {{\varvec{\varepsilon }}}^T \mathbb {C} {{\varvec{\varepsilon }}} +2P +\frac{\mathcal {G}_c}{2\epsilon } \right) + 2\mathcal {G}_c \epsilon \left( {{\varvec{B}}}_I^s \right) ^T {{\varvec{B}}}_J^s \biggl \rbrace \; \mathrm {d}V \end{aligned}$$
(36)

with

$$\begin{aligned} \frac{\mathrm {d}P}{\mathrm {d {{\varvec{\varepsilon }}}}} = qb \langle D_\mathrm {f}-D_\mathrm {c} \rangle ^{b-1} \frac{\mathrm {d}D_\mathrm{f}}{\mathrm {d}\tilde{\sigma }} \frac{\mathrm {d}\tilde{\sigma }}{\mathrm {d} {{\varvec{\varepsilon }}}}. \end{aligned}$$
(37)

The choice of an adequate driving force quantity \(\tilde{\sigma }({{\varvec{\varepsilon }}})\) may depend on the material and also on the type of crack extension expected. For example, if only unidirectional loading is expected, the normal stress component in loading direction, which can be formulated as

$$\begin{aligned} \tilde{\sigma }={{\varvec{ 1 }}}^T\mathbb {C}{{\varvec{\varepsilon }}} \end{aligned}$$
(38)

with \({{\varvec{ 1 }}}=\left( 0,1,0\right) ^T\) for a vertical load direction, is a sufficient driving force quantity for the fatigue mechanism. As a more general driving force quantity in terms of mixed mode loading, the first principal stress should be utilized. This is a consequence of the maximum tangential stress criterion [12], which states that within linear elastic fracture mechanics a crack will extend perpendicular to the direction of the largest tension around the crack tip. For this two-dimensional setting, the first principal can be formulated as

$$\begin{aligned} \tilde{\sigma }={{\varvec{ 1 }}}_+^T \mathbb {C}{{\varvec{\varepsilon }}}+\left[ \left( {{\varvec{ 1 }}}_-^T\mathbb {C}{{\varvec{\varepsilon }}} \right) + \left( {{\varvec{ 1 }}}_{\tau }^T\mathbb {C}{{\varvec{\varepsilon }}} \right) \right] ^{1/2} \end{aligned}$$
(39)

with operators \({{\varvec{ 1 }}}_+ = \left( \tfrac{1}{2},\tfrac{1}{2},0 \right) ^T\), \({{\varvec{ 1 }}}_- = \left( \tfrac{1}{2},-\tfrac{1}{2},0 \right) ^T\) and \({{\varvec{ 1 }}}_{\tau } = \left( 0,0,1 \right) ^T\). Finally, the damping matrix can be derived which reveals zero contributions except for the component

$$\begin{aligned} D_{IJ}^{s\dot{s}}=\frac{\partial F_I^s }{\partial s,_N}= \int _{\mathrm {\Omega }} \frac{H_I H_J}{ \hat{M}} \; \mathrm {d}V. \end{aligned}$$
(40)

4 Numerical examples

In order to illustrate the behavior of the proposed model, different numerical examples and their results are presented in the following. Properties of interest within the field of fatigue crack growth are on one hand the crack growth rates \(\mathrm {d}a/\mathrm {d}N\), and on the another hand, the path a crack follows while extending under load. For both of these properties, different setups where created and simulations where performed. The software FEAP 8.4 was used for these simulations, where the model described in the previous sections was implemented as quadrilateral user element with linear shape functions for both the displacements and the phase field.

4.1 Effect of mean load

As outlined before, the mean load of an applied load time sequence has a very high effect on the fatigue life of a structure. The number of bearable load cycles will be significantly lower for sequences with a higher mean load of the load cycles compared to a sequence with the same amplitude but lower mean load. Within the literature, this effect is referred to as mean stress effect and is characterized in terms of the stress ratio \(R=\sigma _{\text {min}}/\sigma _{\text {max}}\). How this effect is incorporated in our phase field model was discussed in Sect. 2.2. To verify the models response to a variation of the mean load simulations with compact tension (CT) test specimens were performed. This specimen geometry is commonly used for an experimental investigation of the crack growth behavior for quasi static as well as for cyclic loading. The geometry and the dimensions of the used CT-specimen according to the ASME standard [3] is illustrated in Fig. 4. A tensile load in vertical direction was applied within these simulations by a distributed force at the upper hole of the specimen. The contour plot in Fig. 4 shows the evolution of the crack field obtained from a simulation using the proposed phase field model. A crack indicated by areas with \(s=0\) is visible. This crack grows in horizontal direction as a result of the cyclic load. For the following evaluation of crack growth rates, the approximation for the stress intensity factor range \(\Delta K\) proposed in [3] was used, where the crack length is defined as horizontal distance from the crack tip to the center of the holes.

Fig. 4
figure 4

Details of CT-specimen according to [3] with indication of contour plot of phase field obtained from cyclic phase field simulation

Fig. 5
figure 5

Comparison of energy contribution ratio (W\(=\)strain energy, P\(=\) fatigue energy, E\(=\)total energy) of simulations for different stress ratios

Fig. 6
figure 6

Results from simulations with different stress ratios: a crack growth rates (Paris’ law), b curves of crack length vs number of applied load cycles

To evaluate the model behavior regarding mean stress effects several simulations with identical amplitudes and different stress ratios were performed. In Fig. 5, the schematic difference of the six simulations regarding the stress ratio is indicated. As explained by many textbooks (see e.g., [11, 19, 31, 33]), based on experimental data, the crack growth curves (Paris’ law [30])

$$\begin{aligned} \frac{\mathrm {d}a}{\mathrm {d}N}=C \Delta K^m \end{aligned}$$
(41)

with constants C and m are also dependent on the stress Ratio R, where higher R-ratios generally shift the crack growth curve to a level of higher rates. In this regard, a quiet convenient property of the presented phase field model is revealed by a comparison of the relative energy contributions for the state of crack initiation, which is given in Fig. 5. The bar plot shows that the static strain energy delivers a significantly higher contribution to the total energy for stress ratios with a higher mean load. With other words, less accumulated fatigue energy is required to support crack growth at high stress ratios than for lower stress ratios. This behavior is also reflected by the Paris diagram of Fig. 6a, which illustrates the basic characteristic of fatigue crack growth for the six simulations in terms of crack growth rates. The plot consistently indicates higher rates \(\mathrm {d}a/\mathrm {d}N\) for equal stress intensity factor ranges with increasing stress ratio R. Consequently, the evolution of the crack length with respect to the absolute number of load cycles also increases with the stress ratio as the curves in Fig. 6b underline. For these simulations, a value of 0.3 was used for parameter \(\beta \) of the function f(L) in the damage increment Eq. (20). This function controls the threshold stress for local fatigue growth and determines the mean stress effect. For the case of a stress ratio of \(R=0.5\), the factor \(\beta \) was varied. Fig. 7 illustrates that the deviation of the crack growth rates with respect to the basic alternating load decreased significantly by using smaller values for \(\beta \). The choice of the function f(L) within this study is of purely exemplary origin in order to illustrate how an adaption of the proposed phase field model can be accomplished. We do not claim that this kind of function provides a general method for any material. However, the model allows for arbitrary functions for the threshold stress in order to properly fit a certain behavior.

Fig. 7
figure 7

Curves of crack growth rates obtained for simulations with the CT-specimen

Fig. 8
figure 8

Example of crack interaction as test case for the propose model: a simulation setup, b contour plots of the phase field variable s after different numbers of applied load cycles

4.2 Mixed mode loading

Dealing with fracture mechanics, it is a very crucial task to gain knowledge about the path a crack will follow, since, for instance, this information may be included in a fail safe design of a structure. For several cases, this path may be very obvious, as, for example, in the case of an mode-I loading of the CT-specimen of the previous section, but within practical applications variations of conditions occur in terms of load, deformation and crack interaction. This can cause a very complex crack growth behavior. In order to validate the presented model regarding the predicted crack growth direction, the geometry described in Fig. 8a was set up for a phase field simulation. As within this simulation crack, curving is expected, and the criterion of a fatigue driving force has to be generalized. Accordingly, Eq. (39) was used as driving stress measure. The specimen for this simulation is a narrow stripe with sharp notches on both vertical edges. The notches are positioned with a small distance in vertical direction. The lower edge is fixed, and a displacement load in vertical direction is applied on the upper edge, where within this simulation the cyclic load was constructed as purely alternating (\(R=-1\)). Contour plots of the phase field variable s are presented in Fig. 8b. Below the contour plots, the number of applied load cycles to produce the particular state of crack evolution is indicated. These results show that a fatigue crack initiates at about 28000 cycles, and the initial direction of crack extension is purely horizontal. After this period of horizontal crack growth, the cracks start to curve such that the two cracks continue to grow away from each other (\(N\approx 40000\)). With ongoing cycling, the cracks follow this direction for a low number of cycles before both cracks again change their direction of growth in the opposite direction. After another 6000 cycles, the cracks obviously start to grow toward each other. The simulation was stopped after about 55000 cycles. This state is shown in the last contour plot. The question whether this final crack patter is plausible can be addressed by an analytic solution of the crack tip fields within linear elastic fracture mechanics. For a plane deformation of a crack, two basic modes of loading are differentiated. The mode-I case, where a crack is loaded in normal direction with respect to its orientation (e.g., the example for the previous section) and the mode-II case, where the crack is loaded such that the crack surfaces undergo antisymmetric displacements along the crack. A crack will extend straight for a pure mode-I load, and it will grow under an angle of approximately \(-70.5^{\circ }\) (see e.g., [17]) for a pure mode-II load. For an arbitrary load case, proportionality factors for the particular modes, namely \(K_I\) and \(K_{II}\), contribute to the stress and displacement field near the crack tip. Accordingly, a crack that deviates from its straight path undergoes a load with a significant \(K_{II}\) contribution. For two cracks of equal length and orientation in an infinite plate under normal loading, the \(K_{II}\) contribution can be approximated (see [17]) by

$$\begin{aligned} K_{II}=K_0 \frac{a^2}{2d^2}\left[ \sin (4 \theta )-\sin (2 \theta ) \right] . \end{aligned}$$
(42)

In this equation, 2a is the crack length, d is the distance of the cracks and \(\theta \) is the skew angle, which becomes \(\pi /2\) for stacked cracks and zero for collinear cracks. Equation (42) implies that if cracks grow toward each other, the \(K_{II}\) contribution is at first positive until a certain angle \(\theta \) is reached such that \(K_{II}\) changes its sign. This explains the crack paths observed in the simulation. The presented simulation results can also be verified by the crack patterns presented in [43], which were obtained for cyclic experiments and show the same behavior.

5 Conclusion

We presented a phase field fracture model, which in addition to quasi static fracture is also capable to predict fatigue crack growth. The focus of our development lies on relevant loading conditions as various loading amplitudes and arbitrary crack paths. It was shown that the effect of the mean stress is included in the model formulation. Further, we give explanations how a specific material behavior can be incorporated. In a second example, the behavior in terms of the crack growth direction is investigated by means of a mixed mode load case. The obtained crack pattern is confirmed by experimental findings from the literature as well as analytical considerations.