A surface integrity-informed crystal-plasticity based modelling of fatigue crack initiation in aerospace-grade Ti-6Al-4V

Abstract

A Crystal Plasticity Finite Element Method (CPFEM)-based approach is developed to predict fatigue crack initiation sites in aerospace-grade Ti-6Al-4V in an as-machined condition. The model integrates an Electron Backscatter Diffraction (EBSD)-reconstructed microstructure, incorporating prior deformation history, residual stress, and milling-induced roughness. The constitutive formulation employs Armstrong-Frederick kinematics with microstructure-sensitive fatigue indicators. Results show that fatigue initiation is highly localised, driven by surface topology, strain accumulation, and crystallographic variations with depth. Basal and prismatic slip systems dominate deformation, with Extreme Value Statistics (EVS) revealing their susceptibility to fatigue crack initiation. The findings highlight that Schmid factor alone does not fully explain fatigue hotspots, as local stress states, grain interactions, and surface-induced geometric constraints also play a critical role in crack initiation.

Introduction

Titanium and its alloys are often among the chosen materials for critical engineering applications (e.g., rotative parts in jet engines)^1,2 owing to their great set of properties such as a high strength-to-weight ratio, excellent corrosion resistance and biocompatibility. These characteristics make them attractive materials for industrial applications, especially for the aeronautics, astronautics and biomedical fields. Despite their widespread use across various sectors, the aerospace industry remains the primary consumer of titanium alloys. ³ Among them, the α + β Ti-6Al-4V is the most widely used alloy due to its good reproducibility and an outstanding balance of properties, including strength, ductility, and resistant to different modes of fracture, especially fatigue. ⁴ However, high-cycle fatigue (HCF) is still the dominant failure mechanism in rotative Ti-6Al-4V components in jet engines, e.g., fan discs and compressor blades, which can lead to a final fracture with disastrous consequences. This necessitates a comprehensive understanding of the mechanisms governing fatigue crack initiation for optimising manufacturing processes to further enhance the durability of these components.

Experimental studies indicate that high-cycle fatigue behaviour in Ti-6Al-4V is strongly influenced by its microstructure. Lee et al. ⁵ investigated a set of bimodal Ti-6Al-4V alloys manufactured through different fabrication methods and observed substantial differences in fatigue strength due to phase morphology. Wu et al. ⁶ conducted a statistical analysis of experimental data from previous research, suggesting that the microstructure of Ti-6Al-4V alloy plays a major role in resistance to HCF. Their findings ranked microstructural based HCF strength in decreasing order of bimodal, lamellar and equiaxed morphologies. In their study on bimodal microstructures, Wu et al. reported that HCF strength initially increases and then decreases as the volume fraction of primary alpha $({\alpha }_p)$ increases. Additionally, they noted that an increase in the size of the ${\alpha }_p$ grains leads to a reduction in HCF strength.

Beyond the intrinsic properties and microstructural characteristics of Ti-6Al-4V, fatigue crack initiation can also be influenced by post-processing and downstream (i.e., post-forging) manufacturing operations. Machining processes, which are the main method to obtain the final geometries and finishing of aerospace components, can induce surface and subsurface alterations (i.e., induce damage) which determine fatigue crack initiation and propagation behaviours. Liu et al. ⁷ reviewed the effect of machined surface integrity on fatigue performance of metallic components. Their analysis revealed that minimising roughness, inclusions and defects, and increasing compressive residual stress, had beneficial impact on fatigue performance. In the review, Liu et al. also stated that work hardening delays nucleation but accelerates propagation. Sun et al. ⁸ investigated the effects of conventional drilling and helical milling on the fatigue behaviour of Ti-6Al-4V and found that the latter results in longer fatigue performance. The improvement was attributed to lower surface roughness and lower relative severe plastic deformation caused by the helical tool.

Several attempts have been made to predict and model fatigue crack initiation by means of CPFEM. Liu et al. ⁹ examined micro-textured regions of α-grain aggregates using Electron Backscatter Diffraction maps as input and found that crack initiation can result from the accumulation of plastic strain on prismatic ⟨a⟩ slip systems which form dislocation pile-ups and cause stress concentration. Chen et al. ¹⁰ studied the crack initiation mechanism in Ti-6Al-4V under very high-cycle fatigue (VHCF) regime using CPFEM and found that the model effectively captured stress concentration and plastic accumulation, while also highlighting the significant influence of α-colonies on VHCF performance. Dinh et al. ¹¹ modelled the effect of surface roughness on the fatigue behaviour of additively manufactured Ti-6Al-4V, using a fatigue indicator parameter (FIP) and found a strong agreement with experimental data based on the Smith-Watson-Topper variable. Recent CPFEM studies have also extended to model more complex fatigue scenarios, such as fretting fatigue in Ti-6Al-4V using submodelling techniques, ¹² and fatigue crack propagation under texture control in aluminium alloys using coupled crystal plasticity–cohesive zone frameworks. ¹³

This work proposes a methodology for predicting fatigue crack initiation in Ti-6Al-4V, using EBSD maps of the microstructure, and by incorporating the prior deformation history induced by the milling process. The mechanical behaviour of the material is described via a kinematic isotropic hardening crystal plasticity (CP) constitutive model, implemented in the finite element (FE) method. The model accounts for prior accumulated strain (i.e., induced by machining), residual stress, and surface roughness, while fatigue initiation is evaluated using FIPs to provide deeper insight into the complex interplay among the factors influencing fatigue crack initiation.

Methodology

This section outlines the theoretical foundation for developing a micromechanical simulation of Ti-6Al-4V. The methodology encompasses the CP formulation, the FIPs considered in this study, the experimental data acquisition processes, and the integration of material characterisations with the computational modelling.

Crystal plasticity model

High-cycle fatigue crack initiation is affected by microstructural characteristics and prior manufacturing processes. To evaluate the micromechanical response of Ti-6Al-4V, this study implements the CP theory in a FE model. A user-material subroutine (UMAT) based on a code published by Huang ¹⁴ was used to incorporate CP in Abaqus FEA. The UMAT updates the stress, the state dependent variables and provides the Jacobian matrix to the FE solver. This approach builds upon the pioneering kinematic theory of Taylor, ¹⁵ Hill ¹⁶ and Rice ¹⁷ and it is used to describe the deformation at a micro-scale crystal level.

The deformation gradient, $\mathit{F}$ can be decomposed into an elastic and rotational component associated with the stretching of the crystal lattice and rigid body rotation ( ${\mathit{F}}^e$ ), and an inelastic component defined by dislocation slip ( ${\mathit{F}}^p$ ), as follows:

$$\begin{array}{c}\mathit{F}={\mathit{F}}^e{\mathit{F}}^p\end{array}$$

(1)

The velocity gradient can then be calculated from the deformation gradient tensor with the following equation:

$$\begin{array}{c}\mathit{L}=\dot{\mathit{F}}{\mathit{F}}^{-1}=\mathit{D}+\mathbf{\Omega }\end{array}$$

(2)

In Equation (2), the velocity gradient is represented as the sum of the rate of deformation $\mathit{D}$ , which is a symmetric tensor denoting stretch, and an antisymmetric tensor $\mathbf{\Omega }$ denoting the continuum spin. The velocity gradient can be decomposed into an elastic and a plastic components satisfying:

$$\begin{array}{c}{\mathit{L}}_p={\dot{\mathit{F}}}^p{\mathit{F}}^{p-1}={\mathit{D}}^p+{\mathbf{\Omega }}^p=\sum _{\alpha =1}^{N_s}{\dot{\gamma }}^{\alpha }{\mathit{s}}^{\alpha }\otimes {\mathit{n}}^{\alpha }\end{array}$$

(3)

where ${\dot{\gamma }}^{\alpha }$ is the slipping rate in each slip system and $N_s$ represents the total number of slip systems. ${\mathit{s}}^{\alpha }$ and ${\mathit{n}}^{\alpha }$ are directional cosines that represent the slip direction and the normal to the slip plane in the deformed configuration, respectively. Assuming that crystalline slip follows Schmid's law, the resolved shear stress $\tau$ in a slip system $\alpha$ can be calculated as follows:

$$\begin{array}{c}\tau =\mathit{\sigma }:({\mathit{s}}^{\alpha }\otimes {\mathit{n}}^{\alpha })\end{array}$$

(4)

where $\mathit{\sigma }$ is the Cauchy stress tensor, which can be correlated to $\stackrel{\mathrm{\nabla }}{{\mathit{\sigma }}^e}$ , the elastic component of the Jaumann stress rate, with the following relationship:

$$\begin{array}{c}\stackrel{\mathrm{\nabla }}{{\mathit{\sigma }}^e}=\mathit{C}:{\mathit{D}}^e-\mathit{\sigma }(\mathit{I}:{\mathit{D}}^e)\end{array}$$

(5)

$\mathit{C}$ is the elasticity matrix and $\mathit{I}$ is the second order identity tensor. The determination of the stress increments and the Jacobian matrix depend on the shear strain increment associated with each slip system. This analysis utilises a power law relationship for the calculation of slip rates, building upon the foundational work of Peirce, Asaro, and Needleman. ¹⁸

$$\begin{array}{c}{\dot{\gamma }}^{\alpha }={\dot{\gamma }}_0{\left|{\displaystyle \frac{{\tau }^{\alpha }-x^{\alpha }}{g^{\alpha }}}\right|}^{n}sgn({\tau }^{\alpha }-x^{\alpha })\end{array}$$

(6)

where ${\dot{\gamma }}_0$ is the reference strain rate, ${\tau }^{\alpha }$ is the resolved shear stress of the slip system $\alpha$ , $g^{\alpha }$ is the current strength of the slip system which represents resistance to slip, $x^{\alpha }$ denotes back stress and n is the rate sensitivity exponent. The back stress is integrated to the model to account for the Bauschinger effect, that is, to consider the forces that emerge from dislocation pile-ups and their subsequent redistribution upon load reversal during fatigue. For this work, the Armstrong-Frederick model ¹⁹ was selected to represent the nonlinear kinematic hardening effects during fatigue loading:

$$\begin{array}{c}{\dot{x}}^{\alpha }=C{\dot{\gamma }}^{\alpha }-D{x}^{\alpha }|{\dot{\gamma }}^{\alpha }|\end{array}$$

(7)

$C$ is the kinematic hardening modulus and $D$ determines the decay of the back stress up to a maximum saturated stress. The evolution of slip system hardening is given by:

$$\begin{array}{c}{\dot{g}}^{\alpha }=\sum _{\alpha =1}^{N_s}{h}_{\alpha \beta }{\dot{\gamma }}^{\beta }\end{array}$$

(8)

where $h_{\alpha \alpha }$ and $h_{\alpha \beta }$ are the self- and latent-hardening moduli, respectively, and the summation includes all active slip systems. In the present study, the following expression formulated by Peirce et al. ²⁰ was employed:

$$\begin{array}{c}{h}_{\alpha \alpha }=h(\gamma )=h_0\left|{\displaystyle \frac{h_0\gamma }{{\tau }_s-{\tau }_0}}\right|\end{array}$$

(9)

where $h_0$ is the initial hardening modulus, ${\tau }_0$ is the yield strength, ${\tau }_s$ is the shear stress at which large plastic flow initiates (stage I), and $\gamma$ is the cumulative shear strain calculated as follows:

$$\gamma =\sum _{\alpha =1}^{N_s}\underset{0}{\overset{t}{\int }}|{\dot{\gamma }}^{\alpha }|dt$$

The latent hardening is calculated as the product of the self-hardening modulus and a constant $q$ :

$$\begin{array}{c}{h}_{\alpha \beta }=qh(\gamma )\alpha \ne \beta \end{array}$$

(10)

The calibration of the CP parameters was performed using a reduced microstructure with 70 grains, which facilitated the iterative adjustment of the values and was validated against the fatigue results of Kurath. ²¹ In the simulation, 12 slip systems were considered: three basal, three prismatic, and six first-order pyramidal ⟨a⟩ slip systems as described in Table 1.

Table 1.

Slip systems considered in the CPFEM simulation for the HCP crystals of Ti-6Al-4V.

Basal		Prismatic
$B1$	$(0001)⟨1\overline{2}10⟩$	$P1$	$(10\overline{1}0)⟨\overline{1}2\overline{1}0⟩$
$B2$	$(0001)⟨2\overline{1}\overline{1}0⟩$	$P2$	$(0\overline{1}10)⟨2\overline{1}\overline{1}0⟩$
$B3$	$(0001)⟨11\overline{2}0⟩$	$P3$	$(1\overline{1}00)⟨11\overline{2}0⟩$
First-order pyramidal ⟨a⟩
$Py1$	$(10\overline{1}1)⟨\overline{1}2\overline{1}0⟩$	$Py4$	$(\overline{1}011)⟨1\overline{2}10⟩$
$Py2$	$(0\overline{1}11)⟨2\overline{1}\overline{1}0⟩$	$Py5$	$(01\overline{1}1)⟨\overline{2}110⟩$
$Py3$	$(1\overline{1}01)⟨11\overline{2}0⟩$	$Py6$	$(\overline{1}101)⟨\overline{1}\overline{1}20⟩$

The elasticity matrix coefficients were obtained from the work of Mayeur et al., ²² who used generalised plane strain elements for FE simulations. The corresponding elasticity constants were C₁₁ = 162,400 MPa, C₁₂ = 92,000 MPa, C₁₃ = 69,000 MPa, C₃₃ = 180,700 MPa and C₄₄ = 46,700 MPa. The remaining parameters were adjusted based on the studies of Zang et al., ²³ Kapoor et al. ²⁴ and Han et al. ²⁵ Table 2 shows the final calibrated CP parameters.

Table 2.

Values of the calibrated CP parameters.

Parameter	Basal	Prismatic	First-Order Pyramidal	All slip systems
${\dot{\gamma }}_0$				0.001
$n$				50
$h_0[MPa]$	339	235	281
${\tau }_s[MPa]$	242	215	335
${\tau }_0[MPa]$	240	211	330
$C$				14,500
$D$				160

Fatigue indicator parameters

Analysing the formation of cracks due to cyclic loading requires meso-scale parameters to be defined which reflect driving forces for crack generation as a function of state variables and material properties. These parameters, commonly referred to as fatigue indicator parameters (FIPs) in the literature, utilise physical quantities to identify zones in the microstructure that are prone to crack initiation.

The experimental work of Wilkinson et al. ²⁶ and Ahmed, Wilkinson and Roberts ²⁷ using electron channelling contrast imaging demonstrated that persistent slip bands are associated with crack initiation at the single-crystal domain. Their research revealed that crack nucleation occurs along slip band planes. Therefore, the accumulated plastic strain can be considered as a FIP,^28–30 given by:

$$\begin{array}{c}FI{P}_p=\underset{0}{\overset{t}{\int }}\sqrt{{\displaystyle \frac{2}{3}\dot{\mathit{\epsilon }}:\dot{\mathit{\epsilon }}}}dt\end{array}$$

(11)

Equation (11) does not discern the contributions from ratcheting effects and those from reversed cyclic plasticity, as delineated by McDowell and Dunne, ³¹ which may be critical in the context of fatigue. McDowell introduced FIPs that explicitly reflect the roles of directional and reversed plasticity in cyclic fatigue.^32,33 For example, the Fatemi-Socie FIP ³⁴ is intended to represent transgranular slip band formation with the following expression:

$$\begin{array}{c}FI{P}_{FS}=\underset{\alpha \in (1,\dots ,N_s)}{max}\left[{\displaystyle \frac{\mathrm{\Delta }{\gamma }_{p,max}^{\alpha }}{2}\left(1+k⟨{\displaystyle \frac{{\sigma }_{n,max}^{\alpha }}{{\sigma }_0}}⟩\right)}\right]\end{array}$$

(12)

where $\mathrm{\Delta }{\gamma }_{p,max}^{\alpha }$ is the maximum shear strain range over a cycle for slip system $\alpha$ , ${\sigma }_{n,max}^{\alpha }$ is the maximum normal stress on the critical plane, i.e., the plane of maximum shear strain range, k is a material constant controlling the influence of the normal stress, and ${\sigma }_0$ is the macroscale yield strength of the material.

Experimental methods

Samples of an aerospace-grade Ti-6Al-4V material, extracted from a final forged and heat-treated fan disc, were subjected to a finishing milling pass using a face cutter, replicating a typical manufacturing step for aerospace components. The depth of cut was 0.3 mm, and the feed speed at the machined diameter was 737 mm/min. After milling, surface roughness was measured at the centre of the milling path using focus variation microscopy with a Hirox digital optical microscope, applying a cut-off value of 2.5 mm.

Following surface roughness measurements, near surface residual stress was quantified using an Electronic Speckle Pattern Interferometry (ESPI) hole-drilling apparatus, manufactured by Stresstech. The milled Ti-6Al-4V sample was coated with a thin layer of paint to enhance the quality of the surface for more precise surface movement (i.e., caused by residual stress relaxation during hole-drilling) measurements by the ESPI. A 1 mm diameter hole was incrementally drilled into the material, at 30 steps, resulting in stress relaxation due to material removal (i.e., the hole). As the drilling progressed, a monochromatic laser beam, with 532 nm wavelength, was shone on the surface to generate the speckle pattern, such that the hole was in the centre of the image. Progressive step-by-step hole drilling resulted in surface displacement which was then measured by the ESPI. These displacements were analysed to compute the in-plane residual stress components as a function of depth, assuming a Poisson's ratio of 0.342 and a Young's modulus of 113.8 GPa.

Integration of material characterisation and computational modelling

A Ti-6Al-4V sample, taken from the same fan disc (i.e., as those used for residual stress and surface roughness analyses), was prepared for materials characterisations. Following standard metallography preparation, a scanning electron microscopy (SEM) image of the material was taken and is presented in Figure 1(a). EBSD analysis was performed to obtain orientation data from the sample. The EBSD data was acquired using an accelerating voltage of 20 kV over a 149 × 103 μm region with a step size of 0.146 μm. The acquired EBSD data was post-processed using MTEX software, an open-source MATLAB toolbox. A misorientation threshold of 7° was applied to define grain boundaries. Noise reduction was performed by removing grains with fewer than 10 pixels, and smoothing was implemented to minimise the staircase effect. The resulting orientation map with average local crystal orientations is shown in Figure 1(b).

Figure 1.

(a) An SEM image of the Ti-6Al-4V sample, (b) EBSD orientation map of the same region processed by MTEX with mean crystal orientation colour map, (c) the digital microstructure imported to Abaqus FEA, and (d) the 2D digital microstructure made as a part in Abaqus with the imposed surface roughness profile on the top surface.

Geometric and orientation data were exported from the MTEX-created EBSD object. This data served as input for an in-house developed Python script designed to interface with the Abaqus FEA software through its application programming interface (API), to import the microstructure as a 2D part. The script employs a geometric algorithm designed to decrease the vertex count of each polygon representing a crystal within the microstructure. This vertex reduction aims to minimise the requisite elements for meshing the microstructure while preserving the essential morphological features. The digital microstructure is presented in Figure 1(c). Additionally, a separate script was used to modify the digital microstructure by trimming its top boundary to match the experimentally measured roughness profile, as shown in Figure 1(d).

The microstructure was meshed in Abaqus FEA with 202703 plane strain triangular elements. The elements were of quadratic order, with an average shape factor of 0.941 and an average aspect ratio of 1.22. The average edge length was 0.362 μm, which is only slightly larger than the EBSD acquisition step size. This indicates that the mesh resolution is sufficient to preserve the orientation fidelity of the microstructure without requiring additional boundary treatments. Boundary conditions for $x=0$ ( $U_x=U{R}_y=U{R}_z=0$ ) and $y=0$ ( $U_y=U{R}_x=$ $U{R}_z=0$ ) and a cyclic displacement of 0.8% engineering strain ( $R=0$ ) were applied. The simulation was executed for 10 cycles to achieve a steady cyclic state. Littlewood et al. ³⁵ demonstrated that after the initial few cycles of fatigue, both the accumulation of slip and the densities of dislocations progress at nearly steady rates per cycle. Moreover, this approach is consistent with prior CPFEM studies showing that fatigue indicator parameters stabilise within the first few cycles under HCF conditions, with variations of only a few percent beyond the third cycle.^36,37 The material properties were assigned by parsing the input file with a Python script and adding local orientation properties to each element based on the MTEX EBSD map.

To model residual stress, a SIGINI user subroutine was implemented. The SIGINI Fortran code was automatically generated using a Python script that was developed to apply residual stress values taken from the measured residual stress curves as a function of depth, considering the imposed roughness profile. The residual stress components were assigned as σ_xx and as σ_zz, as in-plane and out-of-plane stresses that result from a machining operation.

During machining, the near-surface region of the Ti-6Al-4V component undergoes plastic deformation. The micromechanical modelling of this initial strain was based on a study conducted by Brown et al., ³⁸ who examined different orthogonal cutting conditions and measured the strain by printing a grid using electron beam lithography. Strain values were obtained by comparing SEM images of deformed and undeformed grids. The strain was introduced in the CP framework by developing a Fortran code that iterates in each integration point up to the strain tensor values as a function of the position of the element. This strategy ensures that the strength, the total cumulative shear strain and the kinematic hardening effect are considered from the previous strain history in the fatigue simulation. A flowchart of the computational methodology is depicted in Figure 2, which provides a holistic visual representation of the processes.

Figure 2.

Flowchart of the computational framework employed to build the fatigue simulation.

Results and discussion

Roughness profile and residual stress

Figure 3(a) shows an optical microscopy image of the region where the roughness profile was analysed. Figure 3(b) depicts a 3D topographic map of the centre of the milling path. The roughness profile throughout the red line indicated in the map in Figure 3(b), is plotted in Figure 3(c). The average $R_a$ of the profile was analysed as 12.73 μm.

Figure 3.

(a) Optical microscopy image of the milled Ti-6Al-4V surface, (b) 3D roughness map obtained from the milled surface using a Hirox microscope, (c) roughness profile extracted from the centre of the 3D roughness map, indicated by the red dotted line in (b).

The in-plane residual stress profiles for different components as a function of depth, measured using ESPI hole-drilling method, are presented in Figure 4. Each curve represents the average of three measurements per milling path, with the shaded regions indicating the standard error for all measurements. The results show that the in-plane residual stress components were compressive in both the cutting and transverse directions of the tool. The maximum compressive stress in the cutting direction ( ${\sigma }_{xx}$ ) was −369.9 MPa, while in the transverse direction, it reached −430.9 MPa.

Figure 4.

Profiles of near-surface in-plane residual stress components measured with ESPI hole drilling, as a function of depth, in the as-milled Ti-6Al-4V sample.

To incorporate these experimental findings into the computational model, key methodological considerations must be addressed. The computational method employed to create a digital representation of the microstructure (Figure 1) uses an orientation threshold of 7° to define grain boundaries, as shown in the transformation from Figure 1(a) and (b). However, it is important to note that the EBSD data was obtained from a bulk material region, as severe deformation near the surface led to highly sparse data and unreliable indexing. Consequently, direct measurement of surface roughness through algorithmic reconstruction was not feasible. Instead, the roughness profile was computationally integrated by modifying the digital microstructure to account for surface effects. While this computational strategy enables the incorporation of experimental data into the model, it inherently introduces uncertainties in capturing specific microstructural details from the physical twin and mapping them onto the digital twin for fatigue modelling. These uncertainties may influence the accuracy of fatigue predictions, particularly regarding the interaction between surface integrity and crack initiation.

Stress distribution and accumulated plastic strain

Figure 5(a) shows the predicted Von Mises stress distribution in the microstructure after being subjected to cyclic loading. Von Mises stress ranged from 58 MPa to 1780 MPa, with the highest stresses located at the machined surface of the material. The digitally represented microstructure presented in this study captures local orientation effects, including grain misorientation (i.e., variations inside grains), therefore, the capacity of accommodating strain by the grains results in the final shape of the stress field. Consequently, grain boundary effects are implicitly captured through the spatial variation in crystallographic orientation, which alters the stiffness and slip activity across interfaces. The model predicts regions of high Von Mises stress arising from geometric stress concentrators present in the surface, suggesting that the topological features of the surface profile significantly affect the mechanical response of the material. It can also be observed that the stress field magnitude progressively diminishes in the direction of the subsurface layers. Despite that tendency, the stress distribution remains heterogeneous due to the combined effects of local orientation and morphological characteristics.

Figure 5.

Von Mises stress distribution in the modelled microstructure.

The accumulated plastic strain fatigue indicator parameter distribution after ten cycles is presented in Figure 6. The observed distribution results from the combined effects of machining-induced deformation and the applied cyclic loading. The highest magnitudes in this distribution are near the surface and the $FI{P}_p$ values tend to decrease with depth (i.e., distance from the top surface). This simulation indicates that stress concentrations due to the roughness profile (peaks and troughs) lead to higher $FI{P}_p$ near the surface. Since the $FI{P}_p$ does not differentiate between compressive and tensile strains, regions with high compressive strain may exhibit intense $FI{P}_p$ values, despite being less likely to act as a favourable sites for fatigue crack initiations under cyclic loading conditions. Consequently, a FIP based solely on accumulated plasticity may be insufficient for predicting fatigue crack initiation in machined components with an initial compressive strain profile. Additionally, the distribution in Figure 6 reveals three distinctive bands of non-uniform plasticity accumulation, with higher plastic strain localised near the surface and decreasing towards the bulk.

Figure 6.

(a) FIP distribution based on accumulated plastic strain in the micromechanical model conducted in this work. The model predicts three distinct regions (highlighted), including surface and near-surface zone, a transition zone, and a lower zone into the bulk material, each showcasing different levels of plastic strain accumulation.

Figure 7 shows the evolution of accumulated plastic strain in the full microstructural model, computed as a weighted average of all integration points, along with three additional plots for each of the specific regions highlighted in Figure 6 throughout the fatigue simulation. Initially, a rapid increase in accumulated plastic strain is observed, followed by a plateau, i.e., nearly constant accumulation rate, for all regions. During the first cycle, yielding occurs in the microstructure, potentially leading to an adjustment of grain orientations toward a more stable state, which may result in a reduced rate of plastic strain accumulation in subsequent cycles. At a finer scale, this behaviour results from the combined effects of localised plastic strain accumulation in regions exceeding their yield strength and ratcheting, which is associated with the Bauschinger effect.

Figure 7.

Evolution of the weighted average accumulated plastic strain as a function of time during the fatigue simulation. The results are presented for (a) the full model, and (b), (c) and (d) three distinct regions of the surface and near-surface, a transition zone, and a lower region into the bulk, respectively.

The model suggests that during the loading and unloading phases of the cycle, plastic strain is accumulated at relatively slow rate, while the upper and bottom limits of the cycle are responsible for a discontinuity in the curve. This result indicates that at the highest value of displacement, some regions are undergoing plastic deformation as they surpass the yield surface threshold. When unloading, kinematic hardening effects cause a displacement of the yield surface which also leads to an increase in accumulated plasticity. The non-symmetrical behaviour induced by the Bauschinger effect in specific zones of the model produces a ratcheting effect, which is associated with high $FI{P}_{FS}$ values. Furthermore, the analysis of the slopes in Figure 7 reveals a trend in plastic strain accumulation as fatigue cycles progress. The surface and near-surface region exhibit the highest rate of plastic strain accumulation, with a slope greater than that of the transition and lower regions. Highly localised plastic strain under cyclic loading is a characteristic of persistent slip bands. Experimental studies demonstrated that the presence of persistent slip bands act as a precursors of fatigue crack nucleation.^26,27 Therefore a higher rate of strain accumulation suggests an increased likelihood of fatigue crack initiation. The plastic strain accumulation rate in the surface and near-surface region is approximately 53% higher than the model average, 58% higher than in the transition zone, and 63% higher than in the lower region, which can be considered as a reasonable representation of bulk material behaviour. This disparity underscores the critical role of near-surface plasticity in fatigue damage evolution and is consistent with previous established findings that state that fatigue crack initiation under HCF conditions predominantly occurs at free surfaces due to stress concentration caused by surface irregularities, unconstrained slip bands leading to intrusions and extrusions, residual stresses, and environmental effects.^39–41

Figure 8 shows the weighted cumulative shear strain (WCSS) in each slip system after the fatigue cycles. The basal slip systems collectively possess the highest WCSS sum, contributing approximately 47.8% of the total shear strain, followed closely by prismatic slip systems at 41.2%, and pyramidal slip systems at 11.0%. This distribution indicates a competitive activation of slip systems, with basal and prismatic slip systems contributing significantly to the deformation. Among individual slip systems, Basal 1 shows the highest contribution, followed by Prismatic 2 and Basal 2. The average WCSS per crystallographic slip family further underscores this distribution, with basal slip systems having an average strain of 8.2 × 10⁻³, prismatic slip systems at 7.1 × 10⁻³, and pyramidal slip systems at 9 × 10⁻⁴. This trend suggests that under the given fatigue conditions, both basal and prismatic slip systems are essential contributors to plastic accommodation, while the pyramidal slip systems provide a secondary mechanism. These results are consistent with previous experimental studies on α-Ti by Pourian et al. and S. Hémery et al., which demonstrated that the highest accumulated plastic strain occurs in the basal plane.^42,43 Moreover, Zhang et al. showed that basal slip forms continued slip bands in microtextured regions more frequently compared to other crystallographic families. The simulation results are also aligned with the experimental findings by Bridier et al., ⁴⁴ who investigated the cyclic plasticity activity in a commercial dual-phase α/β Ti-6Al-4V alloy and observed a predominance of slip along basal and prismatic planes. Bridier et al. noted that fatigue crack formation occurs primarily in α grains, initiating from both prismatic and basal slip planes. However, they found that basal cracks tend to form earlier and propagate more rapidly, with the fatal crack often forming on a basal plane. Similarly, Liu et al. ⁹ reported a crack at the end of a long prismatic slip band. Their CPFEM simulation revealed high plastic strain accumulation along this band, which was correlated with fatigue crack initiation.

Figure 8.

Weighted cumulative shear strain in each slip system after fatigue cycles.

Fatemi-Socie fatigue indicator parameter

The $FI{P}_{FS}$ results depicted in Figure 9 indicate that fatigue crack initiation is predominantly influenced by (i) topological stress concentration sites at the surface and (ii) localised elastic and plastic strain accumulation. Additionally, the simulation predicts an inhomogeneous distribution of fatigue crack initiation likelihood within microstructural colonies, primarily due to orientation variations, such as primary $\alpha$ -grains, where the formation of slip bands contributes to localised strain accumulation. The distribution of the critical planes for $FI{P}_{FS}$ per slip system is depicted in Figure 10. The prismatic slip systems possess the highest frequency of occurrences, with Prismatic 2 and Prismatic 1 being the most dominant, followed by Basal 1 and Basal 2. In terms of magnitude, Prismatic 3 has the highest $FI{P}_{FS}$ value at 1.37 × 10⁻², followed by Prismatic 2 at 1.23 × 10⁻². In comparison, the maximum $FI{P}_{FS}$ in the basal slip systems is in Basal 2, with a value of 1.11 × 10⁻². The pyramidal slip systems are less frequently critical planes, and the maximum $FI{P}_{FS}$ is 5.27 × 10⁻⁴ in Pyramidal 5. The distributions in Figure 10 suggest that fatigue crack initiation is more likely to occur in the basal or prismatic crystallographic slip families, consequently, the following analyses will be performed for those slip systems.

Figure 9.

Distribution of the Fatemi-Socie fatigue indicator parameter ( $FI{P}_{FS}$ ).

Figure 10.

Distribution of critical planes as a function of $FI{P}_{FS}$ for basal, prismatic and pyramidal slip systems.

Figure 11 presents a box plot of the $FI{P}_{FS}$ values for the slip systems with critical planes. The prismatic slip systems have higher $FI{P}_{FS}$ medians compared to the basal slip systems, suggesting a more dominant role in fatigue crack initiation. The contribution of extreme outliers to the $FI{P}_{FS}$ is over 90% in Basal 1, Prismatic 1, and Pyramidal 3 slip systems, reinforcing their tendency for localised fatigue initiation. This means that only a small fraction of these critical planes lead to an exceptionally high $FI{P}_{FS}$ , making them the primary fatigue hotspots. In contrast, Prismatic 2 and 3, despite high median values, exhibit lower outlier contributions (31.6% and 39.1%, respectively), indicating a more evenly distributed $FI{P}_{FS}$ throughout the microstructure.

Figure 11.

Box plot of the $FI{P}_{FS}$ as a function of the critical planes associated with basal and prismatic slip systems.

Additional extreme value analysis was conducted on the $FI{P}_{FS}$ outliers, and the results are presented in Figure 12. Given that high-cycle fatigue crack initiation is driven by highly localised fatigue indicator parameters, Extreme Value Statistics (EVS) were applied to assess the distribution tails. ⁴⁵ The generalised Pareto distribution (GPD) was fitted to the upper tail of each slip system's $FI{P}_{FS}$ using the peaks-over-threshold approach, with a threshold set at the 95th percentile. The 95th percentile was chosen following standard practices and confirmed through a sensitivity analysis showing stable GPD shape values across nearby thresholds. This choice ensures that only the most extreme values are considered while maintaining a sufficient number of exceedances for statistical robustness, allowing for a probabilistic characterisation of fatigue-critical sites. The results reveal that Basal 1 and Basal 2 exhibit the heaviest distribution tails, suggesting a higher probability of extreme localised fatigue indicator values. Basal 3 has the lowest shape parameter, indicating a rapid decay of extreme values, which suggests a lower probability of FIP values sufficiently high to initiate fatigue damage. Among the prismatic slip systems, Prismatic 3 shows the highest potential for extreme values. Prismatic 1 and Prismatic 2, show a more moderate tail behaviour, suggesting that while fatigue-prone, they contribute to a broader distribution of plastic strain rather than highly localised extreme values. The heavy-tailed distributions in Basal 1, Basal 2, and Prismatic 3 suggest that these slip systems are more prone to localised extreme values.

Figure 12.

Probability density in the basal and prismatic slip systems as a function of the $FI{P}_{FS}$ . The Generalised Pareto Distribution fit to the tail of the distribution is shown in red and the shape (ξ) and scale parameters of the GPD are included in each plot to quantify the tail behaviour of the distributions.

The distribution of critical planes as a function of Schmid Factor reveals distinct tendencies for slip activation across prismatic and basal systems, as shown in Figure 13. Basal 1 and Basal 3 have the lowest mean Schmid Factors, whereas Prismatic 1 and 2 have the highest. In Basal 2 and the prismatic slip systems, the distribution skews toward higher Schmid Factors, aligning with expectations since critical $FI{P}_{FS}$ planes tend to accumulate more plasticity in those orientations. However, Basal 1 deviates from this trend, suggesting that its slip activation is influenced by additional microstructural factors. Prismatic 2, Prismatic 3, and Basal 2 exhibit the highest $FI{P}_{FS}$ values, with most critical planes aligning with high Schmid Factors. Notably, Prismatic 1, despite having a higher Schmid Factor at its maximum $FI{P}_{FS}$ location than Prismatic 3, shows a lower maximum $FI{P}_{FS}$ value, indicating that slip activity is not dictated by Schmid Factors alone. A similar trend is observed in the basal slip systems, where Basal 1 exhibits the highest Schmid Factor but a lower maximum $FI{P}_{FS}$ compared to Basal 2, which has a lower Schmid Factor but a higher maximum $FI{P}_{FS}$ . These findings reinforce that while Schmid factors provide insight into slip system activation, fatigue crack initiation is governed by multiple factors, including local stress states, neighbouring grain interactions, microstructural morphology, proximity to free surfaces and prior strain history.

Figure 13.

Distribution of critical slip planes as a function of the Schmid Factor.

Conclusion

A novel approach for the analysis of fatigue initiation and its probability distribution was presented for a machined Ti-6Al-4V fan-disc part in its final forged and heat-treated microstructure condition. The digital reconstruction of the microstructural model consisted of a series of steps to integrate experimental data into a CP finite element model where machining induced effects (i.e., accumulated plastic strain, surface roughness and residual stress) were superimposed with EBSD orientation information. The model was subjected to ten fatigue cycles and the following conclusions are drawn from the results:

The surface topology (i.e., roughness) has a significant impact on fatigue initiation due to stress concentration sites, which cause plastic strain accumulation and heterogeneous stress fields at the surface and at near-surface sites.

The total accumulated plastic strain, which is used as a fatigue initiation parameter, is not appropriate to derive initiation sites when considering previous effects from manufacturing process steps on top of cyclic fatigue loading. The $FI{P}_{FS}$ formulation is shown to be adequate to calculate initiation sites and their likelihood since its expression is a function of the maximum shear strain increment among all slip systems in a cycle.

The simulation captures ratcheting and predicts that this phenomenon is one of the key drivers of plastic accumulation which could cause fatigue crack initiation. Furthermore, the velocity in which plasticity is accumulated is substantially higher at the surface and near-surface regions and decreases rapidly with the distance from the surface into the bulk.

Slip system activation in the model is anisotropic and the WCSS is dominated by the basal and prismatic slip systems. The simulation indicates that the pyramidal slip system provides a secondary mechanism for plasticity accommodation.

The $FI{P}_{FS}$ critical planes mostly belong to the basal or prismatic crystallographic slip families. Pyramidal slip systems are less frequently critical planes, and they have lower $FI{P}_{FS}$ values. The contribution of extreme outliers indicates that the $FI{P}_{FS}$ is highly localised in Basal 1 and Prismatic 1. Prismatic 2 and 3 have a more evenly distributed fatigue indicator parameter throughout the microstructure.

Analyses of extreme value statistics shows that Basal 1, Basal 2, and Prismatic 3 have the heaviest-tailed distributions, indicating a higher probability of localised extreme plasticity and potential fatigue crack initiation. The Schmid Factor distributions suggest a strong alignment between Prismatic 1 and 2 with the loading direction, reinforcing their high activation likelihood. However, the lack of direct correlation between Schmid factor distributions and extreme $FI{P}_{FS}$ values highlights that fatigue crack initiation is influenced not only by crystallographic alignment but also by local stress states, grain interactions, and surface-induced geometric constraints.