A phase field modeling approach of cyclic fatigue crack growth

Phase field modeling of fracture has been in the focus of research for over a decade now. The field has gained attention properly due to its benefiting features for the numerical simulations even for complex crack problems. The framework was so far applied to quasi static and dynamic fracture for brittle as well as for ductile materials with isotropic and also with anisotropic fracture resistance. However, fracture due to cyclic mechanical fatigue, which is a very important phenomenon regarding a safe, durable and also economical design of structures, is considered only recently in terms of phase field modeling. While in first phase field models the material’s fracture toughness becomes degraded to simulate fatigue crack growth, we present an alternative method within this work, where the driving force for the fatigue mechanism increases due to cyclic loading. This new contribution is governed by the evolution of fatigue damage, which can be approximated by a linear law, namely the Miner’s rule, for damage accumulation. The proposed model is able to predict nucleation as well as growth of a fatigue crack. Furthermore, by an assessment of crack growth rates obtained from several numerical simulations by a conventional approach for the description of fatigue crack growth, it is shown that the presented model is able to predict realistic behavior.


Introduction
The field of cyclic mechanical fatigue is one of the most important branches of post and current research in engineering. Catastrophic failures of different types of technical structures like trains, airplanes, bridges, turbine rotors, pressure vessels and others, occurred as consequence of the fatigue of materials. Accordingly, this phenomenon must be kept in the focus of present studies. From a material science perspective, the fatigue mechanism principles are mostly understood and explained in several textbooks by e.g. Dowling (2013), Schijve (2009), Haibach (2006. The total fatigue life of a structure is generally divided into two phases. The crack nucleation phase describes the period in which a microcrack is initiated as a consequence of irreversible cyclic slip in grains with proper orientation and the growth of this microcrack to an order where it is visible with the naked eye or at least by utilizing a low rate magnifying glass. The second period contains the propagation of the macrocrack until final failure. For both of these periods different phenomenological prediction frameworks exist. The total life time approach is mainly based on so-called Wöhler curves or S-N curves, respectively. These curves assign a certain number of bearable load cycles to a certain load amplitude. The particular cycle numbers may serve as input for damage accumulation laws, where the linear accumulation hypothesis referred to as the Miner rule (Miner 1945) is widely applied. The most popular method to describe macrocrack growth is Paris' law (Paris and Erdogan 1963). Within this method the cyclic crack growth rate da/dN is described as a function of the cyclic stress intensity factor K , which basically incorporates effects of external load, actual crack length and geometry. It was found that, at least within the region of sub critical crack extension, this function plotted on a double logarithmic scale is a straight line. Accordingly the crack growth rate is approximated by a power law. Extensions of this law were proposed to include effects like mean stress or stress ratio (see e.g. Dowling 2013). Paris' law is also the basis for the most widely used numerical tool for fatigue crack growth simulations namely NASGRO (Forman et al. 2005).
Essential studies attempt to model cyclic fatigue failure using the framework of Continuum Damage Mechanics (CDM), which was comprehensively outlined by Lemaitre (1992) to describe the evolution of a damage variable, were presented by e.g. Marigo (1985) or Chow and Wei (1991). However, within these models, the fatigue life is modeled only for the period until crack nucleation. A numerical model for the simulation of fatigue crack growth was introduced by Fish and Oskay (2005). The damage evolution within this work is modeled by the growth of the void volume fraction, which is governed by extensions of Gurson's failure model. To obtain a crack or its sharp interfaces, respectively, finite elements are deleted once a critical void volume fraction is reached.
As an alternative to conventional crack growth simulation techniques, phase field modeling for fracture has gained attention within the past decade. This method provides a very useful tool for the numerical simulation of problems dealing with sharp interfaces, as it does not require procedures for mesh disconnection, element deletion, and re-meshing for an explicit modeling of the interfaces. The phase at a certain location is incorporated by introducing additional degrees of freedom. In the context of fracture, the framework was applied to dynamic fracture (see e.g. Schlüter et al. 2014), and quasi static brittle fracture in e.g. Kuhn and Müller (2010); Borden et al. (2014); Miehe et al. (2010). Extensions of these model exist for ductile fracture (e.g. Kuhn et al. 2016;Borden et al. 2016) and also effects like anisotropic fracture resistance were studied by e.g. Teichtmeister et al. (2017), Schreiber et al. (2017), Hakim and Karma (2009). Models for fatigue damage were so far presented by Alessi et al. (2017) and Seiler et al. (2018). In these models the phase field is driven by degrading the energy release rate once a damage parameter increases.
In this work, the fatigue mechanisms occurring in a material are interpreted as additional driving force contribution, contribution that originates from cyclically accumulated deformation work. The formulation is an extension of the model proposed by Kuhn and Müller (2010) for brittle fracture and as a linear elastic material model is considered within our work, the range of validity in terms of the number of cycles to failure N F may be restricted to high cycle fatigue, which covers N F ≥ 1000 for many metallic materials. The driving force for the phase field variable s(x, t) incorporates an additional energy contribution accounting for fatigue damage accumulation due to an increasing number of load cycles. The core idea of a phase field fracture model is to introduce an additional field variable to represent cracks. A continuous transition of this field variable s(x, t) within the interval [0, 1] approximates the crack faces. The crack field is 1 as long as the material remains undamaged and drops once fracturing occurs, where the lower limit s = 0 represents cracks. The phase field model from Kuhn and Müller (2010) uses the regularized formulation of a variational model for brittle fracture (Bourdin et al. 2000). As an extension of Griffith's theory, which states that a crack propagates once the energy release rate reaches a certain critical value, it is postulated that the displacement field u and the fracture field s locally minimize the total energy E within a loaded body. (1) In Eq. (1) W is the elastic strain energy density, which is multiplied by a degradation function g(s) to account for the stiffness reduction in cracked regions. The surface energy density contribution accounts for the energy required for the formation of a crack surface. The parameter η shall be chosen 0 < η 1 to avoid numerical problems within the static solution scheme. The elastic strain energy density can be formulated in terms of the infinitesimal strain tensor where ∇u is the spacial gradient of the displacements. In Eq.
(2) the Voigt notation for symmetric tensors is used for the linearized strain tensor ε and the stiffness tensor C . The second energy contribution in Eq. (1) is the fracture energy density which besides a local part, further consists of a nonlocal part, characterized by the spacial gradient of s. The two contributions are multiplied by the critical energy release rate G c . The length scale controls the width of the transition zone between undamaged and broken material. This crack surface contribution regularizes Griffith's theory, which proposes the proportionality of the fracture energy and the crack surface. Assuming a hyperelastic material, the stresses can be derived from the potential ψ by Different types of degradation functions g(s) are proposed in the literature (Bourdin et al. 2000;Kuhn et al. 2015;Borden et al. 2016). However, this function has to fulfill several criteria as g(0) = 0, g(1) = 1 and also g (0) = 0. These requirements are satisfied by the quadratic function which is used in the original regularized formulation proposed by Bourdin et al. (2000). Discussions about a quadratic and a cubic degradation function are given in Kuhn et al. (2015). These functions were shown to shift the stress level at which crack initiation occurs to higher stresses and also to conserve the linear elastic stress strain behavior until crack initiation. A mixed formulation for the degradation function was proposed by Borden et al. (2016). Notice, that for β = 2, g β becomes g b and for β → 0, g β becomes the cubic degradation function introduced by Kuhn et al. (2015). Accordingly, the model behavior is effected by varying the parameter β.
A generalized Ginzburg-Landau equation (Gurtin 1996) is used to obtain the time evolution of the phase field s viȧ where is the Laplace operator. The scalar mobility parameter M in Eq. (8) acts as a viscous control parameter to approach the stationary limit, which is specified by δψ δs = 0. The limit case M → ∞ approximates quasi static conditions.

Phase field formulation for cyclic fatigue
To provide a correct understanding of the model described in the following, it is important to comment on the interpretation of the phase field variable s. In particular the transition zone between the two phases broken and intact (0 < s < 1) requires special consideration. A possible interpretation is to connect s to the fraction of micro voids within a representative volume element, as it was proposed by Voyiadjis and Mozaffari (2013). Here, s is brought in direct connection with the damage variable from the CDM framework (Lemaitre 1992). On the other hand, the phase field variable in Bourdin's original model (Bourdin et al. 2000) did not have a physical interpretation at all, but the scalar phase field s was rather incorporated in the energy formula-tion of a fractured body to enable a regularization of the crack energy. According to this explanation, the meaning of s results from the phase field philosophy itself, which is to approximate sharp interfaces by a smooth transition of an order parameter for numerical reasons.
In the phase field model presented within this work, we follow this second interpretation. The term damage refers to fatigue damage, which will be applied to drive the phase field, but it is not explicitly connected to it. Specifically, s is driven by a proper estimate for the fatigue damage increment dD f occurring in the material within a certain cycle interval.

Phase field fracture model for fatigue failure
As explained in Sect. 2.1, a sequence of local minimizers of the total internal energy, depending on displacement and fracture fields, has to be found in order to obtain the correct crack pattern, or in other words information for crack initiation, kinking or branching, respectively. Suppose a sample is exposed to a monotonically increasing load. In that case the solution reveals an evolution of the crack field, since at some point the surface part of Eq.
(1) will increase as it becomes energetically more favorable to form a crack, which means to decrease s, than to allow for more strain energy. Within the regime of cyclic fatigue the load maximums are rather small to very small compared to quasi static fracture loads. Presuming such a small load and a constant G c , the solution will not yield a decline of s in case of the classical model (Kuhn and Müller 2010), as energy consumed by irreversible processes associated with an increasing state of fatigue is not taken into account. Accordingly, crack growth would be too costly with respect to the total internal energy. A modification of the total energy formulation has been proposed by Alessi et al. (2017) in order to enable crack growth driven by cyclic fatigue. The modification basically applies to the crack energy dV , where the fracture resistance G c is degraded by a function of a plastic strain variable, which is accumulated over an increasing number of load cycles.
Within this work we present an alternative approach for a regularized formulation of the total energy where an additional energy density contribution P is introduced to account for the sum of additional driving forces associated with the mechanism of cyclic mechanical fatigue. The function h(s) degrades P in a similar way as W is degraded by g(s) once a crack propagates. The energy contribution E ac is governed by a piece-wise defined function with the threshold value D fc . According to this definition the amount of additional energy to enter Eq. (9) remains zero until a certain critical value D fc of the accumulated damage D f is reached and increases rapidly (governed by the model parameters q and b) after this limit is exceeded. The actual accumulated fatigue damage is in general described as which is the sum of the previous damage state D 0 and the damage increment dD f as a function of the amplitude σ a , the mean value of the stressσ , the load ratio R, the stain rateε, and the damage history or the previous damage state, respectively. The main task in this context is to find an accurate formula for the increment dD f . Equation (11) indicates that, generally, several effects have to be considered for an estimate of the damage increment. However, the choice has to be made also in consideration of computational effort, which a certain model may cause. Straightforward approaches for the fatigue lifetime estimation of a component, which are commonly used and were also incorporated within approaches for numerical fatigue estimations (see e.g. Bograd et al. 2008) are Miner's rule (Miner 1945) and the fatigue damage model proposed by Chaboche and Lesne (1988). The main difference between these two approaches is that Miner's rule does not consider the effect of the current state of damage D 0 on the increment, which is formulated as with the number of bearable load cycles N F i and the cycle increment dN . Considering Eq. (12) it becomes clear that the damage accumulation is linear and accordingly the influence of load sequence on the The particular critical cycle number N F i can be obtained by an appropriate S-N curve. These curves (see Fig. 1) are plots, obtained within fatigue experiments, of the stress magnitude for cyclic loading depending on the corresponding number of bearable load cycles N F . A straight line is obtained using a double logarithmic scale.
With Eq. (12) and an appropriate S-N curve, the total internal Energy E f of a component that may contain cracks caused by quasi static or by cyclic loading is then with the driving force quantity σ A , which is related to a stress amplitude. This variable will be explained within the next section. Other quantities introduced in Eq. (13) are the fatigue limit A D /2, the knee point cycle number n D and the slope factor k of the S-N curve. The · n denote Macauley brackets with the definition · n = 0 for (·) ≤ 0 (·) n for (·) > 0.
Considering the energy in Eq. (13) the stresses become The first term in Eq. (15) is the usual degraded stress tensor and the second term may be considered an additional stress component accounting for accumulated micro stresses as consequence of the fatigue mechanism.

Cycle resolved simulation scheme
Depending on the magnitude of the load cycles, the fatigue mechanism may require a very large number of cycles before consequences on a macro level like crack nucleation or extension can be observed. Accordingly, an efficient scheme is required to integrate large numbers of cycles with respect to a characteristic cycle number, otherwise simulation times cannot be kept within a reasonable limit. Even for conventional fatigue life approaches certain amounts of similar cycles are summarized to blocks, where the connection of all these blocks represents the whole load history called load collective. Fish and Yu (2002) proposed a scheme called "cycle jump", where the damage at the end of a simulation step including a certain number of cycles D (i+ N i ) is calculated via The damage at the previous step is D i and the particular D i values are obtained by an explicit Euler method. Within this simulation scheme the block size summarizing several cycles, N i is chosen adaptive. The time-cycle transfer scheme presented within this work is quite similar. Suppose a constant block size or number of cycles per time N t , respectively, leads to a certain evolution of damage. By utilizing Eq. (12), the damage of the actual time is The phase field evolution Eq. (8) can be transferred to the cycle domain using the chain rule. Accordingly, the evolution of the phase field can be reformulated in terms of cycles withM denoted the mobility in the cycle domain. The damage will be updated within every simulation step according to the change caused by a chosen cycle number. To ensure good convergence even for a rapid increase of P, adaptive step size adjustment must be implemented.

Numerical implementation
For the evaluation of the enhancement of this phase field model for fatigue failure the quadratic degradation function from Eq. (6)

is chosen for g(s) and h(s).
If volume forces are neglected within the mechanical equilibrium and the modified evolution equation for the phase field Eq. (18) is utilized, the governing set of differential equations for the unknown displacements u and the order parameter s is given by where ∇· is used as divergence operator. Multiplying both equations with the respective virtual quantities δu and δs and integration by parts, the particular weak forms are obtained with the prescribed traction vector t * on ∂ t and where the homogeneous Neumann boundary condition ∇s · n = 0 at ∂ is used. The isoparametric concept is used in order to obtain a discrete description of the field quantities u and s by means of the shape functions H I for the nodes I = 1, . . . , n with and where the denotation( ·) indicates nodal values. By introducing the differential operator matrices for a two dimensional setting and where (·), v indicates partial derivatives ∂(·) ∂v , the values for ε and ∇s can be computed from and where the Voigt notation is used. An analog discretization may be applied for the virtual quantities. Accordingly the contributions to the vector of the internal forces at node I become and In the following, we denote by σ A (ε) any representative component of the stress tensor or property of the stress tensor, chosen as the driving quantity for the fatigue mechanism. This work aims at investigating fundamentals of the integration of fatigue failure within a phase field formulation in terms of mode-I fracture. Accordingly, the normal stress component in loading direction with 1 = (0, 1, 0) T is considered as appropriate stress measure for cyclic tensile loading in vertical direction. In order to find a solution U * = (u * , s * ,ṡ * ) T of the nonlinear system, the Newton-Raphson scheme is applied. The particular components of the current stiffness matrix are found as and the damping matrix is zero, except for the component Using the implicit Euler method for the transient problem is sufficient as the number of cycles is directly proportional to time. However, automatic step size adjustment is applied to ensure a stable computation when the phase field decreases due to rapid increase of the fatigue damage.

Examples and validation
For an experimental characterization of the fatigue crack growth behavior of materials the so-called compact tension (CT) specimen is, among others, frequently used. Accordingly, this specimen geometry was also chosen for the conducted simulations within this study, whose results are presented in the following. The definition of the geometry of the C(T)-specimen is given in the ASTM E 399 standard (ASTM 2009). For the simulations the software FEAP 8.4 was used. The model was implemented into a user element routine for a 4-node quadrilateral element. Figure 2 shows the where each simulation step represents a certain increment of load cycles (see Fig. 3). As only constant amplitude loading was considered the applied force was kept constant after the maximum load was reached. For all numerical experiments a material with Young's modulus E = 2.1 · 10 5 MPa, Poisson's ratio ν = 0.3 and critical energy release rate G c = 2.33 · 10 3 J/m 2 was considered. Further, the parameters A D = 200 MPa and k = 5 were used.

Crack initiation and model characteristics
The main feature of the original phase field formulation for brittle fracture is a crack pattern that locally minimizes the total energy. The solution is allowed to yield cracks and cracks will appear once their formation is beneficial with respect to the total energy. However, as previously explained, the modified phase field model for fatigue fracture has an additional mechanism, which includes a fatigue related driving force. The plots in Fig. 4 show the evaluation of the stress contributions σ e and σ ac of Eq. (15) in vertical direction on a line, which extends in horizontal direction from the notch tip of the C(T)-specimen. Furthermore, the effective stress in vertical direction σ A is indicated. For simplicity all stresses are normalized by the threshold stress A D . A corresponding contour plot of the phase field s within a square section (edge length approximately 1 mm) around the notch tip is indicated in every illustration. Note that for this simulation the initial s field was 1 over the whole specimen. In other words, the specimen was not pre-cracked and the crack initiation period will be illustrated in the following. The plot in Fig. 4a shows the situation after the initial ramp up, i.e. complete load applied. The phase field remains 1 indicating intact material everywhere. The stress at the notch tip is slightly above A D . This state stays unaltered until a number of approximately 1.4 · 10 5 cycles was applied (Fig. 4b). At this point an increase of the σ ac contribution can be recognized, which indicates that the damage at the notch tip has grown to a relevant amount such that the energy density contribution P significantly affects the energy functional. Consequently the phase field decreases and the stresses are degraded. After another 0.3 · 10 5 cycles the phase field at the notch tip is almost zero (Fig. 4c) and σ A increases, while σ e y decreases due to the degradation. This observation illustrates that σ A rather than σ e y is a suitable driving force, as otherwise the fatigue damage would stop to grow. A fatigue crack was initiated after a cycle number of approximately 2.0·10 5 (Fig. 4d). Both stress contributions are zero at the crack tip as they are fully degraded. Note that for the subsequent evolution of the energy density contribution P a higher stress is observed as driving force as for the crack nucleation. This will certainly lead to a smaller number of cycles until the next crack growth increment occurs, compared to the interval for crack initiation. This behavior is in agreement with experimental findings (see e.g. Haibach 2006).
In order to assess the crack growth behavior, another simulation was set up and ran until a fatigue crack with a total length of 15 mm was generated. Figure 5 shows the results obtained from this simulation by means of different illustrations. The plot on the right shows the crack evolution with respect to the number of applied load cycles. For selected N-a states, contour plots, representing the phase field variable s over the specimen, are indicated among the line plot. The crack length is measured from the tip of the V-notch to the actual crack tip (maximum x-coordinate with s = 0). A nonlinear function of the crack length with respect to the number of applied load cycles is obtained from the simulation. The steep gradient at the end of the simulation indicates high crack growth rates. This behavior is qualitatively in line with experimental findings (see e.g. Dowling 2013). An illustration of the particular energy contributions of the total energy (Eq. (9)) with respect to the number of load cycles for the same simulation is shown This explains the small deviation of E Γ from zero at the beginning. However, this evaluation provides good insight into the mechanism of the proposed phase field model. The onset of crack growth is indicated in the plot by the first increase of the crack energy contribution E . It may be recognized, that the energies E ac representing the fatigue energy and E W representing strain energy are small compared to E . Accordingly for the illustration these contributions were factorized by 10. The fatigue energy E ac increases as consequence of the rising fatigue damage induced by cyclic loading. Accordingly, s decreases, which yields an increase of E and a decrease of E ac to compensate the additional crack energy and minimize the total energy. The energy contribution E ac increases subsequently within another area according to the actual crack pattern and the process starts again. The strain energy contribution E W is almost constant for the larger part of the simulation and increases with higher rates only after approximately 2.5 · 10 5 cycles due to higher stresses with increasing crack length.

Evaluation of crack growth rates
The growth behavior of a macro crack, which is exposed to cyclic loads is still mainly described by the a relation between the growth rate da/dN and the cyclic stress intensity factor K , which represents the properties of crack length, loading conditions, and geometry. This law, denoted as Paris' law, was first presented in Paris and Erdogan (1963) and is given by were C and m are material dependent parameters. In a double logarithmic diagram one obtains the power law is utilized with p denoting the p-th data point. For the approximation of the K values associated to the C(T)-specimen we stuck to Srawley (1976), who recommends: where the actual crack length a is considered as the horizontal distance from the center of the holes to the crack tip (see Fig. 2). According to this, for another assessment of the proposed phase field model several simulations with different maximum force values for the introduced mode-I test of the C(T)-specimen were performed. A fatigue crack was generated within every simulation and driven to a maximum length of approximately 6 mm. Figure 6 summarizes the results of these simulations. Within a double logarithmic plot of the crack growth rates with respect to the stress intensity factor all the data points lie close to a straight line. This straight line, which confirms Eq. (39) with C = 2.8 · 10 −7 mm/cycle (MPa √ m) m and m = 2.5 was estimated by means of a least square fit and reveals a coefficient of determination of 0.996. The associated N-a plots are presented in Fig. 6a. According to the different levels of the maximum load for each simulation, large differences are obtained for the number of cycles required to generate a certain crack length. However, all samples fit to the Paris parameters of Fig. 6a, which is indeed a very important property regarding a comparison with experimental findings.

Concluding remarks
A phase field model, which is able to cover the phenomenon of cyclic mechanical fatigue, has been presented as enhancement of the model for brittle  Kuhn and Müller (2010). The driving force for the fatigue mechanism was realized by a new energy contribution accounting for the irreversibility of mechanical fatigue. Numerical simulations show that our model generally yields reasonable results in terms of fatigue crack nucleation and growth and therefore describes the phenomenon. Simulated crack growth rates confirm the Paris law, which represents experimental findings. However, as fatigue of materials is a rather complex field, which is effected by many quantities like load ratio, load sequence, strain rate, temperature, environment and not forgetting the class of the material, further enhancement of the current phase field model has to be done in order to accomplish a comprehensive formulation.