Local Monte Carlo Method for Fatigue Analysis of Coarse-Grained Metals with a Nanograined Surface Layer

: The fatigue resistance of coarse-grained (CG) metals can be greatly improved by introducing a nanograined surface layer. In this study, the Weibull distribution is used to characterize the spatially-random fracture properties of specimens under axial fatigue. For the cylindrical solid specimen, the heterogeneity of element sizes may lead to unfavorable size effects in fatigue damage initiation and evolution process. To alleviate the size effects, a three-dimensional cohesive ﬁnite element method combined with a local Monte Carlo simulation is proposed to analyze fatigue damage evolution of solid metallic specimens. The numerical results for the fatigue life and end displacement of CG specimens are consistent with the experimental data. It is shown that for the specimens after surface mechanical attrition treatment, damage initiates from the subsurface and then extends to the exterior surface, yielding an improvement in the fatigue life. Good agreement is found between the numerical results for the fatigue life of the specimens with the nanograined layer and experimental data, demonstrating the efﬁcacy and accuracy of the proposed method.


Introduction
Since the fatigue resistance of metal is strongly dependent on its surface properties, much attention has been paid to the surface modification [1]. Metals with finer grains usually have greater resistance to fatigue crack initiation. The way that microstructures affect the growth of both small and long fatigue cracks can be captured by changing the fatigue threshold term in the Hartman-Schijve variant of the NASGRO crack growth equation [2]. For small fatigue cracks, both fine-grained metals and coarse-grained (CG) metals have a small fatigue threshold and thus a similar crack growth rate [3,4]. Equivalently, microstructural effects often vanish for small, naturally-occurring cracks [3,4]. For long fatigue cracks, compared with fine-grained metals, CG metals possess a larger fatigue threshold and thus a lower crack growth rate [3,4]. Therefore, metals with a nanograined surface layer and coarse-grained (CG) interior hold promise for applications against fatigue failure. At present, 2 of 15 the surface mechanical attrition treatment (SMAT) technique uses a large quantity of smooth balls to hit the specimen surface so that such a surface layer can be created [5]. It can achieve a nanograined layer (NGL) in various metals including Inconel 600 alloy [6], 304 stainless steel (SS) [7], Fe [8], and Co [9]. Dewo et al. refined the grain size of osteosynthesis plates' surfaces by SMAT so that the plates had a significant improvement in strength [10]. Recently, Sun and Bailey found that much of the improved corrosion behavior of the SMATed 304 SS was attributed to the NGL [11]. The same phenomenon was also observed in the SMATed Pb-Sn alloy [12]. Kargar et al. found that the SMATed Cu specimen exhibited greater hardness [13]. Furthermore, Jin et al. fabricated NGL on biomedical TiNbZrFe alloy by SMAT and found that the micro-hardness of the surface layer was enhanced with treatment duration [14]. Huang et al. found the Ti alloy grains in the NGL generated by SMAT had good thermal stability [15]. The SMAT was employed to generate a gradient structure which exhibited high yield strength and good uniform elongation simultaneously [16,17]. So far, with the emergence of SMATed materials, the improved mechanical properties of such nanostructured metals have been well investigated; the NGL also has a beneficial effect on the fatigue life.
A number of experiments have been conducted to investigate the effects of the NGL on the fatigue behavior of SMATed metallic specimens. This is of importance for the practical applications of nanostructured metals. For instance, Kumar et al. found that the SMATed Ti-6Al-4V specimens treated for 30 min exhibited superior fatigue strength compared to untreated specimens and specimens treated for 60 min [18]. From the viewpoint of engineering, it is important to study the failure mode of fatigue and the approach to improve the corresponding fatigue strength. However, the combination of limited facilities, high cost, and the long period of experimental tests has placed a significant constraint on fatigue study. To meet this challenge, here we perform a systematic numerical investigation and a brief experimental verification. A valid three-dimensional (3D) numerical investigation on the axial fatigue life and damage evolution can provide valuable insights.
In this paper, we will use the cohesive finite element method (CFEM), which has proven to be effective in studying the fatigue failure of structural materials. For instance, Brinckmann and Siegmund adopted it to study the fatigue crack growth in metals successfully [19]. Gong et al. implemented the CFEM to uncover the damage evolution of a double cantilever beam [20]. Ghovanlou et al. applied the CFEM to investigate the fatigue crack growth of the brazed joints and showed that the results agreed well with the experiments [21]. Using the CFEM, Puigvert et al. well predicted the fatigue behavior of adhesively-bonded anchorages [22]. Kim and Yoon employed the CFEM to model the high-cycle fretting fatigue of Al alloy [23]. Furthermore, the CFEM was used to analyze the plasticity-induced crack closure of aluminum alloy and numerical analysis was correlated with the fatigue tests [24]. As a common phenomenon, fatigue test data tend to scatter widely due to internal defects, unavoidable wear, and processing errors. Experiments revealed that test data scattered for nominally identical specimens under the same fatigue load. Therefore, it is necessary to perform probabilistic or statistical analysis to fully understand the fatigue problems. To this end, we will also take a statistical approach in this work.
On the one hand, the damage tolerance method in Ref. [25,26] treats the defects as growing cracks under various load histories in service. On the other hand, statistical distributions can describe the internal grain microstructure or material defects [27]. In recent years, Monte Carlo simulation (MCS) has been widely used to investigate the fracture properties of the material and it is an effective approach to simulate the accumulation of large quantities of random incidents. It uses random variables to establish a simple but effective statistical model according to the statistical law of practical problems. For instance, Bouhamida et al. combined MCS to study the distribution of stress concentration for each fiber [28]. Campanile et al. applied MCS with Weibull distribution to investigate statistical properties of the hull girder section modulus and ultimate strength [29]. Based on the scatter of concrete tensile strength and bond strength, Hartig and Häussler-Combe employed MCS to study the statistics of crack widths [30]. It is quite difficult for the MCS to deal with the high nonlinearity related to the fracture process because of the high computational cost associated with a large number of finite elements [31]. However, it is still an effective method to estimate the uncertainty of structural response. So far, only limited studies have combined MCS with CFEM to study fracture behavior [31,32]. The MCS can be employed to study the fatigue failure process [33,34]. So far, 2D models employed by most of the researchers have not been able to fully capture the complexity in the intrinsic 3D phenomena in the fatigue failure. It is also the intent here that we will undertake a 3D analysis.
Recently, Guo et al. combined the CFEM with the MCS to simulate the fatigue crack propagation of hollow tubular specimens under torsional load successfully [35]. Furthermore, Sun et al. employed this approach to study the fatigue life and damage evolution of hollow tubular specimens subjected to combined axial-torsional loads [36]. In the above two papers, the global MCS (i.e., the fracture properties of all cohesive elements are considered to be spatially random and thus characterized by MCS) is applicable to the hollow tubular specimens. Here, a cylindrical solid is used in our axial-fatigue experiments and simulation; the heterogeneity of the element sizes will cause unfavorable results on fatigue damage initiation and evolution. In real materials, the length scale for variation in the local resistance to fatigue damage initiation and evolution generally scales with the grain size or other microstructural parameter. Indeed, the meshes in our cylindrical specimens have a hub-and-spoke topology and all hexahedron solid elements have the same central angle. This suggests that the "random" fracture properties vary much more rapidly at the interior of the cylinder compared to the exterior surface. Therefore, in order to reduce the element size effects, we performed a local MCS (i.e., the fracture properties of local/partial cohesive elements are considered to be spatially random and thus characterized by MCS, while the other elements are homogenized).
In this study, we will develop the local MCS together with a 3D CFEM model to simulate axial fatigue life and damage evolution. The geometry of the specimen is shown in Figure 1. It is based on our recent fatigue experiments conducted via BOSE ElectroForce 3550 (Bose Corporation ElectroForce Systems Group 10250 Valley View Road, Eden Prairie, MN, USA) under sinusoidal load with the frequency of 10 Hz and stress ratio of 0.1. In Section 2, the CFEM, together with geometric modeling and quantization of the constitutive parameters, will be briefly specified, and the Weibull cumulative distribution function will be introduced to represent the statistical fracture properties. In Section 3, the effects of the random fields and loads on the fatigue life and damage evolution of both the CG and the SMATed specimens will be investigated. The main findings will be summarized in Section 4. study fracture behavior [31,32]. The MCS can be employed to study the fatigue failure process [33,34]. So far, 2D models employed by most of the researchers have not been able to fully capture the complexity in the intrinsic 3D phenomena in the fatigue failure. It is also the intent here that we will undertake a 3D analysis.
Recently, Guo et al. combined the CFEM with the MCS to simulate the fatigue crack propagation of hollow tubular specimens under torsional load successfully [35]. Furthermore, Sun et al. employed this approach to study the fatigue life and damage evolution of hollow tubular specimens subjected to combined axial-torsional loads [36]. In the above two papers, the global MCS (i.e., the fracture properties of all cohesive elements are considered to be spatially random and thus characterized by MCS) is applicable to the hollow tubular specimens. Here, a cylindrical solid is used in our axialfatigue experiments and simulation; the heterogeneity of the element sizes will cause unfavorable results on fatigue damage initiation and evolution. In real materials, the length scale for variation in the local resistance to fatigue damage initiation and evolution generally scales with the grain size or other microstructural parameter. Indeed, the meshes in our cylindrical specimens have a hub-andspoke topology and all hexahedron solid elements have the same central angle. This suggests that the "random" fracture properties vary much more rapidly at the interior of the cylinder compared to the exterior surface. Therefore, in order to reduce the element size effects, we performed a local MCS (i.e., the fracture properties of local/partial cohesive elements are considered to be spatially random and thus characterized by MCS, while the other elements are homogenized).
In this study, we will develop the local MCS together with a 3D CFEM model to simulate axial fatigue life and damage evolution. The geometry of the specimen is shown in Figure 1. It is based on our recent fatigue experiments conducted via BOSE ElectroForce 3550 (Bose Corporation ElectroForce Systems Group 10250 Valley View Road, Eden Prairie, MN, USA) under sinusoidal load with the frequency of 10 Hz and stress ratio of 0.1. In Section 2, the CFEM, together with geometric modeling and quantization of the constitutive parameters, will be briefly specified, and the Weibull cumulative distribution function will be introduced to represent the statistical fracture properties. In Section 3, the effects of the random fields and loads on the fatigue life and damage evolution of both the CG and the SMATed specimens will be investigated. The main findings will be summarized in Section 4.

Cohesive Finite Element Method
The CFEM is formulated via the principle of virtual work [37]. The internal work by the virtual strain and that by the virtual crack opening displacement equal the external work by the virtual displacement. The CFEM has been successfully applied to predict the damage initiation and evolution explicitly. Furthermore, it can model spontaneous multiple crack behavior, so it has been widely used to investigate damage process [38][39][40][41][42][43]. Here, our interest is in the high nonlinearity involving the interaction of multiple microcracks, and thus the CFEM is particularly appealing.
The cohesive law, i.e., the traction-separation relation, can describe the cohesion of the material in the CFEM. Thereby, the damage initiation and evolution together with the subsequent cracking process in the fatigue and fracture failure is a natural outcome of the computational framework for a specimen with or without preexisting cracks. Our experiments showed that the tensile fatigue-

Cohesive Finite Element Method
The CFEM is formulated via the principle of virtual work [37]. The internal work by the virtual strain and that by the virtual crack opening displacement equal the external work by the virtual displacement. The CFEM has been successfully applied to predict the damage initiation and evolution explicitly. Furthermore, it can model spontaneous multiple crack behavior, so it has been widely used to investigate damage process [38][39][40][41][42][43]. Here, our interest is in the high nonlinearity involving the interaction of multiple microcracks, and thus the CFEM is particularly appealing.
The cohesive law, i.e., the traction-separation relation, can describe the cohesion of the material in the CFEM. Thereby, the damage initiation and evolution together with the subsequent cracking process in the fatigue and fracture failure is a natural outcome of the computational framework for a specimen with or without preexisting cracks. Our experiments showed that the tensile fatigue-fractured specimen had predominant mode-I cracks. Therefore, we use the following quadratic stress criterion of damage initiation [32,41]: where is the Macaulay bracket, defined as where T n , T s , and T t represent the normal and two shear tractions, respectively. T 0 n denotes the maximum of the nominal normal traction, and T 0 s (T 0 t ) denotes the maximums of the nominal shear traction, as shown in Figure 2a fractured specimen had predominant mode-I cracks. Therefore, we use the following quadratic stress criterion of damage initiation [32,41]: where   is the Macaulay bracket, defined as where , , and represent the normal and two shear tractions, respectively. denotes the maximum of the nominal normal traction, and ( ) denotes the maximums of the nominal shear traction, as shown in Figure 2a  The damage evolution can be characterized by the scalar stiffness degradation (SDEG), D, representing the overall damage of the cohesive element caused by all active mechanisms. It depends on the effective separation combining , , and as . ( For the bilinear cohesive law in Figure 2, the damage evolves according to where is the maximum effective separation attained during the loading history. and are the characteristic separation where the damage initiates and critical separation where the cohesive element fails, respectively. Equation (4) indicates that D monotonically evolves from 0 to 1 upon further loading after damage initiation [44].

3D CFE Geometric Modeling and Quantification of Constitutive Parameters
The solid dumbbell-shaped specimens are used in our axial fatigue experiments. As illustrated in Figure 3a, the cracks in the fractured specimen propagate from the exterior surface to the interior. In the present paper, therefore, a 3D CFE geometric model is established by using an embedding algorithm. The approach contains the following three main steps: (i) Mesh the circular cross section of the 3D specimen. The initial mesh consists of triangular elements in the core location and quadrilateral elements in the surrounding area; The damage evolution can be characterized by the scalar stiffness degradation (SDEG), D, representing the overall damage of the cohesive element caused by all active mechanisms. It depends on the effective separation u m combining u n , u s , and u t as For the bilinear cohesive law in Figure 2, the damage evolves according to where u max m is the maximum effective separation attained during the loading history. u 0 m and u f m are the characteristic separation where the damage initiates and critical separation where the cohesive element fails, respectively. Equation (4) indicates that D monotonically evolves from 0 to 1 upon further loading after damage initiation [44].

3D CFE Geometric Modeling and Quantification of Constitutive Parameters
The solid dumbbell-shaped specimens are used in our axial fatigue experiments. As illustrated in Figure 3a, the cracks in the fractured specimen propagate from the exterior surface to the interior. In the present paper, therefore, a 3D CFE geometric model is established by using an embedding algorithm. The approach contains the following three main steps: (i) Mesh the circular cross section of the 3D specimen. The initial mesh consists of triangular elements in the core location and quadrilateral elements in the surrounding area; (ii) Generate 3D solid meshes. Read the nodal coordinates and the nodal connectivity of the 2D elements, and then renumber them to establish the index between each new element and the corresponding nodal number. The 2D indices can be extended to 3D so that wedge elements (C3D6) in the core location and hexahedron elements in the surrounding area can be generated; (iii) Insert cohesive elements. Partition every hexahedron element into four wedge elements along the left and right surfaces, as shown in Figure 3b. One node is replaced by eight separate nodes at the same position at the diagonal center, and four 8-noded cohesive elements (COH3D8 in [44]) and four 6-noded cohesive elements (COH3D6) are generated. The right view of this insertion model is also presented in Figure 3b. COH3D8 can be inserted between the frontal and rear surfaces of the wedge elements, COH3D6 and COH3D8 can be inserted between the left and right surfaces of the wedges and between the upper and lower surfaces of the wedge elements, respectively. (ii) Generate 3D solid meshes. Read the nodal coordinates and the nodal connectivity of the 2D elements, and then renumber them to establish the index between each new element and the corresponding nodal number. The 2D indices can be extended to 3D so that wedge elements (C3D6) in the core location and hexahedron elements in the surrounding area can be generated; (iii) Insert cohesive elements. Partition every hexahedron element into four wedge elements along the left and right surfaces, as shown in Figure 3b. One node is replaced by eight separate nodes at the same position at the diagonal center, and four 8-noded cohesive elements (COH3D8 in [44]) and four 6-noded cohesive elements (COH3D6) are generated. The right view of this insertion model is also presented in Figure 3b. COH3D8 can be inserted between the frontal and rear surfaces of the wedge elements, COH3D6 and COH3D8 can be inserted between the left and right surfaces of the wedges and between the upper and lower surfaces of the wedge elements, respectively. To save computational cost, we consider one sixth of the test section in the axial direction and remove fixture part and transitional arcs of the specimen. Figure 4 shows the CG CFE model, while only the highlighted elements are characterized by the MCS. The number of the MCSed and all discretized elements of CG models are listed in Table 1. Note that the MCSed elements refer to the elements where the MCS is performed.  To save computational cost, we consider one sixth of the test section in the axial direction and remove fixture part and transitional arcs of the specimen. Figure 4 shows the CG CFE model, while only the highlighted elements are characterized by the MCS. The number of the MCSed and all discretized elements of CG models are listed in (ii) Generate 3D solid meshes. Read the nodal coordinates and the nodal connectivity of the 2D elements, and then renumber them to establish the index between each new element and the corresponding nodal number. The 2D indices can be extended to 3D so that wedge elements (C3D6) in the core location and hexahedron elements in the surrounding area can be generated; (iii) Insert cohesive elements. Partition every hexahedron element into four wedge elements along the left and right surfaces, as shown in Figure 3b. One node is replaced by eight separate nodes at the same position at the diagonal center, and four 8-noded cohesive elements (COH3D8 in [44]) and four 6-noded cohesive elements (COH3D6) are generated. The right view of this insertion model is also presented in Figure 3b. COH3D8 can be inserted between the frontal and rear surfaces of the wedge elements, COH3D6 and COH3D8 can be inserted between the left and right surfaces of the wedges and between the upper and lower surfaces of the wedge elements, respectively. To save computational cost, we consider one sixth of the test section in the axial direction and remove fixture part and transitional arcs of the specimen. Figure 4 shows the CG CFE model, while only the highlighted elements are characterized by the MCS. The number of the MCSed and all discretized elements of CG models are listed in Table 1. Note that the MCSed elements refer to the elements where the MCS is performed.   For the SMATed specimen, we choose an NGL with a thickness of 20 µm. The modeling approach is similar to that for the CG specimen. Specifically, the circumferential and axial dimensions of NGL elements are the same as those of CG elements, while the thickness of the NGL elements (10 µm) is much thinner than that of the nearby CG elements (~180 µm). The corresponding geometric model of the SMATed specimen is illustrated in Figure 5. The number of the MCSed and discretized elements of SMATed models are listed in Table 1.  For the SMATed specimen, we choose an NGL with a thickness of 20 µm. The modeling approach is similar to that for the CG specimen. Specifically, the circumferential and axial dimensions of NGL elements are the same as those of CG elements, while the thickness of the NGL elements (10 µm) is much thinner than that of the nearby CG elements (~180 µm). The corresponding geometric model of the SMATed specimen is illustrated in Figure 5. The number of the MCSed and discretized elements of SMATed models are listed in Table 1. A bilinear cohesive law has two important parameters, namely the cohesive strength and fracture energy . The fracture energy of the CG 304 SS is obtained from 1 / in terms of its fracture toughness , which is taken as 75 MPa√m [45]. For the NGL, is taken as 36 MPa√m. Young's modulus E and Poisson's ratio v of both the CG and the NGL are taken as 200 GPa and 0.33, respectively [46]. The flow stress of the CG specimen is from the experiments [46], while that of the NGL is estimated from experiments of the SMATed 316L SS [47]. The kinematic hardening model is used to describe the cyclic plastic behavior of the metals. Additionally, the initial stiffness k should be large enough to alleviate the inherent mesh dependence, while too large value may cause numerical ill-conditioning [38,41]. It is determined as 2 10 N/m.

Weibull Random Field
Based on the weakest link model, the Weibull distribution assumes that strength is connected with flaws randomly distributed in the materials. Here, the cohesive strength is modeled by Weibull A bilinear cohesive law has two important parameters, namely the cohesive strength T 0 n and fracture energy G IC . The fracture energy of the CG 304 SS is obtained from G IC = K 2 IC 1 − v 2 /E in terms of its fracture toughness K CG IC , which is taken as 75 MPa √ m [45]. For the NGL, K NGL IC is taken as 36 MPa √ m. Young's modulus E and Poisson's ratio v of both the CG and the NGL are taken as 200 GPa and 0.33, respectively [46]. The flow stress of the CG specimen is from the experiments [46], while that of the NGL is estimated from experiments of the SMATed 316L SS [47]. The kinematic hardening model is used to describe the cyclic plastic behavior of the metals. Additionally, the initial stiffness k should be large enough to alleviate the inherent mesh dependence, while too large value may cause numerical ill-conditioning [38,41]. It is determined as 2 × 10 19 N/m.

Weibull Random Field
Based on the weakest link model, the Weibull distribution assumes that strength is connected with flaws randomly distributed in the materials. Here, the cohesive strength is modeled by Weibull random fields, while the effective separation at failure u f n is taken as a constant. Therefore, the fracture energy is also a random field. The Weibull cumulative distribution function of failure of the material subjected to a stress σ can be formulated into where σ 0 the material parameter, m the Weibull modulus, σ u the threshold stress below which the failure will not occur, V the sample volume, and V 0 the reference one. Its mean and standard variance depend on σ 0 , σ u , and m. Here, we take the Weibull modulus m as 50 [48]. The threshold stress σ u is taken as the yield strength σ y , and the material parameter σ 0 as the ultimate strength σ f [48].

Cohesive Strength Effects on the Fatigue Life for the CG Specimens
As described in Section 2, to relate the mean cohesive strength T As depicted in Figure 6, when T 0 n = 8 σ y , the predicted fatigue life is lower than the experimental results; when T 0 n = 13 σ y , the predicted fatigue life is higher. Therefore, T 0 n is calibrated as 10 σ y , as listed in Table 2. random fields, while the effective separation at failure is taken as a constant. Therefore, the fracture energy is also a random field. The Weibull cumulative distribution function of failure of the material subjected to a stress can be formulated into where the material parameter, m the Weibull modulus, the threshold stress below which the failure will not occur, V the sample volume, and the reference one. Its mean and standard variance depend on , , and m. Here, we take the Weibull modulus m as 50 [48]. The threshold stress is taken as the yield strength , and the material parameter as the ultimate strength [48].

Cohesive Strength Effects on the Fatigue Life for the CG Specimens
As described in Section 2, to relate the mean cohesive strength with yield strength of the CG 304 SS, = 8 , 10 , and 13 , the corresponding , are selected to investigate its effect on the fatigue life under the asymmetric sinusoidal forces, namely, 893.6-8936 N (the corresponding stress amplitude Δσ is 320 MPa), 868.5-8685 N (Δσ = 311 MPa), and 843.3-8433 N (Δσ = 302 MPa). As depicted in Figure 6, when = 8 , the predicted fatigue life is lower than the experimental results; when = 13 , the predicted fatigue life is higher. Therefore, is calibrated as 10 , as listed in Table 2.      By fixing the axial load 8936 N, we investigate the end displacement and damage evolution of the CG specimen under five random fields. As an example, Figure 8 shows the end displacement of the specimen under No. 2 random field. In the early stage, since the material constitutive relation has a stabilization process in the experiment, the experimental value will be smoothly transitioned. After the stabilization, the mean and amplitude of the simulated end displacement are consistent with the experimental results. The damage in the CG specimens is evaluated by the SDEG. All elements in Figure 9 are severely damaged (i.e., SDEG > 0.3). It can be seen that the damage initiates from the exterior surface of the specimen. With the increase in the cycle number, the damage propagates circumferentially and ultimately fails the specimen. Figure 10 shows the damage evolution under the same axial load and the different random fields at the same cycle number (16,000). The severely damaged elements (i.e., SDEG > 0.5) are shown in Figure 10. It shows that the number and the distribution of severely damaged elements vary from each other. In this case, different random fields have various weakened zones, which leads to the different damage degree and distribution during the fatigue loading. The damaged elements are mainly distributed along the circumferential direction.  By fixing the axial load 8936 N, we investigate the end displacement and damage evolution of the CG specimen under five random fields. As an example, Figure 8 shows the end displacement of the specimen under No. 2 random field. In the early stage, since the material constitutive relation has a stabilization process in the experiment, the experimental value will be smoothly transitioned. After the stabilization, the mean and amplitude of the simulated end displacement are consistent with the experimental results. The damage in the CG specimens is evaluated by the SDEG. All elements in Figure 9 are severely damaged (i.e., SDEG > 0.3). It can be seen that the damage initiates from the exterior surface of the specimen. With the increase in the cycle number, the damage propagates circumferentially and ultimately fails the specimen. Figure 10 shows the damage evolution under the same axial load and the different random fields at the same cycle number (16,000). The severely damaged elements (i.e., SDEG > 0.5) are shown in Figure 10. It shows that the number and the distribution of severely damaged elements vary from each other. In this case, different random fields have various weakened zones, which leads to the different damage degree and distribution during the fatigue loading. The damaged elements are mainly distributed along the circumferential direction.
circumferentially and ultimately fails the specimen. Figure 10 shows the damage evolution under the same axial load and the different random fields at the same cycle number (16,000). The severely damaged elements (i.e., SDEG > 0.5) are shown in Figure 10. It shows that the number and the distribution of severely damaged elements vary from each other. In this case, different random fields have various weakened zones, which leads to the different damage degree and distribution during the fatigue loading. The damaged elements are mainly distributed along the circumferential direction.    Figure 11 shows the end displacement of specimens associated with No. 2 random field under the above three axial cyclic loads. It can be seen that both the amplitude and the mean of the end displacement under higher axial stress is larger. This is due to the fact that the higher axial stress generates larger deformation, which results in the larger amplitudes and higher mean of end displacement. On the other hand, when the stress is lower, the amplitude of the displacement is   Figure 11 shows the end displacement of specimens associated with No. 2 random field under the above three axial cyclic loads. It can be seen that both the amplitude and the mean of the end displacement under higher axial stress is larger. This is due to the fact that the higher axial stress generates larger deformation, which results in the larger amplitudes and higher mean of end displacement. On the other hand, when the stress is lower, the amplitude of the displacement is relatively smaller. The elements in Figure 12 are severely damaged (i.e., SDEG > 0.5). The number of  Figure 11 shows the end displacement of specimens associated with No. 2 random field under the above three axial cyclic loads. It can be seen that both the amplitude and the mean of the end displacement under higher axial stress is larger. This is due to the fact that the higher axial stress generates larger deformation, which results in the larger amplitudes and higher mean of end displacement. On the other hand, when the stress is lower, the amplitude of the displacement is relatively smaller. The elements in Figure 12 are severely damaged (i.e., SDEG > 0.5). The number of damaged elements increases with the stress and the damage is more serious under higher stress.

Load Effects on the Fatigue Life and Damage Evolution of CG Specimens
Here, although the current model can simulate the damage process reasonably, the end effect exists, i.e., some severe damage occurs in the end zones due to the limitation on the length of the computational specimens. In the future study, we will further improve the model.

Fatigue Life and Damage Evolution of the SMATed Specimens
The material properties of the nanocrystalline SS are assigned to the elements in the NGL. The other elements have the same material properties of the CG 304 SS, as listed in Table 2. To reduce the element size effects, the local MCS is adopted to all of the elements in the NGL and the outer two layers of CG elements, while the other elements have homogenized properties. In our experiments, the SMATed 304 SS exhibits an enhancement of the fatigue life. Therefore, three elevated force ranges are selected, namely 979.6-9796 N (Δσ = 351 MPa), 942-9420 N (Δσ = 342 MPa), and 893.6-8936 N (Δσ = 320 MPa). According to Weibull distribution, five random fields are generated under each load. Figure 13 shows five random distributions for the cohesive strength of the elements in the NGL with the thickness 20 µm. The mean and standard deviation of the cohesive strength are 5109 MPa and 88 MPa, respectively. As shown in Figure 14, the predicted fatigue life is in a good agreement with the

Fatigue Life and Damage Evolution of the SMATed Specimens
The material properties of the nanocrystalline SS are assigned to the elements in the NGL. The other elements have the same material properties of the CG 304 SS, as listed in Table 2. To reduce the element size effects, the local MCS is adopted to all of the elements in the NGL and the outer two layers of CG elements, while the other elements have homogenized properties. In our experiments, the SMATed 304 SS exhibits an enhancement of the fatigue life. Therefore, three elevated force ranges are selected, namely 979.6-9796 N (Δσ = 351 MPa), 942-9420 N (Δσ = 342 MPa), and 893.6-8936 N (Δσ = 320 MPa). According to Weibull distribution, five random fields are generated under each load. Figure 13 shows five random distributions for the cohesive strength of the elements in the NGL with the thickness 20 µm. The mean and standard deviation of the cohesive strength are 5109 MPa and 88 MPa, respectively. As shown in Figure 14, the predicted fatigue life is in a good agreement with the

Fatigue Life and Damage Evolution of the SMATed Specimens
The material properties of the nanocrystalline SS are assigned to the elements in the NGL. The other elements have the same material properties of the CG 304 SS, as listed in Table 2. To reduce the element size effects, the local MCS is adopted to all of the elements in the NGL and the outer two layers of CG elements, while the other elements have homogenized properties. In our experiments, the SMATed 304 SS exhibits an enhancement of the fatigue life. Therefore, three elevated force ranges are selected, namely 979.6-9796 N (∆σ = 351 MPa), 942-9420 N (∆σ = 342 MPa), and 893.6-8936 N (∆σ = 320 MPa). According to Weibull distribution, five random fields are generated under each load. Figure 13 shows five random distributions for the cohesive strength of the elements in the NGL with the thickness 20 µm. The mean and standard deviation of the cohesive strength are 5109 MPa and 88 MPa, respectively. As shown in Figure 14, the predicted fatigue life is in a good agreement with the experimental results.
The damage in the SMATed specimens is quantified by the SDEG. Figure 15 depicted the damage evolution of the SMATed specimen associated with No. 2 random field under 9796 N. The illustrated elements are severely damaged (i.e., SDEG > 0.3) and the damaged elements initially occur in the CG layer in Figure 15b. With an increase in the cycle number, the elements in the NGL also get damaged gradually, as shown in Figure 15c, and the specimen fails ultimately, which indicates that the damage initiates from the subsurface and extends to the exterior surface. SMAT improves the strength of the exterior surface so that the damage initiation site can transition from the exterior surface to the subsurface of the specimen, which is the intrinsic mechanism for the enhanced fatigue resistance of the SMATed specimens. A similar phenomenon has also been found in axial fatigue tests of other types of surface-treated specimens [49,50]. The number and the distribution of severely damaged elements vary from different random fields. In addition, the higher stress generates more severe damage and also leads to much more damaged elements. resistance of the SMATed specimens. A similar phenomenon has also been found in axial fatigue tests of other types of surface-treated specimens [49,50]. The number and the distribution of severely damaged elements vary from different random fields. In addition, the higher stress generates more severe damage and also leads to much more damaged elements.   resistance of the SMATed specimens. A similar phenomenon has also been found in axial fatigue tests of other types of surface-treated specimens [49,50]. The number and the distribution of severely damaged elements vary from different random fields. In addition, the higher stress generates more severe damage and also leads to much more damaged elements.

Conclusions and Outlooks
The heterogeneity of element sizes in a cylindrical solid may result in size effects in the fatigue process when the conventional global Monte Carlo method is employed. To alleviate the size effects, the local Monte Carlo method together with a 3D cohesive finite element method for solid metallic specimens is proposed to investigate the axial fatigue life and damage evolution. Our model has been verified by comparing the simulated fatigue life and end displacement with relevant experiments. After the surface mechanical attrition treatment, the damage of specimens initiates from the subsurface and extends to the exterior surface, leading to an increase in the fatigue resistance and the fatigue life. Therefore, this study shows that the proposed numerical approach can well predict the axial fatigue failure. Furthermore, it can also be employed to study the effects of other surface treatment on the fatigue life and thus provide guidelines for future experimental investigation.
Currently, fracture mechanics tools coupled with valid small crack growth equations can well compute the life of structural members with complex shapes under variable amplitude load spectra [2,[51][52][53][54][55][56][57]. At present, it is not a trivial task for the numerical model based on the cohesive finite element method and Monte Carlo method to accurately compute the fatigue crack growth history of realistic 3D structures under complex load histories. Therefore, further development of the cohesive finite element method and Monte Carlo method is of great interest.
For a given surface roughness, the probability distribution of lead crack growth rate and hence fatigue life can be determined by that of the local fatigue threshold, and the role of the fatigue threshold is the effect of surface roughness [58,59]. Linking the probability distribution of the local fatigue threshold in Ref. [58,59] with that of the local fracture properties in our numerical model and to establish the relationship between the parameters in the above two probability distributions is stimulating to provide more insights into the fatigue life of metallic materials or structures.
Author Contributions: X.G. conceived this study and made a plan for the simulations; Q.S., T.Y., and C.Z. prepared the numerical framework; Q.S. performed the simulations; X.G., Q.S., T.Y., G.J.W., and X.F. wrote the paper together.

Conclusions and Outlooks
The heterogeneity of element sizes in a cylindrical solid may result in size effects in the fatigue process when the conventional global Monte Carlo method is employed. To alleviate the size effects, the local Monte Carlo method together with a 3D cohesive finite element method for solid metallic specimens is proposed to investigate the axial fatigue life and damage evolution. Our model has been verified by comparing the simulated fatigue life and end displacement with relevant experiments. After the surface mechanical attrition treatment, the damage of specimens initiates from the subsurface and extends to the exterior surface, leading to an increase in the fatigue resistance and the fatigue life. Therefore, this study shows that the proposed numerical approach can well predict the axial fatigue failure. Furthermore, it can also be employed to study the effects of other surface treatment on the fatigue life and thus provide guidelines for future experimental investigation.
Currently, fracture mechanics tools coupled with valid small crack growth equations can well compute the life of structural members with complex shapes under variable amplitude load spectra [2,[51][52][53][54][55][56][57]. At present, it is not a trivial task for the numerical model based on the cohesive finite element method and Monte Carlo method to accurately compute the fatigue crack growth history of realistic 3D structures under complex load histories. Therefore, further development of the cohesive finite element method and Monte Carlo method is of great interest.
For a given surface roughness, the probability distribution of lead crack growth rate and hence fatigue life can be determined by that of the local fatigue threshold, and the role of the fatigue threshold is the effect of surface roughness [58,59]. Linking the probability distribution of the local fatigue threshold in Ref. [58,59] with that of the local fracture properties in our numerical model and to establish the relationship between the parameters in the above two probability distributions is stimulating to provide more insights into the fatigue life of metallic materials or structures.
Author Contributions: X.G. conceived this study and made a plan for the simulations; Q.S., T.Y., and C.Z. prepared the numerical framework; Q.S. performed the simulations; X.G., Q.S., T.Y., G.J.W., and X.F. wrote the paper together.