Modeling and Simulation of Fatigue Crack Initiation Process Based on Field Theory of Multiscale Plasticity (FTMP) Part I: PSB Ladder Formation and Verification

1. Introduction

Research into metal fatigue has a long-standing history, with significant contributions dating back over half a century. Early investigations by Thompson et al. [1] revealed that fatigue cracks initiate within slip bands of a single grain, highlighting the critical role of microscopic dislocation structures in fatigue crack formation. Building on this foundation, Woods et al. [2] emphasized the presence of persistent slip bands (PSBs), characterized by uniformly arranged dislocation walls, often referred to as ladder walls. Polák et al. [3,4,5,6] further explored the relationship between these dislocation substructures and the resultant surface morphology. Essmann et al. [7,8] conducted extensive studies on irreversible deformation associated with the ladder structure of PSBs and its influence on crack initiation at the PSB-matrix interface. They examined the crucial role of vacancies commonly found in fatigued metals, [9,10] including their agglomeration into crack nucleation sites. Subsequent studies along similar lines have continued to advance our understanding of the relationship between PSB ladder structures and the processes of fatigue crack initiation, an inquiry that remains ongoing.

The significance of deformation-induced patterning in plasticity is paramount for understanding the multiscale nature of material mechanics. This importance arises from the dynamic changes in dislocation substructures influenced by long-range stress fields and misorientation, particularly in materials that yield dislocation structures with cellular morphology [11,12]. These dislocation patterns play a crucial role in the storage and release of elastic strain energy, influencing not only the deformation capacity but also the strength and toughness of materials. A notable exception is the ladder structure observed in PSBs, which predominantly facilitates deformation during cyclic loading and typically does not produce long-range stress field due to its composition of edge dislocation dipoles [13]. This PSB ladder structure significantly impacts the mechanisms of fatigue crack initiation, as described above [14,15,16,17,18]. However, the underlying reasons for this patterning and its critical roles remain undefined, partly due to limitations in conventional approaches to plasticity.

While analogy-based methods, such as reaction-diffusion equations grounded in the concept of dissipative structures, can address patterning-related issues, they often provide limited insights beyond mere analogy [19,20]. These methods do not elucidate the mechanisms involved in crack initiation. Recent research, fueled by advances in computational capabilities, has increasingly relied on dislocation dynamics (DD) [21,22,23] and molecular dynamics (MD) simulations [24], employing multiscale methods to model these phenomena. However, these studies often focus narrowly on specific aspects, such as passing stress [13,25], rather than capturing the complete picture that includes the roles of PSB ladders and the processes of crack nucleation e.g., [24,26]

Modeling fatigue presents significant challenges, particularly concerning crack initiation, as it must consider the effects of dislocation substructures, that is, the laddered morphology in PSBs. The foundational work by Repetto and Ortiz [27] on finite element simulations of fatigue crack initiation emphasized the role of vacancy-induced elongation in the PSB region, which leads to surface protrusion as a key factor; however, their study did not account for the laddered morphology explicitly. In contrast, Nakai et al. [28,29,30] reported that cracks primarily evolve from intrusions rather than extrusions. Their research identified a potential critical condition for the transition from intrusions to cracks, based on detailed observations made using atomic force microscopy (AFM) combined with interrupted fatigue tests. According to their findings, extrusions and/or surface reliefs do not develop into cracks. Through their systematic series of i studies on several materials, they successfully identified a potential critical condition for the transition from intrusions to cracks, with further details presented in Section 2.

Focusing on the transition from slip bands to cracks, particularly once this critical condition is reached, may provide a more effective framework and minimize ambiguities in the modeling process. Integrating this critical condition with our Field Theory of Multiscale Plasticity (FTMP)-based approach [31,32,33,34,35] can significantly reduce computational costs, as FTMP does not always require intricate fine mesh divisions to reproduce deformation-induced substructural patterns accurately. This versatile theoretical framework facilitates the study of deformation-induced patterning, which is crucial for controlling mechanical properties across multiscale contexts [31,32].

This study aims to elucidate the roles of the laddered structure in PSBs on the fatigue crack initiation process through FTMP-incorporated crystal plasticity (CP)-based finite element method (FEM) simulations. Part I of the paper presents a preliminary series of simulations utilizing a simplified vacancy model to examine surface roughening behavior and the subsequent grooving processes.

2. Background

Two pivotal research findings have catalyzed the current series of studies. First, Yokoi et al. [36,37] demonstrated that the addition of copper (Cu) to steels can dramatically alter dislocation substructures, transforming them from cellular to planar configurations. This transformation results in finer surface intrusions and extrusions, significantly delaying the initiation of fatigue cracks, as illustrated in Figure 1(a). This finding underscores the crucial role of dislocation substructures in influencing the fatigue crack initiation process.

Second, Nakai et al. [28,29,30]. investigated the transition process from slip bands to cracks through a systematically designed series of interrupted AFM observations. They revealed that surface grooving—an evident manifestation of slip bands—evolves into fatigue cracks when the slip distance (or groove depth) exceeds a critical threshold, which varies among different materials, as shown in Figure 1(b). This finding highlights a key mechanism in dictating the transition from slip banding to fatigue cracking. Notably, surface reliefs (bottom left) do not contribute to crack formation; only grooves serve as precursors to crack nucleation, at least according to their studies. In this context, the effect of Cu addition on the screw core structure of α-Fe has been studied separately through ab initio calculations [38]. Cross slip by screw dislocations is considered critical for dynamic recovery, which drives cell formation. We reported a transition of the core structure from isotropic to fully extended, a change that tends to inhibit cross slip, thereby hindering dynamic recovery and the resulting cell formation. This provides one possible explanation for why the addition of Cu leads to vein-like or planar dislocation substructures, rather than cellular morphology.

By explicitly simulating these critical aspects through deformation analyses, we aim to provide groundbreaking insights that enhance our understanding of the fatigue crack initiation process in practical applications. For the former, the case of copper solid solution yielding a 2D vein structure represents the most fundamental scenario. Here, we confine our focus to the widely observed PSB ladder structure, with other matrix morphologies left for future investigation.

3. Theory and Kinematics

3.1. Brief Description About FTMP

The Field Theory of Multiscale Plasticity (FTMP) [31,32,33,34,35] is a comprehensive framework encompassing three critical aspects: (a) evolution, (b) description, and (c) cooperation, as illustrated in Figure 2 (top center). It primarily addresses “inhomogeneously evolving” deformation fields, which typically manifest as dislocation substructures (top left: simulated examples). Notably, the theory’s ability to describe these aspects of evolution distinguishes it from many others. Furthermore, the interaction formalism facilitates explicit treatments of (c) cooperation across multiple scales.

For the mathematical description of these inhomogeneous fields, differential geometric concepts are predominantly employed, particularly the curvature tensor $R_{klm}^{ \mathrm{...}n}$ and the torsion tensor $S_{kl}^{\mathrm{..}j}$, defined as:

$$\left\{\begin{array}{l}{S}_{kl}^{\mathrm{..}j}={\Gamma }_{[kl]}^j\hfill \\ R_{klm}^{ \mathrm{...}n} =2\left[{\partial }_{[k}{\Gamma }_{l]m}^n+{\Gamma }_{[k\left|p\right|}^n{\Gamma }_{l]m}^p\right]\hfill \end{array}\right.$$

(1)

where ${\Gamma }_{ij}^k$ represents the coefficient of connection. These tensors correspond to the incompatibility tensor (bottom left in Figure 2) and the dislocation density tensor, respectively, drawing on the non-Riemannian plasticity framework advocated by K. Kondo [39,40]. Both can be contracted into second-rank tensors without loss of information by applying the alternating tensor ${\in }_{ijk}$ as follows:

$$\left\{\begin{array}{l}{\alpha }_{ij} ={\displaystyle \frac{1}{2}}{\in }_{ikl}{S}_{kl}^{\mathrm{..}j}\hfill \\ {\eta }_{ij} ={\displaystyle \frac{1}{4 g}}{\in }_{ikl}{\in }_{jmn}{R}_{klm}^{ \mathrm{...}n}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left(g=\det (g_{ij})\right)\hfill \end{array}\right.$$

(2)

with $g_{ij}$ the metric tensor. The incompatibility tensor ${\eta }_{ij}$ is further defined as the double curl of the plastic strain tensor${\epsilon }_{ij}^p$, while the dislocation density tensor ${\alpha }_{ij}$ is derived from the curl of the plastic distortion tensor${\beta }_{ij}^p$.

$$\left\{\begin{array}{l}{\alpha }_{ij} =-{\in }_{ikl}{\partial }_k{\beta }_{lj}^p\hfill \\ {\eta }_{ij} ={\in }_{ikl}{\in }_{jmn}{\partial }_k{\partial }_m{\epsilon }_{ln}^p\hfill \end{array}\right. \mathrm{or} \left\{\begin{array}{l}\alpha =-\nabla \times {\beta }_{‌}^p\hfill \\ \eta =\nabla \times \nabla \times {\epsilon }_{ln}^p\hfill \end{array}\right.$$

(3)

Consequently, strain gradients are intrinsically integrated into the theory. It is important to note that the dislocation density tensor, which represents the first gradient of strain, encompasses the concept of geometrically necessary dislocations (GNDs) [33,35,41].

At the core of FTMP is the concept known as the “flow-evolutionary law,” which serves as a working hypothesis (top right) [34,35]:

$${\eta }_{ij} = \kappa \delta T_{ij}$$

(4)

where $\kappa$ refers to the duality coefficient. This relationship connects the incompatibility tensor ${\eta }_{ij}$ with the fluctuation of the energy-momentum tensor $\delta T_{ij}$, where the fluctuation is defined as the deviation from the spatial average, i.e., $\delta T_{ij} \equiv T_{ij}-〈T_{ij}〉$.

The FTMP framework extensively utilizes the incompatibility tensor, Equation (3)2, among other elements. Its enhanced applications, presented in Figure 3, are categorized into three aspects: an extended definition to four-dimensional (4D) spacetime, the flow-evolutionary hypothesis (Equation (4)), and the interaction formalism applicable across multiple scales. The 4D-extended definition of Equation (3)₂ is expressed as:

$${\eta }_{ij} = {\in }_{iklp}{\in }_{jmnp}{\partial }_k{\partial }_m{\epsilon }_{ln}^p$$

(5)

In this context, the indices indicated by lowercase letters denote four-dimensional spacetime (1, 2, 3 for spatial dimensions and 4 for time), while capital letters are used for spatial components when necessary for clarity. The pure temporal component of the 4D-extended incompatibility tensor ${\eta }_{44}$ results in the spatial trace, represented as:

$${\eta }_{44} = {\in }_{4klp}{\in }_{4mnp}{\partial }_k{\partial }_m{\epsilon }_{ln}^p ={\eta }_{KK} =t{r}^{sp}\eta$$

(6)

Thus, the pure temporal component of Eq.(1) becomes:

$${\eta }_{KK} = \kappa \delta \left(U^e+K\right)$$

(7)

where $U^e$ denotes the elastic strain energy and $K$ represents the kinetic energy. In static conditions, as is often the case, this simplifies to:

$${\eta }_{KK} = \kappa \delta U^e$$

(8)

as presented in the bottom right of Figure 1. Equation (8) visualizes the dynamic interrelationships between excessive strain energy being converted or redistributed into the incompatibility-related degrees of freedom, which drive the field evolutions during the course of elasto-plastic deformation.

In practice, there is no need to directly solve Equation (8). Instead, we can incorporate the incompatibility-based underlying degrees of freedom into the hardening law of the constitutive equation used in CP-FEM simulations, as detailed separately below (also see the bottom right of Figure 2). This integration allows the system to autonomously manage excessive storage of elastic strain energy, leading to significant field evolutions, such as substructure formation. In this context, the evolution of dislocation substructures arises naturally from system-wide accommodations facilitated by incompatibility. If this approach is correct, the necessary dislocation patterning can spontaneously emerge, provided that we properly integrate the incompatibility-related degrees of freedom into the simulation; otherwise, it may not occur. Thus, we anticipate that the laddered patterning in PSBs can be effectively reproduced without the need for any ad hoc models, underscoring the robustness of the FTMP framework.

The divergence free condition for the incompatibility tensor in the 4D space-time reads:

$$di{v}^{4D}\eta = 0\iff {\dot{\eta }}_{44} +div{\eta }_{A4} =0$$

(9)

Here, ${\dot{\eta }}_{44}= {\dot{\eta }}_{KK}$ according to Eq.(3), while the space-temporal mixed component ${\eta }_{4A}$ is given as:

$${\eta }_{A4} = -{\left({\alpha }_{BC} -{\alpha }_{CB}\right)}_{‌}^{·}$$

(10)

Ultimately, we establish the following relationship between the incompatibility rate and the edge dislocation density flux by substituting E.(10) into Eq.(9).

$${\dot{\eta }}_{KK} = div{\left({\alpha }_{BC} -{\alpha }_{CB}\right)}_{‌}^{·}$$

(11)

This relationship allows multiple interpretations regarding pair-wise dislocation creation and annihilation processes, making it useful for modeling various detailed processes related to dislocation fields. The present study utilizes this model to investigate vacancy sources arising from the pair annihilation of edge dislocations in Part II.

3.2. Kinematics of Crystal Plasticity

The FTMP concepts can be integrated with the conventional kinematics framework of crystal plasticity [42,43], serving as a foundational vehicle for their application. We essentially adhere to this framework, beginning with Lee’s elastic-plastic decomposition of the deformation gradient tensor to accurately capture the finite deformation conditions:

$$F = F^{*}\cdot F^p$$

(12)

where $F^{*}$ accounts for elastic deformation, including lattice rotation and $F^p$ serves as its plastic counterpart.

The constitutive equation for elasticity in a rate form, suitable for finite crystal plasticity, is expressed as:

$${({\tau }_{(J)}^{*})}^{○}=C{‌}^e:d^{*}\mathrm{with} {({\tau }_{(J)}^{*})}^{○} ={\sigma }^{○}+\sigma \cdot trd$$

(13)

where ${({\tau }_{(J)}^{*})}^{○}$ is the Jaumann rate of the Kirchhoff Stress tensor, viewed from an observer on the rotating lattice. $\sigma$ expresses the Cauchy stress tensor, $C{‌}^e$ is the elasticity tensor, and $d^{*}$ represents the elastic part of the strain rate tensor $d$. The relation $\dot{J} = trd = tr{d}^{*}$ has been used in the second equation, with $\dot{J}$ denoting the Jacobian rate.

Ultimately, we reach:

$${\tau }_{(J)}^{○} =C^e:d-{\displaystyle \sum _{\alpha =1}^N{R}^{\left(\alpha \right)}{\dot{\gamma }}^{\left(\alpha \right)} }\mathrm{with} \left\{\begin{array}{l}\hspace{0.17em}{R}^{\left(\alpha \right)}=C^e : P^{\left(\alpha \right)} +{\beta }^{\left(\alpha \right)}\hfill \\ \hspace{0.17em}{\beta }^{\left(\alpha \right)} =W^{\left(\alpha \right)}\cdot \sigma -\sigma \cdot W^{\left(\alpha \right)}\hfill \end{array}\right.$$

(14)

Here, $P^{\left(\alpha \right)} ={\left(s^{\left(\alpha \right)} \otimes m^{\left(\alpha \right)}\right)}_{sym}$ and $W^{\left(\alpha \right)}={\left(s^{\left(\alpha \right)} \otimes m^{\left(\alpha \right)}\right)}_{skew}$ are evaluated using the unit vectors for the slip direction $s^{\left(\alpha \right)}$ and slip plane normal $m^{\left(\alpha \right)}$, with $(\alpha )$ referring to the slip systems.

3.3. Constitutive Equation for CP-FEM

We employ the general constitutive equation proposed in [34,35] based on dislocation dynamics, applicable both to FCC and BCC metals in a wide range of strain rate and temperature.

$$\left\{\begin{array}{l}\hspace{0.17em}\hspace{0.17em}{\dot{\gamma }}^{(\alpha )} = {\dot{A}}_{SR}{\tau }_{‌}^{{*}^{(\alpha )}}\left[\left| {\tau }_{‌}^{{*}^{(\alpha )}}\right|B_{SR}\exp {\left(1-{\left|{\displaystyle \frac{〈{\tau }_{‌}^{{‌}^{(\alpha )}}- {\tau }_{Peierls}^{{*}^{(\alpha )}}〉 - {\Omega }^{(\alpha )}}{K^{(\alpha )}}}\right|}^p\right)}^q+C_{SR} \right]\hfill \\ \hspace{0.17em}\hspace{0.17em}{\tau }_{‌}^{{*}^{(\alpha )}}\equiv 〈{\tau }_{‌}^{{‌}^{(\alpha )}}- {\tau }_{Peierls}^{{*}^{(\alpha )}}〉 - {\Omega }^{(\alpha )}\hfill \end{array}\right.$$

(15)

where ${\tau }_{Peierls}^{{*}^{(\alpha )}}$ in the Macaulay bracket ($〈○〉 \equiv \left(○+\left|○\right|\right)/2$) expresses the critical stress for the Peierls mechanism, further given as,

$${\tau }_{Peierls}^{*(\alpha )}={\widehat{\tau }}_{Peierls}^{*}{\left[1-{\left(B_{SR}^{p-1}\ln {\displaystyle \frac{{\dot{\gamma }}_{0P}}{{\dot{\gamma }}_{‌}^{(\alpha )}}}\right)}^{1/q_P}\right]}^{1/p_P}$$

(16)

For FCC metals, we typically set ${\tau }_{Peierls}^{{*}^{(\alpha )}}=0$, while for BCC metals, we can reasonably assume $C_{SR} = 0$. The present study neglects the back stress, i.e., ${\Omega }^{(\alpha )}=0$, for simplicity. The material parameters used in this series of simulations are compiled in Table 1. The evolution of the drag stress is governed by,

$$\left\{\begin{array}{l}\hspace{0.17em}\hspace{0.17em}{\dot{K}}^{(\alpha )}=Q_{\alpha \beta }H\left(\gamma \right)\left|{\dot{\gamma }}^{(\beta )}\right|\hfill \\ \hspace{0.17em}\hspace{0.17em}{Q}_{\alpha \beta }= {\delta }_{\alpha \beta } + f_{\alpha \kappa }{S}_{\kappa \beta }+ F_k({\alpha }_k^{(\beta )};{\eta }_k^{(\beta )}) \hfill \end{array}\right.$$

(17)

where $H\left(\gamma \right)$ represents the referential hardening moduli, $H\left(\gamma \right) =h_0{\left\{\left(h_0 /n{\tau }_0 \right)\gamma +1\right\}}^{n-1}$ with $\gamma ={\displaystyle \sum _{\alpha }{\gamma }^{\left(\alpha \right)}}$, $f_{\alpha \beta }$ is the interaction matrix for pairwise dislocation reactions, and $Q_{\alpha \beta }$ is the hardening ratio that accounts for additional hardening or softening from the straining/loading history. Here, $F_k({\alpha }_k^{(\beta )};{\eta }_k^{(\beta )})$ collectively represents the strain gradient terms for the dislocation density and incompatibility fields. The explicit forms are given as [34,35],

$$\left\{\begin{array}{l}\hspace{0.17em}\hspace{0.17em}F({\alpha }^{(\alpha )})=k_{p_{\alpha }} {\left({\displaystyle \frac{\left|{\alpha }^{(\alpha )}\right|}{b}}\right)}^{1/2}\hfill \\ \hspace{0.17em}\hspace{0.17em}F({\eta }^{(\alpha )})=\mathrm{sgn}({\eta }^{(\alpha )})\cdot k_{p_{\eta }}{\left({\displaystyle \frac{l_{defect}}{b}}\left|{\eta }^{(\alpha )}\right|\right)}^{1/2} \hfill \end{array}\right.$$

(18)

The parameters include $l_{defect}$ that represents the targeted inhomogeneous fields, which coincides with the Burgers vector b in the case of individual dislocations, along with the coefficients $k_{p_{\alpha }}$ and $k_{p_{\eta }}$. For evaluating the slip system-wise dislocation density and incompatibility fields, the following projections are employed.

$$\left\{\begin{array}{l}\hspace{0.17em}\hspace{0.17em}{\alpha }^{(\alpha )} = \left(t_{‌}^{(\alpha )} \otimes s_{‌}^{(\alpha )}+s_{‌}^{(\alpha )} \otimes s_{‌}^{(\alpha )}\right):\alpha \hfill \\ \hspace{0.17em}\hspace{0.17em} {\eta }^{(\alpha )} = \left(t_{‌}^{(\alpha )} \otimes s_{‌}^{(\alpha )}+s_{‌}^{(\alpha )} \otimes s_{‌}^{(\alpha )}+s_{‌}^{(\alpha )} \otimes m_{‌}^{(\alpha )}\right):\eta \hfill \end{array}\right.$$

(19)

where $t_{‌}^{(\alpha )} = s_{‌}^{(\alpha )} \times m_{‌}^{(\alpha )}$ represents the line direction in the case of dislocations. Note that the choice of projection directions can affect the evolved patterning, requiring careful attention, especially regarding the incompatibility term, whereas the dislocation density tensor has a clear physical meaning, specifically representing the edge and screw components via $t_{‌}^{(\alpha )} \otimes s_{‌}^{(\alpha )}$ and $s_{‌}^{(\alpha )} \otimes s_{‌}^{(\alpha )}$ projections, respectively. All the material parameters used in the present study is listed in Table 1.

4. Model and Analytical Condition

The finite element (FE) model used for the preliminary analyses is illustrated in Figure 4, which also outlines the analytical conditions. This model, oriented along a single slip direction $[1 \overline{5} 2]$ is discretized into 24 × 80 × 4 crossed-triangular elements. The evaluation radius for the strain gradient calculation is set to $r_{GR} =1.0\mu m$, within which the average slope of the strain values is determined using the least squares method. Cyclic straining is applied under conditions controlled by the plastic strain amplitude of $\Delta {\epsilon }^p / 2 =0.6\%$.

5. Results and Discussion

5.1. General Features

Figure 5 presents a simulated incompatibility contour on a representative slip system (α), along with its projected counterpart on the primary slip system, as illustrated in the process below. Here, $m^{(\alpha )}$ represents the slip plane normal. The projection onto the primary slip plane, expressed by $m_{primary}^{(\alpha )}$, effectively cancels the incompatibility pattern in the matrix region, leaving only the ladder-like pattern along the slip band. This suggests that the incompatibility-based patterning can essentially be categorized as a fluctuation, resulting in a net value of zero overall, in general. In contrast, the dislocation density contour, represented by the edge component ${\alpha }_{31}$ shown on the rightmost side of Figure 5, naturally exhibits a ladder-like structure even without this projection. This indicates that the resultant dislocation substructure corresponds to the projected incompatibility distribution.

The reproduced PSB accompanied by laddered patterning, as anticipated above (see Section 3.1), underscores the predictive power of the FTMP framework, without any additional ingredients nor specifically-designated intricate mathematical models, moreover, essentially no need elaborates finite element discretization. This reinforces our confidence in modeling seamlessly the fatigue crack initiation processes that such specific dislocation structure can be accurately reproduced without reliance on ad hoc models.

To demonstrate how the ladder-patterned region is softened, the sample is monotonically pulled after interrupting the cyclic straining. Figure 6a compares the samples pulled before the onset of cyclic straining (0 cycles) and after 80 cycles. As expected, the laddered sample exhibits landslide-like localized deformation, while the 0-cycle sample shows relatively uniform elongation. The results of monotonic tension after various straining cycles are presented in Figure 6b, indicating that 15 cycles of straining lead to complete softening of the PSB.

Figure 7 compares the hysteresis loops for the matrix and PSB regions. The PSB region exhibits a plastic strain amplitude three times larger, highlighting strain localization along the laddered band, with this region accounting for approximately 14% of the total volume. Notably, the loop for the PSB region tends to shift toward the tension side to balance system-wide deformation, ensuring symmetric push-pull straining. It is worth mentioning that this concentration of deformation is less than what is typically observed in reality, primarily because the present simulation does not distinctly capture vein structures in the matrix region, which can impede deformation. However, reproducing the vein morphology in the matrix is not within the scope of this study at this stage.

5.2. Specific Features

Let us examine how the reproduced patterns reflect reality. Figure 8 analyzes the laddered region in terms of component-wise dislocation density distribution, while Figure 9 illustrates the growth of their cross-sectional distributions with increasing straining cycles. Here, we display ${\alpha }_{22}$ and ${\alpha }_{31}$ as representatives for the screw and edge components, respectively, while ${\alpha }_{11}$ and ${\alpha }_{32}$ show similar trends, although it is not presented. The following points are clearly demonstrated: first, the ladder morphology comprises both screw and edge dislocation components; second, these components are alternately aligned, closely resembling experimental observations, as also confirmed in Figure 9. Furthermore, there is essentially no misorientation developed throughout the sample, with a deviation of less than 1 degree, consistent with experimental observations. To be noted that incompatibility-induced substructures are accompanied by misorientation developments in general.

On the less desirable side, we observe a slightly larger amplitude for ${\alpha }_{31}$ than for ${\alpha }_{22}$ in Figure 9, indicating that the structure is less capable of carrying deformation, as the screw segments play a key role in this capacity. This also contributes to the previously mentioned smaller strain concentration, alongside the less-developed vein-like morphology in the matrix region. Additionally, another difference remains in the width ratio of the ladder wall to channel regions, which ought to be much smaller, indicating potential room for further improvement.

5.3. Surface Roughening and Groove Formation

The strain concentration described earlier promotes surface roughening around the banded region, ultimately leading to grooving in that area. Figure 10 illustrates the evolving surface profile with increasing straining cycles, highlighting surface undulation up to N=90 cycles. This local deformation creates a sample-wide biased surface slope, including a protruded left edge, as the system seeks overall equilibrium. Notably, the banded region displays a significantly larger slope, indicative of localized strain that contributes to grooving at its right edge. The magnified view on the left captures the early stage of groove formation at the intersection of the PSB and the surface, where a shallow but distinct groove tip is evident at the right edge of the laddered band. These observations clearly demonstrate that spontaneous substructuring alone can effectively reproduce surface roughening and initiate the formation of groove embryos.

5.4. Preliminary Analysis with Vacancy Effect Using a Simplified Model

To anticipate the effect of the formed pattern on the subsequent behavior of the PSB and the resulting development of surface undulations, we introduce a simplified model for vacancy formation. This serves as a prelude to the explicit treatments of vacancy source model coupled with diffusion analysis that will be explored in Part II. Assuming that the gradient of the elastic strain energy attended with the formed laddered pattern can promote vacancy diffusion, combined with the flow-evolutionary law Equation (8), we introduce the following simplified model via the incompatibility term:

$${\eta }_{cv}= {\kappa }_v \delta U_{‌}^e$$

(20)

The idea here is inspired by a systematic series of diffusion analyses describe in Appendix A1, where the diffusion equation for vacancy concentration is solved against an assumed ladder pattern specified by an a priori elastic strain energy distribution. In this diffusion analyses, diffusion is solely enhanced by the gradient of the elastic strain energy, which is now reflected through the incompatibility term $F({\eta }_{cv})$ in the present analysis. We incorporate this term into the original $F(\eta )$ in the hardening ratio $Q_{\alpha \beta }$, (Equation (17)₂) anticipating a positive contribution to groove formation and further growth.

Figure 11 presents a comparison of incompatibility contours (top) and the corresponding elastic strain energy fluctuations (bottom) ($F(\eta )$ and $\delta U^e$) to assess the effect of the simplified vacancy model on the evolution of laddered patterns in the PSB region. In the absence of the vacancy model (right row), the laddered pattern tends to collapse partially along the righthand side edge of PSB as straining cycles progress, as indicated by the incompatibility contours. This ultimately results in a relatively scattered distribution of $\delta U^e$ around the PSB edge on the sample surface. In contrast, incorporating the vacancy model (left row) appears to stabilize the pattern, leading to a concentration of $\delta U^e$. We anticipate that this contribution to elastic strain energy fluctuations will foster the healthy development of surface grooves.

Figure 12 compares the evolution of surface profiles with and without the contribution of ${\eta }_{cv}$. The inclusion of ${\eta }_{cv}$ leads to the formation of a single groove, whereas its absence results in a splitting of the groove peaks, becoming evident after 200 straining cycles. This behavior corresponds to the concentrated and diffused strain energy distributions seen in Figure 11.

To further investigate the role of ${\eta }_{cv}$ in the final surface grooving, we examine the individual variations of the incompatibility terms, $F(\eta )$ and $F({\eta }_{cv})$, with respect to straining cycles in Figure 13. The term $F({\eta }_{cv})$ itself increases parabolically with the number of cycles while exhibiting periodic oscillations, simply reflecting the applied cyclic straining. Adding $F({\eta }_{cv})$ significantly increases the overall variation of $F(\eta )$ in the negative direction (bottom right), while subtracting it leads to a mild decrease in the overall change (bottom left). This upward trend in $F(\eta )$ enhances local deformation, resulting in a converging and sharpened groove, as confirmed by the final surface profile displayed at the top right of Figure 13.

These results clearly demonstrate that the addition of $F({\eta }_{cv})$ positively promotes grooving by significantly sharpening and deepening the groove profile, as anticipated. This ultimately leads to an accumulation of strain energy at the right edge of the PSB-sample surface intersection, whereas, in the absence of this term, the elastic strain energy tends to diffuse. We will utilize the final surface profile (b) for restart analyses to be conducted in Part II of the present study.

5.5. Discussion and Future Scope

In this study, we have successfully demonstrated that the FTMP-incorporated CP-FEM simulation can effectively reproduce the empirically observed laddered morphology and the PSB. Moreover, this simulation adeptly captures the dynamics of surface roughening and the eventual formation of grooves, even in the absence of vacancy formation and diffusion. These findings underscore the significance of the incompatibility tensor-based underlying degrees of freedom, which possess a remarkable capability to reproduce dislocation substructures—features that are often lacking in conventional CP-FEM approaches. This advancement not only helps enhance our understanding of PSB laddered patterning but also opens new avenues for further exploring the complex interplay of microstructural features that dictate the eventual crack initiation process, as will be extended in Part II.

The current simulation can be seamlessly extended to multi-grain scenarios, as demonstrated in Figure 14. This figure showcases three preliminary examples featuring single slip-oriented grains located in (a) the upper half, (b) the surface-facing grain, and (c) the central grain. In each case, we observe similar PSB-ladder formations exclusively in the single-slip oriented grain as observed previously, which intriguingly terminate at the grain boundaries. These results underscore the potential for expanding the present study to polycrystalline simulations, highlighting the versatility and applicability of our approach.

6. Conclusion

In this study, we first demonstrate that PSB ladder substructuring during cyclic straining can be effectively reproduced through simple deformation analyses, showcasing the predictive capability of the FTMP framework. Notably, this achievement does not require any ad hoc models, intricate mathematical formulations, or elaborate finite element discretization. This reinforces our confidence in the practical feasibility of computational modeling for fatigue crack initiation processes and supports our advancement to the next stage of research without relying on complex modeling approaches.

Furthermore, the FTMP-incorporated CP-FEM simulation not only replicates the empirically observed laddered morphology and PSB with satisfactory accuracy, but also adeptly simulates surface roughening and eventual grooving—independent of vacancy formation and diffusion. This highlights the crucial role of laddered morphology in these processes and lays the groundwork for further investigations into the effects of vacancy formation through explicit diffusion analyses in Part II. These advancements stem from the incompatibility tensor-based underlying degrees of freedom, which offer exceptional capabilities in modeling dislocation substructures often overlooked by conventional CP-FEM approaches. Thus, the robustness of the FTMP framework is validated, positioning it as a transformative tool in materials science and paving the way for deeper insights into slip band-fatigue crack transition.