A Comparative Review of Fatigue Modeling Approaches in Metals: From Total-Life Methods to Fracture-Based Models

1. Introduction

Beginning in the mid-19th century, the study of fatigue in metals has become, over the years, a field of major interest for materials scientists and structural engineers, mainly because it represents one of the principal failure mechanisms in metallic engineering structures. Accordingly, early fatigue assessment was largely based on empirical approaches, notably the stress amplitude–life (S–N) curves derived from Wöhler’s pioneering studies [1,2,3] and the strain–life ($ϵ$–N) diagrams emerging from the works of Coffin [4] and Manson [5]. Commonly classified as total-life approaches, these methodologies were originally intended for fatigue life estimation and remain widely employed to predict the service life of metallic components. Since the emergence of these pioneering approaches, fatigue prediction capabilities have advanced considerably, culminating in the development of sophisticated numerical models capable of predicting not only fatigue life but also the progressive degradation of material properties associated with the fatigue damage, including crack initiation and propagation processes that ultimately lead to complete structural failure.

Approaches exhibiting these characteristics are commonly referred to as fracture-based fatigue models. Over the years, several classes of models encompassing such features have been developed. Their diversity and breadth make it impractical to provide a comprehensive treatment of all existing formulations within a single work. Nevertheless, since the introduction of the energy-based fracture theory proposed by Griffith in 1921 [6], three classes of models have gained sufficient relevance to deserve a comparative review, namely those grounded in frameworks of fracture mechanics [7,8] and continuum damage mechanics [9,10], and phase-field formulations [11,12]. Although originally developed to describe fracture under monotonic loading conditions, many formulations belonging to these classes have subsequently been extended to account for fatigue loading scenarios. Accordingly, in addition to the total-life approaches, this work presents a review of relevant fracture formulations within these classes of models and their extension to fatigue, followed by a general comparison among the presented approaches. Furthermore, the finite element method (FEM) and the extended finite element method (XFEM) are considered to frame the discussion from a numerical perspective, as they are among the most widely employed discretization strategies for modelling fracture and fatigue.

The remainder of this paper is organized as follows. Section 2 reviews the classical total-life approaches, namely the stress–life and strain–life methods. Section 3 presents fracture-based modelling approaches grounded in fracture mechanics, continuum damage mechanics, and phase-field formulations, together with their extensions to fatigue. A comparative assessment of the reviewed approaches is provided in Section 4 and their advantages and limitations are comparatively summarized in a table. Finally, Section 5 highlights the main features of each class of model discussed in this work.

2. Total-life Approaches

Total-life approaches constitute classical methods employed in fatigue analysis and design. With a focus on failure rather than the evolution of fatigue cracks, they are fundamentally characterized by life curves described in terms of cyclic stress or strain ranges, which correspond to the so-called stress- and strain-based methods, respectively. These curves are often obtained through standard fatigue tests conducted on polished coupon specimens under either constant amplitude stress- or strain-controlled regimes [13,14]. As a results, the number of cycles to nucleate and propagate a dominant crack until the specimen complete failure is measured. A general description of the main features of each method is presented in the following.

2.1. Stress-life Method

Grounded on the premise of global elastic response of the material under sub-critical stresses, the stress-life method has been traditionally applied in high-cycle fatigue (HCF) analysis, i.e., in cases in which the fatigue failure occurs after a high number of load cycles, often more than ${10}^4$. This analysis is commonly performed using Wöhler curves obtained by the fitting of S-N experimental points. Accordingly, different approaches have been proposed to mathematically describe these curves [15,16]. Among them, the equation introduced by Basquin [17] has become particularly popular due to its simplicity. In the logarithmic form, it reads

$$log\left(N_f\right)=log\left(C\right)-k·log\left({\sigma }_a\right),$$

(1)

where ${\sigma }_a$ and $N_f$ represent the stress amplitude and number of cycle to failure, whereas C and m are model parameters fitted through the least square method. The stress amplitude, in turn, may be defined as

$${\sigma }_a=\left|{\displaystyle \frac{{\sigma }_{max}-{\sigma }_{min}}{2}}\right|,$$

(2)

with ${\sigma }_{max}$ and ${\sigma }_{min}$ being the maximum and minimum stresses, respectively.

Figure 1 schematically illustrates the main features of typical S-N curves resulting from such a fitting procedure. As can be seen, these curves exhibit distinct patterns for ferrous and non-ferrous metallic material. In the case of ferrous alloys, they are often designed with a knee point followed by a plateau called fatigue limit. Beyond this point, a structure submitted to cyclic loads is assumed to withstand an infinite number of cycles, as far as the corresponding stress amplitude remains below this limit. In real applications, however, it is difficult to clearly define this knee point due to the scatter inherent to the fatigue phenomenon, as shown in the graphic by the range in which finite and infinite lives are overlapped. Despite that, it is a common practice to adopt values between $\sim {10}^5$ and $2\times {10}^6$ cycles for this point [18,19], avoiding – whenever possible – stress amplitudes close to the fatigue limit in infinite-life designs. Contrarily to ferrous metals, non-ferrous alloys do not exhibit a clear plateau, as illustrated in the graphic. Cutoff values are usually set for infinite-life designs though. For aluminum grades, for instance, values between ${10}^7$ and $5\times {10}^8$ are typically assumed for this purpose [19,20].

S-N curves can be determined not only for coupon specimens but also for mechanical component and engineering structures in general. However, large-scale fatigue tests are extremely costly and, as a result, the use of specimen fatigue data – often assumed as material properties – obtained under controlled conditions is preferred whenever possible. Thus, when the design and life assessment of actual components are concerned, material S-N curves are modified in an attempt to adapt this specimen-level data to the part in-service conditions. This is commonly done by modifying the material nominal fatigue limit by the so-called Marin reduction factors to account for the influence of aspects such as surface condition, part size, load type, reliability, among others [21]. In its general form, the fatigue limit of the actual component is estimated through

$$S_e=k_a{k}_b{k}_c{k}_d{k}_e{k}_f{S}_e^{\prime },$$

(3)

where $S_e^{\prime }$ is the material nominal fatigue limit and $k_a$, $k_b$, $k_c$, $k_d$, $k_e$, $k_f$ represent surface, size, loading, temperature, reliability and miscellaneous-effect factors.

In addition, the effect of stress concentrations induced by geometrical discontinuities such as grooves, holes and notches is also considered in this type of analysis. For this purpose, fatigue notch factors are applied to the nominal stress calculated using elasticity-based formulae to determine a "fatigue-effective" stresses [22,23], i.e.,

$${\sigma }_{max}=K_f·{\sigma }_{\mathrm{avg}}.$$

(4)

According to Peterson [23], for instance, this factor is given by

$$K_f=1+q(K_t-1)$$

(5)

where $K_t$ is the classical stress concentration factor used in elasticity problems and q the so-called notch sensitivity. Typically depending on geometrical features, material properties and loading condition (e.g., axial, bending and torsional), q is defined between 0 and 1, assuming values close to the unity for large notch radius. Accordingly, in many practical cases, full notch sensitivity is assumed (i.e., $q=1$ and $K_f=K_t$) as a conservative measure, given the common scatter of the data used in the q estimation. In addition to the classical approaches, alternative methods such as those based on stress gradients [24], highly-stressed-volume [25] and theory of distances [26] have also been proposed to incorporate the notch effect in these fatigue analyses. Figure 2 schematically illustrates the notch effect on the stress field of an axially loaded notched specimen and how the stress concentration and notch factors are employed to estimate the peak and maximum stress, ${\sigma }_{peak}$ and ${\sigma }_{max}$, from the nominal counterpart, ${\sigma }_{avg}$.

Standardized fatigue tests are typically conducted under a fully reversed load regime such that, the stress ratio

$$R={\displaystyle \frac{{\sigma }_{min}}{{\sigma }_{max}}}$$

(6)

is equal to -1. In addition, a mean stress

$${\sigma }_m={\displaystyle \frac{{\sigma }_{max}+{\sigma }_{min}}{2}}$$

(7)

equal to zero naturally follows from this load condition. However, load regimes leading to non-zero mean stresses are often found in real world applications. Since such mean stresses may have a significant influence in the material fatigue response, many different approaches have been proposed to account for their effects in fatigue analysis [27,28,29]. In practical terms, these approaches provide a way to calculate an equivalent stress amplitude at $R=-1$, given a stress amplitude and a mean stress at $R\ne -1$, which may be then confronted with the S-N data from fully reversed tests. Table 1 summarizes some of the most relevant mean stress correction models and their dependence on specific material parameters. Among them, the classical Geber [30], Goodman [31], Soderberg [33] and Morrow [34] models are also graphically shown in Figure 3.

Alternatively to the models presented in Table 1, Haigh [41] proposed an approach in which the mean stress correction is carried out according to four different stress ratio regimes (i.e., $R>1$, $-\infty \le R\le 0$, $0\le R\le 0.5$ and $R\ge 0.5$), as detailed in the FKM-Guideline [42].

Real fatigue load histories often exhibit varying amplitudes. Accordingly, numerous damage accumulation models have been proposed to handle the effect of cyclic load blocks of different stress amplitudes in fatigue calculations [43]. These models rely on damage parameters that cumulatively compute the contribution of each block to the overall material fatigue damage, acting ultimately as remaining life estimators. Among the existing models, the Palmgren-Miner linear accumulation rule

$$D={\displaystyle \sum _{i=1}^k}{\displaystyle \frac{N_{c_i}}{N_{f_i}}}$$

(8)

is widely used, mainly due to its simplicity [44,45]. Here $N_{c_i}$ is the number of cycles applied at a certain stress amplitude, ${\sigma }_{a_i}$, and $N_{f_i}$ the corresponding number of cycles to failure. D, in turn, represents the accumulated damage, which indicates failure when it reaches a critical value – typically equal to 1. Despite its wide practical use, this rule assumes that the damage accumulation is independent of the load block application order. As a consequence, it fails to capture non-linear load sequences and overload effects [46,47]. The appropriate use of damage accumulation rules require the adoption of an efficient cycle counting procedure, particularly for the cases of highly irregular load histories. Accordingly, the rain-flow method is the most widely used approach for this purpose [46]. It is based on the comparison of adjacent stress ranges such that a cycle is counted whenever the second range is greater than the first. Figure 4 shows an application example of the rain-flow counting procedure together with the corresponding percentage damage accumulated in each block of cyclic load.

The analytical solution of elasticity problems is usually limited to specific, idealized cases. For this reason, standard linear elastic finite element analyses (FEA) are commonly employed to computed the stress field, particularly when more complex geometries are involved. Hence, a stress norm such as von Mises is used to handle the multiaxial states arising from the FEA. Finally, the stress-life method is then applied as a post-processing operation to perform the fatigue analysis. The results from this type of analysis are usually represented either as a map of the number cycles to failure or the accumulated damage due to the application of predefined load blocks. In addition, a map of the safety factor, n, – conceptually obtained by substituting ${\sigma }_{a,R}$ and ${\sigma }_m$ by $n{\sigma }_{a,R}$ and $n{\sigma }_m$ in the equations of Table 1 – can also be plotted for design and verification purposes. Figure 5 illustrates the results of stress-life fatigue analyses carried out in conjunction with the linear elastic FEA.

2.2. Strain-life Method

In contrast to the stress-life approach, the strain-life method explicitly accounts for plasticity, providing as a result, a more robust foundation to handle fatigue problems in the low cycle domain – where macroscopic plastic strains become a relevant factor. As such, it constitutes a local approach in which the fatigue life – typically expressed by the number of reversals, $2N_c$, – is related to the total strain range, $\Delta \epsilon$, through the well-known Coffin-Manson relationship [4,5], which is given by

$$\Delta \epsilon ={\displaystyle \frac{{\sigma }_F^{\prime }}{E}}{\left(2N_c\right)}^b+{\epsilon }_F^{\prime }{\left(2N_c\right)}^c,$$

(9)

where E is the material Young’s modulus, ${\sigma }_F^{\prime }$ and ${\epsilon }_F^{\prime }$ are the fatigue strength and ductility coefficients and b and c the fatigue strength and ductility exponents. In addition, the first and second terms on the right side of this equation represent, respectively, the elastic and plastic strain ranges (i.e., $\Delta {\epsilon }_e$ and $\Delta {\epsilon }_p$). As an alternative to the classical Coffin-Manson approach, numerous strain-based methods have been proposed as fatigue life estimators. Examples include the modified universal slopes [51], uniform materials [52], Roessle-Fatemi’s hardness [53], modified Mitchell [54] and medians methods [55]. For a more comprehensive survey on them, the reader is referred to [56,57].

The fatigue constants of Equation (9) are usually determined by fitting experimental data from fully reversed strain-controlled fatigue tests. Stable stress-strain hysteresis loops – commonly recorded near the half of the specimen fatigue life – are employed for this purpose [58]. Figure 6 schematically represents a typical hysteresis loop used in this procedure, highlighting the stress range, $\Delta \sigma$, and the total, elastic, and plastic strain ranges – $\Delta \epsilon$, $\Delta {\epsilon }_e$ and $\Delta {\epsilon }_p$ – associated with this type of curve.

The fatigue life estimation resulting from Equation (9) and the respective experimental data are usually represented in terms of the number of reversal, $2N_c$ and the total strain amplitude, $\Delta \epsilon /2$, as shown in Figure 7.

Despite the usual calibration through data from fully reversed tests, Equation (9) can be modified to account for the mean stress effect in fatigue analyses involving strain ratios, $R_{\epsilon }$, different from -1. This can be done by correcting the number of cycles, $N_c$, in this expression as follows

$$N_{c{R}_{\epsilon }}=N_c{\left[{\displaystyle \frac{{\sigma }_a}{f\left({\sigma }_a,{\sigma }_m\right)}}\right]}^{-1/b},$$

(10)

where $f\left({\sigma }_a,{\sigma }_m\right)$ can be any mean stress correction model depending on the stress amplitude and the mean stress such as those presented in Table 1 [58].

Analogous to the stress-life approach, damage accumulation rules and the rain-flow cycle counting procedure are commonly employed to handle variable amplitude load scenarios in fatigue analyses using the strain-life method. However, special care should be taken when analyzing highly irregular load histories in the low cycle regime, where plastic strains may play a critical role in the fatigue failure. This is particularly important because plasticity is a path-dependent phenomenon and Equation (9) assumes independent, stabilized cycles for the fatigue life estimation. Furthermore, an additional source of error may be introduced if the Palmgreen-Miner rule is adopted as the damage accumulation model, since it does not distinguish the sequence of the applied load blocks, as previously discussed. Consequently, the application of the strain-life method in these scenarios may lead to significant inaccuracies in fatigue life estimations.

As a local approach, the application of the strain-life method at the component level requires the determination of the total strain at the part critical locations – task that is hardly feasible through analytical solutions, particularly when dealing with intricate geometries. Thus, nonlinear FEA, incorporating a suitable plasticity material model (e.g., with isotropic or kinematic hardening), may be employed to obtain stable hysteresis stress-strain loops at the component critical points. These loops can then be used in a post-processing step to estimate the fatigue life through the Coffin-Manson relationship. For complex, lengthy load histories, however, plasticity models may lead to computationally highly expensive simulations, hampering their application in such scenarios. Alternatively, linear elastic FEA have been employed in conjunction with notch correction factors to estimate the elasto-plastic stresses and strains in the part hot spots [59,60].

3. Modeling Approaches for Fracture and Their Extension to Fatigue

Unlike the total-life methods, fracture models – when applied in the context of fatigue analysis – not only serve as life indicators but also provide a description of the crack growth process. Classical examples of them include the approaches based on the fracture mechanics (FM). Beyond modeling the crack propagation, some classes of fracture models – such as those grounded on the continuum damage mechanics (CDM) and phase-field (PF) method – are conceived to capture the entire degradation process, from the material intact condition up to its final failure. Since they three represent some of the most relevant classes of fracture model, a general overview of the principles underlying each of them is presented in the following.

3.1. Models Based on the Fracture Mechanics

Relying on the premise that engineering components are inherently flawed, classical fracture models formulated within the framework of the FM are often referred to as defect-tolerant approaches. Accordingly, these models focus on the growth of a pre-existing crack such that, in the context of fatigue, the corresponding life is defined as the number of cycles required to propagate this crack from its initial size up to a certain critical dimension.

Crack growth models within the framework of FM are generally classified into main categories: the linear elastic fracture mechanics (LEFM) and the elastic-plastic frature mechanics (EPFM). The primary distinction between them lies in the extent of the plastic zone ahead of the crack tip. In LEFM, this zone is assumed to be negligible, allowing linear elastic assumptions to hold. In contrast, EPFM accounts for significant plastic deformation in the crack-tip region.

Under the LEFM assumption, the stress singularity at the crack tip is handled by the scalar-valued intensity factor,

$$K=\beta \sigma \sqrt{\pi a},$$

(11)

where $\sigma$ represents the nominal stress and a the crack size, whereas $\beta$ is a dimensionless correction factor depending on the component geometry, size, shape of the crack and loading type. Polynomials and curves defining $\beta$ for basic configurations can be found in various handbooks [62,63,64,65].

The expression in Equation (11) serves as a general formulation, but it is important to note that crack propagation can occur in three fundamental modes, or a combination of them (see Figure 8), each associated with a distinct type of loading and a corresponding stress intensity factor (SIF):

• Mode I: Opening (tensile) mode, characterized by $K_I$;
• Mode II: Sliding (in-plane shear) mode, characterized by $K_{II}$;
• Mode III: Tearing (anti-plane shear) mode, characterized by $K_{III}$.

In general terms, K provides a measure of the crack severity, such that – under monotonic loading conditions – when it reaches a critical value related with the material fracture toughness, crack propagation occurs suddenly. However, for cyclic loading, the cracks may propagate at K values well below the material fracture toughness. Thus, in fatigue scenarios, the crack growth – typically expressed in terms of the rate per cycles, $da/dN$, – depends on the SIF range, $\Delta K$, as first postulated by Paris et al. [66] through the following empirical law

$${\displaystyle \frac{da}{d{N}_c}}=C\Delta K^m,$$

(12)

where C and m are scaling constants. The SIF range, in turn, is given by

$$\Delta K=K_{max}-K_{min},$$

(13)

with $K_{max}$ and $K_{min}$ being the SIF corresponding to the maximum and minimum stresses (i.e., ${\sigma }_{max}$ and ${\sigma }_{min}$) within a given load cycle, respectively. Figure 9 illustrates, in a log-log scale, a typical crack growth curve which is usually divided into three distinct regions, namely: nucleation (I), stable propagation (II) and unstable propagation (III). Accordingly, crack propagation only starts for $\Delta K$ larger than the $\Delta K_{th}$ threshold shown in region I. After the nucleation region, the crack follows a linear growth rate described by the Paris law, which is only applicable to region II. Finally, a sudden and unstable propagation occurs in the region III.

Despite its conceptual importance, the applicability of the Paris law is rather limited because it is restricted to $\Delta K$ values in the region II of Figure 8 at constant amplitudes. In addition, the constant C and m are influenced by the stress ratio, R. Other important effects such as crack closure and environmental influences are also neglected. To address these limitation, several crack propagation models considering different factors have been proposed [67,68,69,70] and [71,72,73,74] within the LEFM and EPFM frameworks, respectively. Among them, the so-called NASGRO equation [75],

$${\displaystyle \frac{da}{dN}}=C{\left[{\displaystyle \frac{\left(1-f\right)}{\left(1-R\right)}}·\Delta K\right]}^m·{\displaystyle \frac{{\left(1-{\displaystyle \frac{\Delta K_{th}}{\Delta K}}\right)}^p}{{\left(1-{\displaystyle \frac{K_{max}}{K_c}}\right)}^q}},$$

(14)

has become a widespread model since it provides a comprehensive description of the cracking process, capturing the crack-growth behaviour across all three phases of Figure 9, as well as crack-closure effects and different loading modes. This generality, however, comes at the cost of increased complexity, requiring the calibration of up to 11 model parameters. In Equation (14), $K_c$ represent the critical stress intensity at which fracture occurs, f a crack opening function, whereas p and q are empirical constants.

Once the crack-growth model is set, the fatigue life can be determined by its integration from the initial to the final crack length. Since, the SIF range, $\Delta K$, is usually a function of the crack advance, and also the structure geometry through the correction factor $\beta$ in Equation (11), the analytical integration of these models is rarely possible so that numerical methods are commonly used for this purpose. Thus, the integration procedure is incrementally performed in a way that, in each step, the crack is advanced by an incremental length, and then the corresponding number of cycles is computed using the propagation model. Within a numerical framework, this procedure is typically broken down into the following workflow [76,77]:

• Determination of the displacement field and the maximum and minimum stress norms ahead of the crack front along the current load cycle;
• Evaluation of the minimum and maximum stress intensity factors;
• Identification of the direction for the crack propagation;
• Representation of the updated crack.

This process is recursively repeated until $K_{max}$ reaches the threshold $K_c$ and then the number of cycles to failure is obtained through the summation of the respective fractions at each step. As the most adopted numerical methods in these cases, the FEM and its variation XFEM are framed in more detail in rest of this section.

3.1.1. Determination of the Stress Intensity Factor

When considering a FEM setting, different methods may be employed to compute the SIF after determining the displacement, strain and stress fields of the structure. Among them, the stress, displacement and energy release rate methods are of common use. For mode I, the stress method prediction is given by

$$K_I=\underset{r\to 0}{lim}{\sigma }_{yy}(r,\theta =0)\sqrt{2\pi r},$$

(15)

where r is the radial distance from the crack tip and ${\sigma }_{yy}$ the stress normal to the crack plane. Similarly, the displacement approach yields

$$K_I=\underset{r\to 0}{lim}{\displaystyle \frac{E^{\prime }}{4}}\sqrt{{\displaystyle \frac{2\pi }{r}}}{u}_y,$$

(16)

where $u_y$ is the displacement perpendicular to the crack and $E^{\prime }$ the effective elastic modulus which is given by

$$\begin{array}{c}\hfill E^{\prime }=\left\{\begin{array}{cc}E\hfill & \mathrm{for}\mathrm{plane}\mathrm{stress}\hfill \\ {\displaystyle \frac{E}{1-{\nu }^2}}\hfill & \mathrm{for}\mathrm{plane}\mathrm{strain},\hfill \end{array}\right.\end{array}$$

(17)

with E being the Young’s modulus and $\nu$ the Poisson’s ratio.

In the common practice, $K_I$ values are computed at mesh nodes as function of r and then extrapolated to $r=0$. Accordingly, the inability of standard finite elements to accurately capture the stresses field in the vicinity of the crack tip makes the stress method highly sensitivity to the mesh arrangement, which often renders it as the least accurate approach. Although originally applied for mode I by Chan et al. [78], the stress and displacement methods are also extendable to modes II and III upon the consideration of the appropriate stress components and mode-specific displacement field equations.

Energy release rate methods, in turn, are based on the notion that the potential energy of a body decreases per unit of crack extension, as originally proposed by Griffith [6]. Accordingly, the energy release rate, G, may be computed by means of a path-independent integral as proposed by Rice [79]. For two dimension problems this integral reads

$$J={\int }_{\mathcal{C}}Wdy-\mathbf{t}·{\displaystyle \frac{\partial \mathbf{u}}{\partial x}}\mathrm{d}s,$$

(18)

where $\mathcal{C}$ is an arbitrary curve surrounding the crack tip, W the strain energy density, $\mathbf{t}$ the traction vector associated with the outward normal to $\mathcal{C}$, $\mathbf{u}$ the displacement vector and $\mathrm{d}s$ is an incremental arc length along $\mathcal{C}$. Following this, the mode I SIF can be obtained by

$$K_I=\sqrt{G{E}^{\prime }}=\sqrt{J{E}^{\prime }},$$

(19)

as demonstrated by Chan et al. [78]. Here J is computed through the numerical integration of Equation (18). It is important to note that Equation (19) is only valid under the LEFM assumption, i.e., under small-scale yielding at the crack tip – condition in which the plastic energy dissipation is negligible.

Despite its improved accuracy, this methods is rarely directly applicable in three-dimension problems because the integration of Equation (18) must be numerically carried out over a surface which is often difficult to define. As an alternative, the divergence theorem is used to convert it into a volume domain integral, which is computed through Gaussian quadrature [80]. In addition, virtual extension crack methods have also been proposed to overcome this limitation [81,82,83].

The standard J integral can provide the total energy release rate, G, for all fracture modes. However, while it is effective for pure mode I, II, or III, it does not inherently separate the individual mode contributions in mixed-mode fracture. To overcome this limitation, approaches such as the virtual crack closure technique (VCCT) [84] and the interaction integral (also known as the M integral) [85] have been developed to independently compute the energy release rates and, consequently, the stress intensity factors associated with each fracture mode.

3.1.2. Determination of the Crack Growth Direction

The direction of the crack propagation is reasonably well-defined when dealing with benchmark cases under pure load modes. Under mixed modes, however, the propagation direction is not known a priori. Since the FM does not inherently account for the crack path, an additional criterion must be incorporated to address this limitation [47]. Based on the notion that crack advances in the direction of least resistance, numerous approaches have been proposed for this purpose. Examples include the criteria of maximum tangential stress [86], maximum energy release rate [87], minimum strain energy [88], maximum dilatational strain energy density [89] and the minimum accumulated elastic stain energy in the crack tip vicinity [90]. In the FEM context, the maximum tangential stress criterion is often used in fatigue analysis as it offers an explicit approximated solution for the crack growth as a function of $K_I$ and $K_{II}$, which is ease to implement.

3.1.3. Representation of the Crack Growth

Numerical models grounded in the FM framework are often categorized as discontinuity-based approaches, since the topology of a crack is represented as a sharp discontinuity over the domain (see Figure 10). This discontinuous representation is arguably challenging from the numerical standpoint, particularly when using standard finite elements.

In practical terms, after the computation of the SIF and the determination of the propagation direction, a predefined crack increment, $\Delta a$, is represented in the numerical model and the corresponding number of cycles computed through the integration of the selected propagation law. In particular, to represent this crack increment in a standard FEM setting, the model mesh must be accordingly updated. Since the SIF usually depends on the crack current length, this update is commonly followed by the SIF recalculation in the subsequent propagation step. Variations of two main approaches are typically employed to perform the mesh update: one based on locally refined fixed meshes and the other on remeshing techniques.

In the locally refined approach, a global mesh is generated with a sufficiently high level of refinement in the region where the crack is expected to propagate. Then techniques such as node debonding [82], element removal at the crack front [91] or even cohesive zone elements (zero-thickness elements following a traction-separation law) [92] are employed to represent the crack advance, $\Delta a$, in each propagation step. While effective to solve simple benchmark problems, the mesh generation procedure in this approach requires a priori knowledge of the crack path to ensure accurate results, making its application in more complex cases particularly challenging. Furthermore, it suffers from mesh dependency, as the crack path is constrained to follow the mesh, regardless of the refinement level. Nevertheless, this approach is frequently used in ductile fracture simulations to avoid the high computational cost and accuracy loss associated with the mapping of the plastic strain history between different meshes during remeshing.

To address the limitations of the discretization scheme based on locally refined fixed meshes, remeshing techniques have been employed. Accordingly, the remeshing procedure can be performed at the global level [93,94,95], when the entire domain is remeshed, and at the local one [96,97], cases in which a new mesh is generated only in the region around the crack. Although global remeshing techniques may be required in cases involving complex crack paths that the other alternatives are unable to capture, local approaches are often preferred due to their lower computational cost. As illustrated in Figure 11, the mesh generation procedure in the local approaches commonly rely on the four following steps: (i) removal of the elements in the vicinity of the crack; (ii) advance of the crack increment, $\Delta a$; (iii) creation of quarter-point elements [98,99] in a rosetta pattern at the crack tip and; (iv) remeshing of the area where the elements have been previously removed.

Following a similar rational, Colombo and Gliogle [76] have recently proposed a remeshing methodology in which only the elements intersected by the crack increment, $\Delta a$, are removed. The resulting open region is then remeshed using a triangulation algorithm that duplicates the nodes along the crack faces, which are subsequently debonded (see, Figure 12). This technique optimizes the remeshing procedure by reducing its computational cost to a very minimum, while maintaining the crack advance mesh-independent. Nevertheless, the triangular elements created at the crack tip can be distorted in many cases, leading to inaccurate solutions for the SIF – particularly, when using the stress- and displacement-based methods. For this reason, the authors recommend the use of this methodology in conjunction with a submodelling approach [100] in which a specialized mesh employing the quarter-point technique is generated in the surrounding of the crack front to accurately computed the SIF.

Supported by the advances in computational power achieved over recent decades, the study of three-dimensional fatigue problems has gained increasing traction within the FM community. Although based on some of the same principles underlying the two-dimensional approaches previously discussed, the application of remeshing techniques to three dimensional problems may present significantly greater complexity. As a consequence, this topic has evolved into a research field of its own and will not be addressed in detail here. Thus, for a broader overview of existing three-dimensional remeshing techniques, the reader is referred to [101].

3.1.4. Crack Modeling Using the Extended Finite Element Method

As mentioned earlier, FM simulations using the standard FEM usually rely on either models with fixed meshes that are locally refined to properly represent the crack or on remeshing techniques. Although these approaches can handle a wide range FM problems, each has notable shortcomings: the first, requires the a priori knowledge of the crack path and is plagued by mesh dependency, while the second can lead to significant computational cost in large-scale fracture problems, as at least part of the mesh must be recreated at each crack increment. To address these limitations, Belytschko and Black [102] proposed a mesh-independent method that requires little remeshing. Then this approach was further developed by Moës, Dolbow and Belytschko [103] to eliminate the need for mesh recreation. This method eventually became widely known as extended finite element method (XFEM). A comprehensive review of it can be found in [104]. In this approach, the standard finite element formulation is enriched by additional nodal degrees of freedom to enable the representation of discontinuities with non-conforming meshes, as shown in Figure 13. As such, the XFEM has been successfully applied not only to solve crack growth problems but also to the modeling of material interfaces, as those found in two-phase flow and fluid-structure interaction cases [104].

In practical terms, only the nodes of the elements intersected by the discontinuity are enriched to represent the displacement, strain and stress fields in this region. Accordingly, in the context of crack growth problems, the displacement field can be approximated by

$${\mathbf{u}}^h\left(\mathbf{x}\right)=\underset{\begin{array}{c}\mathrm{standard}\mathrm{FE}\\ \mathrm{approximation}\end{array}}{\underbrace{\sum _{i\in I}{N}_i\left(\mathbf{x}\right){\mathbf{u}}_i}}+\underset{\mathrm{enrichment}}{\underbrace{\sum _{i\in I^{*}}{N}_i^{*}\left(\mathbf{x}\right)·\left[\chi \left(\mathbf{x}\right)-\chi \left({\mathbf{x}}_i\right)\right]{\mathbf{q}}_i}}$$

(20)

where $N_i\left(\mathbf{x}\right)$ are the standard FEM shape functions and ${\mathbf{u}}_i$ the corresponding nodal degrees of freedom. $N_i^{*}\left(\mathbf{x}\right)$ and $q_i$, in turn, are shape functions fulfilling the partition unity property – often chosen to be equal to $N_i\left(\mathbf{x}\right)$ – and the associated additional degrees of freedom of the enriched nodes, which belong to the subset $I^{*}\subset I$ of all nodes. Finally, $\chi \left(\mathbf{x}\right)$ is the enriching function while $\chi \left({\mathbf{x}}_i\right)$ acts as a shifting function that deactivates ${\mathbf{q}}_i$ at the respective nodes, ensuring displacement compatibility between adjacent elements. Under the LEFM assumption, this enriching function can be expressed in polar coordinates as follows [102]

$$\chi (r,\theta )=\left\{\sqrt{r}sin\left({\displaystyle \frac{\theta }{2}}\right),\sqrt{r}cos\left({\displaystyle \frac{\theta }{2}}\right),\sqrt{r}sin\left({\displaystyle \frac{\theta }{2}}\right)sin\left(\theta \right),\sqrt{r}cos\left({\displaystyle \frac{\theta }{2}}\right)sin\left(\theta \right)\right\},$$

(21)

where r is the radial distance from the crack tip, while $\theta$ represents the angular coordinate, with $\theta$ = 0 defined along the direction of the crack extension. Alternative enriching functions for the cases involving significant plasticity occurring at the crack tip have been introduced by Elguedj et al. [106]. Complementary to these enrichment functions, Moës et al. [103] have proposed an additional enrichment based on the Heaviside function to effectively represent long, severely curved or three-dimensional cracks – eliminating the need for even minimal remeshing, as required in the original approach conceived by Belytschko and Black [102].

Upon the determination of the displacement, strain and stress fields, the SIF is usually extracted in a XFEM framework using the domain form of the interaction integral approach [102,103,106,107,108]. Energy release methods are generally preferred over the displacement-based approach to this end, because obtaining the displacement field at the crack front is less straightforward in a XFEM setting. The crack advance direction, in turn, is typically determined through the maximum tangential stress criterion [102,103,108,109].

Finally, a mathematical formulation is required to implicitly represent the crack due to the mesh non-conformity characteristic of the XFEM. The level set method is commonly employed for this purpose. Accordingly, two level set curves are used to describe a crack, $\mathsf{\Gamma }\left(t\right)$, of arbitrary shape.

The first is a signed distance function defining the crack geometry and position. It is typically denoted by

$$\varphi (\mathbf{x},t)=\pm mi{n}_{{\mathbf{x}}_{\mathsf{\Gamma }}\in \mathsf{\Gamma }\left(t\right)}\parallel \mathbf{x}-{\mathbf{x}}_{\mathsf{\Gamma }}\parallel$$

(22)

where ${\mathbf{x}}_{\mathsf{\Gamma }}$ is the closest point on the crack interface $\mathsf{\Gamma }\left(t\right)$ to $\mathbf{x}$. The sign of the minimum distance depends on which side of the interface $\mathbf{x}$ lies.

The second is also a signed distance function, used to locate the crack front. It is orthogonal to $\mathsf{\Gamma }\left(t\right)$ and equal to zero at the crack fronts. Thus, it reads

$$\gamma (\mathbf{x},t)=(\mathbf{x}-{\mathbf{x}}_i)·{\widehat{\mathbf{t}}}_i$$

(23)

where ${\widehat{\mathbf{t}}}_i$ is a unit vector tangent to the crack tip ${\mathbf{x}}_i$ at the time t.

With the two level set functions defined, the crack geometry at time t can be implicitly described as follows

$$\mathsf{\Gamma }\left(t\right)=\left\{\mathbf{x}:\varphi (\mathbf{x},t)=0\wedge \gamma (\mathbf{x},t)\le 0\right\},$$

(24)

such that $\mathsf{\Gamma }\left(t\right)$ is incrementally updated according to the crack growth length, $\Delta a$, and its propagation direction. The fast marching method proposed by Sukumar et al. [107] may be employed for updating planar cracks in three dimension problems. An overview of different approaches used to perform this update is provided in [110]. Figure 14 schematically illustrates the crack representation procedure using the level set method.

3.2. Models Based on the Continuum Damage Mechanics

Emerging from the broader field of continuum mechanics, the CDM has been primarily conceived to describe the effects of distributed microdefects, such as microcracks and microvoids, on the progressive deterioration of materials. As such, it provides a continuum-level framework that captures the evolution of these diffuse microdefects up to the initiation of macrocracks and their impact on the material response – typically manifested as a reduction in one or more material properties, including strength, stiffness, toughness or residual life. Figure 15 schematically illustrates the domain originally covered by the CDM theory, emphasizing its complementarity with the FM.

Building on the popularization of numerical methods, the CDM domain of applicability has been gradually extended to model not only the diffuse damage phase but also the localization stage, characterized by macrocrack propagation leading to the material complete failure. One of the earliest contributions in this direction was made by Lemaitre [112] who discussed the notion of structural crack in the context of the local approach to fracture. This approach assumes that an elementary volume of material governed by a certain constitutive response (e.g., elastic, plastic, viscoplastic, etc.) may be damaged to a critical state that triggers the initiation of a local macrocrack. Accordingly, the successive fulfillment of this initiation condition in adjacent elementary volumes gives rise to a structural crack.

Despite its conceptual importance, the local CDM theory leads to ill-posed boundary value problems (BVP) when modeling the strain-softening behavior typically observed in damaged materials. This ill-posedness creates difficulties in solving these problems with computational methods such as FEM. In fact, Bažant and Belytschko [113] proved that, in a local strain-softening formulation, damage localizes in a region of zero volume – collapsing to a point in 1D problems, to a line in 2D and to a surface in 3D. They further demonstrated that the corresponding dissipated energy tends to zero under increasing mesh refinement, which lacks physical meaning.

To address these inconsistencies, regularization techniques – most notably non-local formulations [114,115] and gradient-enhanced [116,117] approaches – have been employed, ensuring that the dissipated energy remains finite and consistent with the fracture energy of Griffith’s theory. In these regularized settings, a crack is no longer represented as a sharp discontinuity but rather through a strain-softening response that is either spatially averaged over a finite neighborhood (nonlocal integral) or regularized via higher-order differential operators (gradient damage), both delimited by a characteristic length scale, $l_c$, (see Figure 16). This localized region is commonly interpreted as the fracture process zone, while the crack itself may be idealized as the locus where the damage reaches a certain critical value, $d_c$ [118]. Here it is important to point out that, in this text, the lowercase letter d is used in regularized formulations to denote the damage variable associated with the initiation and also the propagation of macrocracks leading to complete material failure. Conversely, the uppercase D is employed in classical approaches to describe the diffuse damage phase, which extends only up to the macrocrack onset.

Since these approaches introduce a length-scale parameter at the continuum level, in a FEM setting, the discretization scheme needs to have enough resolution to resolve the strain localization occurring in the region delimited by this parameter. This often requires highly refined meshes in this zone, leading to computationally costly FE models.

As an alternative, discrete-level regularization has been widely adopted through smeared crack formulations. Based on the crack band theory (CBT) introduced by Bažant and Oh [119], these approaches smear the fracture process zone across the intersecting elements. Accordingly, the effective bandwidth of the localization zone is linked to the characteristic length of these element, while consistency with FM predictions is ensured by enforcing energy equivalence between the constitutive softening law and Griffith’s fracture energy.

This direct link between softening laws and fracture energy highlights that the post-peak degradation response plays a central role in regularized formulations. As a consequence, the choice of a softening law tailored to the material behavior is a key requirement for consistent energy dissipation. For instance, in the case of semi-brittle and brittle materials, two widely adopted choices are the linear and exponential softening curves, as schematically shown in Figure 17. These laws define how stress decays after the peak load and are typically calibrated such that, for a given characteristic length, the energy dissipated within the fracture process zone is consistent with the fracture energy.

When properly calibrated, both continuum-regularized and smeared models are able to reproduce the well-known Bažant size effect. This capability stems from their energetic consistency with fracture mechanics and from the strength criterion naturally embedded in the constitutive law. Accordingly, this effect can be expressed as

$$S_N=B{S}_u{\left(1+{\displaystyle \frac{L}{L_0}}\right)}^{-1/2},$$

(25)

where $S_N$ denotes the nominal strength, $S_u$ the material ultimate strength, L the characteristic structural size, and $L_0$ a reference structural size. The parameter B is a dimensionless constant that depends on the structure geometry and loading conditions.

Equation (25) evidences that, for small characteristic lengths ($L\ll L_0$), the nominal strength approaches the material ultimate strength, and failure follows the predictions of the limit analysis. Conversely, for large structures ($L\gg L_0$), a significant reduction in the nominal strength is observed, with the collapse governed by the material’s fracture toughness, as predicted by the LEFM. For intermediate structural sizes, L, the fracture behavior is influenced by both mechanisms. Figure 18 schematically describes the Bažant size effect in a log–log scale.

Finally, the numerical reproduction of Grégoire’s three-point bending tests [121] is presented in Figure 19 to illustrate the capability of smeared crack models to capture the Bažant size effect.

3.2.1. Early Approaches and Basic Concepts

The foundation of the CDM was laid by Kachanov [123] in 1958 through his pioneering work on creep in metals. In this context, he introduced a scalar damage variable $\varphi$ to represent the growth of voids and the corresponding load-carrying net area of tensile bars over time. Accordingly, $\varphi$ was defined to range from 1 in the material undamaged state to 0 in the fully damaged condition. Although Kachanov did not provide a detailed discussion on the physical description of $\varphi$, this variable can be interpreted as the ratio of the load-carrying net area to the corresponding nominal one, i.e.,

$$\varphi ={\displaystyle \frac{A}{A_n}}.$$

(26)

This concept was later further developed by Rabotnov [124], who introduced the alternative form

$$D=1-\varphi ,$$

(27)

for the damage variable, which varies from 0 to 1 instead and accounts for the reduction in the load-carrying net area. In addition, he incorporated to the original theory a coupling that accounts for the influence of the damage state on the creep rate, enabling the description of the tertiary creep observed in experiments under constant pressure.

The change of load-carrying net area described by these damage variables naturally leads to a magnification of the stresses in the corresponding area. This gives rise to the fundamental CDM concept of effective stress, which is expressed as

$${\mathsf{\sigma }}_0={\displaystyle \frac{\mathsf{\sigma }}{\varphi }}={\displaystyle \frac{\mathsf{\sigma }}{1-D}}.$$

(28)

This effective stress concept is usually accompanied by a corresponding hypothesis of equivalent strains, which ensures a consistent relationship between the effective (undamaged) and real (damaged) spaces – i.e., same strain, $\epsilon$, in both configurations –, as schematically represented in Figure 20.

Following these pioneering works, the CDM became a very active research field from the early 1970s, with mathematical formulations being proposed to model not only creep damage, but also elastic [126,127,128,129,130,131], elastic-plastic [120,126,129,132,133,134,135,136,137,138,139,140], fatigue damage [120,126,131,132,141,142,143,144,145,146,147,148,149,150,151,152] and corrosion [126,129,132,136,153], to name a few. The classification presented in Table 2 summarizes some of these damage mechanisms and their associated microstructural characteristics, as originally compiled by Murakami [154] and extended here to incorporate other relevant as well as more recent works.

Most of the early formulations developed to model damage degradation processes employed scalar variables and assumed isotropic evolution laws [123,134,168,169,170]. These models adequately idealize damage processes characterized by randomly oriented microdefects or by homogeneous distributions of spherical voids. However, they fail when applied to describe microdefects with significant geometrical directionality or spherical voids with markedly orientation in their spacial distribution. To address these limitations, subsequent developments introduced tensor-valued damage variables to represent anisotropic damage states, as proposed by Murakami and Ohno [171] through the second rank symmetric damage tensor

$$\mathsf{\Omega }={\displaystyle \sum _{i=1}^3}{\mathsf{\Omega }}_i{\mathbf{n}}_i\otimes {\mathbf{n}}_i,$$

(29)

where ${\mathsf{\Omega }}_i$ and ${\mathbf{n}}_i$ are the principal values and the respective principal directions of $\mathsf{\Omega }$, respectively. Accordingly, the corresponding effective stress tensor is given by

$${\mathsf{\sigma }}_0={\displaystyle \frac{1}{2}}\left[{\left(\mathit{I}-\mathsf{\Omega }\right)}^{-1}\mathsf{\sigma }+\mathsf{\sigma }{\left(\mathit{I}-\mathsf{\Omega }\right)}^{-1}\right].$$

(30)

Many other second rank [172,173,174,175,176] as well as higher rank [177,178] tensorial models have been proposed to address these directional effects. In addition to these tensor-based approaches, some models relying on vetorial variables were also proposed, primarily to describe damage states characterized by in-plane microdefects distributions [128,179,180,181,182].

3.2.2. Particular Models

Among the various approaches that followed the initial concepts on creep problems, two major formulations emerged in the context of ductile fracture and became the backbone of subsequent CDM developments. The first, introduced by Gurson [134] and later extended by Tvergaard and Needleman [183], is based on micromechanical considerations. It employs the void volume fraction f as a physically measurable scalar damage variable (void volume fraction), which evolves from $f=0$ in the undamaged state to a critical value $f_c$ at the void coalescence onset. Within this framework, the yield surface is modified to account for the influence of void growth on the material degradation, assuming the form

$$\mathbb{F}\left(\mathsf{\sigma },f\right)={\left({\displaystyle \frac{{\sigma }_{eq}}{S_y}}\right)}^2+2q_{1}fcosh\left({\displaystyle \frac{3q_2{\sigma }_h}{2S_y}}\right)-\left(1+q_3{f}^2\right)=0,$$

(31)

where ${\sigma }_{eq}$ denotes the von Mises equivalent stress, ${\sigma }_h$ the mean (hydrostatic) stress and $S_y$ the current yield stress of the matrix material. The parameters $q_1$, $q_2$, and $q_3$, in turn, were introduced by Tvergaard [184] to improve the agreement with results of numerical cell models. In this formulation, the material elastic properties remain unaffected, since the damage variable, f, governing the void evolution, is coupled solely with the yield surface.

Later, the modified volume fraction

$$f^{*}=\left\{\begin{array}{cc}f,\hfill & \mathrm{for}f\le f_c,\hfill \\ f_c+{\displaystyle \frac{f_u^{*}-f_c}{f_f-f_c}}(f-f_c),\hfill & \mathrm{for}f>f_c,\hfill \end{array}\right.$$

(32)

was introduced by Tvergaard and Needleman [183] to account for the void coalescence phase and final failure. Here $f_f$ is the void volume fraction at the final fracture and $f_u^{*}=1/q_1$ represents the ultimate void volume fraction.

The second approach, due to Lemaitre [135,185], relies on a thermodynamically consistent phenomenological framework in which damage is coupled with plasticity through the effective stress concept described in Equation (28). In this formulation, a scalar damage variable $D\in \left[0,1\right]$ is defined to represent the progressive reduction of the load-carrying net area, while the yield function is given by

$$\mathbb{F}\left({\mathsf{\sigma }}_0,\kappa \right)={\sigma }_{eq}-S_y\left(\kappa \right)=0,$$

(33)

where ${\sigma }_{eq}$ represent the von Mises equivalent stress, $S_y\left(\kappa \right)$ the current yield stress and $\kappa$ the isotropic hardening variables.

Following the general principles of standard constitutive models based on internal variables, the energy density release rate (thermodynamic force conjugate to damage) is defined as

$$Y=-{\displaystyle \frac{\partial {\mathsf{\Psi }}_e}{\partial D}},$$

(34)

where ${\mathsf{\Psi }}_e$ represent the elastic strain energy density. This ensure compliance with the Clausius–Duhem inequality and guarantees non-negative dissipation, i.e.,

$$Y\dot{D}\ge 0.$$

(35)

Alongside this requirement, the following damage evolution law

$$\dot{D}={\left({\displaystyle \frac{Y}{S}}\right)}^s{\dot{\overline{\epsilon }}}_p$$

(36)

is postulated. Here S and s are the damage strength and exponent (material parameters), whereas ${\dot{\overline{\epsilon }}}_p$ is the accumulated plastic strain rate. Finally, the damage internal variable affects the material elastic property such that

$$E=\left(1-D\right)·E_0,$$

(37)

where E and $E_0$ represent the elastic modulus in the damage and undamaged states, respectively.

• Other relevant plastic-damage models

Building on this thermodynamic foundation, Simo and Ju [137,138] developed strain- and stress-based CDM formulations that provided a rigorous variational framework for integrating plasticity and damage within a FEM setting. Their work showed that strain- and stress-driven descriptions could be derived within a single thermodynamically consistent structure, while also introducing consistent tangent operators essential for robust finite element implementation. Other important contributions include the energy-based coupled elastoplastic-damage theory by Ju [186], which reinforced the thermodynamic structure of CDM, the plastic–damage model by Lubliner et al. [139] and the coupled thermo-plastic–damage formulation by Luccioni et al. [140], both of which extended Lemaitre-type models to quasi-brittle materials.

• Pure-damage models

In contrast to these plastic-damage formulations, Krajcinovic [128,181] demonstrated that Lemaitre-type models can be particularized, by neglecting plasticity-related effects, to describe the response of semi-brittle and brittle materials, such as concrete, rocks and glass. More specifically, he showed that the theory remains thermodynamically consistent when the plasticity-related terms are eliminated and the inelastic response becomes primarily governed by distributed microcracks rather than plastic deformation. This led to constitutive models in which the damage variable directly controls the material stiffness degradation and crack initiation.

Following this line, Oliver et al. [130] proposed an isotropic damage model for concrete that, while building on the variational structure of Simo and Ju [137,138] and incorporating aspects of Lubliner’s formulation [139], also neglects plasticity effects. In their approach, concrete degradation and crack localization are represented exclusively through the evolution of the damage internal variable, further consolidating the role of pure-damage formulations as an alternative to coupled plastic–damage approaches for modeling the response of semi-brittle and brittle materials.

3.2.3. Extension to Fatigue

One of the earliest fatigue models formulated within the CDM framework was proposed by Chaboche in 1974 [141]. Originally derived for uniaxial stress states, this approach was later extended to multiaxial loading conditions [187]. In this formulation, the damage evolution, expressed in terms of the number of cycles,

$$dD=f\left({\sigma }_{max},{\sigma }_m,D\right)dN,$$

(38)

was defined to describe the progressive deterioration of the material resulting from the application of cyclic loads. Here $f\left({\sigma }_{max},{\sigma }_m,D\right)$ is a function depending on the maximum stress, ${\sigma }_{max}$, the mean stress, ${\sigma }_m$ and the previously accumulated damage, D. The integration of Equation (38) yields a non-linear damage accumulation law, which can be incorporated into the CDM framework either through total-life approaches or the effective stress concept of Equation (28). The respective constants of integration, in turn, are obtained from conventional material data, including S-N curves.

Following a similar line, Lemaitre and Plumtree [163] proposed a formulation to model creep-fatigue interactions. Their approach introduced strain-based damage evolution equations to describe pure fatigue processes (time-independent), fatigue-creep (time-dependent), and cases in which both mechanisms act either simultaneously, i.e., in a combined fatigue-sustained loading setting, or in sequence, as in separate phases of fatigue and creep.

• Chaboche-type thermodynamically consistent pure-damage models for fatigue

Although grounded in the CDM principles, these two classical approaches were not explicitly derived within a rigorous thermodynamically consistent framework. This step was later taken by Paas et al. [143], Peerling [188] and Xiao et al. [145] (the latter focusing on the high cycle domain), who formulated rate-form damage evolution equations, with the cycle-wise increment obtained through time integration over a cycle period, under the assumption that damage remains constant within each cycle. As such, damage is primarily defined as an internal variable, making these models naturally compatible with the principles of the thermodynamics of irreversible processes. More recently, Fu et al. [152] presented a formulation that follows the same rational but employs a tensorial damage internal variable to capture the anisotropy effects of additively manufactured Ti-6Al-4V coupon specimens in the high cycle and very high cycle regimes.

• Lemaitre-type plastic-damage models for fatigue

Based on pure-damage formulations, these approaches do not account for plasticity effects that are characteristic of metal fatigue in the low cycle domain. Thus, fatigue in this regime is often addressed by plastic-damage formulations such as the one presented by Dufailly and Lemaitre [144], who adapted the Lemaitre’s type model to describe very low cycle thermomechanical fatigue. The same damage law was later associated with viscoplasticity using kinematic hardening to describe creep-fatigue interaction in the low cycle domain [167]. Building on this foundation, Lemaitre [146] proposed in 1999 a two scale plastic-damage formulation to also model fatigue in the high cycle domain. In this approach, plasticity and damage are accounted for at the microscale while the mesoescale is considered to be elastic. More recently, Liu et al. [150] extended this formulation by incorporating a probabilistic method to account for the scatter typically found in the HCF life of metals.

Table 3 summarizes the damage evolution laws of some of the fatigue formulations previously mentioned.

• Chaboche-type plastic-damage models for fatigue

Another class of relevant fatigue plastic-damage models builds on plasticity formulations, most often employing kinematic hardening, coupled with the Chaboche-type damage evolution law. In these approaches, the coupling is established through the effective stress concept, as described by Lemaitre and Chaboche [131]. Models of this class have already been applied to evaluate fatigue damage in both low and high cycle regimes under a variety of strain controlled load paths [189], as well as to assess fatigue in weld joints considering the inherent porosity [190], metallic bolted joints in the high cycle domain [191] and the influence of additive manufacturing (AM) on the fatigue performance of aerospace commonly-used alloys [192]. More recently, such models have also been integrated with machine learning (ML) approaches to analyze fatigue in the AM context [149,151].

• Residual strength-based isotropic thermo-plastic-damage model for fatigue

As discussed earlier, Chaboche- and Lamaitre-type fatigue models rely on damage evolution laws that are either phenomenologically defined in term of cycle increments or derived from rate equations that are subsequently integrated over a cycle, yielding cycle-wise incremental expression. Framing fatigue modeling from a different perspective, Oller and co-workers [148] proposed a thermodynamically consistent fatigue model – based on the isotropic thermo-plastic-damage formulation presented by Luccioni et al. [140] – which does not introduce an explicitly cycle-wise damage accumulation law but instead retains a rate-form structure. The central idea is that fatigue in metals manifests as a progressive reduction in the material strength. This strength loss then induces the material inelastic response, which can be interpreted as the coalescence of microdefects into macrodefects (voids and cracks) and their subsequent propagation. Within this framework, a fatigue state function depending on the number of cycles, $N_c$, stress ratio, R, maximum stress, ${\sigma }_{max}$, and temperature $\theta$, is introduced to describe the evolution of the residual strength which is linked to the yield and damage surfaces, governed by standard plastic and damage internal variables. Accordingly, the yield surface in this formulation may be equivalently expressed as

$${\mathbb{F}}^P\left({\sigma }_{ij},{\mathcal{K}}^P\right)=f^P\left({\sigma }_{ij}\right)-{\mathcal{K}}^P({\sigma }_{ij},{\alpha }^p)·f_{red}(N_c,R,{\sigma }_{max},\theta )=0,$$

(39)

and

$${\mathbb{F}}^P\left({\sigma }_{ij},{\mathcal{K}}^P\right)={\displaystyle \frac{f^P\left({\sigma }_{ij}\right)}{f_{red}(N_c,R,{\sigma }_{max},\theta )}}-{\mathcal{K}}^P({\sigma }_{ij},{\alpha }^p)=0.$$

(40)

Analogously, the damage surface can by written as

$${\mathbb{F}}^D\left({\sigma }_{ij},{\mathcal{K}}^D\right)=f^D\left({\sigma }_{ij}\right)-{\mathcal{K}}^D({\sigma }_{ij},d)·f_{red}(N_c,R,{\sigma }_{max},\theta )=0,$$

(41)

and

$${\mathbb{F}}^D\left({\sigma }_{ij},{\mathcal{K}}^D\right)={\displaystyle \frac{f^D\left({\sigma }_{ij}\right)}{f_{red}(N_c,R,{\sigma }_{max},\theta )}}-{\mathcal{K}}^D({\sigma }_{ij},d)=0,$$

(42)

where $f^P\left({\sigma }_{ij}\right)$ and $f^D\left({\sigma }_{ij}\right)$ denote stress norms related with the plasticity and damage processes, ${\mathcal{K}}^P({\sigma }_{ij},{\alpha }^p)$ and ${\mathcal{K}}^D({\sigma }_{ij},d)$ the corresponding thresholds functions, ${\alpha }^p$ and d the plastic and damage internal variables. As previously mention, d here is expressed in lowercase letter to emphasize that the damage internal variable in this model describes the macrocrack initiation and propagation leading to the material complete failure, i.e., the damage localization phase. In addition to that, $f_{red}(N_c,R,{\sigma }_{max},\theta )$ in Equations (39)–(42) represents the fatigue state function. Benefiting from the extensive experimental data available in the literature, this function is postulated in a way that can be fully calibrated from the information provided by S-N data.

Following this, the coupling between damage and plasticity processes is enforced through the total dissipated energy. Accordingly, the normalized energy release rate reads

$$\dot{q}={\dot{\mathsf{\Xi }}}_m·\mathfrak{R}\left({\sigma }_{ij}\right)=\left({\dot{\mathsf{\Xi }}}^P+{\dot{\mathsf{\Xi }}}^D\right)\left[{\displaystyle \frac{r\left({\sigma }_{ij}\right)}{g_f}}+{\displaystyle \frac{1-r\left({\sigma }_{ij}\right)}{g_c}}\right],$$

(43)

where

$$r\left({\sigma }_{ij}\right)=\left\{\begin{array}{cc}1,\hfill & \mathrm{for}\mathrm{pure}\mathrm{tension}\hfill \\ 0,\hfill & \mathrm{for}\mathrm{pure}\mathrm{compression}.\hfill \end{array}\right.$$

(44)

In addition, ${\dot{\mathsf{\Xi }}}^P$ and ${\dot{\mathsf{\Xi }}}^D$ represent the rate of plastic and damage energy dissipation, whereas $g_f$ and $g_c$ define the maximum energy that can be dissipated in tension and compression at a given material point at the end of the inelastic process.

This formulation has already been applied to analyze fatigue in hydraulic cylinders [193] and has also been particularized, under the general assumptions of infinitesimal strains and negligible plasticity effects, to model a stiffness-based rapid fatigue test [194], HCF of railway structures [195], an automotive suspension component [196], as well as thermomechanical fatigue in the same domain [197]. Under this simplification, the underlying approach assumes the form of the pure-damage model proposed by Oliver et al. [130] which is then modified by the fatigue state function. The explicit damage-integration scheme enabled by Oliver’s model, combined with the use of cycle-jump algorithms [199,200], makes this formulation highly efficient for HCF simulations when embedded within a FEM framework. More recently, this constitutive model has also been extended to incorporate manufaturing-induced effects in the HCF modeling of metals [198].

3.3. Phase-field Models

Dating back to the late 1990s, phase-field (PF) models for fracture originated in the context of brittle materials (elastic solids) and have since received increasing attention due to their ability to naturally capture crack initiation, propagation, and complex patterns such as merging and branching, without the need for any ad-hoc techniques. In these approaches, an additional scalar field – the phase field – is introduced to represent the crack in a regularized manner within a finite region of the domain. Analogously to the internal variable used in regularized damage models, the order parameter associated with the phase field describes the transition between the intact and fully broken states of the material in a continuous way. In addition to this characteristic, PF models inherently include a gradient term in their formulation and therefore exhibit many features in common with gradient-enhanced damage models [11].

Since the early formulations, PF models for fracture have been developed in parallel and largely independently within the physics and mechanics communities. With particular focus on dynamic fracture problems, PF approaches emerging from the physics community (e.g., [201,202,203]) are grounded in the Landau–Ginzburg phase transition theory [204]. However, these models do not explicitly include the fracture energy in their formulation and are, therefore, not directly linked to Griffith’s theory. Moreover, they do not explicitly incorporate a characteristic length-scale parameter governing the size of the fracture process zone.

On the other hand, PF models developed within mechanics community stem from the variational formulation of brittle fracture proposed by Francfort and Marigo [205]. In this framework, the processes of crack initiation, propagation, and branching are formulated as a minimization problem of the total energy functional

$$\mathbb{E}(\mathbf{u},\mathsf{\Gamma })={\int }_{\mathsf{\Omega }}{\mathsf{\Psi }}_0\left(\mathsf{\epsilon }\right)\mathrm{d}V+{\int }_{\mathsf{\Gamma }}{G}_c\mathrm{d}S,$$

(45)

where ${\mathsf{\Psi }}_0$ denotes the strain energy density of the undamaged solid occupying the domain $\mathsf{\Omega }$, $G_c$ is the material fracture energy, $\mathsf{\Gamma }$ represents the crack surface, and $\mathbf{u}$ is the displacement field. From a computational standpoint, however, the direct evaluation of the surface integral over $\mathsf{\Gamma }$ is not feasible, since the crack geometry and its evolution are generally unknown a priori. To overcome this limitation, Bourdin [206] introduced a regularized formulation of the total energy functional, which in its general form can be expressed as

$${\mathsf{\Pi }}_{\ell }={\int }_{\mathsf{\Omega }}g\left(\varphi \right){\mathsf{\Psi }}_0\left(\mathsf{\epsilon }\right)\mathrm{d}V+{\int }_{\mathsf{\Omega }}{G}_c\gamma (\varphi ,\ell )\mathrm{d}V,$$

(46)

where $\varphi$ is the phase-field order parameter – taking the value 0 for the intact material and 1 for the fully broken state –, ℓ denotes the length-scale parameter that governs the fracture process zone width, $g\left(\varphi \right)$ is a degradation function and $\gamma (\varphi ,\ell )$ represents the so-called crack surface density function. Alternatively, PF formulations grounded on thermodynamic principles [207,208] and employing a higher-order regularization [209] have also been proposed within the mechanics community. Figure 21 schematically compares the sharp crack discontinuity with its corresponding PF approximation.

Following the regularized energy functional of Equation (46), the constitutive stress response is obtained by differentiation of the strain energy density with respect to the strain. This yields the stress tensor

$$\mathsf{\sigma }={\displaystyle \frac{\partial \mathsf{\Psi }\left(\mathsf{\epsilon }\right)}{\partial \mathsf{\epsilon }}}.$$

(47)

For linear elasticity, the elastic strain energy density is

$${\mathsf{\Psi }}_e\left(\mathsf{\epsilon }\right)=g\left(\varphi \right){\mathsf{\Psi }}_0\left(\mathsf{\epsilon }\right)=g\left(\varphi \right){\displaystyle \frac{1}{2}}\mathsf{\epsilon }:{\mathbb{C}}_0:\mathsf{\epsilon },$$

(48)

where ${\mathbb{C}}_0$ represents the constitutive tensor of the material in the undamaged state. Thus, the constitutive law of Equation (47) reduces to

$$\mathsf{\sigma }=g\left(\varphi \right)\left({\mathbb{C}}_0:\mathsf{\epsilon }\right),$$

(49)

so that the elastic stiffness progressively degrades as the phase-field order parameter evolves from $\varphi =0$ (undamaged) to $\varphi =1$ (fully broken).

When it comes to the crack surface density and degradation functions, it can be demonstrated that the regularized formulation converges to the total energy functional of Equation (45) (i.e., Griffith’s solution) as $\ell \to 0^{+}$ [210], provided that these functions are appropriately defined. Among the various constitutive choices available for the crack surface density function, the formulations proposed by Ambrosio and Tortorelli (AT) [211] – originally developed as an approximation of the Mumford–Shah functional for image segmentation [212] – are the most commonly used. In their general form, they read

$${\int }_{\mathsf{\Gamma }}{G}_c\mathrm{d}S\approx {\int }_{\mathsf{\Omega }}\gamma \left(\varphi \right)\mathrm{d}V={\int }_{\mathsf{\Omega }}{\displaystyle \frac{1}{4c_w\ell }}\left(w\left(\varphi \right)+{\ell }^2{|\nabla \varphi |}^2\right)\mathrm{d}V,$$

(50)

where $w\left(\varphi \right)={\varphi }^2$ with $c_w=1/2$ corresponds to the so-called standard or AT2 model, while $w\left(\varphi \right)=\varphi$ with $c_w=2/3$ defines the AT1. Accordingly, the following quadratic degradation function is used with both

$$g\left(\varphi \right)={(1-\varphi )}^2.$$

(51)

For a homogeneous 1D response ($\nabla \varphi =0$), these surface density functions lead to distinct stress–strain behaviors. The AT1 model exhibits a linear elastic behavior up to a critical stress, whereas the AT2 model shows an earlier deviation from linearity. In addition, the AT1 model displays a steeper softening decay, as shown in Figure 22. As a result, the corresponding critical stresses are

$${\sigma }_c^{\mathrm{AT}1}=\sqrt{{\displaystyle \frac{3E{G}_c}{8\ell }}},$$

(52)

and

$${\sigma }_c^{\mathrm{AT}2}=\sqrt{{\displaystyle \frac{27E{G}_c}{256\ell }}}.$$

(53)

Equations (51) and (52) show that, as $\ell \to 0^{+}$, the critical stress tends to infinity, which is consistent with the $\mathsf{\Gamma }$-convergence argument. On the other hand, by considering $\ell >0$, a finite material strength – i.e., an upper-bounded critical stress – is introduced. In this sense, ℓ acts as an intrinsic material property (i.e., internal length scale) alongside the Young’s modulus, E, and the fracture energy, $G_c$. As a direct consequence, these PF models combine stress- and toughness-based fracture criteria, enabling them to naturally capture size effects, as illustrated in Figure 23 for the example of the single-edge notched specimen presented in [213].

Up to this point, the PF formulation presented here assumes an isotropic degradation of the strain energy density, meaning that the crack-driving mechanism does not distinguish between tension and compression. This assumption, however, may lead to nonphysical crack growth under compressive loading, since compressed regions may still contribute to the degradation of the material. Moreover, the contact between crack faces under closure conditions is not naturally enforced in this type of formulation.

To overcome these limitations, several tension–compression splits of the strain energy density have been proposed, giving rise to the so-called anisotropic PF fracture models. One of the earliest split schemes was introduced by Amor et al. [214], who decomposed ${\mathsf{\Psi }}_0$ into volumetric and deviatoric contributions. In this approach, the positive (active) and negative (inactive) parts of the strain energy density are defined as

$$\begin{array}{cc}\hfill {\mathsf{\Psi }}_0^{+}\left(\mathsf{\epsilon }\right)& ={\displaystyle \frac{1}{2}}\left(\lambda +{\displaystyle \frac{2}{3}}\mu \right){\langle tr\left(\mathsf{\epsilon }\right)\rangle }_{+}^2+\mu \left({\mathsf{\epsilon }}^{\mathrm{dev}}:{\mathsf{\epsilon }}^{\mathrm{dev}}\right),\hfill \\ \hfill {\mathsf{\Psi }}_0^{-}\left(\mathsf{\epsilon }\right)& ={\displaystyle \frac{1}{2}}\left(\lambda +{\displaystyle \frac{2}{3}}\mu \right){\langle tr\left(\mathsf{\epsilon }\right)\rangle }_{-}^2,\hfill \end{array}$$

(54)

where $\lambda$ and $\mu$ are the Lamé constants of linear elasticity and ${\langle •\rangle }_{\pm }$ denote the Macaulay brackets, defined as ${\langle \square \rangle }_{\pm }=\frac{1}{2}(\square \pm |\square \left|\right)$. In addition, the deviatoric strain is given by

$${\mathsf{\epsilon }}^{\mathrm{dev}}=\mathsf{\epsilon }-{\displaystyle \frac{1}{3}}tr\left(\mathsf{\epsilon }\right)\mathbf{I},$$

(55)

with $\mathbf{I}$ being the identity matrix. While this split provides physically meaningful results for most loading conditions, it still permits spurious crack growth under states in which all three principal strains are negative.

As an alternative to this approach, Miehe [208] proposed a tension–compression split based on the spectral decomposition of the elastic strain tensor. Accordingly, the positive and negative parts of the strain energy density read

$${\mathsf{\Psi }}_0^{\pm }\left(\mathsf{\epsilon }\right)={\displaystyle \frac{1}{2}}\lambda {\langle tr\left(\mathsf{\epsilon }\right)\rangle }_{\pm }^2+\mu tr\left({\mathsf{\epsilon }}_{\pm }^2\right),$$

(56)

where the strain tensor is spectrally decomposed as

$${\mathsf{\epsilon }}_{\pm }={\displaystyle \sum _{i=1}^3}{\langle {\epsilon }_i\rangle }_{\pm }{\mathbf{n}}_i\otimes {\mathbf{n}}_i.$$

(57)

Here ${\epsilon }_i$ and ${\mathbf{n}}_i$ denote the principal strains and their associated eigenvectors, respectively. It is noteworthy that this split leads to a strongly nonlinear stress–strain relationship, typically resulting in a higher computational cost when compared with the volumetric–deviatoric split of Amor et al. [214].

Another tension–compression split employed in PF fracture models is the no-tension scheme proposed by Freddi et al. [215]. In this formulation, the positive and negative parts of the strain energy density render

$${\mathsf{\Psi }}_0^{\pm }\left(\mathsf{\epsilon }\right)={\displaystyle \frac{1}{2}}\lambda {tr}^2\left({\mathsf{\epsilon }}_{\pm }\right)+\mu tr\left({\mathsf{\epsilon }}_{\pm }^2\right),$$

(58)

where the split of the strain tensor is obtained through the projection on its symmetric positive and negative parts,

$${\mathsf{\epsilon }}_{\pm }={\mathrm{sym}}_{\pm }\left(\mathsf{\epsilon }\right).$$

(59)

Here ${sym}_{\pm }(•)$ denote the positive and negative projections of a symmetric tensor, ensuring that only tensile stress states contribute to the phase-field evolution.

Upon the tension-compression split choice, the elastic strain energy density in the anisotropic formulations is commonly rewritten such that only the positive part is affected by the degradation function $g\left(\varphi \right)$. Thus, it reads

$${\mathsf{\Psi }}^e(\mathsf{\epsilon },\varphi )=g\left(\varphi \right){\mathsf{\Psi }}_0^{+}\left(\mathsf{\epsilon }\right)+{\mathsf{\Psi }}_0^{-}\left(\mathsf{\epsilon }\right).$$

(60)

Following this, the stress tensor is given by

$$\mathsf{\sigma }=g\left(\varphi \right)·{\displaystyle \frac{\partial {\mathsf{\Psi }}_0^{+}\left(\mathsf{\epsilon }\right)}{\partial \mathsf{\epsilon }}}+{\displaystyle \frac{\partial {\mathsf{\Psi }}_0^{-}\left(\mathsf{\epsilon }\right)}{\partial \mathsf{\epsilon }}}.$$

(61)

Figure 24 illustrates how these formulations effectively separate tensile (active) and compressive (inactive) contributions under fatigue loading.

The solution of a PF boundary value fracture problem follows from the minimization of the total energy functional of Equation (46). Thus, by applying the variational principle, one may write the balance equations in the differential form as follows

$$\begin{array}{cc}\hfill \nabla ·\mathsf{\sigma }& =0\hfill \\ \hfill \nabla ·\mathsf{\xi }-\omega & =0\hfill \end{array}\mathrm{in}\mathsf{\Omega },$$

(62)

with the corresponding natural boundary conditions

$$\begin{array}{cc}\hfill \mathsf{\sigma }·\mathbf{n}& =\mathbf{t}\hfill \\ \hfill \mathsf{\xi }·\mathbf{n}& =0\hfill \end{array}\mathrm{on}\partial \mathsf{\Omega },$$

(63)

where $\omega$ and $\mathsf{\xi }$ are microstress quantities work-conjugate to the phase-field order parameter $\varphi$ and its gradient $\nabla \varphi$, respectively. In addition, $\mathbf{t}$ represents the traction and $\mathbf{n}$ the outward unit normal to the boundary surface, $\partial \mathsf{\Omega }$.

The balance equation given by Equation (62)₂ governs the phase-field evolution. As devised by Miehe [208], this evolution law may be expressed as

$${\displaystyle \frac{G_c}{2c_w\ell }}\left[{\displaystyle \frac{w^{\prime }\left(\varphi \right)}{2}}-{\ell }^2{\nabla }^2\varphi \right]+g^{\prime }\left(\varphi \right)\mathcal{H}=0,$$

(64)

with

$$\mathcal{H}(\mathbf{x},t)=\underset{\tau \in [0,t]}{max}{\mathsf{\Psi }}_0^{+}(\mathbf{x},\tau )$$

(65)

being a history-field variable that accounts for irreversibility in the crack phase-field evolution, i.e., in cyclic and loading-unloading scenarios.

Within a FEM framework, two main numerical strategies have been traditionally employed to solve the system composed by Equations (62)₁ and (62)₂. In the so-called monolithic approach, the displacement and phase-field sub-problems are solved simultaneously, while in the staggered scheme, a sequential solution procedure is adopted [213].

The monolithic approach is known for its unconditional stability and the ability to use larger time steps, which makes it computationally efficient. However, since the total potential energy is non-convex with respect to both $\mathbf{u}$ and $\varphi$, convergence difficulties may arise. In contrast, the energy functional is convex in $\mathbf{u}$ for a fixed $\varphi$ and vice versa, which makes the staggered approach robust in terms of convergence but less efficient due to its conditional stability. Recently, it has been found that the quasi-Newton method Broyden-Fletcher-Goldfarb-Shanno (BFGS) is able to reconcile the efficiency of the monolithic approach with the convergence robustness characteristic of staggered schemes.

Efficient solution schemes are particularly important when extending PF formulations to fatigue, since resolving long cyclic loading histories typically leads to a significant increase in computational cost. In this context, some of the main developments of PF models for fatigue are discussed in the following [213].

3.3.1. Extension to Fatigue

PF models for fatigue are relatively recent, with most developments emerging over the last decade. Despite this, a substantial research effort has been devoted to the topic due to the engineering relevance of fatigue-induced fracture. Most formulations proposed in this period can be broadly classified into two model classes [217].

In the first class, the material fracture energy (toughness), $G_c$, is progressively reduced as a result of the application of cyclic loads. In practical terms, this is achieved by introducing a fatigue degradation function, $f\left(\overline{\alpha }\right)$, into the phase-field evolution law of Equation (64), so that it becomes

$${\displaystyle \frac{f\left(\overline{\alpha }\right)G_c}{2c_w\ell }}\left[{\displaystyle \frac{w^{\prime }\left(\varphi \right)}{2}}-{\ell }^2{\nabla }^2\varphi \right]+g^{\prime }\left(\varphi \right)\mathcal{H}=0,$$

(66)

where $\overline{\alpha }$ denotes a fatigue history variable which can be taken as the accumulation of any scalar quantity that can describe the cyclic history of the material. Here it is important to mention that $f\left(\overline{\alpha }\right)$ is a dimensionless function and varies from 1 to 0 as it progressively degrades the material toughness.

In contrast, the second class of models relies on an additional energy term, ${\mathcal{H}}_F\left(\overline{\alpha }\right)$, which adds the fatigue contribution to the crack driving force. The phase-field evolution law is then written in this case as

$${\displaystyle \frac{G_c}{2c_w\ell }}\left[{\displaystyle \frac{w^{\prime }\left(\varphi \right)}{2}}-{\ell }^2{\nabla }^2\varphi \right]+g^{\prime }\left(\varphi \right)\mathcal{H}+{\tilde{g}}^{\prime }\left(\varphi \right){\mathcal{H}}_F\left(\overline{\alpha }\right)=0,$$

(67)

where ${\tilde{g}}^{\prime }\left(\varphi \right)$ is a degradation function (potentially chosen equal to $g^{\prime }\left(\varphi \right)$) acting on the additional energy term. Unlike $f\left(\overline{\alpha }\right)$, the term ${\mathcal{H}}_F\left(\overline{\alpha }\right)$ is an energy quantity that has no upper bound and its order magnitude is defined upon parametric calibration.

• First class of models: Based on a fatigue degradation function

The first class of models originates from the work by Alessi et al. [218], who proposed – for the one-dimensional setting – four different fatigue functions using the accumulated strain, $\overline{\mathsf{\epsilon }}$, as fatigue history variable. Later, Carrara et al. [47] extended this framework to handle three dimensional problems and reformulated the fatigue functions such that they may take any scalar quantity representative of the material cyclic history (chosen in their work as the strain energy density). Accordingly, these functions read

$$f\left(\overline{\alpha }\left(t\right)\right)=\left\{\begin{array}{cc}1\hfill & \mathrm{if}\overline{\alpha }\left(t\right)\le {\alpha }_T\hfill \\ {\left({\displaystyle \frac{2{\alpha }_T}{\overline{\alpha }\left(t\right)+{\alpha }_T}}\right)}^2\hfill & \mathrm{if}\overline{\alpha }\left(t\right)>{\alpha }_T,\hfill \end{array}\right.\left(\mathrm{Asymptotic}\right)$$

(68)

and

$$f\left(\overline{\alpha }\left(t\right)\right)=\left\{\begin{array}{cc}1\hfill & \mathrm{if}\overline{\alpha }\left(t\right)\le {\alpha }_T\hfill \\ {\left[1-{\kappa }_{f}log\left({\displaystyle \frac{\overline{\alpha }\left(t\right)}{{\alpha }_T}}\right)\right]}^2\hfill & \mathrm{if}{\alpha }_T\le \overline{\alpha }\left(t\right)\le {\alpha }_T{10}^{1/{\kappa }_f}\hfill \\ 0\hfill & \mathrm{if}\overline{\alpha }\left(t\right)\ge {\alpha }_T{10}^{1/{\kappa }_f},\hfill \end{array}\right.\left(\mathrm{Logarithmic}\right)$$

(69)

where ${\alpha }_T$ is a fatigue threshold parameter and ${\kappa }_f$ governs the slope of the logarithmic degradation function. The essential difference between them is that the asymptotic function (68) decreases smoothly toward zero for $\overline{\alpha }\left(t\right)\to \infty$, whereas the logarithmic function (69) reaches zero at a finite value of $\overline{\alpha }\left(t\right)$, as illustrated in Figure 25.

Carrara et al. [47] further proposed two accumulation strategies for the fatigue history variable $\overline{\alpha }\left(t\right)$. The first, defined as

$$\overline{\alpha }(\mathbf{x},t)={\int }_0^{t}H\left(\alpha \dot{\alpha }\right)\left|\dot{\alpha }\right|\mathrm{d}\tau ,$$

(70)

is independent of the mean load. The second, which accounts for mean load effects, is expressed as

$$\overline{\alpha }(\mathbf{x},t)={\displaystyle \frac{1}{{\alpha }_n}}{\int }_0^{t}H\left(\alpha \dot{\alpha }\right)\alpha \dot{\alpha }d\tau ,$$

(71)

where $H\left(\alpha \dot{\alpha }\right)$ is the Heaviside function, defined as

$$H\left(\alpha \dot{\alpha }\right)=\left\{\begin{array}{cc}1\hfill & \mathrm{during}\mathrm{loading}(\alpha \dot{\alpha }>0)\hfill \\ 0\hfill & \mathrm{during}\mathrm{unloading}(\alpha \dot{\alpha }\le 0).\hfill \end{array}\right.$$

(72)

The parameter ${\alpha }_n$ in Equation (71) is a normalization constant introduced to ensure dimensional consistency.

As evidenced by the Heaviside function in Equation (72), the fatigue history variable accumulates only during the loading phase in these approaches, which leads to nonphysical accumulation under monotonic loads. To overcome this limitation, Seleš et al. [219] proposed an alternative strategy in which the accumulation occurs in the unloading phase. However, it has been found that this method may introduce artificial increases in the fatigue history variable particularly behind the crack tip, because material points in this region may experience local unloading as the crack propagates.

More recently, Golahmar et al. [216] proposed a reversal-based accumulation strategy, in which the fatigue history variable increases only once per cycle, at each load reversal (see Figure 26). This approach not only eliminates the nonphysical accumulation issues of previous methods but also significantly improves the computational efficiency, as the fatigue variable is updated only once per cycle rather than at every intermediate time step of the loading or unloading phases. In this formulation, the increment in fatigue history variable is expressed as

$$\Delta \overline{\alpha }={\left({\displaystyle \frac{{\alpha }_{max}}{{\alpha }_n}}\right)}^n{\left({\displaystyle \frac{1-R}{2}}\right)}^{2\gamma n}H\left(\underset{\tau \in [0,t]}{max}{\alpha }_{max}{\left({\displaystyle \frac{1-R}{2}}\right)}^{2\gamma }-{\alpha }_e\right)$$

(73)

where the exponent n and the term ${({\alpha }_{max}/{\alpha }_n)}^n$ are introduced to provide the model with flexibility to properly capture the slope of the resulting S-N curves. As in Carrara’s model, the parameter ${\alpha }_n$ is included to ensure dimensional consistency. The constant ${\alpha }_e$ is a fatigue threshold introduced in the model to allow the representation of the material endurance limit below which fatigue does not occur. The Heaviside function H is equal to one for positive arguments and zero otherwise, ensuring that the fatigue variable does not increase when the driving force is below ${\alpha }_e$. When the strain energy density is adopted as the fatigue variable, one may estimate this constant as ${\alpha }_e=S_e^2/\left(2E\right)$. The parameter, $\gamma \in [0,1]$ is the Walker mean stress parameter, characterizing the material’s sensitivity to the stress ratio. Because this approaches incorporates the stress ratio R in its formulation, it explicitly accounts for mean stress effects.

Following a different rationale, several accumulation strategies have been formulated directly in the cycle domain, rather than in the physical time domain, so that the fatigue process is described continuously with respect to the cycle variable instead of relying on each loading and unloading phase. As a consequence, these approaches employ a representative (often constant) loading cycle – characterized by a stress amplitude and stress ratio – rather than a fully time-varying load history, which considerably improves computational efficiency. On the other hand, since the accumulation variables are based on empirical total-life concept, these models typically require experimental data such as S–N curves or Paris-law parameters as input.

A representative example of cycle-dependent accumulation strategy is the Miner-rule-based formulation adopted by Seiler [220], which is expressed as

$$\overline{\alpha }\left(N\right)={\int }_0^{N_c}{\displaystyle \frac{1}{N_f}}\mathrm{d}N,$$

(74)

where $N_f$ denotes the number of cycles to failure associated with the current loading state.

Alternatively to the Miner-type approach, Grossmann–Ponemon et al. [221] presented an energy-based accumulation strategy of the form

$$\overline{\alpha }\left(N\right)=1-{\displaystyle \frac{1}{{\alpha }_n}}{\int }_0^{N_c^{*}}{〈1-{\displaystyle \frac{{\mathsf{\Psi }}_{min}^e}{g\left(\varphi \right){\mathsf{\Psi }}_{vol,+}^e\left(N\right)}}〉}_{+}\mathrm{d}N,$$

(75)

which stems from the lifetime formulation originally introduced by Mesgarnejad et al. [222]. Here ${\mathsf{\Psi }}_{min}^e$ denotes the minimum elastic strain energy attained within a representative load cycle, ${\mathsf{\Psi }}_{vol,+}^e\left(N\right)$ is the tensile volumetric elastic strain energy at cycle N, and ${\alpha }_n$ is a normalization constant ensuring dimensional consistency. In addition, the upper integration limit is defined as

$$N_c^{*}=min\left\{N_c;N\mathrm{such}\mathrm{that}\varphi =0.5\wedge \alpha ={\alpha }_1\right\},$$

(76)

so that the accumulated fatigue variable $\alpha$ is evaluated only up to the point where the phase-field order parameter, $\varphi$, reaches 0.5, which the authors use as the criterion for the onset of failure. Finally, the Macaulay bracket ${\langle ·\rangle }_{+}$ ensures that accumulation takes place only when the fatigue driving quantity, $g\left(\varphi \right){\mathsf{\Psi }}_{vol,+}^e\left(N\right)$, exceeds ${\mathsf{\Psi }}_{min}^e$.

The PF fatigue formulations presented so far are all grounded in elastic material models. With the exception of the approach proposed by Seiler [220], which relies on the Neuber rule combined with the local strain method to capture the transition between the LCF and HCF regimes, these models remain applicable only within the high cycle domain.

Nevertheless, PF formulations have also been extended to model fatigue in the low cycle regime, which is typically characterized by significant plastic strain development near the crack tip. To account for this effect, the crack driving force term, $g^{\prime }\left(\varphi \right)\mathcal{H}$ in Equation (66) is usually supplemented by the plastic energy contribution, ${\mathsf{\Psi }}_p$.

Following this approach, Seleš et al. [219] proposed a PF fatigue model that includes an elastic-plastic constitutive law with isotropic and kinematic Chaboche-type hardening. Besides capturing the low cycle regime, this model degenerates to an elastic response under small stress amplitudes, thus also capturing the high cycle fatigue essential features. In addition, it is able to naturally reproduce the mean stress effect. In their work, the accumulation of the fatigue history variable is chosen to depend only on the tensile elastic strain energy, ${\mathsf{\Psi }}_e^{+}$.

Another formulation in which the crack-driving force is supplemented by the plastic energy was proposed by Ulloa et al [223]. Accordingly, their model features gradient-enhanced isotropic and multi-surface kinematic hardening, softening, and an explicit ratchetting strain variable. In addition, they employed Carrara’s logarithmic expressions as fatigue function and – unlike in Seleš et al. [219] – the fatigue history variable accumulates according to the combined contribution of the tensile elastic and plastic energy densities, i.e., ${\mathsf{\Psi }}_e^{+}$ and ${\mathsf{\Psi }}_p$.

Different from these two previous approaches, the model proposed by Khallil et al. [224] does not include the plastic energy in the crack-driving force term. Instead, this formulation extends Carrara’s model by implementing a Chaboche-type constitutive law with isotropic and nonlinear kinematic hardening to describe the material response. In this case, the accumulation of the fatigue history depends on the maximum temporal tensile elastic and plastic energy densities, ${\mathsf{\Psi }}_{e,max}^{+}\left(t\right)$ and ${\mathsf{\Psi }}_p$.

• Second class of models: Based on an additional energy term

As mentioned earlier, when it comes to the second class of models, the energy-based term, ${\mathcal{H}}_F\left(\overline{\alpha }\right)$ added to the crack-driving force represents the central mechanism by which fatigue effects are incorporated. According to this, several formulation have been proposed for ${\mathcal{H}}_F\left(\overline{\alpha }\right)$, each coupled with a specific accumulation model for the fatigue history variable. For example, in the work by Shreiber et al. [225], this term is defined as

$${\mathcal{H}}_F\left(\overline{\alpha }\right)=\left\{\begin{array}{cc}0\hfill & \mathrm{for}{\mathcal{H}}_F<{\mathcal{H}}_{F_1}\hfill \\ {\kappa }_{\mathcal{H}}{\langle {\mathcal{H}}_F-{\mathcal{H}}_{F_1}\rangle }^b\hfill & \mathrm{for}{\mathcal{H}}_F\ge {\mathcal{H}}_{F_1},\hfill \end{array}\right.$$

(77)

where ${\mathcal{H}}_{F_1}$ acts as a fatigue threshold, ${\kappa }_{\mathcal{H}}$ is a material scaling parameter, b controls the nonlinearity of the fatigue contribution and $\langle ·\rangle$ denotes the Macaulay bracket. Since this approach relies on S-N curves and the representative-cycle concept, the corresponding fatigue accumulation variable is expressed in the cycle domain as

$$\overline{\alpha }\left(N\right)={\int }_0^{N_c}{\displaystyle \frac{1}{N_f}}{\left[{\displaystyle \frac{\tilde{\sigma }\left(N\right)(1-L)}{A_f{(1-L)}^{\beta }}}\right]}^{{\kappa }_{\mathcal{H}}}\mathrm{d}N,$$

(78)

where $N_f$ is the number of cycles to failure under the current representative loading, $\tilde{\sigma }\left(N\right)$ is the equivalent stress amplitude, L the load ratio and $A_f$ and $\beta$ material constants derived from S–N data.

As Shreiber’s model is grounded in an elastic constitutive law, it is applicable only to HCF problems. A more comprehensive formulation in this class was presented by Haveroth et al. [226] who incorporated an elasto-plastic constitutive framework including non-isothermal processes, rate dependence and inertia effects in a thermodynamically consistent manner. In this formulation, the additional energy term takes the form

$$H\left(\overline{\alpha }\left(t\right)\right)=-{\displaystyle \frac{1}{\ell }}\overline{\alpha }\left(t\right),$$

(79)

where ℓ is the internal length scale and $\overline{\alpha }\left(t\right)$ represents the fatigue history variable. Accordingly, this variable is defined in the physical time domain as

$$\overline{\alpha }\left(t\right)={\int }_0^t{\displaystyle \frac{-\tilde{a}(1-\varphi )\parallel \dot{\mathsf{\epsilon }}\parallel }{\ell \theta }}\varphi \mathrm{d}\tau ,$$

(80)

where $\tilde{a}$ is a material parameter controlling the rate of fatigue damage, $\theta$ is the absolute temperature and $\parallel \dot{\mathsf{\epsilon }}\parallel$ is the norm of the strain-rate tensor.

• Other models Some PF fatigue models do not fit clearly into either the first or the second class. One example is the ductile PF model proposed by Aygun et al. [227] which employs the Armstrong-Frederick elasto-plastic constitutive law to describe the material cyclic response without introducing an explicit fatigue history variable or an additional energy-based fatigue term. In this formulation, the fatigue crack growth is governed solely by the accumulation of the plastic strain energy density, ${\mathsf{\Psi }}_p$, which is included in the crack-driving force term of the PF evolution law. Owing to this characteristics, this model is destined to LCF problems, where significant plastic strains develops at the crack tip.

Another example outside the two main classes is the formulation proposed by Lo et al. [228], who introduced a viscous term, $\eta \dot{\varphi }$, on the right-hand side of the evolution law in Equation (64) to account for fatigue crack growth. In this approach, the Paris law is used as model input and the load is applied statically so that no cycle-wise simulation is carried out.

• Simulation acceleration strategies

PF fatigue formulations are usually computationally demanding, particular when elasto-plastic constitutive laws are included in the model. This poses a major challenge for performing cycle-by-cycle simulations in the HCF regime, especially for components of engineering-relevant scales. As an alternative to this drawback, PF models formulated in terms of representative (often constant) loads have been introduced, as previously mentioned [220,225].

In addition to representative-load approaches, several acceleration strategies based on numerical considerations have been developed. Kristensen et al. [229], for example, proposed a modified Newton method to reduce the number of matrix factorizations. They further introduced a cycle-jump strategy, in which only the load cycles that effectively contribute to fatigue evolution are considered. This constant-load accumulation approach significantly reduces the number of load increments required to perform the simulations.

Another acceleration approach was presented by Yang and Shen [230]. Building on the time-scale homogenization theory, they introduced an efficient scheme for Carrara’s model in which fatigue degradation is accelerated through an adaptive time-stepping algorithm.

Figure 27(a) schematically illustrates the load-application procedure in models based on a representative load, in contrast to those employing a cycle-jump schemes, while Figure 27(b) shows how the fatigue history variable accumulates throughout the simulation in both scenarios.

4. Comparative Assessment of the Reviewed Approaches

As classical methods for fatigue analysis, the total-life approaches presented in Section 2 are often employed in practical applications due to their simplicity. Relying on empirical data such as Wöhler and strain-number of reversal curves, these methods act, however, as purely life-estimation tools and do not provide information about the progressive degradation of the material properties and the corresponding crack propagation process.

In contrast, the classes of fracture models discussed in Section 3 encompass formulations designed to describe at least the crack-growth phase, as in the case of FM, or even the entire material degradation process – from its intact to fully broken state – as captured by some CDM approaches and PF models. Accordingly, to illustrate their crack-representation capabilities, numerical results of the well-known Arrea and Ingraffea beam test [231] are comparatively shown in Figure 28 (a), (b), and (c), respectively, for an FM XFEM-based model, a CDM formulation using mixed finite elements, and a PF model, as reported in [232].

As discussed in Section 3.3, PF models derived from the variational formulation originally proposed by Francfort and Marigo [205] and later regularized by Bourdin [206] rely on a physically sound theoretical framework, as they $\mathsf{\Gamma }$-converge to Griffith’s theory of fracture when $\ell \to 0^{+}$. At the same time, a strength-based failure response can be incorporated in these formulations through the appropriate choice of the degradation function and the crack surface density function (e.g., AT1 and AT2 models), as demonstrated by Kristensen et al. [229].

Regularized CDM models – often formulated within a thermodynamically consistent framework – can, in turn, be made compatible with FM by tailoring the post-peak softening behavior such that the energy dissipated within the fracture process zone matches the material fracture energy. In these approaches, the strength criterion is naturally embedded in the constitutive response.

The ability of both model families to consistently reconcile the energy-based criterion of FM with a strength-based failure response enables them to capture the well-known size effect, as discussed in Section 3.2 and Section 3.3.

From a numerical standpoint, FM models require the definition of an ad hoc criterion to determine the crack-propagation direction under mixed-mode loading. Since, in such cases, the crack path is generally not known a priori, its accurate representation often demands either the use of remeshing procedures or XFEM formulations. Moreover, the extension of these models to three dimensional problems is usually not straightforward. As a result, FM approaches exhibit a limited capability to model phenomena such as multiple cracking, branching, and merging under complex geometrical configurations or loading conditions.

When applied to fatigue problems, the fact that FM approaches focus solely on crack growth poses a significant limitation for life-estimation purposes, since the initiation phase preceding the macrocrack nucleation may account for up to 90% of a component total fatigue life [217], particularly in the high cycle regime. This positions CDM and PF models as more suitable candidates for life-estimation analyses of this type, as they naturally capture both the initiation and propagation phases of fatigue damage. In addition, they are able to seamlessly describe the complex crack phenomena previously mentioned.

Nevertheless, FM-based approaches generally exhibit a reduced computational cost when compared to both PF and CDM counterparts. In particular, PF models typically require highly refined meshes in the vicinity of the crack (on the order of 4–10 finite elements across the regularization bandwidth) to accurately resolve the length scale ℓ, which is introduced through a continuum-level regularization scheme, as discussed in Section 3.3. This often leads to prohibitive computational costs, especially for high cycle fatigue simulations in large-scale three dimensional geometries.

In this sense, CDM smeared approaches emerge as more computationally affordable alternatives, since the regularization is introduced at the discrete level rather than at the continuum scale, typically requiring only 1–2 elements across the localization band to ensure objective energy dissipation. However, when standard finite element formulations are employed, mesh-objective results are not generally guaranteed. This limitation can be overcome through the use of mixed finite element formulations, in which an independent strain field is introduced and its continuity is enforced [232]. Although this strategy entails an additional computational cost, it remains considerably more affordable than phase-field formulations in most practical applications.

Table 4 summarizes the comparative assessment discussed in this section. Total-life approaches are not included in this comparison, as they are not suited to describe crack propagation phenomena.

5. Conclusions

In this work, different approaches for modelling fatigue in metals have been reviewed. First, total-life methods were presented as the classical approaches for fatigue life estimation. Although widely employed in practical applications due to their simplicity, these methods do not provide a description of crack propagation and therefore act primarily as life-estimation tools.

Fracture mechanics-based models, often referred to as defect-tolerant approaches, were subsequently reviewed. These formulations are primarily conceived to describe the crack-growth process through propagation laws driven by fracture parameters. As defect-tolerant approaches, they assume the existence of a pre-existing crack or defect and therefore do not describe the progressive degradation process leading from the intact material state to crack initiation. Instead, their domain of applicability begins once a crack has already nucleated and becomes amenable to fracture-mechanics analysis.

More comprehensive descriptions of fatigue degradation can be achieved through continuum damage mechanics and phase-field formulations, which are capable of representing the progressive deterioration of the material from its intact state to complete failure, including both crack initiation and propagation. Phase-field models, derived from a variational framework that $\mathsf{\Gamma }$-converges to Griffith’s theory of fracture, provide a physically sound and versatile description of complex fracture phenomena. However, their computational cost is generally high, particularly in high cycle fatigue analyses of large-scale three-dimensional structures. In this context, continuum damage models emerge as computationally more affordable alternatives. When properly regularized and calibrated, they also retain consistency with fracture mechanics, making them particularly attractive for engineering applications involving large-scale structures.