# Towards a comprehensive regularized continuum damage mechanics model for fatigue crack growth

O. Voreux, S. Feld-Payet, P. Kanouté, S. Kruch

► To cite this version:

O. Voreux, S. Feld-Payet, P. Kanouté, S. Kruch. Towards a comprehensive regularized continuum damage mechanics model for fatigue crack growth. Engineering Fracture Mechanics, 2025, 328, pp.111494. 10.1016/j.engfracmech.2025.111494. hal-05455919

HAL Id: hal-05455919

https://hal.science/hal-05455919v1

Submitted on 13 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

ELSEVIER

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

# Towards a comprehensive regularized continuum damage mechanics model for fatigue crack growth

O. Voreux<sup>ID*</sup>, S. Feld-Payet<sup>ID</sup>, P. Kanouté, S. Kruch

DMAS, ONERA, Université Paris-Saclay, 29, avenue de la Division Leclerc, Châtillon, 92320, France

## ARTICLE INFO

**Keywords:**

Fatigue crack growth

Continuum damage model

Implicit gradient formulation

Finite element method

Continuous-discontinuous approach

Mesh adaptation

## ABSTRACT

Efficient crack propagation models are required to ensure that critical parts of an aircraft continue to perform their function even in the presence of cracks. One of the main difficulties is to model the different types of damage that can occur in the same formalism given the variety of loading paths during in-service conditions that these parts are subjected to. To meet this challenge, incremental models are of particular interest, especially those based on *Continuum Damage Mechanics* (CDM). However, in order to obtain mesh-independent results, these models require the use of regularization methods. In addition, to check the ability of the model to catch the crack growth kinetics when the model is intrinsically continuous, a damage-to-crack transition is required in order to monitor the crack size. All these issues need to be addressed before any comparison can be made between numerical results from finite element calculations and those from a conventional *Linear Elastic Fracture Mechanics* (LEFM) approach. This paper proposes a possible solution by combining a regularized time-incremental model that takes into account two different types of damage (here fatigue and ductile damage) with a continuous-discontinuous strategy to insert an easy-to-track discrete crack. The procedure is investigated on a nickel-based superalloy subjected to fatigue loading conditions at elevated temperatures. The high-temperature cyclic behavior of the material is modeled within the framework of unified elastic-viscoplasticity. The discrete crack, within structural calculations, is represented by successive crack-path tracking and remeshing steps. Finite element calculations illustrate the ability of the proposed strategy to recover features identical to those obtained with conventional LEFM methods.

## 1. Introduction

### 1.1. General context

In the aerospace industry, accurately predicting the response of critical components, such as the Nickel-based superalloy high-pressure turbine disks of aircraft engines, is essential to ensure the safety of the passengers in service while controlling both development and maintenance costs. To this end, numerical models are established to predict the response to a number of specific loading spectra that are representative of in-service conditions. These include regular operational loading cases, which may be represented by cyclic loadings, as well as punctual severe loading conditions corresponding to aborted landings or so-called go-around procedures. In these cases, the material may lose its load-carrying capacity due to progressive microstructural internal damage development. This may eventually lead to crack initiation and propagation. Once a crack has initiated, it is necessary to ensure that

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: olivier.voreux@onera.fr (O. Voreux).</footnote>
Nomenclature

<table><tr><td>ε̇<sub>0</sub>, K, n</td><td>Over-stress parameters</td></tr><tr><td>ℓ̇<sub>c</sub></td><td>Internal length scale</td></tr><tr><td>η</td><td>Damage weighting parameter</td></tr><tr><td>λ</td><td>Viscoplastic multiplier</td></tr><tr><td>Ω<sup>*</sup></td><td>Thermodynamic dissipation potential</td></tr><tr><td>ζ</td><td>Non-local variable</td></tr><tr><td>ψ</td><td>Thermodynamic state potential</td></tr><tr><td>ρ, ρ<sub>0</sub></td><td>Current, initial density</td></tr><tr><td>α<sup>k</sup></td><td>Conjugated force to X<sup>k</sup></td></tr><tr><td>σ</td><td>Cauchy stress tensor</td></tr><tr><td>ε, ε<sup>e</sup>, ε<sup>in</sup></td><td>Total, elastic or inelastic strain tensors</td></tr><tr><td>n</td><td>Normal to the yield surface</td></tr><tr><td>X<sup>k</sup></td><td>k-th component of the back-stress tensor</td></tr><tr><td>ζ</td><td>Effective variable (strongly coupled to damage)</td></tr><tr><td>E</td><td>4-th order isotropic elastic tensor</td></tr><tr><td>H</td><td>4-th order elliptic tensor for equivalent stress</td></tr><tr><td>φ<sub>c</sub></td><td>Pressure sensitivity parameter</td></tr><tr><td>C<sup>k</sup>, γ<sup>k</sup></td><td>Kinematic hardening constants</td></tr><tr><td>D<sub>f</sub>, D<sub>v</sub>, D<sub>t</sub></td><td>Fatigue, ductile or total damage</td></tr><tr><td>D<sub>v<sub>0</sub></sub></td><td>Initial porosity</td></tr><tr><td>f(D)</td><td>Damage-related degradation function</td></tr><tr><td>p</td><td>Accumulated viscoplasticity</td></tr><tr><td>Q<sup>i</sup></td><td>Isotropic hardening constants</td></tr><tr><td>R<sup>i</sup></td><td>i-th component of the drag stress</td></tr><tr><td>r<sup>i</sup></td><td>Conjugated force to R<sup>i</sup></td></tr><tr><td>S<sub>f</sub>, m<sub>f</sub>, Y<sub>f<sub>0</sub></sub>, β<sub>f</sub></td><td>Fatigue damage parameters</td></tr><tr><td>CDM</td><td>Continuum Damage Mechanics</td></tr><tr><td>FCG</td><td>Fatigue crack growth</td></tr><tr><td>FE</td><td>Finite element</td></tr><tr><td>FPZ</td><td>Fracture process zone</td></tr><tr><td>HCF</td><td>High-cycle fatigue</td></tr><tr><td>LCF</td><td>Low-cycle fatigue</td></tr><tr><td>LEFM</td><td>Linear Elastic Fracture Mechanics</td></tr><tr><td>SEN-T</td><td>Single edge notched-tension</td></tr><tr><td>SIF</td><td>Stress intensity factor</td></tr></table>

the part will continue to perform its structural functions between two inspections, despite possible crack propagation. This is the principle of fail-safe design with the so-called *damage tolerant approach* [1].

A crack can notably propagate due to alternating cyclic stresses. It is then natural to describe the fatigue damage process by considering the number of load cycles *N*, which becomes the temporal measure for fatigue. This is indeed the approach taken by Paris & Erdogan [2] with the so-called Paris fatigue crack growth law: the corresponding equation expresses the derivative of the crack length with respect to the load cycles *N* as a function of the stress intensity factor (SIF) variation to characterize crack propagation under cyclic loading. Such a modeling approach remains to date the most popular *Fracture Mechanics*-based method for fatigue lifespan analysis. It is also suitable for large cracks relative to the microstructure, thus overriding the microstructural characteristics and possible crack growth dependency associated with short cracks.

In order to have a better description of the whole degradation process that may occur during in-service conditions, it would be interesting to consider the combined effects of cyclic loading with possibly complex or random loading conditions (e.g. for dwell-fatigue). This is easier if the different damage mechanisms are described by the same formalism. In this sense, punctual or random loading is more appropriately described by time-incremental formalism. Moreover, thermo-mechanical fatigue problems (strong coupling) could also be easily addressed using a time-incremental formalism. This motivates the development of time-incremental models with different damage mechanisms for fatigue crack growth (FCG) and the verification that results equivalent to those given by the Paris law are obtained when periodic cyclic conditions are involved.
## 1.2. Time-incremental models for fatigue crack growth

There are three main types of time-incremental models for FCG. The first relies on the geometrical characteristics of the plastic blunting at the crack front to predict the crack growth rate [3–5]. These models give satisfactory results in agreement with experimental data for a large number of cases, but seem to reach their limits for complex three-dimensional problems. The second type of time-incremental models considers cohesive zone models with their traction-separation law [6,7]. Although these models seem very promising, their ability to reproduce complex crack paths with mixed mode failure remains an open problem. Finally, the third type of time-incremental model is the CDM-based model. Unlike the previously mentioned models, this type of model has the advantage of being able to describe damage growth in the bulk material. This could open the door to complex 3D mixed mode crack propagation. However, it has been mainly considered for fatigue lifetime estimation up to crack initiation and quite rarely for crack propagation. This is mainly due to the need to have solutions for additional issues related to the numerical process of fracture modeling.

Indeed, when considering crack propagation, it is more useful to consider a coupled damage model (i.e. damage affects the constitutive behavior) to account for damage-induced stress redistribution than an uncoupled one. The problem is that coupled damage models exhibit a pathological dependence on the spatial discretization, with strain localization in the smallest band that can be captured by the mesh [8] and a singularity of the damage rate at the crack front [9]. This well-known mesh size and orientation sensitivity is related to the strain-softening behavior of the models and arises mathematically from the loss of uniqueness of the solution to the initial boundary value problem. To solve this issue, several regularization methods have been proposed, with notably the integral approach [10], the implicit gradient formulation [11,12], the enhanced energy [13] or the phase-field approach [14] (to name a few). Only some of them have been applied to fatigue models.

Examples of regularized CDM-based time-incremental models for FCG include the model proposed by Maniar et al. [15]. The authors choose to apply a nonlocal convolution integral formulation to a unified viscoplastic material model to solve the localization issue associated with the local formulation. The strategy is evaluated in finite element (FE) calculations performed on standardized specimens under low-cycle fatigue (LCF) conditions. Yin et al. [16] regularized a cyclic viscoplastic damage model for thermomechanical fatigue with an implicit gradient function. The damage evolution variable is controlled through a non-local internal variable related to the viscoplastic activity. This model has also been used to simulate low-cycle thermomechanical fatigue. Finally, Khalil et al. [17] proposed a dual phase-field model to deal with fatigue crack growth in metals. The authors used an elastic and plastic strain-based fatigue damage evolution.

In all cases, these regularization methods result in a solution that converges with the mesh size at the cost of finer discretization in the strain localization band. This mesh refinement can be a daunting additional cost if the entire band is meshed at the start of the calculation (assuming the crack path is known). To counterbalance this aspect, a mesh refinement strategy, as proposed in [18] and used in [19], can reduce the computational cost by automatically refining only the region where non-linear phenomena occur.

## 1.3. From a continuous damage model to a crack growth rate measure

When studying crack growth, it is important to note that regularizing the continuum damage problem is only the first step. In fact, continuous damage models such as those mentioned above can only provide qualitative results regarding the position of the crack and its evolution. In order to quantitatively study the evolution of the crack growth rate and to be able to relate it, for instance, to the Paris crack growth expression, it is necessary to be able to determine the position of the crack front.

**Principle.** To determine the position of the crack front, the idea is to first define where the crack is located: so-called *crack path tracking algorithms* are generally used to determine a hypersurface (i.e. a line in 2D or a surface in 3D) that serves as a crack support (see for example [20,21] for a review). Then, an additional criterion is used to estimate the position of the crack front along this crack support. The most commonly used criterion is based on a critical damage value. Once the crack front is defined, it is possible to either stay with a continuous model to avoid increasing the computational cost or to introduce a true discontinuity to better capture the crack kinematics.

**With a continuous damage model only.** Several authors choose to estimate the crack front position with a regularized continuous damage model without any discrete crack insertion. For example, Mesgarnejad et al. [22] have proposed a phase-field approach dedicated to cyclic fracture and were able to recover the Paris law under simple mode I loading. Similarly, Xue et al. [23] have recently proposed a phase-field damage model for creep-fatigue prediction using Fast Fourier Transform (instead of the Finite Element Method, which is the framework considered in the present paper). The work of Xue et al. [23] also recovers the Paris law for pure fatigue in ductile materials. Creep-fatigue problems are also addressed, although not through a damage variable. Finally, it should be noted that the measurement of a crack front position from a continuous damage model has also been performed using other regularization methods than phase field: Baruah and Singh [24] used an implicit gradient regularization for an elastic-plastic damage model dedicated to fatigue crack growth modeling under LCF conditions. In the latter case, a classical $da/dN - \Delta K$ representation is proposed to link the model provision to experimental results.

All these works seem to use a critical damage value to determine the crack front position. Unfortunately, none of them clearly specify how the crack support is determined.
With a continuous–discontinuous approach. When a crack is inserted, the closer the crack representation is to physics observations, the more obvious it is to measure the position of the front. In fact, according to [21], it is possible to distinguish three types of strategies for representing a crack:

* Removal of an uncontrolled damaged area that no longer transmits energy (e.g., using the element deletion technique);

* Constraining degradation in a controlled damaged area around the crack to dissipate energy while modeling gradual crack opening (e.g., with controlled element deletion technique [25] or embedded elements [26]);

* limiting crack opening at an interface, which may be combined with a cohesive zone model to dissipate energy (e.g., mesh adaptation or X-FEM [27]).

However, only a few of these have been implemented for fatigue problems with CDM-based models: mainly mesh refinement and element deletion techniques.

Mesh adaptation has been used by Khoei et al. [28] for crack propagation in ductile material under cyclic loading based on a coupled damage model regularized with a gradient-enhanced formulation. Cyclic calculations have been performed for structural problems, but mainly for very low-cycle fatigue (up to 10 cycles) and not for low-cycle fatigue (LCF) conditions. Due to the considered framework with a limited number of loading cycles, a link to the Paris law could not be made. Mesh adaptation provides the most accurate physical description of the crack geometry. However, the significant change in discretization requires field transfer and possibly some additional operations to recover equilibrium (just as X-FEM would). This may explain why previous works have preferred the element deletion technique to represent the discrete crack.

The principle of the element deletion technique is to remove an element from the problem when its material properties are degraded enough (i.e. usually when a critical value is reached). This technique has been used by Peerlings et al. [9] with a gradient implicit fatigue damage model. This technique has also been used with a standard fatigue damage model (i.e. without any regularization) by Hamon et al. [29]. In the latter case, to reduce the problem of mesh dependence, the mesh size was explicitly accounted for through the value of the critical energy. However, there is still a dependence on the mesh orientation. Nevertheless, the authors were able to obtain a fairly accurate estimate of the FCG curves for two aircraft alloys on relatively simple 2D problems. It should be noted that, in general, there is no guarantee with the element deletion technique that the resulting crack path will be continuous. Also, the width of the crack is not necessarily limited to one element (e.g., see [9]) unless a crack path tracking algorithm has been used beforehand and the element deletion has been limited to the found crack path.

## 1.4. Objective

The objective of this paper is to propose a regularized CDM-based time-incremental model for fatigue crack growth and a continuous–discontinuous strategy that allows to:

1. recover the main features of linear elastic fracture mechanics fatigue crack growth;

2. easily account for different damage mechanisms associated with potentially complex loading schemes.

This work is therefore in line with that proposed by Hamon et al. [29], which to our knowledge is the only one with two damage variables: one for static damage mechanisms (ductile fracture modes) at high $\Delta K$ values, and the other for cyclic damage processes (quasi-cleavage) from low to high $\Delta K$. In this work, however, we propose an implicit gradient-regularized formulation to overcome the mesh dependence issue associated with stress softening. Furthermore, we use a mesh adaptation strategy based on discretization error estimation to limit computational cost while ensuring a sufficiently fine discretization of the strain localization band.

To provide a realistic description of the crack and better capture its kinematics, a continuous–discontinuous strategy is used to build crack increments based on the damage distribution. For crack path tracking, we have chosen the *Marching Ridges* algorithm [19], which ensures the continuity of the crack path. For crack representation, mesh adaptation is also used, similar to Khoei et al. [28]. However, in this work we consider more cycles to be able to recover the Paris crack propagation law within the standard *Linear Elastic Fracture Mechanics* (LEFM) framework for long fatigue crack.

Note that the modeling framework is generic enough to consider different metallic materials (including nickel-based superalloys) and either LCF or HCF<fnref n="1" /> conditions. In fact, we propose an elastic–viscoplastic damage model based on the unified viscoplasticity framework [30] and CDM concepts [31] while all the considered numerical methods are independent of the constitutive model.

## 1.5. Layout

This paper is organized as follows. Section 2 describes the proposed multi-mechanism elastic–viscoplastic damage model. First, the thermodynamic formulation of the local damage model for cyclic loading is established within the framework of generalized standard materials. Then, the possibility of considering other types of damage mechanisms is illustrated with the addition of quasi-static ductile damage. Finally, the model is regularized by an implicit gradient approach to solve the spurious mesh dependence effect associated with damage-induced stress softening.

<footnote n="1">¹ Note that although the model is able to capture crack initiation, only the propagation of a (long) fatigue crack is considered in this paper.</footnote>
The continuous-to-discontinuous strategy is presented in Section 3. The first component is the optimization of the mesh size during the computation to properly capture the damage evolution while controlling the computational cost. The second component is the insertion of crack increments to more accurately represent and track crack propagation. In Section 4, the implementation of the non-local elastic-viscoplastic damage model is validated by verifying that a near-zero non-local length gives results similar to those of a local model. Then, the mesh insensitivity is also validated for monotonic loading up to crack initiation. Finally, the proposed model combined with the continuous-discontinuous strategy is applied in FE fatigue crack propagation applications to recover a Paris law. The sensitivity of the model to mesh size, applied boundary conditions, load ratio and load frequency is analyzed. Note that for all these tests, material parameters corresponding to a Nickel-based superalloy are considered and the complete multi-field problem is solved using a monolithic scheme and implicit backward Euler time integration with the Z-set software. More details about the implementation are provided in the appendices.

## 2. Time-incremental formulation of a unified elastic-viscoplastic non-local damage model

### 2.1. Thermodynamic formulation of the viscoplasticity-damage model for cyclic loading

In this section, the authors propose an isotropic multi-mechanism elastic-viscoplastic model featuring a strong coupling with damage for cyclic loading. Although this model has been initially developed for Nickel-based superalloys [32], it is general enough to be applied to other types of metallic alloys.

#### 2.1.1. Strain partition

Assuming here a small strain framework and no thermal effects for the sake of simplicity, the total strain tensor $\epsilon$ can be decomposed in a classical way into a recoverable elastic strain and an irreversible inelastic strain, denoted respectively $\epsilon^e$ and $\epsilon^{in}$:

$$ \epsilon = \frac{1}{2} (\nabla \bar{u} + (\nabla \bar{u})^T) = \epsilon^e + \epsilon^{in} \quad (1) $$

where $\bar{u}$ is the displacement vector.

#### 2.1.2. State couplings

For the present approach, the effect of fatigue damage, denoted $D_f$, is taken into account through the introduction of a scalar softening function $f(D_f)$. This damage function should be at least twice differentiable, should take the value 1 for the undamaged state ($D_f = 0$) and tend to 0 when damage is close to its critical value $D_c$. In this study, an exponential form for the degradation function has been chosen, following the work of [33–35]:

$$ f(D_f) = \exp(-\eta D_f) \quad (2) $$

where $\eta$ is a material parameter enabling to adjust the weighting influence of damage on the degradation function. If $\eta \neq 0$, there is a full coupling between damage and material behavior, whereas the case with $\eta = 0$ yields a non-coupling case, i.e. no behavior-damage coupling and no damage-induced stress-softening.

Let us note that it could have been possible to define three different degradation functions ($f_e$, $f_k$, $f_l$) to adapt the effect of damage on the elastic and both kinematic and isotropic hardenings [32]. However, in the present case, since there are no specific data enabling to differentiate the effects of damage, we simply consider a unique function $f(D_f)$, viz.:

$$ f_e(D_f) = f_k(D_f) = f_l(D_f) = f(D_f) \quad (3) $$

#### 2.1.3. Relation between primal and dual variables

**Free energy definition.** The proposed model is established within the Generalized Standard Materials framework hence assuming that the material behavior can be defined with both state $\psi$ and dissipation $\Omega^*$ potentials [31,36,37]. The model enables to take into account viscoplastic behavior with both kinematic and isotropic hardening mechanisms. Following the work of [31] the specific Helmholtz free energy $\psi$ (or the Helmholtz free energy $\bar{\psi}$) can be split into elastic and inelastic contributions:

$$ \psi(\epsilon^e, \alpha^k, r^i, D_f) = \frac{1}{\rho_0} \bar{\psi}(\epsilon^e, D_f) + \sum_k \bar{\psi}_{in}^k(\alpha^k, D_f) + \sum_l \bar{\psi}_{in}^l(r^i, D_f) \quad (4) $$

where $\rho_0$ is the initial density, $\alpha^k$ and $r^i$ are the strain-like variables associated with the different kinematic or isotropic hardenings, respectively. In Eq. (4), the different contributions to $\bar{\psi}$ are defined by:

$$ \bar{\psi}_e(\epsilon^e) = \frac{1}{2}(\epsilon - \epsilon^{in}) : f(D_f) \mathbb{E} : (\epsilon - \epsilon^{in}) = \frac{1}{2}\epsilon^e : f(D_f) \mathbb{E} : \epsilon^e \quad (5a) $$

$$ \bar{\psi}_{in}^k(\alpha^k) = \frac{1}{2}\alpha^k : f(D_f) \mathbb{C}^k : \alpha^k \quad (5b) $$

$$ \bar{\psi}_{in}^l(r^i) = \frac{1}{2} b^l f(D_f) Q^l r^i \quad (5c) $$

In Eq. (5a), $\mathbb{E}$ designates the elastic stiffness defined by:

$$ \mathbb{E} = \lambda_e \mathbf{1} \otimes \mathbf{1} + 2\mu_e \mathbb{I} = \frac{\nu E}{(1 + \nu)(1 - 2\nu)} \mathbf{1} \otimes \mathbf{1} + \frac{E}{2(1 + \nu)} \mathbb{I} \quad (6) $$
where 1 and I refer to the 2nd- and 4th-order unity tensors, respectively.

Notations $\mathbb{C}^k$ in Eq. (5b) and $Q^i$ in Eq. (5c) designate the moduli of each kinematic and isotropic hardening mechanism, respectively while $\gamma^k$ and $b^i$ are material parameters for kinematic and isotropic hardening rules, respectively. We assume a pure deviatoric kinematic strain-hardening in Eq. (11) with:

$$ \mathbb{C}^k = \frac{2}{3} \mathbb{C}^k \mathbb{P}^{dev} \qquad (7) $$

with $\mathbb{C}^k$ the modulus for the k-component of the back-stress and $\mathbb{P}^{dev}$ the 4th-order deviatoric projection tensor defined thanks to the 4th-order volumetric projection tensor $\mathbb{P}^{vol}$ such as:

$$ \mathbb{P}^{vol} = \frac{1}{3} (\mathbf{1} \otimes \mathbf{1}) \qquad (8) $$

$$ \mathbb{P}^{dev} = \mathbb{I} - \frac{1}{3} (\mathbf{1} \otimes \mathbf{1}) = \mathbb{I} - \mathbb{P}^{vol} \qquad (9) $$

These projectors are used to decompose any second-order tensor into its volumetric and deviatoric parts, respectively.

The coupling of the elastic-viscoplastic behavior with damage is realized through the diminution of the apparent moduli $\mathbb{E}$, $\mathbb{C}^k$ and $Q^i$ in Eq. (5a) to Eq. (5c) thanks to the softening function $f(D_f)$.

**Derivatives of the free energy.** From the expression of the thermodynamic state potential, $\psi$, by applying the Coleman-Noll relations [38], we get the state laws relating strain-like variables ($\epsilon^e$, $\alpha^k$, $r^i$, $D_f$) to their stress-like counterparts ($\sigma$, $X^k$, $R^i$, $Y_f$ respectively):

$$ \sigma = \frac{\partial \bar{\psi}}{\partial \epsilon^e} = f(D_f) \mathbb{E} : (\epsilon - \epsilon^{in}) = f(D_f) \mathbb{E} : \epsilon^e \qquad (10) $$

$$ X^k = \frac{\partial \bar{\psi}}{\partial \alpha^k} = f(D_f) \mathbb{C}^k : \alpha^k \qquad (11) $$

$$ R^i = \frac{\partial \bar{\psi}}{\partial r^i} = b^i f(D_f) Q^i r^i \qquad (12) $$

$$ Y_f = -\frac{\partial \bar{\psi}}{\partial D_f} = -\frac{1}{2} f'(D_f) \left( \epsilon^e : \mathbb{E} : \epsilon^e + \sum_k \alpha^k : \mathbb{C}^k : \alpha^k + \sum_i b^i Q^i r^i \right) \qquad (13) $$

The notation $\sigma$ denotes the Cauchy stress tensor, $X^k$ the k-th contribution to the back stress (kinematic hardening), $R^i$ the i-th contribution to the drag stress (isotropic hardening) and $Y_f$ the thermodynamic force associated with fatigue damage corresponding to the so-called energy release rate. This last term involves the derivative of the degradation function $f(D_f)$ which, according to Eq. (2), can be written as:

$$ f'(D_f) = -\eta \exp(-\eta D_f) \quad \leftrightarrow \quad f'(D_f) = -\eta f(D_f) \qquad (14) $$

Based on Eqs. (13) and (14), the energy release rate can simply be written as:

$$ Y_f = -\eta \left( -\frac{1}{2} \epsilon^e : f(D_f) \mathbb{E} : \epsilon^e - \frac{1}{2} \sum_k \alpha^k : f(D_f) \mathbb{C}^k : \alpha^k - \frac{1}{2} \sum_i b^i f(D_f) Q^i r^i \right) = \eta \bar{\psi} \qquad (15) $$

### 2.1.4. Pseudo-dissipation potential

The evolution laws for the internal variables are derived from a convex dissipative pseudo-potential $\Omega^*$. This pseudo-potential is the sum, on the one hand, of a viscoplastic contribution denoted $\Omega_{vp}^*$ and, on the other hand, of a fatigue damage contribution denoted $\Omega_{D_f}^*$:

$$ \Omega^* = \Omega_{vp}^* (\tilde{\sigma}, \tilde{X}^k, \tilde{R}^i; \alpha^k, r^i) + \Omega_{D_f}^* (\epsilon^e, \alpha^k, r^i, D_f) \qquad (16) $$

The stress-like variables are written here with the notation $\cdot$ to indicate an effective quantity (i.e. coupled to damage) according to the principle of strain equivalence [39] and are thus defined as follows:

$$ \tilde{\sigma}(D_f) = \frac{1}{f(D_f)} \sigma \qquad (17) $$

$$ \tilde{X}^k(D_f) = \frac{1}{f(D_f)} X^k \qquad (18) $$

$$ \tilde{R}^i(D_f) = \frac{R^i}{f(D_f)} \qquad (19) $$

Firstly, for the viscoplastic contribution to the pseudo-dissipation potential $\Omega_{vp}^*$, we consider a simplified form of the hyperbolic sine flow rule first introduced in the work of [40] with specific particularization given in [41,42]:

$$ \Omega_{vp}^*(\tilde{\sigma}, \tilde{X}^k, \tilde{R}^i; \alpha^k, r^i) = \dot{\epsilon}_0 \frac{K}{n+1} \cosh \left\langle \frac{F_p(\tilde{\sigma}, \tilde{X}^k, \tilde{R}^i; \alpha^k, r^i)}{K} \right\rangle_+^{n+1} \qquad (20) $$
with

$$ \frac{\partial \bar{\sigma}_{eq}}{\partial \bar{\mathbf{X}}^k} = -\frac{1}{2} \times \frac{3}{2} \times 2 \times \frac{\left( \bar{\sigma}^{dev} - \sum_k \bar{\mathbf{X}}^k \right)}{\bar{\sigma}_{eq}} = -\bar{n} \quad (31) $$

Using Eqs. (11) and (18) to express differently the second term between brackets in Eq. (30), we finally get:

$$ \dot{\alpha}^k = \frac{\dot{\lambda}}{f(D_f)} (\bar{n} - \gamma^k \alpha^k) \quad (32) $$

Derivation with regards to $R^i$. Similarly, the derivation of $\Omega_{vp}^*$ with regards to $R^i$ provides the evolution law for $r_i$:

$$ \dot{r}^i = -\frac{\partial \Omega_{vp}^*}{\partial R^i} = -\frac{\partial \Omega_{vp}^*}{\partial F_p} \frac{\partial F_p}{\partial \bar{R}^i} \frac{\partial \bar{R}^i}{\partial R^i} = -\dot{\lambda} \left( \frac{\partial F_y}{\partial \bar{R}^i} + \frac{1}{Q^i} \bar{R}^i \right) \frac{\partial \bar{R}^i}{\partial R^i} \quad (33) $$

Using Eq. (22) together with Eq. (12) to express differently the term between brackets and Eq. (19) to determine $\frac{\partial \bar{R}^i}{\partial R^i}$ leads to:

$$ \dot{r}^i = \frac{\dot{\lambda}}{f(D_f)} (1 - b^i r^i) \quad (34) $$

Derivation with regards to $Y_f$. Finally, the evolution of fatigue damage is given by the derivative of the pseudo-potential related to fatigue damage $\Omega_{D_f}^*$:

$$ \dot{D}_f = -\frac{\partial \Omega^*}{\partial(-Y_f)} = \frac{\partial \Omega^*}{\partial Y_f} = \frac{\partial \Omega_{D_f}^*}{\partial Y_f} \quad (35) $$

According to Eq. (25), when $(Y_f - Y_{f0}) > 0$:

$$ \dot{D}_f = \frac{\partial \Omega_{D_f}^*}{\partial Y_f} = \dot{p} f(D_f) \left\langle \frac{Y_f - Y_{f0}}{S_f} \right\rangle_+^{m_f} \exp(\beta_f D_f) \quad (36) $$

It is possible to express the fatigue damage rate $\dot{D}_f$ as a function of the equivalent viscoplastic multiplier $\dot{\lambda}$ by relating $\dot{\lambda}$ to the accumulated viscoplastic strain rate $\dot{p}$ thanks to the norm of the inelastic strain rate tensor, viz:

$$ \dot{p} = \sqrt{\frac{2}{3} \dot{\epsilon}^{in} : \dot{\epsilon}^{in}} \quad (37) $$

Using Eq. (28), it is then possible to express $\dot{\epsilon}^{in}$ as a function of $\dot{\lambda}$:

$$ \dot{p} = \frac{\dot{\lambda}}{f(D_f)} \sqrt{\frac{2}{3} \bar{n} : \bar{n}} = \frac{\dot{\lambda}}{f(D_f)} \quad (38) $$

since:

$$ \sqrt{\frac{2}{3} \bar{n} : \bar{n}} = \sqrt{\frac{2}{3} \times \frac{3}{2} \frac{(\bar{\sigma}^{dev} - \sum_k \bar{\mathbf{X}}^k)}{\bar{\sigma}_{eq}} : \frac{3}{2} \frac{(\bar{\sigma}^{dev} - \sum_k \bar{\mathbf{X}}^k)}{\bar{\sigma}_{eq}}} = 1 \quad (39) $$

Consequently:

$$ \dot{\lambda} = \dot{p} f(D_f) \quad (40) $$

Replacing $\dot{p} f(D_f)$ by $\dot{\lambda}$ in Eq. (36), we finally get the following time evolution of fatigue damage:

$$ \dot{D}_f = \dot{\lambda} \left\langle \frac{Y_f - Y_{f0}}{S_f} \right\rangle_+^{m_f} \exp(\beta_f D_f) \quad (41) $$

## 2.2. Extension to ductile damage modeling

One of the main advantages of a time-incremental formulation for the fatigue model is the possibility to easily consider other degradation mechanisms. As an example, we propose in this section to consider ductile damage.

### 2.2.1. Ductile damage definition

During plastic flow in metallic alloys, so-called ductile damage associated with the growth of porosity in the bulk may occur for relatively high plastic strain levels [44]. For fatigue crack growth problems in metallic alloys, it is well admitted that plasticity occurs at the crack-tip and governs the crack advance [1,3]. To take into account this additional damage phenomenon due to void growth into the model, we propose to introduce a new damage variable $D_v$ that can be defined based on the variation of the density $\rho$:

$$ \rho = (1 - D_v) \rho_0 \quad \leftrightarrow \quad D_v = 1 - \frac{\rho}{\rho_0} \quad (42) $$
### 2.2.2. Evolution equation

Derivation of Eq. (42) leads to:

$$ \dot{D}_v = - \frac{\dot{\rho}}{\rho_0} \qquad (43) $$

Neglecting the elastic strain-induced volume variation, the mass conservation principle enables to obtain the following expression for the density rate $\dot{\rho}$:

$$ \dot{\rho} = -\rho \operatorname{tr} \dot{\epsilon}^{in} \qquad (44) $$

Replacing $\dot{\rho}$ in Eq. (43) and expressing the ratio $\frac{\rho}{\rho_0}$ as a function of $D_v$ thanks to Eq. (42), we get:

$$ \dot{D}_v = \frac{\rho \operatorname{tr} \dot{\epsilon}^{in}}{\rho_0} = (1 - D_v) \operatorname{tr} \dot{\epsilon}^{in} \qquad (45) $$

Note that since the physical basis of the development is the same, the evolution equation given in Eq. (45) is identical to the cavity growth equation in the Gurson [45] or Rousselier [46] models.

### 2.2.3. Ductile damage coupling

Ductile damage directly impacts the viscoplastic flow through the introduction of a 4th-order tensor $\mathbb{H}$ in the yield function initially defined in Eq. (22), following [44,47,48], so that:

$$ F_y = \sigma_v = \sqrt{\frac{3}{2} \left( \tilde{\sigma} - \sum_k \tilde{X}^k \right) : \mathbb{H} : \left( \tilde{\sigma} - \sum_k \tilde{X}^k \right) - \sum_i \tilde{R}^i - R_0} \qquad (46) $$

We choose an elliptic form for the tensor $\mathbb{H}$, as already used by [47] for porous materials:

$$ \mathbb{H} = [\mathbb{P}^{dev} + 2 \varphi_c g_D (D_v) \mathbb{P}^{vol}] \qquad (47) $$

In this equation, $\varphi_c$ is a material parameter to modulate the pressure sensitivity: e.g. with $\varphi_c = 0$, the elliptic yield criterion Eq. (46) reduces to the standard pressure-independent von Mises criterion given in Eq. (23). Notation $g_D$ designates a function of the ductile damage $D_v$ that enables to take into account the effect of an initial porosity $D_{v0}$ and makes the criterion pressure-dependent [49]. Here, it is simply chosen to be equal to the following sum:

$$ g_D (D_v) = D_v + D_{v0} \qquad (48) $$

The introduction of the $\mathbb{H}$ tensor in the yield criterion results in a change of the effective normal to the yield surface $\tilde{n}$ previously defined in Eq. (27) which now reads:

$$ \tilde{n} = \frac{3}{2} \frac{\mathbb{H} : \left( \tilde{\sigma}^{dev} - \sum_k \tilde{X}^k \right)}{\tilde{\sigma}_{eq}} \qquad (49) $$

The ductile damage coupling also impacts the cumulative viscoplastic strain rate which now reads:

$$ \dot{p} = \sqrt{\frac{2}{3} \dot{\epsilon}^{in} : \mathbb{H}^{-1} : \dot{\epsilon}^{in}} \qquad (50) $$

### 2.2.4. Damage interactions

The various damage processes can be summed up to give a linear cumulative damage $D_t$, as proposed by [50], viz.:

$$ D_t = D_f + D_v \quad \text{and} \quad \dot{D}_t = \dot{D}_f + \dot{D}_v \qquad (51) $$

This is the simplest type of interaction model, formulated here under the assumption of an equally weighted linear combination of all damage mechanisms [51]. This type of model is usually associated with the assumption that failure occurs when the sum of fatigue and ductile damage equals a critical value. This summation results in the direct coupling between damage and both elasticity and (visco)plasticity through the softening function given in Eq. (2). Indeed, the softening function becomes a function of the total damage $D_t$ that impacts all the effective variables through Eq. (17) to Eq. (19). Consequently, the relation between the cumulative viscoplastic strain rate and the viscoplastic multiplier is also impacted and now reads:

$$ \dot{p} = \frac{\dot{\lambda}}{f(D_t)} \qquad (52) $$

## 2.3. Implicit gradient enhancement

Since damage is fully coupled to the material behavior, material softening is expected to occur. The mathematical problem then becomes ill-posed since the governing partial differential equations lose their ellipticity [52]. As a consequence, the solution becomes mesh dependent. To overcome this well-known problem, we consider a non-local implicit gradient approach [11].
### 2.3.1. Additional equations

The implicit gradient enhancement is applied to the total damage variable $D_t$ which includes all possible damage mechanisms (i.e. fatigue and ductile damage, see Eq. (51)). Note that the chosen variable is not bounded thanks to the exponential form used in the degradation function $f(D_t)$ given in Eq. (2), which is an essential requirement [12]. Moreover, from practical considerations, the choice for $D_t$ as the non-local variable rather than $\epsilon_{eq}^{in}$ or $p$, which is also very common [53–55], stems mainly from the need to dissociate the cyclic material response from damage-related mechanisms. In this way, the elastic–viscoplastic constitutive equations of the material can be calibrated according to the hysteresis loops and damage-related parameters can then be calibrated according to crack growth data. Such an aspect is discussed in detail in [32].

A non-local total damage variable $\overline{D}_t$ is thus defined from the local total damage $D_t$:

$$ \overline{D}_t - c \nabla^2 \overline{D}_t = D_t \quad \text{in } \Omega \qquad (53) $$

In this equation, $\Omega$ is the domain occupied by the body of interest, $\nabla^2$ is the Laplacian operator and $c = \frac{1}{2}\ell_c^2$ is a non-local parameter that is homogeneous to the square of a length (hence the notion of so-called internal characteristic length) [9]. This non-local parameter $c$ indirectly controls the width of the localization band [55]. It is worth noting that if $c = 0$ mm² the non-local model recovers its local form and there is no effect of the non-local extension on the damage model.

Besides, the non-local variable should verify the following homogeneous Neumann-type boundary condition:

$$ \overline{\nabla} \overline{D}_t \cdot \vec{n}_d = 0 \quad \text{on } \partial\Omega \qquad (54) $$

where $\overline{\nabla}$ is the gradient operator and $\vec{n}_d$ is the unit normal to the boundary $\partial\Omega$.

In order to ensure a positive dissipation along the whole loading path, a decreasing evolution for $\overline{D}_t$ must be precluded. To this end, healing effects should be excluded, which can be, in an analytical and incremental manner, be achieved by setting:

$$ d\overline{D}_t \ge 0 \quad \leftrightarrow \quad \overline{D}_t(t) = \max[\overline{D}_t(\tau)], \quad 0 \le \tau \le t \qquad (55) $$

### 2.3.2. Constitutive equations modifications

To benefit from the regularization effect of the implicit gradient enhancement, we replace the local total damage by its non local counterpart in all the equations (except for its definition as the sum of the local fatigue and ductile damages). The resulting equations for the so-called non-local elastic–viscoplastic damage model are summarized in Table 1 where they can be compared with their previous local version.

Both local and non-local damage models have been implemented in the Z-set finite element code [56] for structural calculations (see the forthcoming Section 4).

## 3. Continuous–discontinuous strategy

One of the aims of this paper is to present the application of a continuous–discontinuous strategy for fatigue crack propagation (see Fig. 1). As the strategy has already been applied within a different framework, namely monotonic loading [19,57], the main components are only briefly described here. Note that all the numerical developments have been implemented in the Z-set FE code [56].

### 3.1. Generic mesh-adaptation procedure for continuous media

In order to maximize the accuracy/cost ratio, we choose to use the automatic mesh adaptation strategy proposed in [18] and used in [19]. In addition, such a procedure prevents excessive element distortion in large deformation and damage problems through mesh refinement and FE shape optimization.

#### 3.1.1. Error estimation

The chosen strategy relies on an a posteriori error estimation (step 4a. of Fig. 1) based on a recovery method, as proposed by Zienkiewicz and Zhu [58,59]. This method is based on the construction of a more continuous field that serves as an estimation of the exact solution. A discretization error $\epsilon_f^e$ is obtained by considering the difference with the finite element solution. Following the work of [60,61], Feld-Payet [18] proposed to construct an error indicator for any gradient-like scalar field $f$ (i.e. a field stored at the integration points such as the damage field) for any element $\Omega_e$ or the whole structure $\Omega$ in the following way:

$$ \epsilon_f^e = 100 \times \sqrt{\frac{\int_{\Omega_e} |(f^* - f) \times (df^* - df)| d\Omega}{|\int_{\Omega_e} (f^* \times df^*) d\Omega|}} \quad \text{and} \quad \epsilon_f = 100 \times \sqrt{\frac{\int_{\Omega} |(f^* - f) \times (df^* - df)| d\Omega}{|\int_{\Omega} (f^* \times df^*) d\Omega|}} \qquad (56) $$

In Eq. (56), $f^*$ denotes the smoothed field obtained with the Superconvergent Patch Recovery technique [58,59] (step 3 of Fig. 1). Choice is made here that $f = \overline{D}$ and thus $f^* = \overline{D}^*$.

At the end of each converged increment, if the global error $\epsilon_{\overline{D}}$ exceeds a user-defined authorized value, mesh adaptation is triggered.
<table><caption>**Table 1**

Summary of the constitutive equations for both local and non-local elastic-viscoplastic damage models.</caption><thead><tr><th>Local damage model</th><th>Non-local damage model</th></tr></thead><tbody><tr><td>σ̄ = f(D<sub>t</sub>)<sup>-1</sup> σ</td><td>σ̄ = f(D̄<sub>t</sub>)<sup>-1</sup> σ</td></tr><tr><td>σ = f(D<sub>t</sub>) E : ε<sup>e</sup></td><td>σ = f(D̄<sub>t</sub>) E : ε<sup>e</sup></td></tr><tr><td>X̂<sup>k</sup> = f(D<sub>t</sub>)<sup>-1</sup> X̂<sup>k</sup></td><td>X̂<sup>k</sup> = f(D̄<sub>t</sub>)<sup>-1</sup> X̂<sup>k</sup></td></tr><tr><td>X<sup>k</sup> = f(D<sub>t</sub>) Ĉ<sup>k</sup> : α<sup>k</sup></td><td>X<sup>k</sup> = f(D̄<sub>t</sub>) Ĉ<sup>k</sup> : α<sup>k</sup></td></tr><tr><td>R̂<sup>i</sup> = f(D<sub>t</sub>)<sup>-1</sup> R<sup>i</sup></td><td>R̂<sup>i</sup> = f(D̄<sub>t</sub>)<sup>-1</sup> R<sup>i</sup></td></tr><tr><td>R<sup>i</sup> = b<sup>i</sup> f(D<sub>t</sub>) Q<sup>i</sup> r<sup>i</sup></td><td>R<sup>i</sup> = b<sup>i</sup> f(D̄<sub>t</sub>) Q<sup>i</sup> r<sup>i</sup></td></tr><tr><td>T<sub>y</sub> = σ̄<sub>eq</sub> - Σ<sub>k</sub> R̂<sup>i</sup> - R<sub>0</sub></td><td>σ̄<sub>eq</sub> = √(σ̄<sup>dev</sup> - Σ<sub>k</sub> X̂<sup>k</sup>) : H : (σ̄<sup>dev</sup> - Σ<sub>k</sub> X̂<sup>k</sup>)</td></tr><tr><td>h̄ = <sup>3</sup>&frasl;<sub>2</sub> (H : (δ<sup>dev</sup> - Σ<sub>k</sub> X̂<sup>k</sup>) / δ<sub>st</sub>)</td><td>h̄ = <sup>3</sup>&frasl;<sub>2</sub> (δ<sup>dev</sup> - Σ<sub>k</sub> X̂<sup>k</sup>)</td></tr><tr><td>H = [P<sup>dev</sup> + 2 φ<sub>c</sub> g<sub>D</sub> (D<sub>t</sub>) P<sup>vol</sup>]</td><td>H = [P<sup>dev</sup> + 2 φ<sub>c</sub> g<sub>D</sub> (D̄<sub>t</sub>) P<sup>vol</sup>]</td></tr><tr><td>ε<sup>in</sup> = <sup>λ</sup>&frasl;<sub>f(D<sub>t</sub>)</sub> h̄</td><td>ε<sup>in</sup> = <sup>λ</sup>&frasl;<sub>f(D<sub>t</sub>)</sub> h̄</td></tr><tr><td>α<sup>k</sup> = <sup>λ</sup>&frasl;<sub>f(D<sub>t</sub>)</sub> (h̄ - γ<sup>k</sup> α<sup>k</sup>)</td><td>α<sup>k</sup> = <sup>λ</sup>&frasl;<sub>f(D<sub>t</sub>)</sub> (h̄ - γ<sup>k</sup> α<sup>k</sup>)</td></tr><tr><td>r<sup>i</sup> = <sup>λ</sup>&frasl;<sub>f(D<sub>t</sub>)</sub> (1 - b<sup>i</sup> r<sup>i</sup>)</td><td>r<sup>i</sup> = <sup>λ</sup>&frasl;<sub>f(D<sub>t</sub>)</sub> (1 - b<sup>i</sup> r<sup>i</sup>)</td></tr><tr><td>λ = ε<sub>0</sub> sinh (∫<sub>0</sub><sup>λ/κ</sup> d<sub>+</sub>) = p̂ f(D<sub>t</sub>)</td><td>λ = ε<sub>0</sub> sinh (∫<sub>0</sub><sup>λ/κ</sup> d<sub>+</sub>) = p̂ f(D̄<sub>t</sub>)</td></tr><tr><td>D̂<sub>f</sub> = λ (∫<sub>0</sub><sup>Y<sub>f</sub>-Y<sub>0</sub>/S<sub>f</sub></sup> d<sub>+</sub>)<sup>m<sub>f</sub></sup> e<sup>(β<sub>f</sub> D<sub>t</sub>)</sup></td><td>D̂<sub>f</sub> = λ (∫<sub>0</sub><sup>Y<sub>f</sub>-Y<sub>0</sub>/S<sub>f</sub></sup> d<sub>+</sub>)<sup>m<sub>f</sub></sup> e<sup>(β<sub>f</sub> D̄<sub>t</sub>)</sup></td></tr><tr><td>D̂<sub>v</sub> = (1 - D<sub>t</sub>) <sup>λ</sup>&frasl;<sub>f(D<sub>t</sub>)</sub> tr h̄</td><td>D̂<sub>v</sub> = (1 - D̄<sub>t</sub>) <sup>λ</sup>&frasl;<sub>f(D<sub>t</sub>)</sub> tr h̄</td></tr><tr><td>D̂<sub>t</sub> = D̂<sub>f</sub> + D̂<sub>v</sub></td><td></td></tr></tbody></table>

### 3.1.2. Size map computation

The local error estimation $\epsilon_D^e$ can be used to construct the new size map to reach the expected global error $\epsilon_D^-$ on the new mesh according to the resmashing criterion proposed in [62] that aims at imposing equal distribution of the global error. This criterion enables to compute the corresponding element sizes for the new elements $h_e'$:

$$h_e' = h_e \times \frac{\epsilon_D^{1/q}}{\epsilon_D^{e/2/(2q+d)} \left[ \sum \Omega_e \frac{\epsilon_D^e}{D} \right]^{1/2q}} \quad (57)$$

where $h_e$ is the current element size, $d$ is the problem dimension and $q$ is a scalar value related to the convergence rate (it is equal to 1 if the displacement field is discretized with linear shape functions for the displacement field and 2 if the displacement field shape functions are quadratic). This operation corresponds to step 5 in Fig. 1.

Note that the user always has the possibility to control the computation cost by imposing a minimum allowable element size $h_{min}$.

### 3.1.3. Mesh generation

In this work, the new (optimized) mesh is obtained (step 6 of Fig. 1) thanks to a fully automatic mesh generator that has been developed at INRIA by the Distene work-group, namely MeshGems [63]. Since this mesh generator requires a size map defined at the nodes, a list of optimal sizes for the new mesh is computed by averaging the sizes for the new elements $h_e'$:

$$h_i' = \frac{1}{n_e} \sum_{\Omega_e \subset \Omega_i} h_e' \quad (58)$$

In this equation, $\Omega_i$ designates the union of the $n_e$ elements connected to node $i$.

### 3.1.4. Fields transfer

Once the new mesh is generated, it is necessary to transfer the mechanical fields onto it in order to continue computation (step 7 of Fig. 1). Different transfer operators are used depending on the type of variables. Variables stored at the nodes are transferred by
<figure><img src="image_5.png" /><figcaption>**Fig. 1.** Flowchart for the mesh adaptation and crack propagation strategy. For each time increment, if the finite element computation has converged, a smoothed field is computed to enable error estimation and crack path tracking. If the criterion for crack increment insertion or mesh adaptation is met, then a size map is computed in order to guide the new mesh generation. The fields are then transferred onto the new mesh and a re-equilibrium step helps retrieve equilibrium to continue computation for the next time increment.</figcaption></figure>

direct interpolation from the nodal values in the previous discretization using the element's shape functions. Variables stored at the integration points are transferred thanks to a moving least squares technique using a limited number of data, namely 13 integration points for 6-node triangles [19]. See [18] for more details on the implementation.

### 3.1.5. Retrieving equilibrium

Since the transferred fields are no longer consistent, this may cause some convergence issues for the next time step. It is then useful to enable time step division, as in [18,64]. This solution may be sufficient as long as the transfer-related perturbations remain relatively minor. Since they can be more severe in the case of crack increment insertion, we use the re-equilibrium procedure proposed by El Ouazani Tuhami [57,65] to help retrieve equilibrium more efficiently (i.e. with a reduced number of time steps). It consists in introducing an additional step with the same boundary conditions as the increment leading to mesh adaptation. For this re-equilibrium step (i.e. step 8 of Fig. 1), the problem is solved with the full constitutive behavior to enable adjustment of both plasticity and damage.

## 3.2. Crack propagation

Mesh adaptation is also used in this study to insert crack increments when necessary. To determine the crack increments, we choose to use the continuous-discontinuous strategy of [18,19]. The idea is to use the regularized damage distribution to construct the crack path.

### 3.2.1. Crack support determination

The first step is to determine the support of the crack path thanks to a crack path tracking algorithm. Among the different possibilities, we use the so-called *Marching Ridges* algorithm, described in [19] and illustrated in Fig. 2, for its versatility and
<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Principle of the search for the crack direction thanks to the marching ridges algorithm. The scalar field with a ridge representing material degradation *f* and its gradient vectors are represented in orange and the polar grid and its orthoradial vectors are represented in blue. The resulting scalar product *g* averaged by 10° sectors is in green. Its sign changes from positive to negative when the considered angle *θ* corresponds to the crack.</figcaption></figure>

accuracy. It consists in projecting the smoothed damage field $\bar{D}^*$ (see Section 3.1.1) onto the orthoradial vectors $\bar{e}_\theta$ of a polar grid (centered at the last known point on the crack) to find the next points on the ridge. These ridge points correspond to the angular direction where the considered scalar product changes sign from positive to negative. For more robustness, the results are averaged by $\Delta\theta = 10^\circ$ sectors which gives a very satisfactory resolution. The scalar products are evaluated at a maximum distance that enables to have the desired angular accuracy based on the spatial distribution of the integration points and that is at least larger than the desired crack increment length.

### 3.2.2. Crack increment construction and insertion

The crack increment length $\ell_{inc}$ is fixed and crack insertion is triggered when damage has exceeded a critical value $\bar{D}^{crit}$ along the entire crack increment. The choice of this length should be a compromise: it should not be too long to minimize the difference of mechanical state before and after crack increment insertion, but it should not be too short to avoid too frequent remeshing steps and field transfers. It is advised to choose a length slightly larger than the damage localization bandwidth to minimize the risk of spurious crack path oscillations. Such a value has been set for the whole study, see Table 4. Note that the sensitivity analysis from [57] suggested that longer crack increments could be considered. However, the final choice should be guided by a test-calculation comparison.

Once the support of the crack and the length of the crack increment have been defined, we have an explicit representation of the crack increment surface.

### 3.2.3. Cracked mesh generation

From the explicit definition of the crack surface, there are several possible ways (X-FEM, E-FEM, controlled erosion, ...) to take into account the effect of the crack in the computation when crack insertion is required. Since mesh adaptation is already used to maximize the accuracy/cost ratio, it makes sense for this work to also use mesh adaptation to represent the crack.

To propagate the crack in the mesh, we resort to the mesh intersection algorithm developed by [66]. This algorithm uses the crack surface to provide a contour mesh of the cracked structure. The *MeshGems* software [63] is then used to generate the entire mesh. During this operation, the mesh generator can still take into account any size map instructions resulting from the error estimation, as described in Section 3.1, for an optimal cracked mesh.

## 4. Applications

In this section, we discuss the results of some numerical applications to evaluate and validate our proposal. We consider different test cases with a material that may be used for critical aircraft engine components (e.g. high-pressure turbine disks in Nickel-based superalloy).

All FE calculations are performed with the implicit Z-set solver [56]. In addition, all these test cases share the same FE formulation and numerical parameters which are described at the beginning in Section 4.1 and Section 4.2, respectively. We focus first on the validation of the constitutive model in Section 4.3 and then, in Section 4.4, on the validation of its combination with the continuous-discontinuous strategy in LCF conditions to ensure reasonable computational cost. The final Section 4.6 discusses the performance of the proposed approach.
<table><caption>**Table 2**
Nominal chemical composition (in wt%) of AD730™ Nickel-based superalloy.</caption><thead><tr><th>Element</th><th>C</th><th>Ni</th><th>Cr</th><th>Mo</th><th>W</th><th>Al</th><th>Co</th><th>Ti</th><th>Nb</th><th>Fe</th><th>B</th><th>Zr</th></tr></thead><tbody><tr><th>Content</th><td>0.01</td><td>Bal.</td><td>15.46</td><td>2.82</td><td>2.59</td><td>2.22</td><td>8.27</td><td>3.48</td><td>1.11</td><td>4.16</td><td>0.0094</td><td>0.034</td></tr></tbody></table>

<table><caption>**Table 3**
Material parameters for the cyclic constitutive model of AD730™ at 550 and 700°C. For industrial confidentiality reasons, specific cyclic material parameters have been made dimensionless with respect to the yield stress $R_0$.</caption><thead><tr><th></th><th>E/R<sub>0</sub><br>(GPa)</th><th>ν<br>(-)</th><th>ε'<sub>0</sub><br>(s<sup>-1</sup>)</th><th>K/R<sub>0</sub><br>(MPa)</th><th>n<br>(-)</th><th>S<sub>f</sub><br>(MPa)</th><th>m<sub>f</sub><br>(-)</th><th>β<sub>f</sub><br>(-)</th><th>Y<sub>f</sub><br>(MPa)</th><th>φ<sub>c</sub><br>(-)</th><th>D<sub>v0</sub><br>(-)</th><th>η<br>(-)</th></tr></thead><tbody><tr><td>550°C</td><td>280.6</td><td>0.3</td><td>2.6e-8</td><td>0.01</td><td>0.48</td><td>2.0</td><td>1.0</td><td>1.0</td><td>0.0</td><td>0.2</td><td>1.e-5</td><td>4.55</td></tr><tr><td>700°C</td><td>540.6</td><td>0.3</td><td>7.8e-8</td><td>0.06</td><td>0.35</td><td>2.0</td><td>1.0</td><td>1.0</td><td>0.0</td><td>0.2</td><td>1.e-5</td><td>4.55</td></tr><tr><th></th><th>R<sub>0</sub>/R<sub>0</sub><br>(MPa)</th><th>Q<sup>1</sup>/R<sub>0</sub><br>(MPa)</th><th>b<sup>1</sup><br>(-)</th><th>Q<sup>2</sup>/R<sub>0</sub><br>(MPa)</th><th>b<sup>2</sup><br>(-)</th><th>C<sup>1</sup>/R<sub>0</sub><br>(MPa)</th><th>γ<sup>1</sup><br>(-)</th><th>C<sup>2</sup>/R<sub>0</sub><br>(MPa)</th><th>γ<sup>2</sup><br>(-)</th><th>C<sup>3</sup>/R<sub>0</sub><br>(MPa)</th><th>γ<sup>3</sup><br>(-)</th><th>ε'<sub>c</sub><br>(mm)</th></tr><tr><td>550°C</td><td>1.0</td><td>-0.60</td><td>1.5</td><td>-0.09</td><td>2</td><td>925.4</td><td>8.7e+3</td><td>238.8</td><td>1.e+3</td><td>0.8</td><td>3.e+1</td><td>0.03</td></tr><tr><td>700°C</td><td>1.0</td><td>-0.97</td><td>3.5</td><td>0.0</td><td>0.0</td><td>1937.5</td><td>6.5e+3</td><td>437.5</td><td>1.e+3</td><td>34.4</td><td>3.e+1</td><td>0.03</td></tr></tbody></table>

## 4.1. About the finite element formulation

The structures are meshed with 6-node quadratic triangular elements with reduced integration points. To overcome the well-known volumetric locking problem that occurs with isochoric plastic flow [67], we use the p2p1p1p1 mixed non-local finite element formulation proposed by [53]: we consider a quadratic interpolation (p2) for the displacement $\bar{u}$ and a linear interpolation (p1) for the non-local damage $\bar{D}_t$ and the two additional degrees of freedom, i.e. the assumed pressure $P$ and the assumed dilatation $\theta$.

The mechanical problem is solved using a monolithic scheme. Although computationally expensive, such a procedure ensures a strong coupling between damage and material behavior without inconsistencies in the solution [54]. Furthermore, the consistent tangent operator can be estimated from the implicit resolution schemes considered for both global equilibrium and local constitutive equations (see Appendix C). This choice allows unconditionally stable and relatively large resolution time steps, which is of first importance for cyclic computations.

## 4.2. Numerical parameters

**Constitutive model.** We consider the recently developed Nickel-based superalloy AD730™ [68]. Its chemical composition is given in the following Table 2. The applied heat treatments to control the fine-grained microstructure are detailed in the following reference from the authors [32]. The resulting average grain size is approximately 8–10 µm.

Corresponding parameters associated with the cyclic material response of AD730™ are presented in Table 3. This choice is based in particular on the comprehensive study of the cyclic behavior of AD730™ carried out in [32], which proposes a parameter calibration for the elastic–viscoplastic constitutive equations (i.e. without damage). Two temperatures are considered, 550 and 700 °C, which are common target temperatures for high-pressure rotating aircraft components subjected to cyclic loading. On the one hand, the first applications have been studied at 700 °C, a temperature at which material strength is lower, hence resulting in shorter fatigue lives and faster computations. On the other hand, a change to 550 °C is considered in the last Section 4.6 in order to evidence the validity of the model for the second target temperature.

Using experimental data to calibrate the fatigue and ductile crack growth parameters of a regularized CDM-based model is still a challenge to this day. Indeed, these parameters characterize the strain localization stage and standard calibration methods relying on homogeneous fields are no longer appropriate. Therefore, this study proposes to consider an arbitrary set of damage-related parameters that could serve as an initial guess for an optimization loop in a subsequent study. In particular, these damage parameters have been chosen to ensure a rapid damage evolution and relatively fast crack propagation compatible with LCF conditions. For example, the magnitude of the internal length $\ell_c$, introduced by the implicit gradient method, was chosen to be sufficiently small to avoid slowing down the evolution of damage. Indeed, the characteristic length impacts the width of the localization band and thus the dissipated energy [55]: the smaller the internal length, the smaller the dissipated energy. However, we have ensured that it is sufficiently large so that the resulting localization band width is greater than three finite elements at the onset of failure for mesh convergence issues, following guidelines from [55].

**Mesh adaptation and crack path tracking related parameters.** As specified in Section 3.1, mesh adaptation techniques are considered in this work to insert a discrete crack increment and to optimize the computational cost. The associated parameters considered for all the simulations in Section 4.4 are given in Table 4.

In the end, the computational cost depends largely on the choice of the minimum element size $h_{min}$, which is not trivial, especially in cyclic situations, and also rarely discussed.

It is somehow admitted that the size of the finite elements in the fracture process zone should be small enough to properly discretize the monotonic and cyclic plastic zones [69]. To estimate the appropriate element size in the vicinity of the crack
<table><caption>**Table 4**
Parameters for the mesh adaptation and crack insertion strategy.</caption><thead><tr><th>Parameters</th><th>ε<sub>D</sub><br>(-)</th><th>h<sub>min</sub><br>(mm)</th><th>Δθ<br>(°)</th><th>D<sup>crit</sup><br>(-)</th><th>ℓ<sub>inc</sub><br>(mm)</th></tr></thead><tbody><tr><td>Value</td><td>0.01</td><td>0.010</td><td>10</td><td>0.95</td><td>0.06</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 3.** Comparison of the local and non-local responses for a cyclic loading with the parameters for AD730™ at 550 °C on a single FE: the local and non-local damage models ($ℓ_c$ = 0.01 mm) lead to the same curves.</figcaption></figure>

front, [69,70] proposed the following criterion:

$$h_{min} = \frac{1}{10} r_{mpz} = \frac{1}{10} \left[ \frac{\alpha}{\pi} \left( \frac{K_{max}}{\sigma_y} \right)^2 \right] \qquad (59)$$

where $r_{mpz}$ is the radius of the monotonic plastic zone (in Irwin's sense [71]), $\alpha$ is a constraint factor ($\alpha$ = 1 for plane stress, and $\alpha$ = $(1 - 2\nu)^2$ for plane strain conditions) and $K_{max}$ is the maximum applied level reached during the FCG test. Although it is clear that the size of the monotonic plastic zone $r_{mpz}$ depends on the maximum load level, a compromise has to be found for the entire simulation.

For example, for the SEN-T specimen geometry subjected to fatigue loading described later on in Section 4.4, the estimated minimum element size is derived from two limit situations corresponding to two different crack lengths: $a$ = 1.0 mm ($K_{max}$ = 25.6 MPa$\sqrt{m}$) and $a$ = 8.0 mm ($K_{max}$ = 179.4 MPa$\sqrt{m}$). The criterion described by Eq. (59) leads to the corresponding evaluations for the minimum element size: $h_{min}$ ≈ 4.0 µm and $h_{min}$ ≈ 186.0 µm (at 550 °C), respectively. A compromise can be made by choosing an intermediate refinement from $h_{min}$ = 10.0 up to 20.0 µm for the whole simulation in order to minimize the computational cost while ensuring the reliability of the numerical results.

### 4.3. Validation of the constitutive damage model

This subsection is dedicated to the numerical validation of the continuous model presented in Section 2.1, without any mesh adaptation or crack insertion.

#### 4.3.1. Validation of the incremental fatigue damage model

The first objective is to validate the numerical implementation of the local and non-local incremental models with respect to fatigue damage. This is done by simulating a simple cyclic uniaxial tension–compression test with the parameters for AD730™ at 550 °C thanks to a single FE.<fnref n="2" /> To verify that the time evolution is similar for both local and non-local models, the non-local length was set to a value close to zero (i.e. $ℓ_c$ = 0.01 mm).

The results presented in Fig. 3 show that the numerical implementation of the gradient-enhanced damage model ensures a good convergence towards the limit case of vanishing characteristic length $ℓ_c$ → 0 mm (i.e. when the non-local model degenerates into a local one), as expected. In both cases, due to the strong coupling between damage and elastic–viscoplastic constitutive equations, the stiffness decreases with increasing damage, leading to a stress-softening cyclic response, as visible in Fig. 3(a). Furthermore, Fig. 3(b) shows that the evolution of fatigue damage is exactly the same for both models.

<footnote n="2">² Note that the implicit non-local formulation allows to perform the simulation of a non-local model on a single FE as Eq. (53) only requires data that can be obtained thanks to the product of the form functions and nodal values of an element, not the values in the surrounding area as for an integral formulation.</footnote>
<figure><img src="image_8.png" /><figcaption>Fig. 4. Geometry of the double-notched specimen subjected to monotonic tensile loading with corresponding dimensions and associated boundary conditions.</figcaption></figure>

### 4.3.2. Validation of the ductile damage

Once the fatigue damage evolution and its impact on material stiffness have been validated with a single FE, the proper spatial regularization and the impact of ductile damage still need to be validated. This is done all at once thanks to a second validation test with monotonic loading. A 2D tensile test on a double-notched specimen, originally proposed by [72] and presented in Fig. 4, is considered. Plane strain conditions are assumed and a displacement loading rate of $\dot{u}_y = 1 \times 10^{-2} \text{ mm s}^{-1}$ is applied. The material parameters are those listed in Table 3 for 700°C, except for the internal length $\ell_c$ set to 0.4 mm. This is necessary to be able to consider larger element sizes and still have at least three elements within the localization band.

To evaluate the problem regularization, three meshes with average element sizes $h_e$ of 0.3 mm, 0.2 mm and 0.1 mm are considered. The results for the local and non-local models are compared in terms of total damage distribution Fig. 5 and global force-displacement curves Fig. 6.

For the local model, the total damage fields in Fig. 5(a)–5(b) tend to localize in the smallest band that can be represented with the chosen spatial discretization: consequently, the smaller the element size, the smaller the localization band. Thus, the smaller the element size, the smaller the dissipated energy. This can be confirmed by observing the areas under the global curves in Fig. 6(a), which decrease with element size. This dependence of the results on the element size shows that the local model does not adequately represent the physics of the damage and fracture process.

In contrast, for the non-local model, the size of the localization band does not depend on the element size, Fig. 5(c)–5(d). In addition, all global force-displacement curves in Fig. 6(b) tend to converge to a unique solution with mesh refinement.

These results thus illustrate that the proposed regularized model satisfactorily enables to overcome the mesh dependency problems encountered with the local model for monotonic loadings.

## 4.4. Validation of the model with damage-to-crack transition for fatigue crack growth problems

Since the continuous incremental constitutive model is validated, we can now proceed to the application of the entire numerical strategy – from regularized damage evaluation up to crack increment insertion using crack path tracking, mesh adaptation and field transfer – to a fatigue crack growth problem. Our goal is to show that it is possible to obtain results that are similar to those found in the literature on fatigue crack growth problems, especially within the framework of LEFM.

### 4.4.1. Problem geometry and load case

In this study, we consider a Single Edge Notched tension (SEN-T) specimen, which is a fairly standard choice for fatigue crack growth analysis [32]. Only the region of interest of the SEN-T specimen is modeled. The associated geometry and boundary conditions are shown in Fig. 7. An initial pre-crack of 0.6 mm has been inserted into the initial geometry to consider sufficiently long crack propagation within the LEFM regime. The area in the vicinity of the crack front is meshed with the finest FE size, according to Eq. (59), to properly capture the non-linear phenomena occurring near the singularity.

A uniform stress load is applied continuously and cyclically to the upper boundary (in the y-direction) with a 1 Hz frequency, a load ratio of $R_\sigma = 0.05$ and a maximum stress at peak of $\sigma_\infty = 400 \text{ MPa}$. The displacements along the x-direction on both the upper and lower edges are free. A node is fixed to prevent the rigid body motion. The specimen is modeled with plane strain conditions (the specimen thickness is $t = 4.0 \text{ mm}$).
<figure><img src="image_9.png" /><figcaption>**Fig. 5.** Total damage fields obtained for the tensile test on the double notched specimen. For the local model, at the end of the simulations, the damage is localized in a band that extends over the areas associated with one or two integration points. For the non-local model, the damage localization bandwidth for $u_y = 0.1689$ mm is similar for all mesh sizes.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Force-displacement curves for the double-notched specimen with three mesh sizes.</figcaption></figure>

### 4.4.2. Parallel with the Paris regime in FCG

**Methodology to draw an analogy with LEFM.** One of our goals is to demonstrate the ability of the model to reproduce the Paris regime for a sufficiently long fatigue crack. For LCF conditions, the well-known Paris law relates the crack propagation to the range of the SIF through a relation $da/dN \sim (\Delta K)^m$ where $m$ is an exponent depending on the material and testing conditions (temperature,
<figure><img src="image_11.png" /><figcaption>Fig. 7. Left: Geometry, boundary conditions and mesh of the considered SEN-T specimen subjected to cyclic loading in the 2D case. Right: Close-up view of the FE mesh around the notch and the pre-crack.</figcaption></figure>

frequency, environment, etc.) [2]. However, with a CDM-based model, such global information is not directly accessible in the FE calculation. The methodology for estimating these quantities is detailed here.

First, even if the model is established within a time-incremental formalism, with a sinusoidal loading in the present case, the conversion from time to a number of loading cycles N is trivial. Besides, the position of the crack front is easily obtained thanks to the mesh adaptation strategy and leads directly to the estimation of the crack length a. Finally, the evaluation of the SIF range ΔK requires the knowledge of the so-called Murakami's formula [73] for the SEN-T specimen. According to this formula, a polynomial function relates the loading parameters to the geometric properties of the specimen:

$$ \Delta K = \frac{F_{max} - F_{min}}{t \cdot w} \sqrt{\pi a} f(\phi) = (1 - R_{\sigma}) \sigma_{\infty} \sqrt{\pi a} f(\phi) \quad \text{with} \quad \phi = a/w \\ = (1 - R_{\sigma}) \sigma_{\infty} \sqrt{\pi a} (1.122 - 0.231 \phi + 10.550 \phi^2 - 21.710 \phi^3 + 30.382 \phi^4) \quad (60) $$

where a is the actual crack length, σ∞ the far-field applied stress, Rσ the load ratio, t and w the thickness and width of the specimen, respectively.

It is then possible to estimate the current crack lengths after successive remeshing steps together with Eq. (60) to obtain an equivalent ΔK.

**Results analysis.** The fatigue crack growth simulation for the considered SEN-T specimen with the material parameters for AD730™ at 700 °C resulted in 86 crack increment insertions to obtain a final crack length of approximately 2.01 mm.

Contour plots of the non-local damage around the crack front are reported in Fig. 8(a) & 8(b) for two different crack lengths (i.e. a = 0.71 mm and a = 2.01 mm). Damage is localized in a band with a width that is controlled by the non-local length scale $\ell_c$ and that is larger than one element. This tends to indicate an efficient spatial regularization. The width of the localization band is quasi-constant as the crack grows, although a slight widening of this zone can be observed in the wake of the crack for the longest crack (i.e. a = 2.01 mm). This is due to the coarsening of the mesh in areas that are no longer active in order to save computational cost. This coarsening should have little or no effect, since the active areas with non-linear phenomena are normally concentrated ahead of the crack front (where the mesh is fine), while the areas near the open lips behind the crack front correspond to very low energy regions. This argument is supported by Fig. 8(d) & 8(e) showing the contour plots of the von Mises equivalent stress field. Note that among the white zones indicating very low energy areas, those immediately around the crack lips are the result of stress softening in the damage localization band. The crack therefore propagates in a near zero von Mises region where material degradation is complete. In addition, ahead of the crack front, the von Mises stress distribution seems conform to the theoretical repartition. Comparison of the distributions for the two crack lengths shows that the stress magnitude increases with the crack length, as expected.

The resulting FCG curve of the simulation is shown in Fig. 8(c). A remarkable result is the quasi-linear trend of the FCG curve with increasing ΔK for a sufficiently long crack (i.e. a ≥ 500 µm). This indicates a stable propagation regime that can be assimilated to a Paris regime. The resulting Paris' exponent here is equal to m = 4.90, which is a slightly larger value than usual ones that can be explained by the chosen set of severe damage parameters (see Section 4.2), hence resulting in faster crack growth rates and shorter lifetimes. Nevertheless, the order of magnitude is equivalent to common values for which m lies in the range 1.8–4.0 [1].

Finally, note that the proposed strategy could allow fatigue crack propagation to continue. However, going beyond a = 2.01 mm for this particular test would no longer allow the assumption of homogeneous boundary conditions (see the distribution of the von Mises equivalent stress concentration in Fig. 8(e)).
<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Overall fatigue crack growth curve for the FE calculation with some mechanical fields during a standard fatigue loading at 700°C (sinus 1 Hz, $\sigma$ = 400 MPa, $R_g$ = 0.05).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 9. Distribution of the non-local damage (left) and von Mises equivalent stress (right) for the three different mesh sizes at the maximum load of the 700th cycle (view at contours) - 700 °C.</figcaption></figure>

## 4.5. Sensitivity of the damage-to-crack numerical strategy

Before discussing fatigue crack growth results obtained by the proposed numerical strategy, it is worth checking the sensitivity of the damage-to-crack transition method to its numerical parameters. First, the mesh independence needs to be valid for fatigue crack propagation prior to any further cyclic calculations. Then, the sensitivity of the continuous-discontinuous strategy to the length of the crack increment $\ell_{inc}$ and to the critical non-local damage $\bar{D}_{crit}$ triggering crack insertion need to be assessed.

### 4.5.1. Mesh convergence analysis for fatigue crack growth

The gradient enhanced damage model has been shown in Section 4.3.2 to provide mesh insensitive results under monotonic loading case up to crack initiation. We now seek to validate our proposal to combine it with a continuous-discontinuous strategy to describe the crack growth process under cyclic loading conditions thanks to an additional mesh convergence analysis. For this purpose, the same SEN-T fatigue crack growth problem has been simulated with three minimum element sizes (only in the FPZ where the non-linear phenomena are localized), namely $h_{min} = 0.015$, 0.01 and 0.005 mm. The loading and numerical parameters all remain the same as in Section 4.4.2.

The results for the three minimum mesh sizes $h_{min}$ can be compared in terms of distribution of non-local damage and von Mises stress fields in Figs. 9 & 10 (on the left and right images, respectively), for two cycles: the 700th and the 850th. For the 700th cycle, the non-local damage shows comparable iso-contours. The coarser mesh in Fig. 9(a) shows a slight delay compared to the two finest discretizations in Fig. 9(b) & 9(c), probably due to a higher error discretization. The same applies for the equivalent stress profiles.

For the 850th cycle in Fig. 10, the crack is longer and the difference in terms of crack length is a little bit more noticeable: the two finest mesh sizes (0.01 and 0.005 mm in Fig. 10(b)–10(c)) have almost the same crack length while the delay of the coarser mesh, in Fig. 10(a), is more pronounced.

These qualitative observations are confirmed by the more quantitative evaluation of the crack length evolution with the number of loading cycles reported in Fig. 11(a): the crack lengths are very close for the three different minimum sizes with a growing gap with the coarser mesh curve, as can be expected with a larger error discretization. Most importantly, it can be observed that, for the two smallest element sizes, the difference in the number of loading cycles to reach the same crack length is about 10 cycles, which is negligible compared to the number of loading cycles that are generally considered in LCF conditions. The results then appear to effectively converge upon mesh refinement.

Since one of the aims of this study is to demonstrate the ability to obtain results comparable to standard LEFM analyses, the FCG curves are presented in a $da/dN$ vs. $\Delta K$ plot in Fig. 11(b). One can easily see in this plot an equivalent quasi-linear evolution of the FCG rates with increasing ranges of SIF for the three minimum mesh sizes. The fact that the corresponding Paris lines have almost the same slope confirms that similar crack growth rates are obtained for the three mesh sizes. This again tends to confirm the efficiency of the proposed numerical approach for the simulation of FCG in LCF conditions.
<figure><img src="image_14.png" /><figcaption>Fig. 10. Distribution of the non-local damage (left) and von Mises equivalent stress (right) for the three different mesh sizes at the maximum load of the 850th cycle (view at contours) - 700°C.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 11. Mesh convergence analysis on a SEN-T specimen under fatigue loading conditions at 700°C (sinus 1 Hz, $\sigma_{\infty} = 400$ MPa, $R_{\sigma} = 0.05$, plane strain conditions). Each point corresponds to a single mesh adaptation procedure associated with a crack increment insertion, since the crack length remains constant between two remeshing steps.</figcaption></figure>

Note that the slight difference with the coarser mesh is not even visible in the Paris curves since only the points corresponding to the linear Paris regime are taken into account. The last points correspond to the final stage with an unstable crack growth, discussed hereafter in Section 4.6.1.

### 4.5.2. Sensitivity of the numerical strategy to the critical damage

The first numerical parameter controlling the crack growth kinetics in our modeling approach is the critical non-local damage $\bar{D}_{crit}$ (see Table 4). The *Marching Ridges algorithm* is used to search for the crack orientation thanks to the local maximum of the degradation variable $\bar{D}$, see Section 3.2. Once the crack direction has been found, the insertion is triggered only if damage has reached its critical value $\bar{D}_{crit}$ over the specified length $\ell_{inc}$ of the crack increment (a sensitivity analysis to this second parameter is discussed hereafter in Section 4.5.3). This section is dedicated to the sensitivity of the numerical continuous-discontinuous strategy to the value of the critical damage $\bar{D}_{crit}$.
<figure><img src="image_16.png" /><figcaption>**Fig. 12.** Sensitivity of the numerical damage-to-crack strategy to the critical non-local damage $\bar{D}_{crit}$ on the FCG on a SEN-T specimen under pure fatigue loading conditions at 550 °C (sinus 1 Hz, $R_{\sigma}$ = 0.05, plane strain conditions).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 13.** Sensitivity of the numerical damage-to-crack strategy to the length of the crack increment $\ell_{inc}$ on the FCG on a SEN-T specimen under pure fatigue loading conditions at 550 °C (sinus 1 Hz, $R_{\sigma}$ = 0.05, plane strain conditions).</figcaption></figure>

The Fig. 12(a) plots the $a-N$ results for two distinct values of the critical damage, $\bar{D}_{crit} \in [0.8; 0.9]$. A similar trend can be seen in each curve, suggesting that the fatigue lifetime is identical regardless the value of critical damage. Moreover, the $da/dN - \Delta K$ plots shown in Fig. 12(b) demonstrate identical Paris lines, thus confirming the statement regarding equivalent fatigue crack growth kinetics. The value of the critical damage $\bar{D}_{crit}$ here only triggers crack insertion: a lower value indicates an earlier insertion of the crack increment. As can be seen from the exponential form of the damage-related softening function given in Eq. (2), there is only a slight effect between the two considered (but large enough) values for $\bar{D}_{crit}$. Nevertheless, a sufficiently large value for this numerical parameter must be considered in order to insert a crack into a zone that is almost completely damaged and experiencing maximum stress softening, where only a small amount of energy needs to be dissipated during the crack growth process. This aspect is important for the FE solver to easily retrieve the equilibrium.

### 4.5.3. Sensitivity of the numerical strategy to the crack increment size

The second numerical parameter controlling the crack growth kinetics in the proposed modeling approach is the size of the crack increment $\ell_{inc}$. Once the crack direction has been found and the non-local damage has reached its critical value over the length of the crack increment, conform remeshing is triggered and a crack increment of length $\ell_{inc}$ is inserted. The aim of this section is to study how sensitive the damage-to-crack strategy is to the size of the crack increment $\ell_{inc}$.

The Fig. 13 shows the insensitivity of the numerical method to the length of the crack increment $\ell_{inc}$ on the FCG predictions. One can see in the Fig. 13(a) the crack length $a$ against the number of loading cycles $N$ for different values of the crack increment length. Using three different values for $\ell_{inc} \in [10; 15; 20]$ µm but equivalent minimum element size $h_{min} = 10$ µm, one can notice that the fatigue lifetime is equivalent. The numerical parameter $\ell_{inc}$ yields similar results, thus not influencing the results for FCG regime. This result was expected, given that the numerical tools had already been evaluated for ductile fracture prediction in steel
<figure><img src="image_18.png" /><figcaption>Fig. 14. Effect of the prescribed boundary conditions on the FCG on a SEN-T specimen under pure fatigue loading conditions at 550 °C (sinus 1 Hz, R<sub>σ/μ<sub>y</sub></sub> = 0.05, plane strain conditions).</figcaption></figure>

in the work of [57]. Using an equivalent damage-to-crack transition strategy and various crack increment sizes, the authors of that study obtained similar crack paths and crack growth kinetics for cup-cone failure.

As for the estimated crack growth kinetics, the Fig. 13(b) gives the corresponding da/dN-ΔK plots for the 3 values of crack increments. Similar crack growth kinetics can be observed for the three different configurations, as evidenced by identical Paris' lines.

## 4.6. Complementary analyses and discussion

Apart from the Paris-like crack growth kinetics observed in Section 4.4.2 validating the proposed FCG modeling approach and the robustness of the numerical damage-to-crack transition strategy discussed in Section 4.5, there are other features of the LEFM theory that can be recovered by the proposed approach (in particular, the effect of the boundary condition nature, the load ratio and loading frequency).

### 4.6.1. Sensitivity to applied boundary conditions

The relatively rapid increase in the crack growth rate observed after its linear evolution can be attributed to the fact that the loading is stress-controlled (i.e. force-controlled), as in the vast majority of standard LEFM models. Such an asymptotic trend is generally attributed to the intrinsic fracture toughness of the material, which is reached when the crack is long enough to cause sudden failure. The observed trend hence seems reasonable regarding standard FCG problems and is expected to depend on the loading mode.

To demonstrate this, we compare the effects of stress- and displacement-controlled conditions in the case of the SEN-T specimen. In order to avoid non-homogeneous boundary conditions In addition, to evaluate the capability of the cyclic damage model under different loading conditions, we consider the material parameters for AD730™ at 550 °C (see Table 3) for which material stiffness, cyclic hardening and viscous behavior are different to those at 700 °C [32,68].

The evolution of the crack length as a function of the loading cycles is shown in Fig. 14(a), for both loading cases. In the case of stress-controlled loading, the FCG rates are globally higher. This may be explained by the constant increase of the SIF. This is not the case for displacement-controlled loading which may allow a deceleration of the crack advance and thus a longer fatigue lifetime.

In Fig. 14(b), the FCG curves are plotted in a da/dN vs. ΔK frame. It can be seen that the stress-controlled case has a higher slope in the Paris regime, which is consistent with the previously observed faster crack growth rate. In this Fig. 14(b), note that the asymptotic tendency of the curve in the case of stress-controlled loading should not be confused with the real fracture toughness ΔK<sub>c</sub> of the material: the unstable crack growth is due to the nature of the boundary condition which, due to the increasing stress level, causes significant damage growth in the FPZ before the fracture toughness limit is reached.

### 4.6.2. Link with Irwin's theory

In close connection with the stress field, the process zone ahead of the crack-tip is subjected to plastic yielding. The existence of such a plastic zone at the level of the singularity in an elastic-plastic strain-hardening material has been extensively studied in the last decades [71]. A common estimate of its size is attributed to Irwin's theory [71] (especially for perfect plasticity):

$$r_{mpz} = \frac{1}{6\pi} \left( \frac{K_{max}}{\sigma_y} \right)^2 \qquad (61)$$
<figure><img src="image_19.png" /><figcaption>**Fig. 15.** Evolution of the size of the plastic zone at the tip of an advancing fatigue crack at maximum load (550 °C, sinus 1 Hz, $\sigma_{\infty}$ = 400 MPa, $R_{\sigma}$ = 0.05, plane strain conditions).</figcaption></figure>

with $K_{max}$ the maximum SIF (in mode I), $\sigma_y$ the yield strength and $r_{mpz}$ the plane strain monotonic plastic zone size at the crack-tip.

Such a plastic zone is known to be responsible (among other mechanisms) for crack closure effects as well as for the modification of both crack kinetics and orientation. It is also well known that crack growth is largely influenced by the effect of this plastic zone at the front of the crack. The latter, as it grows, leaves a wake of yielded material, also referred to as the *plastic wake*.

In order to study the evolution of a sufficiently significant plastic zone with the proposed non-local model, an additional test on a SEN-T specimen is considered. Since the calculation discussed in Section 4.4.2 at 700 °C has yielded to inhomogeneous boundary conditions (see the distribution of the von Mises equivalent stress concentration in Fig. 8(e)), material parameters at 550 °C are still considered here. In fact, the material behavior of AD730™ at 700 °C proved to be less stiff [32]. The loading case at 550 °C therefore appears to be less detrimental in terms of (visco)plastic and damage activity ahead of the crack front.

One can easily observe from Fig. 15(a) to 15(e) that the size of the (visco)plastic zone increases with the crack length. As the crack advances, the range of stress intensity factors $\Delta K$ also increases from 22.7 to 40.1 MPa$\sqrt{m}$. In addition, the “butterfly” shape of the plastic zone is also noticeable, as expected for pure fatigue loading under plane strain conditions [74]. According to [75], the cyclic plastic zone size is assumed to be equal to the radius of the zone starting from the crack front to the boundary of the plastic zone along the crack direction. In Fig. 15(a) to 15(e), the theoretical radius of the plastic zone is represented by an orange circle, while the calculated viscoplastic activity is represented at integration points by the “butterfly-like” pattern. It can be observed that the numerical and theoretical radii for the viscoplastic zones are in fairly good agreement according to Irwin’s theory [71] as long as the crack grows.

### 4.6.3. Effect of load ratio on FCG

In this section, the effect of load ratio is investigated. Simulations on a SEN-T specimen are still performed with material parameters of AD730™ at 550 °C, see Table 3. The same stress-controlled loading case was simulated with two different positive loading ratios $R_{\sigma}$, namely 0.05 and 0.12. The range of these load ratios is relatively narrow in order to still perform ECG computation with a limited number of loading cycles while avoiding convergence problems due to fast damage growth. Note that negative load ratios are not considered to avoid mechanisms associated with crack lip closure (not considered in this study). The stress amplitude is kept unchanged and is equal to $\Delta\sigma/2 = 190$ MPa. Consequently, the minimum and maximum stress levels evolve ($\sigma_{max}$ = 400 and 430 MPa, respectively), and since the load ratio is increased, the mean stress follows the same trend.

The evolution of the crack length as a function of the cycle number in Fig. 16(a) shows that a higher load ratio results in a lower fatigue lifetime, as expected. Indeed, since the mean stress is higher at high load ratios, the material is more likely to develop plastic strain and damage. This is confirmed by the slight upward shift of the fatigue crack growth curve for the higher loading ratio Fig. 16(b), with the slope of the Paris lines remaining almost similar.

### 4.6.4. Sensitivity to applied loading frequencies

The effect of the loading frequency on the FCG is examined here. Simulations are run with material parameters of AD730™ at 550 °C, with $\ell_c$ = 0.15 mm and $h_{min}$ = 0.05 mm (different from previous analyses). The loading is stress-controlled with an applied stress of $\sigma_{\infty}$ = 400 MPa at the upper boundary. Three loading frequencies are prescribed, namely 0.1, 1.0 and 5.0 Hz.

The FCG curves in Fig. 17(a) show that there is no significant difference in the slopes of the $da/dN - \Delta K$ curves for the three loading frequencies in the Paris regime. The crack growth rates seem to be almost equivalent, as expected for pure fatigue [76].

However, the evolution of the crack length with the number of cycles in Fig. 17(b) indicates an effect of the loading frequency on the fatigue lifetime for this particular test case. Indeed, there is a difference in the asymptotic limit of the number of cycles: the lower the loading frequency, the shorter the fatigue lifetime. In the present model, the frequency effect is more pronounced than
<figure><img src="image_20.png" /><figcaption>Fig. 16. Load ratio effect on the FCG on a SEN-T specimen under pure fatigue loading conditions at 550 °C (sinus 1 Hz, Δσ/2 = 190 MPa, plane strain conditions).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 17. Effect of loading frequency on the FCG on a SEN-T specimen under pure fatigue loading conditions at 550 °C (sinus 0.1, 1.0 & 5.0 Hz, σ∞ = 400 MPa, Rσ = 0.05, plane strain conditions).</figcaption></figure>

reported in the literature for FCG [76], but this can be explained by the fact that the test is stress-controlled and that, at some point, crack propagation becomes unstable with increasing SIF, as discussed in Section 4.6.1. Crack evolution is then highly dependent on the stresses, which can be affected by frequency changes since the model is intrinsically viscoplastic. Indeed, using Eqs. (22) and (24) of the material model, one gets the viscous stress σv:

$$ \sigma_v = K \arcsin h \left( \frac{\lambda}{\dot{\epsilon}_0} \right)^{1/n} \quad \text{if } F_y > 0 \qquad (62) $$

The contribution of this viscous stress to both the total stress level and associated viscoplasticity will therefore depend on the loading frequency and will affect the material response.

To illustrate the difference in material response with the change in loading frequency, the von Mises stress fields for the three loading frequencies are shown in Fig. 18 for the 600th cycle, when the slower loading frequency already led to a larger crack (see Fig. 17(b)). It can be seen that the slower loading frequency leads to higher stress levels over larger areas.

As a result, the increment of the non-local damage ΔDt is larger for the slowest loading frequency, as can be seen in Fig. 19(a). This explains why slower loading frequencies lead to globally shorter fatigue lifetimes in this stress-controlled case. Similar observations can be made when comparing both Figs. 18(b)–18(c) and 19(b)–19(c) with increasing loading frequencies.

Apart from these numerical results, it should be mentioned that for another Nickel-based superalloy, namely DA Inco718, similar results of shorter fatigue lifetime with decreasing loading frequency in pure fatigue at 550 °C have also been evidenced by Fessler et al. [76]. This finally supports the numerical results given in Fig. 17 for AD730™ at 550 °C.

Finally, the time-dependent sensitivity of the model is here associated with the unified viscoplastic framework, as no creep mechanisms have been introduced in this preliminary work. In the case of creep-fatigue interaction or predominant time-dependent effects, a dedicated creep damage variable may be required in accordance with experimental results.
<figure><img src="image_22.png" /><figcaption>Fig. 18. von Mises equivalent stress at peak level of the 600th cycle for the three different loading frequencies.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 19. Increment of non-local damage at peak level of the 600th cycle for the three different loading frequencies.</figcaption></figure>

## 5. Conclusions and future scope

In this paper, a non-local time-incremental fatigue damage model for predicting fatigue crack growth has been established. First, a local unified elastic-viscoplastic material model, which is fully coupled to damage mechanisms, has been proposed based on the thermodynamics of irreversible processes together with the concepts of Continuum Damage Mechanics. This model allows the simulation of complex loading schemes thanks to two types of damage (i.e. fatigue and ductile). Then, an implicit gradient extension of this local model has been proposed to solve the spurious spatial localization associated with the stress softening response. The last part of the proposed strategy is the application of a continuous-to-discontinuous transition using mesh adaptation to accurately capture the crack growth kinetics.

After validating some basic features of the model on simple cases, the entire modeling strategy was evaluated for LCF crack propagation on a Nickel-based superalloy with 2D computations. Macroscopic results in terms of fatigue crack growth curves $a = f(N)$ and crack growth kinetics $da/dN - \Delta K$ have been shown to be similar to those obtained from Linear Elastic Fracture Mechanics theory. In particular, the proposed CDM-based approach enabled to recover the Paris crack growth regime. Although the continuous-discontinuous strategy could be applied with the element deletion technique, the use of mesh adaptation enabled to obtain results that are mesh independent during crack propagation. Moreover, it has been demonstrated that the numerical strategy yields identical results with respect to the parameters (namely $\ell_{inc}$ and $\bar{D}_{crit}$) associated with the damage-to-crack transition step. Finally, the study of the effect of the type of boundary conditions (i.e. stress or displacement controlled) has enabled to better understand the results of the following sensitivity analyses related to some extrinsic loading parameters (i.e. load ratio and loading frequency). In all the considered analyses, the proposed modeling has given satisfactory results with respect to the theoretical LEFM results.

Now that the concept has been validated with 2D computations, further studies will focus on calibrating the damage model based on experimental tests on AD730™ in order to reproduce the real LCF crack growth kinetics. In addition, more complex loading schemes with possible overloads will be considered to take full advantage of the proposed time-incremental multi-damage model. Additional damage mechanisms could also be considered, e.g. to take into account creep and oxidation, as they tend to develop near the crack front at high temperatures. Note that the 3D extension should not pose any technical difficulty, as all the tools are already 3D compatible (see, for example, the works from [19,57] for ductile fracture prediction using the Z-set FE solver). However, it may require some parallelization work to reduce the computational cost. In this regard, recent developments have also been achieved in the Z-set solver (see e.g. the work from [77]).

Finally, thermo-mechanical fatigue problems could also be studied using the proposed time incremental approach, since temperature changes could be intrinsically taken into account at each successive time step.
CRediT authorship contribution statement

O. Voreux: Writing – review & editing, Writing – original draft, Visualization, Methodology, Investigation, Formal analysis, Conceptualization. S. Feld-Payet: Writing – review & editing, Writing – original draft, Visualization, Supervision, Methodology, Conceptualization. P. Kanouté: Writing – review & editing, Supervision. S. Kruch: Writing – review & editing, Supervision.

Declaration of generative AI in scientific writing

During the preparation of this work the authors used the DeepL/write tool in order to check for the readability of their writing. After using this tool/service, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work has been carried out within a PhD project co-founded by ONERA - The French Aerospace Lab and SafranTech - SAFRAN group. The latter is gratefully acknowledged for supporting this research project.

Appendix A. Material integration scheme

In Section 2.1 a multi-mechanism elastic-viscoplastic damage model has been established within a consistent thermodynamic framework associated with CDM. In order to solve the constitutive evolution equations governing the system, a relevant numerical integration algorithm should be considered. Indeed, the whole set of evolution equations Table 1 governing the non-linear and damageable behavior of the material needs to be integrated with respect to time during the whole loading. Such an integration scheme occurs at the element level for each integration (Gauss) point. The numerical implementation of the present model is achieved within the Z-set implicit solver [56]. Time incremental equations are integrated following the fully implicit (backward Euler) scheme [78], ensuring a good precision as well as a good stability under large time steps.

From a numerical point of view and assuming the use of the standard FE method, the purpose of the constitutive model is to compute the stress tensor σ at the end of the load increment accordingly to the strain tensor increment Δε associated with the prescribed displacement of the structure, during a time step Δt. Moreover, the non-local extension presented in the Section 2.3 associated with the strong coupling assumption between non-local damage and constitutive equations result in the simultaneous resolution of both equilibrium equation associated with strain increment Δε and diffusion-like equation associated with non-local total damage increment ΔD̅. As a consequence, the non-local constitutive model relates increments of gradient variables ΔV̇ = {Δε, ΔD̅} to their flux counterparts ΔV̇ = {Δσ, ΔD̅}.

In the following, for the sake of readability, we omit the index “t” for both total local and non-local damages. The FE resolution scheme available in the Z-set code [56] allows for the fast determination of the consistent tangent operator needed for the construction of the elementary stiffness matrix. To that extent, a strain-based formulation should be used in order to consistently derive the material tangents. Thus, considering a strain-driven problem, that is, controlled by the total strain tensor ε, the time interval is given by [t_n, t_{n+1}] whereas the time increment is Δt. On the one hand, the values of the state variables at t_n are known {ε_n, ε_n^e, ε_n^i, α_n^k, r_n^i, λ_n, D_n, D̅_n}. On the other hand, the prescribed variables are given at current time t_{n+1} = t_n + Δt by {ε_{n+Δt}, D̅_{n+Δt}}.

For the strain-driven problem, the set of internal variables is chosen as V̇ = {ε^e, α^k, r^i, λ, D}. These independent variables need to be updated at the end of the time step, that is at t_n + Δt. Through time discretization of the evolution equations given in Table 1, we obtain the corresponding residual equation system to be solved:

$$
\mathcal{R} = \{ \mathcal{R}_{\epsilon^e}, \mathcal{R}_{\alpha^k}, \mathcal{R}_{r^i}, \mathcal{R}_{\lambda}, \mathcal{R}_D \}^T = 0 \quad (A.1)
$$

where $\mathcal{R}$ is the vector containing the material residuals.

For the sake of clearness, we omit the index for quantities at time t_n + Δt which gives in an incremental form:

$$
\mathcal{R}_{\epsilon^e} = \Delta\epsilon^e - \Delta\epsilon + \frac{\Delta\lambda}{f(\bar{D})} \bar{n} = 0 \qquad (A.2a)
$$

$$
\mathcal{R}_{\alpha^k} = \Delta\alpha^k - \frac{\Delta\lambda}{f(\bar{D})} (\bar{n} - \gamma^k \alpha^k) = 0 \qquad (A.2b)
$$

$$
\mathcal{R}_{r^i} = \Delta r^i - \frac{\Delta \lambda}{f(\bar{D})} (1 - b^i r^i) = 0 \qquad (A.2c)
$$

$$
\mathcal{R}_{\lambda} = \Delta\lambda - \Delta t \dot{\epsilon}_0 \sinh \left\langle \frac{F}{K} \right\rangle_{+}^{n} = 0 \qquad (A.2d)
$$
<caption>**Table A.5**
Numerical integration algorithm associated with the strain-driven problem.</caption>

1. Given $\{\epsilon_n, \epsilon_n^e, \epsilon_n^{im}, \alpha_n^k, r'_n, \lambda_n, D_n, \bar{D}_n\}$ at $t_n$ and $\{\epsilon, \bar{D}\}$ at $t_{n+1}$
2. Evaluate the trial stress: $\sigma^{tr}$
3. Set the initial guess of the variables: $\mathcal{V}_{int}^0 = \{\epsilon_n^e, \alpha_n^k, r'_n, \lambda_n, D_n\}$
4. Formulate the residual equation system: $\mathbf{R} |_{\mathcal{V}_{int}^0 = \mathcal{V}_{int}^0}$ and check the convergence: if $|| \mathbf{R} |_{\mathcal{V}_{int}^0 = \mathcal{V}_{int}^0} || < \epsilon_{tol}$ then EXIT
5. Calculate the algorithmic tangent operator: $\mathbf{J} |_{\mathcal{V}_{int}^0 = \mathcal{V}_{int}^0}$
6. Solve the linearized system: $\mathbf{R} |_{\mathcal{V}_{int}^0 = \mathcal{V}_{int}^0} + \mathbf{J} |_{\mathcal{V}_{int}^0 = \mathcal{V}_{int}^0} \delta\mathbf{\Delta}\mathcal{V}_{int} = 0$
7. Update the unknowns: $\mathcal{V}_{int}^0 \leftarrow \mathcal{V}_{int}^0 + \delta\mathbf{\Delta}\mathcal{V}_{int}$ and go to 3.

$$ \mathbf{R}_D = \Delta D - \Delta \lambda \left\langle \frac{Y_f - Y_{f0}}{S_f} \right\rangle_+^{m_f} e^{(\beta_f \bar{D}_t)} - (1 - \bar{D}_t) \frac{\Delta \lambda}{f(\bar{D})} \text{tr} \bar{n} = 0 \quad (A.2e) $$

By solving Eq. (A.1), the value of $\mathcal{V}_{int}$ at time $t_n + \Delta t$ can be obtained. As Eq. (A.1) is a set of non-linear equations, the Newton-Raphson method is used to derive the numerical solution. For that purpose, the residual equation system is linearized into:

$$ \mathbf{R} |_{\mathcal{V}_{int} = \mathcal{V}_{int}^0} + \mathbf{J} |_{\mathcal{V}_{int} = \mathcal{V}_{int}^0} \delta\mathbf{\Delta}\mathcal{V}_{int} = 0 \quad (A.3) $$

where $\mathcal{V}_{int}^0$ is the initial value of the variables at each iteration step and $\mathbf{J}$ is the algorithmic tangent operator, also known as material Jacobian matrix. The matrix $\mathbf{J}$, although cumbersome to calculate is highly necessary. For the sake of brevity, the values of each term are given in another section, see Appendix B.

Once the residual equation system Eq. (A.1) is solved, the independent internal variables $\mathcal{V}_{int} = \{\epsilon^e, \alpha^k, r^i, \lambda, D\}$ at time $t_{n+1}$ can be determined. Finally, thermodynamics associated variables can be updated at the end of the time step following state relations recalled in Table 1.

As a matter of synthesis, the numerical integration algorithm for a strain-driven problem is summarized in Table A.5.

## Appendix B. Material Jacobian matrix

Last but not at least, the exact evaluation of the Jacobian matrix $\mathbf{J}$ in Eq. (A.3) allows for a fast identification of the consistent tangent operator of the material associated with the time discretization [78,79]. Indeed, the non-linear equations at the global FE level are usually treated by Newton's method. Thus, associated with the need for implementing this material model in a general finite element code, the tangent moduli that are consistent with the integration algorithm should also be calculated, apart from the updated stress and local damage (seen as a *source term*). Note that all these moduli form the global material tangent. Such moduli, if well evaluated, ensure a quadratic convergence of the global resolution scheme [78].

The algorithmic Jacobian matrix $\mathbf{J}$ is formally expressed as:

$$ \mathbf{J} = \frac{\partial \mathbf{R}}{\partial \Delta \mathcal{V}_{int}} = \begin{bmatrix} \frac{\partial R_{\epsilon^e}}{\Delta \epsilon^e} & \frac{\partial R_{\epsilon^e}}{\Delta \alpha^k} & \frac{\partial R_{\epsilon^e}}{\Delta r^i} & \frac{\partial R_{\epsilon^e}}{\Delta \lambda} & \frac{\partial R_{\epsilon^e}}{\Delta D} \\ \frac{\partial R_{\alpha^k}}{\Delta \epsilon^e} & \frac{\partial R_{\alpha^k}}{\Delta \alpha^k} & \frac{\partial R_{\alpha^k}}{\Delta r^i} & \frac{\partial R_{\alpha^k}}{\Delta \lambda} & \frac{\partial R_{\alpha^k}}{\Delta D} \\ \frac{\partial R_{r'_n}}{\Delta \epsilon^e} & \frac{\partial R_{r'_n}}{\Delta \alpha^k} & \frac{\partial R_{r'_n}}{\Delta r^i} & \frac{\partial R_{r'_n}}{\Delta \lambda} & \frac{\partial R_{r'_n}}{\Delta D} \\ \frac{\partial R_{\lambda^e}}{\Delta \epsilon^e} & \frac{\partial R_{\lambda^e}}{\Delta \alpha^k} & \frac{\partial R_{\lambda^e}}{\Delta r^i} & \frac{\partial R_{\lambda^e}}{\Delta \lambda} & \frac{\partial R_{\lambda^e}}{\Delta D} \\ \frac{\partial R_{\bar{D}}}{\Delta \epsilon^e} & \frac{\partial R_{\bar{D}}}{\Delta \alpha^k} & \frac{\partial R_{\bar{D}}}{\Delta r^i} & \frac{\partial R_{\bar{D}}}{\Delta \lambda} & \frac{\partial R_{\bar{D}}}{\Delta D} \end{bmatrix} \quad (B.1) $$

where each residual is given in Eq. (A.2) in the last Appendix A. Details for the determination of each term of the Jacobian matrix $\mathbf{J}$ Eq. (B.1) can be found in [32].

## Appendix C. Consistent tangent moduli

An efficient scheme for calculating the consistent tangent moduli is introduced here. First, according to Simo and Taylor [79], the tangent moduli consistent with the algorithm and necessary for the elementary stiffnesses are defined by:

$$ C_{uu} = \frac{\partial \sigma (\epsilon, \mathcal{V}_{int, n}, \bar{D}_n, \bar{D}_{n+1})}{\partial \epsilon} \quad (C.1a) $$

$$ C_{u\bar{D}} = \frac{\partial \sigma (\epsilon, \mathcal{V}_{int, n}, \bar{D}_n, \bar{D}_{n+1})}{\partial \bar{D}} \quad (C.1b) $$

$$ C_{\bar{D}u} = \frac{\partial D (\epsilon, \mathcal{V}_{int, n}, \bar{D}_n, \bar{D}_{n+1})}{\partial \epsilon} \quad (C.1c) $$
$$C_{\overline{DD}} = \frac{\partial D(\epsilon, \mathbf{V}_{int,n}, \overline{D}_n, \overline{D}_{n+1})}{\partial \overline{D}} \qquad (C.1d)$$

where $C_{uu}$ has the dimension of a 4th-order tensor, $C_{u\overline{D}}$ and $C_{\overline{Du}}$ that of a 2nd-order tensor and $C_{\overline{DD}}$ that of a scalar. The full coupling between damage and behavior leads to a non-symmetric tangent modulus $C_{\sigma\epsilon}$. Moreover, the non-local extension of the local damage model also leads to a non-symmetric material tangent matrix $\mathbb{K}^{\text{mat}}$ defined as:

$$\mathbb{K}^{\text{mat}} = \frac{\partial \Delta \mathbf{V}_{\text{OUT}}}{\partial \Delta \mathbf{V}_{\text{IN}}} = \begin{bmatrix} \frac{\partial \Delta \sigma}{\partial \Delta \epsilon} & \frac{\partial \Delta \sigma}{\partial \Delta \overline{D}} \\ \frac{\partial \Delta \overline{D}}{\partial \Delta \epsilon} & \frac{\partial \Delta \overline{D}}{\partial \Delta \overline{D}} \end{bmatrix} = \begin{bmatrix} C_{uu} & C_{u\overline{D}} \\ C_{\overline{Du}} & C_{\overline{D}\overline{D}} \end{bmatrix} \qquad (C.2)$$

For any given values of strain $\epsilon$ and non-local damage $\overline{D}$, stress tensor $\sigma$ and total damage $D$ are in fact determined by the residual equation system Eq. (A.1), thus implying that $\mathbf{R} = 0$ holds for any $\epsilon$ and $\overline{D}$. By calculating the total derivatives of $\mathbf{R}$ in terms of $\epsilon$ and $\overline{D}$, we obtain the following results:

$$\left[ \frac{\mathrm{d}\mathbf{R}}{\mathrm{d}\epsilon} \frac{\mathrm{d}\mathbf{R}}{\mathrm{d}\overline{D}} \right] = \left\{ \frac{\mathrm{d}\mathbf{R}_{\epsilon^e}}{\mathrm{d}\epsilon} \frac{\mathrm{d}\mathbf{R}_{\epsilon^e}}{\mathrm{d}\overline{D}}, \frac{\mathrm{d}\mathbf{R}_{\alpha^k}}{\mathrm{d}\epsilon} \frac{\mathrm{d}\mathbf{R}_{\alpha^k}}{\mathrm{d}\overline{D}}, \frac{\mathrm{d}\mathbf{R}_{\overline{D}^e}}{\mathrm{d}\epsilon} \frac{\mathrm{d}\mathbf{R}_{\overline{D}^e}}{\mathrm{d}\overline{D}}, \frac{\mathrm{d}\mathbf{R}_{\overline{D}^k}}{\mathrm{d}\epsilon} \frac{\mathrm{d}\mathbf{R}_{\overline{D}^k}}{\mathrm{d}\overline{D}}, \frac{\mathrm{d}\mathbf{R}_{\lambda^e}}{\mathrm{d}\epsilon} \frac{\mathrm{d}\mathbf{R}_{\lambda^e}}{\mathrm{d}\overline{D}}, \frac{\mathrm{d}\mathbf{R}_{\overline{D}^k}}{\mathrm{d}\epsilon} \frac{\mathrm{d}\mathbf{R}_{\overline{D}^k}}{\mathrm{d}\overline{D}} \right\}^T = 0 \qquad (C.3)$$

Such result implies that:

$$\mathbf{J} \cdot \mathbf{Y} + \mathbb{B} = 0 \qquad (C.4)$$

where $\mathbf{J}$ is the material's Jacobian matrix (or algorithmic tangent operator) given in Eq. (B.1) and gained at convergence of the local iterative scheme for material integration, and $\mathbb{B}$ corresponds to:

$$\mathbb{B} = \left[ \frac{\partial \mathbf{R}}{\partial \mathbf{V}_{\text{IN}}} \right] = \left[ \frac{\partial \mathbf{R}_{\epsilon^e}}{\partial \epsilon} \frac{\partial \mathbf{R}_{\epsilon^e}}{\partial \overline{D}}, \frac{\partial \mathbf{R}_{\alpha^k}}{\partial \epsilon} \frac{\partial \mathbf{R}_{\alpha^k}}{\partial \overline{D}}, \frac{\partial \mathbf{R}_{\overline{D}^e}}{\partial \epsilon} \frac{\partial \mathbf{R}_{\overline{D}^e}}{\partial \overline{D}}, \frac{\partial \mathbf{R}_{\overline{D}^k}}{\partial \epsilon} \frac{\partial \mathbf{R}_{\overline{D}^k}}{\partial \overline{D}}, \frac{\partial \mathbf{R}_{\lambda^e}}{\partial \epsilon} \frac{\partial \mathbf{R}_{\lambda^e}}{\partial \overline{D}}, \frac{\partial \mathbf{R}_{\overline{D}^k}}{\partial \epsilon} \frac{\partial \mathbf{R}_{\overline{D}^k}}{\partial \overline{D}} \right]^T \qquad (C.5)$$

$$\mathbf{Y} = \left[ \frac{\partial \mathbf{V}_{\text{int}}}{\partial \mathbf{V}_{\text{IN}}} \right] = \left[ \frac{\partial \epsilon^e}{\partial \epsilon} \frac{\partial \epsilon^e}{\partial \overline{D}}, \frac{\partial \alpha^k}{\partial \epsilon} \frac{\partial \alpha^k}{\partial \overline{D}}, \frac{\partial \overline{D}^e}{\partial \epsilon} \frac{\partial \overline{D}^e}{\partial \overline{D}}, \frac{\partial \lambda^e}{\partial \epsilon} \frac{\partial \lambda^e}{\partial \overline{D}}, \frac{\partial \overline{D}}{\partial \epsilon} \frac{\partial \overline{D}}{\partial \overline{D}} \right]^T \qquad (C.6)$$

Finally, we notice here that consistent tangent moduli given in Eq. (C.1) can be deduced from the resolution of the linear system Eq. (C.4) These moduli consist on some blocks of $\mathbf{Y}$ once the system is solved. Some mathematical operations are finally needed to gain the proper tangent moduli as the problem is strain-driven, thus written in terms of strain fields, whereas the consistent tangent moduli Eqs. (C.1a) and (C.1b) need to be written in terms of the Cauchy stress tensor.

The numerical computation of the tangent moduli directly comes from the resolution of the Jacobian matrix Eq. (B.1). On the one hand, a small variation of the input variables $\delta \mathbf{V}_{\text{IN}}$ directly results in a change of the equilibrium through the state variables $\mathbf{V}_{\text{int}}$. Using the chain rule and tensorial calculus, the problem can be expressed in terms of the derivatives of the problem's variables, viz.:

$$\delta \mathbf{R} = \frac{\partial \mathbf{R}}{\partial \Delta \mathbf{V}_{\text{IN}}} \cdot \delta \Delta \mathbf{V}_{\text{IN}} + \frac{\partial \mathbf{R}}{\partial \Delta \mathbf{V}_{\text{int}}} \cdot \delta \Delta \mathbf{V}_{\text{int}} = \frac{\partial \mathbf{R}}{\partial \Delta \mathbf{V}_{\text{IN}}} \cdot \delta \Delta \mathbf{V}_{\text{IN}} + \mathbf{J} \cdot \delta \Delta \mathbf{V}_{\text{int}} = \mathbf{0} \qquad (C.7)$$

with $\mathbf{J}$ given by Eq. (B.1) so that:

$$\delta \Delta \mathbf{V}_{\text{int}} = -\mathbf{J}^{-1} \cdot \frac{\partial \mathbf{R}}{\partial \Delta \mathbf{V}_{\text{IN}}} \cdot \delta \Delta \mathbf{V}_{\text{IN}} \qquad (C.8)$$

On the other hand, the output variables are expressed as functions of the internal variables only. Using Eq. (C.8), one gets:

$$\delta \mathbf{V}_{\text{OUT}} = \frac{\partial \Delta \mathbf{V}_{\text{OUT}}}{\partial \Delta \mathbf{V}_{\text{int}}} \cdot \delta \Delta \mathbf{V}_{\text{int}} + \frac{\partial \Delta \mathbf{V}_{\text{OUT}}}{\partial \Delta \mathbf{V}_{\text{IN}}} \cdot \delta \Delta \mathbf{V}_{\text{IN}} \qquad (C.9)$$

which, combining Eq. (C.8) with Eq. (C.9) leads to:

$$\delta \mathbf{V}_{\text{OUT}} = \left\{ \frac{\partial \Delta \mathbf{V}_{\text{OUT}}}{\partial \Delta \mathbf{V}_{\text{int}}} \left[ -\mathbf{J}^{-1} \cdot \frac{\partial \mathbf{R}}{\partial \Delta \mathbf{V}_{\text{IN}}} \right] + \frac{\partial \Delta \mathbf{V}_{\text{OUT}}}{\partial \Delta \mathbf{V}_{\text{IN}}} \right\} \cdot \delta \Delta \mathbf{V}_{\text{IN}} \qquad (C.10)$$

Looking at both Eqs. (C.2) and (C.10), the consistent tangent matrix is equal to:

$$\mathbb{K}^{\text{mat}} = \frac{\partial \Delta \mathbf{V}_{\text{OUT}}}{\partial \Delta \mathbf{V}_{\text{IN}}} = \left\{ \frac{\partial \Delta \mathbf{V}_{\text{OUT}}}{\partial \Delta \mathbf{V}_{\text{int}}} \left[ -\mathbf{J}^{-1} \cdot \frac{\partial \mathbf{R}}{\partial \Delta \mathbf{V}_{\text{IN}}} \right] + \frac{\partial \Delta \mathbf{V}_{\text{OUT}}}{\partial \Delta \mathbf{V}_{\text{IN}}} \right\} \qquad (C.11)$$

Details for the determination of each term of the matrices $\mathbb{B}$ in Eq. (C.5) and $\mathbf{Y}$ in Eq. (C.6), useful for the determination of the tangent matrix $\mathbb{K}^{\text{mat}}$ in Eq. (C.11), can be found in [32].

## Data availability

The data that has been used is confidential.
References

[1] Jones R. Fatigue crack growth and damage tolerance. Fatigue Fract Eng Mater Struct 2014;37:463-83. http://dx.doi.org/10.1111/ffe.12155, URL: https://onlinelibrary.wiley.com/doi/abs/10.1111/ffe.12155.

[2] Paris P, Erdogan F. A critical analysis of crack propagation laws. J Basic Eng 1963;85:528-33. http://dx.doi.org/10.1115/1.3656900.

[3] Pommier S, Risbet M. Time-derivative equations for fatigue crack growth in metals. Int J Fract 2005;131:79-106. http://dx.doi.org/10.1007/s10704-004-3633-9, URL: https://link.springer.com/article/10.1007/s10704-004-3633-9.

[4] Pommier S. Development of an incremental model for fatigue crack growth predictions. Aerosp Lab 2015;9:1-12. http://dx.doi.org/10.12762/2015.AL09-11.

[5] Lu Z, Liu Y. Small time scale fatigue crack growth analysis. Int J Fatigue 2010;32:1306-21. http://dx.doi.org/10.1016/j.ijfatigue.2010.01.010, URL: http://www.sciencedirect.com/science/article/pii/S014211231000023X.

[6] Bouvard JL, Chaboche JL, Feyel F, Gallerneau F. A cohesive zone model for fatigue and creep-fatigue crack growth in single crystal superalloys. Int J Fatigue 2009;31:868-79. http://dx.doi.org/10.1016/j.ijfatigue.2008.11.002, URL: http://www.sciencedirect.com/science/article/pii/S0142112308002521.

[7] Nijin IS, Banerjee A. An incremental elastic-plastic triaxiality dependent fatigue model. Int J Fract 2021;227:39-55. http://dx.doi.org/10.1007/s10704-020-00494-8.

[8] Besson J, Berdin C, Bugat S, Desmorat R, Feyel F, Forest S, et al. Local approach to fracture. 1st ed. Paris, France: Les Presses de l'Ecole des Mines; 2004.

[9] Peerlings RHJ, Brekelmans WaM, Borst Rd, Geers MGD. Gradient-enhanced damage modelling of high-cycle fatigue. Internat J Numer Methods Engrg 2000;49:1547-69. http://dx.doi.org/10.1002/1097-0207(20001230)49:12<1547::AID-NME16>3.0.CO;2-D, URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/1097-0207%2820001230%2949%3A12%3C1547%3A%3AAID-NME16%3E3.0.CO%3B2-D.

[10] Pijaudier-Cotot G, Bažant ZP. Nonlocal damage theory. J Eng Mech 1987;113:1512-33. http://dx.doi.org/10.1061/(ASCE)0733-9399(1987)113:10(1512), URL: http://ascelibrary.org/doi/10.1061/%28ASCE%290733-9399%281987%29113%3A10%281512%29.

[11] Peerlings RH, de Borst R, Brekelmans W, de Vree J. Gradient enhanced damage for quasi-brittle materials. Internat J Numer Methods Engrg 1996;39:3391-403. http://dx.doi.org/10.1002/(SICI)1097-0207(19961015)39:19<3391::AID-NME7>3.0.CO;2-D, URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/%28SICI%291097-0207%2819961015%2939%3A19%3C3391%3A%3AAID-NME7%3E3.0.CO%3B2-D.

[12] Bažant ZP, Jirásek M. Nonlocal integral formulations of plasticity and damage: Survey of progress. J Eng Mech 2002;128:1119-49. http://dx.doi.org/10.1061/(ASCE)0733-9399(2002)128:11(1119), URL: http://ascelibrary.org/doi/10.1061/%28ASCE%290733-9399%282002%29128%3A11%281119%29.

[13] Lorentz E, Cano V. A nonlocal formulation applied to ductile damage. Lat Am J Solids Struct 2005;2:17-28, URL: https://lajss.org/index.php/LAJSS/article/view/65.

[14] Francfort GA, Marigo JJ. Revisiting brittle fracture as an energy minimization problem. J Mech Phys Solids 1998;46:1319-42. http://dx.doi.org/10.1016/S0022-5096(98)00034-9, URL: https://www.sciencedirect.com/science/article/pii/S0022509698000349.

[15] Maniar Y, Métais B, Kuczynska M, Kabachiev A, Binkele P, Schmauder S. Nonlocal continuum damage mechanics approach in the finite element simulation of lead-free solder joints. In: 2017 18th international conference on thermal, mechanical and multi-physics simulation and experiments in microelectronics and microsystems. 2017, p. 1-8. http://dx.doi.org/10.1109/EuroSimE.2017.7926254.

[16] Yin B, Zreid I, Zhao D, Ahmed R, Lin G, Kaliske M. Thermomechanical fatigue life prediction of metallic materials by a gradient-enhanced viscoplastic damage approach. Internat J Numer Methods Engrg 2022;123:2042-75. http://dx.doi.org/10.1002/nme.6926, URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.6926, eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.6926.

[17] Khalil Z, Elghazouli AY, Martínez-Pañeda E. A generalised phase field model for fatigue crack growth in elastic-plastic solids with an efficient monolithic solver. Comput Methods Appl Mech Engrg 2022;388:114286. http://dx.doi.org/10.1016/j.cma.2021.114286, URL: https://www.sciencedirect.com/science/article/pii/S0045782521005867.

[18] Feld-Payet S. Amorçage et propagation de fissures dans les milieux ductiles non locaux [Ph.D. thesis], Paris, France: Ecole Nationale Supérieure des Mines de Paris; 2010, URL: https://pastel.archives-ouvertes.fr/pastel-00583011/document.

[19] Feld-Payet S, Chiaruttini V, Besson J, Feyel F. A new marching ridges algorithm for crack path tracking in regularized media. Int J Solids Struct 2015;71:57-69. http://dx.doi.org/10.1016/j.ijsolstr.2015.04.043, URL: http://www.sciencedirect.com/science/article/pii/S0020768315002656.

[20] Besson J, Bleyer J, Feld-Payet S, Gourgues-Lorenzon A-F, Hannard F, Helfer T, et al. MEALOR II damage mechanics and local approach to fracture. 1st ed. Zenodo; 2023, http://dx.doi.org/10.5281/zenodo.10125169.

[21] Feld-Payet S. Damage-to-crack transition. In: Numerical methods for strong nonlinearities in mechanics. John Wiley & Sons, Ltd; 2024, p. 293-361, URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/9781394340507.ch7.

[22] Mesgarnejad A, Imanian A, Karma A. Phase-field models for fatigue crack growth. Theor Appl Fract Mech 2019;103:10228. http://dx.doi.org/10.1016/j.tafmec.2019.102282, URL: https://www.sciencedirect.com/science/article/pii/S0167844218306712.

[23] Xue F, Cheng T-L, Lei Y, Wen Y-H. Phase-field modeling of crack growth under coupled creep-fatigue. Int J Fatigue 2024;108577. http://dx.doi.org/10.1016/j.ijfatigue.2024.108577, URL: https://www.sciencedirect.com/science/article/pii/S0142112324004365.

[24] Baruah S, Singh IV. Low-cycle fatigue simulation of ductile materials using elasto-plastic gradient damage approach. Int J Mech Sci 2024;276:109370. http://dx.doi.org/10.1016/j.ijmecsci.2024.109370, URL: https://www.sciencedirect.com/science/article/pii/S0020740324004120.

[25] Yun K, Wang Z, Chang M, Liu J, Kim T-J, Son N, et al. A computational methodology for simulating quasi-brittle fracture problems. Comput Struct 2019;215:65-79. http://dx.doi.org/10.1016/j.compstruc.2019.02.003, URL: https://www.sciencedirect.com/science/article/pii/S004579491831513X.

[26] Oliver J. Modelling strong discontinuities in solid mechanics via strain softening constitutive equations. Part 1: Fundamentals. Internat J Numer Methods Engrg 1996;39:3575-600. http://dx.doi.org/10.1002/(SICI)1097-0207(19961115)39:21<3575::AID-NME65>3.0.CO;2-E, URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/%28SICI%291097-0207%2819961115%2939%3A21%3C3575%3A%3AAID-NME65%3E3.0.CO%3B2-E.

[27] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Internat J Numer Methods Engrg 1999;45:601-20. http://dx.doi.org/10.1002/(SICI)1097-0207(19990620)45:5<601::AID-NME598>3.0.CO;2-S, URL: http://onlinelibrary.wiley.com/doi/10.1002/(SICI)1097-0207(19990620)45:5%3C601::AID-NME598%3E3.0.CO;2-S/abstract.

[28] Khoei AR, Eghbalian M, Azadi H, Saffar H. Numerical simulation of ductile crack growth under cyclic and dynamic loading with a damage-viscoplasticity model. Eng Fract Mech 2013;99:169-90. http://dx.doi.org/10.1016/j.engfracmech.2012.12.009, URL: http://www.sciencedirect.com/science/article/pii/S0013794412004705.

[29] Hamon F, Henaff G, Halm D, Gueguen M, Billaudreau T. A damage model for fatigue crack propagation from moderate to high DeltaK levels. Fatigue Fract Eng Mater Struct 2012;35:160-72. http://dx.doi.org/10.1111/j.1460-2695.2011.01603.x, URL: http://onlinelibrary.wiley.com/doi/10.1111/j.1460-2695.2011.01603.x/abstract.

[30] Chaboche J-L. Constitutive equations for cyclic plasticity and cyclic viscoplasticity. Int J Plast 1989;5:247-302. http://dx.doi.org/10.1016/0749-6419(89)90015-6, URL: http://www.sciencedirect.com/science/article/pii/0749641989900156.

[31] Lemaitre J, Chaboche J-L, Benallal A, Desmorat R. Mécanique des matériaux solides. In: Sciences sup. 3ème ed. Paris, France: Dunod; 2009.

[32] Voreux O. High temperature fatigue crack growth modeling in nickel-based superalloys using a local approach to fracture [Ph.D. thesis], Paris, France: Université Paris-Saclay; 2022, URL: tel.archives-ouvertes.fr/tel-03714432.

[33] Murakami S, Liu Y. Mesh-dependence in local approach to creep fracture. Int J Damage Mech 1995;4:230-50. http://dx.doi.org/10.1177/105678959500400303.

[34] Dimitrijević B, Hackl K. A method for gradient enhancement of continuum damage models. Tech Mech 2008;28:43-52.
[35] Kaminski M, Kanouté P, Kruch S, Busso EP, Chaboche J-L. A high temperature fatigue damage model for single crystal superalloys. Mater High Temp 2016;33:412-24. http://dx.doi.org/10.1080/09603409.2016.1193995.

[36] Halphen B, Nguyen QS. Sur les matériaux standards généralisés. J Mécanique 1975;14:39-63.

[37] Germain P, Nguyen QS, Suquet P. Continuum thermodynamics. J Appl Mech 1983;50:1010-20. http://dx.doi.org/10.1115/1.3167184.

[38] Coleman BD, Gurtin ME. Thermodynamics with internal state variables. J Chem Phys 1967;47:597-613. http://dx.doi.org/10.1063/1.1711937, URL: https://aip.scitation.org/doi/abs/10.1063/1.1711937.

[39] Lemaître J, Chaboche J-L. Aspect phénoménologique de la rupture par endommagement. J Mécanique Appliquée 1978;2:317-65, URL: https://trid.trb.org/view/1048150.

[40] Sellars CM, McTegart WJ. On the mechanism of hot deformation. Acta Metall 1966;14:1136-8.

[41] Chaboche J-L. A review of some plasticity and viscoplasticity constitutive theories. Int J Plast 2008;24:1642-93. http://dx.doi.org/10.1016/j.ijplas.2008.03.009, URL: http://www.sciencedirect.com/science/article/pii/S0749641908000582.

[42] Pétry C. Caractérisation et modélisation du comportement et de l'endommagement d'alliages métalliques sur une grande plage de température [Ph.D. thesis], Cachan, France: École normale supérieure de Cachan - ENS Cachan; 2006, URL: https://tel.archives-ouvertes.fr/tel-00133696.

[43] Saanouni K, Forster C, Hatira FB. On the anelastic flow with damage. Int J Damage Mech 1994;3:140-69. http://dx.doi.org/10.1177/105678959400300203.

[44] Besson J. Continuum models of ductile fracture: A review. Int J Damage Mech 2010;19:3-52. http://dx.doi.org/10.1177/1056789509103482.

[45] Gurson AL. Continuum theory of ductile rupture by void nucleation and growth: Part I—Yield criteria and flow rules for porous ductile media. J Eng Mater Technol 1977;99:2-15. http://dx.doi.org/10.1115/1.3443401.

[46] Rousselier G. Ductile fracture models and their potential in local approach of fracture. Nucl Eng Des 1987;105:97-111. http://dx.doi.org/10.1016/0029-5493(87)90234-2, URL: http://www.sciencedirect.com/science/article/pii/S0029549387902342.

[47] Besson J, Guillemer-Neel C. An extension of the Green and Gurson models to kinematic hardening. Mech Mater 2003;35:1-18. http://dx.doi.org/10.1016/S0167-6636(02)00169-2, URL: http://www.sciencedirect.com/science/article/pii/S0167663602001692.

[48] Kaminski M. Modélisation de l'endommagement en fatigue des superalliages monocrystallins pour aubes de turbine en zone de concentration de contrainte [Ph.D. thesis], Paris, France: Ecole Nationale Supérieure des Mines de Paris; 2007, URL: http://www.theses.fr/2007ENMNP1504.

[49] Chaboche J-L, Boudifa M, Saanouni K. A CDM approach of ductile damage with plastic compressibility. Int J Fract 2006;137:51-75. http://dx.doi.org/10.1007/s10704-005-3061-5, URL: https://link.springer.com/article/10.1007/s10704-005-3061-5.

[50] Taira S. Lifetime of structures subjected to varying load and temperature. In: Hoff NJ, editor. Creep in structures, IUTAM symposia. Berlin, Heidelberg: Springer; 1962, p. 96-124. http://dx.doi.org/10.1007/978-3-642-86014-0-6.

[51] Chaboche J-L, Gallerneau F. An overview of the damage approach of durability modelling at elevated temperature. Fatigue Fract Eng Mater Struct 2001;24:405-18. http://dx.doi.org/10.1046/j.1460-2695.2001.00415.x, URL: http://onlinelibrary.wiley.com/doi/abs/.

[52] Peerlings RHJ, de Borst R, Brekelmans WAM, Geers MGD. Localisation issues in local and nonlocal continuum approaches to fracture. Eur J Mech A Solids 2002;21:175-89. http://dx.doi.org/10.1016/S0997-7538(02)01211-1, URL: http://www.sciencedirect.com/science/article/pii/S0997753802012111.

[53] Feld-Payet S, Besson J, Feyel F. Finite element analysis of damage in ductile structures using a nonlocal model combined with a three-field formulation. Int J Damage Mech 2011;20:655-80. http://dx.doi.org/10.1177/1056789511405935, URL: https://journals.sagepub.com/doi/abs/10.1177/1056789511405935.

[54] Broumand P, Khoei AR. The extended finite element method for large deformation ductile fracture problems with a non-local damage-plasticity model. Eng Fract Mech 2013;112:113:97-125. http://dx.doi.org/10.1016/j.engfracmech.2013.10.002, URL: http://www.sciencedirect.com/science/article/pii/S0013794413003202.

[55] El Ouazani Tuhami A, Feld-Payet S, Quilici S, Osipov N, Besson J. A two characteristic length nonlocal GTN model: Application to cup-cone and slant fracture. Mech Mater 2022;171:104350. http://dx.doi.org/10.1016/j.mechmat.2022.104350, URL: https://www.sciencedirect.com/science/article/pii/S0167663622001272.

[56] ONERA, Mines ParisTech. Z-set software - Non-linear material & structure analysis suite. 2021, URL: http://www.zset-software.com/.

[57] El Ouazani Tuhami A, Feld-Payet S, Quilici S, Osipov N, Besson J. Simulations of complex crack paths using a robust and cost-efficient continuous-discontinuous approach. Int J Solids Struct 2023;274:112199. http://dx.doi.org/10.1016/j.ijsolstr.2023.112199, URL: https://www.sciencedirect.com/science/article/pii/S0020768323000963.

[58] Zienkiewicz OC, Zhu ZJ. The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique. Internat J Numer Methods Engrg 1992;33:1331-64. http://dx.doi.org/10.1002/nme.1620330702, URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620330702.

[59] Zienkiewicz OC, Zhu ZJ. The superconvergent patch recovery and a posteriori error estimates. Part 2: Error estimates and adaptivity. Internat J Numer Methods Engrg 1992;33:1365-82. http://dx.doi.org/10.1002/nme.1620330703, URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620330703.

[60] Perić D, Yu J, Owen D. On error estimates and adaptivity in elastoplastic solids: Applications to the numerical simulation of strain localization in classical and Cosserat continua. Internat J Numer Methods Engrg 1994;37:1351-79. http://dx.doi.org/10.1002/nme.1620370806, URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620370806.

[61] Andrade Pires F, de Souza Neto E, Owen D. On the finite element prediction of damage growth and fracture initiation in finitely deforming ductile materials. Comput Methods Appl Mech Engrg 2004;193:5223-56. http://dx.doi.org/10.1016/j.cma.2004.01.038, URL: http://www.sciencedirect.com/science/article/pii/S0045782504002737.

[62] Ladevèze P, Pelle J-P. La maîtrise du calcul en mécanique linéaire et non linéaire. In: Etudes en mécanique des matériaux & des structures. 1st ed. Paris: Hermès Science - Lavoisier; 2001.

[63] Frey P, George P-L. Mesh generation: application to finite elements. 2nd ed. Wiley, ISTE Ltd. & John Wiley & Sons; 2008, URL: https://onlinelibrary.wiley.com/doi/book/10.1002/9780470611116.

[64] Boromand B, Zienkiewicz OC. Recovery procedures in error estimation and adaptivity. Part II: Adaptivity in nonlinear problems of elasto-plasticity behaviour. Comput Methods Appl Mech Engrg 1999;176:127-46. http://dx.doi.org/10.1016/S0045-7825(98)00333-8, URL: http://www.sciencedirect.com/science/article/pii/S0045782598003338.

[65] El Ouazani Tuhami A. Simulation of ductile failure over long distances with a non-local GTN model [Ph.D. thesis], Paris, France: PSL Research University; 2022, URL: https://pastel.archives-ouvertes.fr/tel-03869005.

[66] Chiaruttini V, Riolo V, Feyel F. Advanced remeshing techniques for complex 3D crack propagation. In: Proceedings of the 13th international conference on fracture. 2013, p. 1-9.

[67] Oñate E, Rojek J, Taylor RL, Zienkiewicz OC. Finite calculus formulation for incompressible solids using linear triangles and tetrahedra. Internat J Numer Methods Engrg 2004;59:1473-500. http://dx.doi.org/10.1002/nme.922, URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.922.

[68] Devaux A, Picqué B, Gervais M, Georges E, Poulain T, Héritier P. AD730 - A new nickel-based superalloy for high temperature engine rotative parts. In: Superalloys 2012-12th international symposium on superalloys. 2012, p. 911-9.

[69] McClung RC, Sehitoglu H. On the finite element analysis of fatigue crack closure - 1. Basic modeling issues. Eng Fract Mech 1989;33:237-52. http://dx.doi.org/10.1016/0013-7944(89)90027-1, URL: https://www.sciencedirect.com/science/article/pii/0013794489900271.

[70] Dougherty JD, Padovan J, Srivatsan TS. Fatigue crack propagation and closure behavior of modified 1070 steel: Finite element study. Eng Fract Mech 1997;56:189-212. http://dx.doi.org/10.1016/S0013-7944(96)00104-X, URL: https://www.sciencedirect.com/science/article/pii/S001379449600104X.
[71] Irwin G. Plastic zone near a crack and fracture toughness. In: Proceedings of the 7th sagamore ordnance materials research conference, vol. 4, New York: Syracuse University; 1960, p. 63-78.

[72] Brokken D. Numerical modelling of ductile fracture in blanking [Ph.D. thesis], Eindhoven, The Netherlands: Technische Universiteit Eindhoven; 1999, URL: https://research.tue.nl/en/publications/numerical-modelling-of-ductile-fracture-in-blanking.

[73] Murakami Y, Aoki S. Stress intensity factors handbook. Pergamon; 1987.

[74] Zhao Y, Hu D, Zhang M, Dai W, Zhang W. In situ measurements for plastic zone ahead of crack tip and continuous strain variation under cyclic loading using digital image correlation method. Metals 2020;10:273. http://dx.doi.org/10.3390/met10020273, URL: https://www.mdpi.com/2075-4701/10/2/273, number: 2 Publisher: Multidisciplinary Digital Publishing Institute.

[75] McClung RC. Crack closure and plastic zone sizes in fatigue. Fatigue Fract Eng Mater Struct 1991;14:455-68. http://dx.doi.org/10.1111/j.1460-2695.1991.tb00674.x, URL: https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1460-2695.1991.tb00674.x.

[76] Fessler E, Andrieu E, Bonnand V, Chiaruttini V, Pierret S. Relation between crack growth behaviour and crack front morphology under hold-time conditions in DA Inconel 718. Int J Fatigue 2017;96:17-27. http://dx.doi.org/10.1016/j.ijfatigue.2016.11.015, URL: http://www.sciencedirect.com/science/article/pii/S0142112316303784.

[77] Rannou J, Bovet C. Domain decomposition methods and acceleration techniques for the phase field fracture staggered solver. Internat J Numer Methods Engrg 2024;125:e7544. http://dx.doi.org/10.1002/nme.7544, URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.7544. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.7544.

[78] Besson J, Cailletaud G, Chaboche J-L, Forest S. Non-linear mechanics of materials. 2nd ed. Dordrecht: Springer; 2010, http://dx.doi.org/10.1007/978-90-481-3356-7.

[79] Simo JC, Taylor RL. Consistent tangent operators for rate-independent elastoplasticity. Comput Methods Appl Mech Engrg 1985;48:101-18. http://dx.doi.org/10.1016/0045-7825(85)90070-2, URL: http://www.sciencedirect.com/science/article/pii/0045782585900702.