A Comparison of Two Methods Modeling High-Temperature Fatigue Crack Initiation in Ferrite–Pearlite Steel

Abstract

High-temperature fatigue tests are carried out to investigate the fatigue crack initiation behavior of ferrite–pearlite steel. Two approaches to modeling fatigue crack initiation under high temperatures in ferrite–pearlite steel are developed and compared in this study. In the first approach, basal energy is introduced as the criterion for crack initiation. Various initial basal energy values are assigned to each potential location. The increments in basal energy after each cycle can be calculated. In the second approach, a crystal plasticity model is developed, combined with the cyclic hardening/softening effect and three fatigue indicator parameters (two of them are adopted from previous works and one is developed in this work). Results indicate that both approaches are capable of predicting the fatigue crack initiation location and the corresponding life. The first approach requires a lower computational cost but has many limitations, and the second approach is suitable for calculating mesoscale stress and strain distributions but with a higher computational cost.

1. Introduction

The fatigue cracking of ferrite–pearlite steel has been widely investigated by different approaches since it is one of the most common materials in industry. However, the high-temperature fatigue crack initiation of ferrite–pearlite steel is still not fully predictable due to its complexity and inhomogeneity in microstructure [1]. For ferrite–pearlite steel, a fatigue crack always tends to initiate and propagate in ferrite grains, and crack propagation rates decrease rapidly at ferrite–pearlite grain boundaries [2]. In the pearlite phase, no cracking has been observed in either room-temperature or high-temperature fatigue tests [3].

Several different approaches have been proposed for predicting fatigue crack initiation at the mesoscale level. For instance, Hoshide et al. [4,5] proposed a methodology to simulate the fatigue crack initiation and final failure of Ti alloy with different microstructures. Compared with the traditional finite element method, the simulation process is faster, with less computational complexity, and is suitable for whole fatigue life simulation, i.e., from crack initiation to final collapse. However, this method is not able to calculate the stress and strain distribution inside each grain.

In the past decade, the crystal plasticity finite element method (CPFEM) has been a powerful method to investigate material deformation at the mesoscale level. The early-stage investigation of the crystal plastic deformation mechanism was completed by Taylor et al. [6,7,8] and the mathematical description was proposed by Hill et al. [9,10]. Huang [11] assumed that plastic deformation is caused by crystallographic dislocation slip only, which means that other types of plastic deformation are not considered. Based on this, Li et al. [12] proposed a new constitutive model considering twinning and detwinning mechanisms to model the uniaxial ratcheting of polycrystalline extruded magnesium (Mg) alloys. Yue et al. [13] conducted cyclic loading simulation using crystal plasticity theory developed by Asaro et al. [14]. The result indicated that a fatigue crack may initiate at the deformation-localized position, which means that CPFEM has the potential to be implemented for fatigue crack investigation.

The past decade has witnessed the rapid development of CPFEM. Sweeney et al. [15] carried out fatigue simulation based on the crystal plasticity constitutive model developed by Dunne et al. [16]; statistically stored dislocations, temperature, and other parameters were involved in this constitutive model. Experimental and simulation results were in good agreement, indicating that CPFEM with appropriate constitutive equations is suitable for predicting fatigue crack initiation precisely. Zhang et al. [17] investigated the fatigue life of a superalloy under symmetrical cyclic loading at 650 °C using CPFEM, assuming that fatigue failure takes place when the standard deviation of the strains reaches a critical value, which can be determined by the cyclic experiments with only one strain amplitude.

Crystal structure is an important aspect in CPFEM because resolved stress and slip deformation depends on it strongly. A number of CPFEM constitutive models have been developed with different crystal structures. The material under investigation in this work mainly consists of ferrite with a body-centered cubic (BCC) structure. Theoretically, there is only one set of slip systems in face-centered cubic (FCC) structures, i.e., $\left\{111\right\}〈110〉$, but three sets of slip systems in BCC structures, i.e., $\left\{110\right\}〈111〉$, $\left\{112\right\}〈111〉$ and $\left\{123\right\}〈111〉$, which makes it more complicated to calculate the slip deformation of BCC structures than FCC. Several simplified models have been proposed; for example, Ashton et al. [18] adopted $\left\{110\right\}〈111〉$ slip systems and the other two sets were neglected, since it has been shown by Franciosi et al. [19] that $\left\{110\right\}〈111〉$ were the principal slip systems in BCC iron. Sweeney et al. [15] adopted two types of slip systems, i.e., $\left\{110\right\}〈111〉$, $\left\{112\right\}〈111〉$. Most recently, Du et al. [20] performed uniaxial micro-tensile tests on single-crystal ferrite specimens with three different orientations; results indicated that only $\left\{110\right\}〈111〉$ and $\left\{112\right\}〈111〉$ slip systems were activated and the contributions made by these two systems for slip deformation were equal. The critical resolved shear stress (CRSS) values of the two active slip systems were close to each other, i.e., ${\mathrm{CRSS}}_{\left\{110\right\}〈111〉}=(1.0\pm 0.1)\times {\mathrm{CRSS}}_{\left\{112\right\}〈111〉}$.

The most intuitive way to simulate fatigue crack initiation should be modeling cracks directly, which can be achieved by different methods, such as extended finite element methods (XFEM) and cohesive element methods, etc. Recently, Zhang et al. [21] conducted CPFEM simulations using XFEM and a self-defined damage criterion to predict crack propagation in a single-crystal nickel-based superalloy; the crack path suited the experiment results very well. However, modeling multiple cracks’ initiation and propagation using the finite element method often causes severe convergence difficulties because of the contact and stress field singularity caused by crack tips; thus, an indicator to show whether an element is cracked was proposed to avoid the difficulties, which is called the fatigue indicator parameter (FIP). It has been demonstrated that FIPs are related to grain scale fatigue crack formation and microstructurally small fatigue crack growth by Mcdowell et al. [22] and can be used as surrogate measures of the driving force for fatigue crack initiation [23], which means that FIPs are capable of predicting the fatigue crack initiation position and life. A wide variety of FIPs have been proposed and implemented in CPFEM simulation; for instance, Dunne et al. [24] adopted the accumulated plastic strain as the FIP, and predicted persistent slip bands accurately. Hallberg et al. [25] conducted a series of simulations to investigate the influence of surface roughness and microstructure heterogeneities on fatigue crack initiation, and employed two additional modified FIPs to show the crack initiation position. However, the FIP needs to be defined based on the constitutive model, and each constitutive model should correspond to one or more appropriate FIPs. Based on previous works, two FIPs are adopted in this work to calculate crack initiation life and crack density, and a modified FIP is developed.

As discussed previously, CPFEM has been implemented for cyclic deformation simulation, and the cyclic hardening/softening effect has been introduced into CPFEM by several works. For example, Lin et al. [26] investigated a nickel-based superalloy with a modified crystal plasticity constitutive model and the stress relaxation behavior was calculated. Sajjad et al. [27] adopted the Chaboche model to simulate the cyclic hardening effect in martensitic high-nitrogen stainless steel. The simulation was carried out for dozens of cycles, and the ratcheting effect was captured well by this model. However, most cyclic hardening/softening simulations using CPFEM were conducted for relatively low cycles (less than 1000 cycles), and a CPFEM model to predict the cyclic hardening/softening effect under relatively high cycle loading needs to be developed.

In this paper, high-temperature fatigue experiments are conducted to investigate the fatigue crack initiation (Section 2). After this, two approaches to fatigue simulation are conducted. In the first approach of this work (Section 3.1), simulation is carried out by assigning suitable crack initiation energy (basal energy) for each grain, which is determined by experimental observation and will increase after each cycle. In the second approach (Section 3.2), a crystal plasticity constitutive model including the cyclic hardening/softening effect, kinematic hardening effect, and fatigue indicator parameters (FIPs) is established, and the simulations are carried out using ABAQUS with UMAT code. The geometric model used in both approaches is generated based on the microstructure of the specimen, and a Python script is written to generate the model. Both 2D and 3D models are adopted in the second approach. Considering computational efficiency, representative volume element (RVE) models are used. Here, 3D models are used to calculate the averaged stress–strain response and the evolution of peak–valley stress to demonstrate the validity of the constitutive model because of the increase in computing time [22]. Two FIPs are adopted from previous works, and another FIP is developed. Finally, a comparison is made to show the pros and cons of the two approaches.

2. Materials and Methods

2.1. Material

The material investigated in this work is low-carbon structure steel treated with annealing. The chemical composition and mechanical properties are listed in

Table 1. Chemical composition of ferrite–pearlite steel as percentage by weight (wt%).

C	Si	Mn	Cr	S	P
$0.2$	$0.23$	$0.51$	$0.013$	$0.011$	$0.024$

and

Table 2. Mechanical properties of ferrite–pearlite steel at different temperatures.

Name of Parameters	Units	Room Temperature	High Temperature
Temperature	°C	20	500
Yield strength	MPa	$335.5$	$316.1$
Tensile strength	MPa	$471.98$	$412.37$
Elongation	%	32	45
Elastic modulus	MPa	$2.28\times {10}^5$	$2.05\times {10}^5$
Poisson’s ratio		$0.303$	$0.301$
Fracture energy density	KJ/m${}^2$	2	2
Shear modulus	GPa	94	81

. Specimens are cut with an annular notch to make the fatigue crack more concentrated, as shown in Figure 1, which will facilitate fatigue crack observation.

2.2. Methods

Fatigue tests are conducted using the MTS Landmark system under strain control with a strain ratio of R = −1 by an axial extensometer with 25 mm gauge, and the temperature is set to 500 °C. Since the cross-sectional area is not invariant in the gauge section, the strain value read from the extensometer is not the actual strain at the notch position, shown in Figure 2. In this work, the average strain amplitude is set to 0.12% (controlled by the extensometer), and the actual strain at the notch position is 0.263%, which is measured by a strain gauge. All experimental parameters are shown in

Table 3. Experimental conditions.

Parameters	Units	Values
Temperature	°C	500
Strain rate	s${}^{-1}$	$0.0024$
Strain ratio		$-1$
Strain amplitude		$0.12\%$

.

The microstructure of the surface is obtained by Nital etching, shown in Figure 3. It can be seen that the pearlite distribution is significantly influenced by the extrusion direction. The “banded” distribution of pearlite is observed along the extrusion direction. The fatigue tests are carried out three times, and dozens of areas are observed. Based on the experimental results shown in Figure 4, cracks initiate at the ferrite–pearlite grain boundary (F-PGB) and ferrite–ferrite grain boundary (F-FGB) and a persistent slip band (PSB) can be observed; these are marked by the polygon, ellipse, and rectangle, respectively. Meanwhile, no cracks initiate in pearlite. The total fatigue life is 20,097, and almost all fatigue cracks initiate in 3000 cycles, which is nearly 15% of the whole life. After this, cracks begin to propagate and only a few new cracks initiate. Since the main purpose of this work is to compare the following two approaches predicting fatigue crack initiation, only the first 3000 cycles are investigated.

3. Results

3.1. The First Approach for Crack Initiation Prediction

The fundamentals of these two approaches are different. In the first approach, the program is written based on experimental observation without involving any types of deformation; in other words, no continuum mechanics are involved.

3.1.1. Generation of Geometric Model

The geometric model was generated based on Voronoi tessellation [28]. As shown in Figure 5, red polygons represent pearlite and the rest are ferrite. Based on Figure 3, the “banded” distribution of pearlite is noticed, and some ranges are defined to determine whether a polygon is pearlite. For example, as shown in Figure 5, if the center point coordinate of a polygon in the horizontal direction is between (0.02 mm, 0.05 mm, or 0.09 mm, 0.11 mm, etc.), it is considered to be pearlite. The arbitrary crystal orientations are assigned to each ferrite grain. The geometric model is restricted to 0.32 × 0.96 mm${}^2$ containing 1350 grains.

3.1.2. Crack Initiation Model

Based on the experimental observation shown in Figure 4, PSBsin ferrite grains, F-FGBs, and F-PGBs are considered to be crack initiation locations, while other locations, such as the pearlite phase, will not initiate cracks. Basal energy is assigned to each grain and grain boundary. Since no pearlite crack is observed in the experimental results, the basal energy of the pearlite phase is set to zero. The main equations introduced in this approach are adopted from previous work done by Wang et al. [1]. The increment in basal energy after each cycle is calculated by Equations (1) and (2).

$$\mathsf{\Delta }{E}_{PSB}=\left\{\begin{array}{cc}\frac{\pi \left(1-\mu \right)d_i{\left({\tau }_{mn}-{\tau }_c\right)}^2\lambda }{2G{W}_c}\hfill & \hfill \mathrm{if}{\tau }_{mn}\ge {\tau }_c\hfill \\ 0\hfill & \hfill \mathrm{if}{\tau }_{mn}<{\tau }_c\hfill \end{array}\right.$$

(1)

$$\mathsf{\Delta }{E}_{GB}=\frac{D_{ij}}{d_j}{k}_1\lambda$$

(2)

$\mathsf{\Delta }{E}_{PSB}$ is the increment in basal energy for PSB and $\mathsf{\Delta }{E}_{GB}$ is for the grain boundary. Here, G is the shear modulus, $\mu$ is the Poisson’s ratio, $W_c$ is the fracture energy density, ${\tau }_c$ is the critical resolved shear stress (CRSS), $\lambda$ is the stress relaxation coefficient, which equals 1 in this paper (following previous work [1]), $d_i$ is the slip band length of the $i_{th}$ grains, ${\tau }_{mn}$ is the shear stress on the associated slip plane, as shown in Figure 6, $D_{ij}$ is the $i_{th}$ edge length of the $j_{th}$ grain boundary, and $k_1$ is a parameter that can be calculated by Equation (3).

$$k_1=f\left(\sigma ,T\right)=C_1·\sqrt{\frac{1}{2}\left[{\left({\sigma }_1-{\sigma }_2\right)}^2+{\left({\sigma }_2-{\sigma }_3\right)}^2+{\left({\sigma }_3-{\sigma }_1\right)}^2\right]}·T$$

(3)

where T is temperature (500 °C in this work), and $C_1$ is the coefficient to scale the equivalent stress. ${\tau }_{mn}$ is defined by:

$${\tau }_{mn}={\gamma }_{m\theta }{\gamma }_{n\theta }{\sigma }_{\theta }+{\gamma }_{mz}{\gamma }_{nz}{\sigma }_z+{\tau }_{\theta z}\left({\gamma }_{n\theta }{\gamma }_{mz}+{\gamma }_{m\theta }{\gamma }_{nz}\right)$$

(4)

where ${\gamma }_{ij}$ are the components of the transformation tensor, which maps a stress tensor from the global system ($z,\theta ,r$) to the crystal system.

When the basal energy increases to 1, a crack is initiated. The Cauchy stress tensor is assigned to each grain based on Equation (4) and only one slip direction is considered for each ferrite grain, although, theoretically, there are three sets of slip systems in ferrite grain. The slip direction and stress assignment are elaborated in Figure 6. The magenta dotted lines represent the slip direction, which is the projection of the slip direction vector on the geometric model. In this model, the lines are also the PSBs. The stress components in Equation (4) are according to the coordinate system in Figure 6, z direction is the loading direction, and r is the normal direction of the sample surface (point out of sample surface). The shear stress in each grain is calculated by Equation (4). Therefore, the shear stress ${\tau }_{mn}$ can be calculated; details can be found in [1].

3.1.3. Computational Model Parameters’ Determination

The FEM is used only to calculate stress components prior to the fatigue simulation process in the first approach. Therefore, the stress field after fatigue tests is not able to be obtained by this method. There are several parameters that need to be determined and calibrated in this model. ${\sigma }_{\theta }$, ${\sigma }_Z$, and ${\tau }_{\theta Z}$ are calculated based on the finite element method, and the slip band length $d_i$ is assumed to be the length of the magenta dotted lines in each grain; the shear modulus G and Poisson’s ratio $\mu$ are listed in Table 2. Since the shear modulus is 81 GPa, ${\tau }_c$ is set to 81 MPa, which is $0.1\%$ of the shear modulus, and the fracture energy density $W_c$ is 2 KJ/m${}^2$, which is adopted from previous work done by Hoshide et al. [29]. Fracture energy density is calibrated during the simulation procedure. Based on Figure 4, cracks are more likely to initiate in PSBs, and the same phenomenon was observed in the previous work [1]. The crack initiation is influenced by a number of mechanisms (dislocation pile-up, second phase particles, etc.). The basal energy criterion is not able to reflect all the mechanics. Therefore, a range of initial values is defined, and the basal energies are random but within this range. The initial basal energy of PSBs in ferrite grains, F-FGBs, and F-PGBs is set to different values to fit the experimental results. All the parameters of the first approach are listed in

Table 4. Parameters in the first approach.

Parameters	Units	Values
G	GPa	81
$\mu$		0.3
${\tau }_c$	MPa	81
$W_c$	KJ/m${}^2$	2
$C_1$		$1.3\times {10}^{-10}$
Initial basal energy of PSB		0.3∼0.5
Initial basal energy of F-FGB		0.4∼0.6
Initial basal energy of F-PGB		0.45∼0.65

.

3.1.4. Results of the First Approach

The simulation results of the first 3000 cycles are shown in Figure 7, and the fatigue life is 19,500 cycles, which is close to the experimental result (20,097 cycles). In order to observe fatigue crack initiation more clearly, the specimen without Nital etching is used and the results are shown in Figure 8. Based on previous work done by Wang et al. [1], before 15% cycles of the whole life, the number of cracks increases rapidly, and it is called stage I. After stage I, cracks almost stop initiating. Since the main purpose of this paper is to predict fatigue crack initiation and the total fatigue life of the specimen is 20,097 cycles, only the first 3000 cycles are investigated. In this work, crack density is calculated using the cracked area divided by the total area. Based on Figure 7 and Figure 8, the crack density (at 3000th cycle) of the experiment and the first approach is 6.09% and 5.92%, respectively, and the evolution is shown in Figure 9. The crack density of the experiment is obtained by calculating the grey level of Figure 8. The lighter pixels are eliminated, and the darker pixels are considered to be cracks. An in-house MATLAB code is written to perform the calculation; in the meantime, the averaged crack width is also calculated. After this, the crack density of the experimental results (as shown in Figure 7) is equal to the total area of cracks (the small black rectangles) divided by the whole observation area, where the width of the rectangles is defined based on previous calculation.

3.2. The Second Approach for Crack Initiation Prediction

The second approach is based on the crystal plasticity finite element method (CPFEM). Unlike the first approach, CPFEM is developed from the crystal plasticity constitutive model; the stress and strain fields can be calculated, and details are discussed in this section.

3.2.1. Generation of Microstructures

A geometric model is critical for CPFEM simulation because the stress and strain distribution rely on it significantly. There are both two-dimensional (2D) and three-dimensional (3D) models adopted in many papers; for example, Berisha et al. [3] adopted a 2D model to calculate the stress distribution in pearlite grains. McDowell et al. [22] constructed both 2D and 3D models to investigate the evolution of the fatigue indicator parameter. In their work, 2D models are used not only to calculate the averaged stress–strain response and cyclic hardening/softening effect, but also the precise positions of crack initiation, i.e., the high FIP positions. Meanwhile, 3D models are used to calculate the averaged stress–strain response to support the design of the CPFEM model, because the increase in computing time when moving from 2D to 3D microstructures is often several orders of magnitude, and such computations are highly intensive, even in 2D [22]. Recently, Stopka et al. [30] exploited a statistical technique to predict maximum FIPs in different materials with 3D models, and the effect of the finite element and the size of the statistical volume element (SVE) are investigated.

The 2D geometric model is generated in the same way as in the first approach, by a Python script, and it is meshed using the open-source software NEPER [31]. To consider the calculation efficiency, three different representative volume element (RVE) models are generated, containing 115 grains, 247 grains, and 421 grains, respectively, as shown in Figure 10. The calculation would be faster with fewer grains. After this, a Python script is written to set orientations and other material parameters for each grain in the model. At this moment, the orientations are also arbitrary.

3.2.2. Constitutive Equations

Based on Huang [11] and Kang et al. [32], a crystal plasticity constitutive model considering the cyclic hardening/softening effect is developed.

The overall deformation gradient tensor can be subdivided into elastic and plastic parts:

$$F=F^e{F}^p$$

(5)

Assuming that plastic deformation is attributed to slip only, we have

$$L_p=\sum _{\alpha =1}^n\dot{{\gamma }^{\alpha }}{s}^{\alpha }\otimes m^{\alpha }$$

(6)

where $L_p$ is the plastic part of the strain rate tensor, $\alpha$ stands for the slip system, and n is the total slip systems activated in the current model. s is the slip direction vector and m is the normal vector of the slip plane. The power-law relation proposed by Hutchinson [33] is used in this work, modified by adding a back stress term, a slip resistance term [32], and a cyclic hardening/softening term:

$$\dot{{\gamma }^{\alpha }}={〈\frac{|{\tau }^{\alpha }-{\chi }^{\alpha }|-Q^{\alpha }+C^{\alpha }}{N}〉}^{N_d}sign({\tau }^{\alpha }-{\chi }^{\alpha })$$

(7)

where $C^{\alpha }$, ${\tau }^{\alpha }$, ${\chi }^{\alpha }$, and $Q^{\alpha }$ are the cyclic hardening/softening term, resolved shear stress, back stress, and isotropic deformation resistance of slip system $\alpha$, respectively, N and $N_d$ are material parameters. 〈〉 is the Macauley’s brackets:

$$〈X〉=\left\{\begin{array}{cc}X\hfill & \hfill \mathrm{if}X>0\hfill \\ 0\hfill & \hfill \mathrm{if}X\le 0\hfill \end{array}\right.$$

(8)

The $C^{\alpha }$ term is written as:

$$C^{\alpha }=S+K·\frac{e^{B·D·{\gamma }^{\alpha }}-1}{B}$$

(9)

where S is the initial cyclic hardening/softening parameter, and D is the cyclic parameter, which will be explained in the next section. K and B are parameters controlling cyclic hardening/softening, and ${\gamma }^{\alpha }$ is the cumulative slip strain of the $\alpha$_th slip system. The evolution of back stress $\chi$ and deformation resistance Q is characterized by

$$\dot{Q^{\alpha }}=\sum _{\beta =1}^n{H}_{\alpha \beta }|\dot{{\gamma }^{\beta }}|$$

(10)

$$\dot{{\chi }^{\alpha }}=A\dot{{\gamma }^{\alpha }}-A_d{\chi }^{\alpha }|\dot{{\gamma }^{\alpha }}|$$

(11)

$$H_{\alpha \beta }=qH+(1-q)H{\delta }_{\alpha \beta }$$

(12)

where $H_{\alpha \beta }$ is the hardening moduli, and q is a constant representing the ratio of latent over self-hardening moduli. A and $A_d$ are direct hardening and dynamic recovery coefficients, controlling the evolution of back stress. The hardening law proposed by Bassani and Wu [34] is adopted and simplified in this work:

$$H_{\alpha \alpha }=\left(h_0-h_s\right)sec{h}^2\left[\frac{(h_o-h_s){\gamma }^{\left(\alpha \right)}}{{\tau }_s-{\tau }_0}\right]+h_s$$

(13)

$$H_{\alpha \beta }=q·H_{\alpha \alpha }$$

(14)

where ${\tau }_s$, ${\tau }_0$, and $h_0$ are the stage I stress, yield stress, and the initial hardening modulus, respectively.

The microstructure of pearlite is more complicated than that of ferrite, which makes it difficult to predict plastic deformation using the finite element method. A direct modeling of pearlite was carried out by Berisha et al. [3] using a spectral solver developed by Roters et al. [35]. There are some different approaches based on homogenization techniques; for example, Ashton et al. [18] used the same constitutive model for both the pearlite and ferrite phase, with identical parameters except for the critical resolved shear stress; the simulation results showed good agreement with the experimental results. Based on the fatigue test results, no fatigue crack in the pearlite phase was observed, which made it possible to treat the pearlite phase as a homogeneous structure. Following Ashton et al. [18], the elastic constants of the ferrite phase are set to 70% of the pearlite phase, and other parameters are same as in the pearlite phase.

In this paper, two different FIPs are adopted based on previous works [25]. The first one is the accumulated equivalent plastic strain, which is defined as

$$FI{P}_{ϵ}=\int {\left(\frac{2}{3}{L}^P:L^P\right)}^{1/2}dt$$

(15)

The second FIP is defined as

$$FI{P}_E=E=\int \sum _{i=1}^n{\tau }_{\alpha }\dot{{\gamma }_{\alpha }}dt$$

(16)

where E stands for the overall stored energy in the current structure without considering mechanical dissipation, and ${\tau }_{\alpha }$ is the resolved shear stress for the ${\alpha }_{th}$ slip system, as mentioned before.

3.2.3. Parameter Calibration

There are over a dozen parameters that need to be calibrated in this approach. The fatigue tests conducted in this work are under 500 °C, the elastic modulus at this temperature is 0.899 times that at room temperature (listed in Table 2), the elastic constants at room temperature of ferrite ($C_{11}$, $C_{12}$, $C_{44}$) can be obtained from [3], and the constants at 500 °C can be calculated. There are 12 parameters left to be determined; each parameter has an impact on the simulation results, which makes it difficult to calibrate all these parameters together.

Several works have been done for calibrating CPFEM parameters. For example, Bandyopadhyay et al. [36] developed an algorithm by combining MATLAB and ABAQUS to automatically calibrate parameters. Du et al. [20] determined some parameters by analyzing the model first, and then made some assumptions for the remaining parameters, and determined all parameters by varying them slightly to fit the experimental results. Since the CPFEM model developed in this work is applied to multiple cycles, parameter calibration needs to be done for all simulation cycles; thus, it is more convenient to make assumptions for some parameters at the beginning and then calibrate them based on the results. The material constants N and $N_d$ in Equation (7) control the slip rate predominantly; therefore, these two parameters should be settled first. Following previous work [37], an initial guess of N = 20 MPa and $N_d=50$ is made. The rest of the parameters are also set to initial values, as discussed in the following paragraph. Since N and $N_d$ only affect the flow rule of the material, after running several simulations, a stress–strain curve similar to the experimental result is obtained, and the N and $N_d$ are adjusted to 22.8 MPa and $18.8$, respectively. In the following calibration process, these two parameters will not change.

As mentioned in the Introduction, there are three sets of slip systems in ferrite, but a recent investigation showed that only $\left\{110\right\}〈111〉$ and $\left\{112\right\}〈111〉$ were activated during deformation, and the mechanical properties of these two systems were similar. Therefore, the activated slip systems are set to $\left\{110\right\}〈111〉$ and $\left\{112\right\}〈111〉$ with identical mechanical properties. In Equations (9) and (11), the back stress term and the cyclic term are affected by the slip rate and cumulative slip strain. To determine the parameters in both terms, a simulation without a cyclic term and back stress term (i.e., $C^{\alpha }=0$ and ${\chi }^{\alpha }=0$ in Equation (7)) is conducted and the average cumulative strain of 24 slip systems is measured in each increment during the whole process. Results indicate that the average cumulative strain increases $1\times {10}^{-4}$ after each cycle. In this work, only the first 3000 cycles of the experiment are investigated based on the analysis of the first approach. Considering computational efficiency, it is impractical to calculate 3000 real cycles using CPFEM; hence, the cyclic parameter D is introduced, which can scale the cumulative strain in each increment, making it possible that each cycle of CPFEM simulation represents 200 cycles of actual tests. The cyclic parameter D in Equation (9) should be set to 200 to scale the cumulative strain; however, the value of 200 times ${\gamma }^{\alpha }$ after each cycle is still too small. For the convenience of the parameter calibration procedure, D is defined as the cyclic parameter times 100; in this case, D is set to $2\times {10}^4$. The initial cyclic hardening/softening parameter S in Equation (9) is defined randomly but at a relatively small value to avoid a zero or negative value of $C^{\alpha }$; K is the factor controlling the value of $C^{\alpha }$—in other words, $C^{\alpha }$ increases as K increases. B controls the concavity and convexity of $C^{\alpha }$, which is determined to be a negative value based on the peak and valley stress evolution of the specimen shown in Figure 11. The precise values of K and B are calibrated by fitting the simulation results to the experimental results.

$h_0$ and $h_s$ in Equation (13) control the hardening effect; previous work [37] has shown that a simplified hardening law is effective to simulate the hardening effect. In this work, $h_0$ is set equal to $h_s$; therefore, only one parameter needs to be calibrated in Equation (13), which can be done by fitting the tensile stress–strain curve of the simulation results to the experimental results. An assumption is made that the interaction parameter $q=1.4$ [3,20].

The parameters in the back stress term, i.e., A and $A_d$ in Equation (11), cannot be determined directly based on the slip rate. Since the back stress is not the main purpose of this work, the initial values of A and $A_d$ are set to 100 and 10, respectively. The parameter Q in Equation (7) represents the initial slip resistance, affecting the slip rate jointly with other parameters in Equation (7). At this time, all parameters except Q, A, and $A_d$ are already determined; a few simulations are carried out to calibrate the three parameters left. All calibrated parameters in this work are listed in

Table 5. Calibrated parameters in the second approach; $C_{11}^f$, $C_{12}^f$, $C_{44}^f$ are elastic constants of ferrite and $C_{11}^p$, $C_{12}^p$, $C_{44}^p$ are elastic constants of pearlite.

Parameters	Values	Units
$C_{11}^f$, $C_{12}^f$, $C_{44}^f$	$210.36,122.07,106.31$	GPa
$C_{11}^p$, $C_{12}^p$, $C_{44}^p$	$300.51,174.39,151.87$	GPa
A	$102.8$	MPa
$A_d$	$7.8$
q	$1.4$
Q	$108.55$	MPa
S	$0.35$	MPa
K	$0.032$	MPa
B	$-0.001$
D	20,000
$h_0$	$4.0$	MPa
$h_s$	$4.0$	MPa
N	$22.8$(2D) 110(3D)	MPa
$N_d$	$18.8$

.

3.2.4. Simulation and Results

Firstly, a series of simulations are carried out to investigate the influence of RVE size, and three models are chosen, shown in Figure 10. Periodic boundary conditions are assigned to each side of the RVE models and cyclic strain with 0.263% amplitude is applied, corresponding to the experimental conditions; the results are shown in Figure 12. In this approach, stress–strain curves are calculated by multiplying the stress and strain at each integration point by the volume of the point and summing them, and then dividing this by the total volume. The integration volume is obtained from the ABAQUS field output file, and a Python script is written to perform the calculation. In a 2D model, “volume” means “area”. Based on the results, the overall stress–strain curves of the three RVE models are slightly different from each other. The stress–strain results of the models with 115 grains and 421 grains are very close, but a slight deviation is observed in the 247-grain model. As discussed in the previous section, the orientation of each grain is random. In order to investigate the influence of model size, the same random seed is assigned to the three models using NumPy [38]. Therefore, the orientation of the same parts in the three models is identical, but the rest of the parts are different. In other words, all three models contain part of (or whole) 115 grains, and the authors can only guarantee that the orientations of 115 grains in the three models are identical. As a result, the stress–strain responses of the three models are different, but the difference is acceptable. Under the consideration of computational efficiency, the RVE model containing 115 grains is adopted. After the size of the RVE model is determined, another series of simulations are carried out to investigate the influence of mesh size. The RVE model is meshed by four different Relative Characteristic Lengths (defined in [31]) containing 3600 elements, 8100 elements, 14,161 elements, and 22,201 elements, respectively. The simulation results are shown in Figure 13, and the model containing 3600 elements is satisfactory. Finally, the simulation is conducted three times to evaluate the influence of grain orientation; three random orientation distributions are assigned to the model. The different orientations are achieved by assigning different random seeds to the three models (numpy.random.seed [38]).The results are shown in Figure 14, which indicates that the influence of orientation for the stress–strain curve is not significant.

In this paper, CPFEM simulation is conducted for 16 cycles, each representing 200 cycles in the fatigue experiment.

Table 6. Correspondence of cycles between experiment and simulation.

Experiment	1	200	400	600	800	1000	1200	1400	1600	1800	2000	2200	2400	2600	2800	3000
Simulation	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16

lists the correspondence of cycles between experiment and simulation.

The hysteresis of the simulation and experiment is shown in Figure 15; the simulation results of the first, 6th, 11th, and 16th cycle and the corresponding experimental results are shown here, representing the first, 1000th, 2000th, and 3000th cycle. The evolution of peak and valley stress is shown in Figure 11. In Figure 15, the back stress and flow stress calculated by CPFEM fit the experimental results very well, indicating that the back stress term (Equation (11)) and the simplified hardening term (Equation (13)) are valid. The evolution of the peak and valley stress values is also in good agreement with the experiment, which means that the constitutive model developed here is suitable for predicting the cyclic hardening/softening effect. It is noteworthy that the cyclic hardening/softening effect is introduced by reducing the strength of the crystal, explained in Equations (7) and (9). To demonstrate the validity of the two equations, another simulation is carried out with D changed from 20,000 to 30,000, while other parameters remain unchanged. The results indicate that each cycle in CPFEM is able to represent 300 cycles in the experiment, as shown in Figure 11 and Figure 16.

Crack density and crack initiation life are determined by FIPs. The distributions of the two FIPs ($FI{P}_{ϵ}$ and $FI{P}_E$) are similar, as shown in Figure 17 and Figure 18. Areas with relatively high FIP values are more likely to be the crack initiation area. Both FIPs in the pearlite phase are set to zero because no crack is observed in it. Different crack initiation positions are observed, which are F-PGBs (marked with polygons), F-FGBs (marked with circles), and PSBs (marked with rectangles). Attempts have been made to determine the crack initiation life by setting critical values for the two FIPs, i.e., once the FIP value of an integration point reaches the critical value, it is cracked. To demonstrate the validity of the FIPs, critical values are set to $0.26$ for $FI{P}_E$ and $48.0$ for $FI{P}_{ϵ}$ to ensure that the crack density is nearly 6% at 3000 cycles, as shown in Figure 9. The evolution of crack density calculated from both FIPs is not in accordance with the experimental result, which is an expected phenomenon, because both FIPs are calculated based on the slip rate (defined in Equations (6), (15) and (16)), which is obtained from Equation (7). In Equation (7), the slip rate is positive when the value in the brackets is greater than zero, N and $N_d$ are invariant, $C^{\alpha }$ increases to the same value for all grains, and ${\tau }^{\alpha }$ is calculated based on the Schmid factor, which changes slightly during simulation. Therefore, as long as the slip rate is positive at an integration point, the values of $FI{P}_{ϵ}$ and $FI{P}_E$ will continue to increase, even if its stress is relatively low (less than critical resolved shear stress). This is incorrect because stress is also an important aspect of crack initiation, which means that it is difficult to determine the fatigue crack initiation life based on the previous two FIPs in this approach.

In order to predict the crack initiation life, another FIP is proposed by modifying $FI{P}_{ϵ}$, explained as follows:

$$FI{P}_{eff}=\int {\left(\frac{2}{3}{L}_{eff}^P:L_{eff}^P\right)}^{1/2}dt$$

(17)

$$L_{eff}^P=\sum _{\alpha =1}^n\left(\dot{{\gamma }^{\alpha }}{s}^{\alpha }\otimes m^{\alpha }{\tau }_{eff}^{\alpha }\right)$$

(18)

$${\tau }_{eff}^{\alpha }=\left\{\begin{array}{cc}{\left(\frac{{\tau }^{\alpha }}{{\tau }_c^{\alpha }}\right)}^i\hfill & \mathrm{if}{\tau }^{\alpha }>{\tau }_c^{\alpha }\hfill \\ 0\hfill & \mathrm{if}{\tau }^{\alpha }\le 0\hfill \end{array}\right.$$

(19)

where $\tau$ is the resolved shear stress and ${\tau }_c$ is the critical resolved shear stress; i is a constant controlling the influence of the resolved shear stress on $FI{P}_{eff}$.

In this FIP, the accumulated equivalent plastic strain is modified. An assumption is made that the resolved shear stress and slip deformation affect fatigue crack initiation together. In this CPFEM model, according to Equation (7), slip deformation occurs under not only high resolved shear stresses (greater than critical resolved shear stress) but also relatively low stresses (less than critical resolved shear stress); therefore, the slip deformation under low stresses should not be considered in the FIP. To this end, the effective resolved shear stress ${\tau }_{eff}$ is invoked in Equation (19), and the effective plastic strain rate tensor is written as Equation (18), which indicates that the contribution of the slip rate to fatigue cracking is different due to the resolved shear stress.

In the first approach, the critical resolved shear stress is set to 0.1% of the shear modulus (in Equation (1) ${\tau }_c$ = 81 MPa and G = 81 GPa), which is a conservative assumption, because only one slip system is adopted in each grain, and the critical resolved shear stress must be set to a relatively low value in case some activated slip system is neglected (according to Equation (1)). In this approach, the critical resolved shear stress is set to $121.5$ MPa, which is 1.5 times that in the first approach. As shown in Equation (7), the slip rate is calculated by some parameters. Theoretically, crystals will not slip if the shear stress does not reach the critical value. However, Equation (7) cannot map this situation accurately, which means that some slip deformation under low shear stress exists, and a number of CPFEM models have the same drawback [11]. This phenomenon may not influence the stress–strain response (as shown in this work, the results are in good agreement with experiments), but it affects the FIP values because the FIP strongly depends on slip deformation. Therefore, in order to eliminate the contribution of the “invalid slip deformation”, a larger estimation of the critical resolved shear stress is made, which is 1.5 times that in the first approach, and it gives a good result for crack density prediction. The distribution and evolution of $FI{P}_{eff}$ are shown in Figure 19 (pearlite phase is also neglected). As with the previous two FIPs, $FI{P}_{eff}$ is also suitable for predicting crack initiation positions. To calculate the fatigue crack initiation life, the critical value of $FI{P}_{eff}$ is set to $0.011$ to ensure that the crack density equals 5.9% at the end of the simulation. In order to investigate the influence of orientation on crack density, another orientation distribution is assigned to the model; the critical value of $FI{P}_{eff}$ is also determined following the same step. The crack density evolution is shown in Figure 9 and the critical values are $0.011$ and $0.0115$ for the two orientation distributions, respectively.

A phenomenon is observed in the magnification box of Figure 9 wherein parts of the FIP curves are flat; the explanation is as follows. According to Equation (19), if the stress is relatively low, the increment in the FIP is zero. As shown in Figure 20, the corresponding stress and strain of the flat parts are relatively low, which are near $-200$ MPa to 150 MPa of stress and $-0.0015$ to $0.001$ of strain; therefore, the increment in the FIP at these parts should be zero. The threshold values of stress and strain are not symmetric, which may be caused by hardening and back stress effects.

3.2.5. 3D Geometric Models

Here, 3D geometric models are also generated and meshed by open-source software NEPER [31]. Three models are generated with 200 grains and 300 grains (shown in Figure 21), containing 3375 elements and 4913 elements, respectively, and the element type is C3D8 (3D 8-node brick element). As discussed in previous work [22], the computations are rather intensive in 2D CPFEM simulations, and when moving to 3D models, the computing time is different by several orders of magnitude. For instance, the smallest 2D model in this work contains 115 grains and 3600 CPS4 elements; a simple calculation shows that if a cubic model represents the same area, it should contain approximately 1200 grains and 216000 C3D8 elements, and the computing time will increase hugely. Therefore, we use 2D models as a simplification procedure to calculate the distribution of FIPs and other details. Here, 3D models with relatively coarse meshes and fewer grains are adopted to calculate the averaged stress–strain response to demonstrate the validity of the CPFEM model. In this work, the computational time of Model (A) is 14 times that of the 115-grain 2D model.

In Equation (7), the shear stress affects the slip rate significantly. With the element type changing from CPS4 to C3D8, the stress will also change. Therefore, the parameters in Equation (7) need to be recalibrated. By fitting the simulation to the experimental result, N is changed from $22.8$ to 110, and all other parameters remain unchanged. The increase in N is expected because the 2D models are simplified structures using plain stress elements to represent three-dimensional structures, and stress components ${\sigma }_{xz}$, ${\sigma }_{yz}$, and ${\sigma }_{zz}$ are ignored. Since the slip rate in Equation (7) is calculated by the shear stress (determined by Equation (20)), the slip rate calculated by the 2D models is smaller than that obtained by the 3D models. Thus, the associated parameters need to be recalibrated.

$${\tau }^{\alpha }=m^{\alpha }·\sigma ·s^{\alpha }$$

(20)

As with 2D models, periodic boundary conditions are assigned and cyclic strain with 0.263% amplitude is applied, corresponding to the experimental conditions. Two models with different grains (Model (A) and Model (C) in Figure 21) are calculated first to investigate the influence of RVE size, and the results are shown in Figure 11. Comparing the experiment and 2D models, the evolution of the stress peak and valley is in good agreement. After this, four models containing 200 grains but with various orientations are calculated to demonstrate the influence of orientations. The results are similar, and two of them are shown in Figure 22. The difference between the two orientations is acceptable, which means that the influence of orientation on the models with 200 grains is negligible. Finally, another model containing 200 grains but with different grain structures (Model (B) in Figure 21) is calculated, as shown in Figure 11. The structures of Model (A) and Model (B) are different, but both curves fit the experiment very well, indicating that the 3D RVE model with 200 grains is sufficient to represent the material.

4. Discussion

The main objective of this paper was to compare the two different approaches to fatigue crack initiation simulation on ferrite–pearlite steel. According to the experimental results, almost all cracks initiate before 3000 cycles, which is approximately 15% of the whole life and is called stage I; the rest of the fatigue life is not discussed in this work.

In the first approach, a 2D geometric model is established based on the microstructure of the specimen. The model is developed based on the works of Hoshide et al. [4] and Wang et al. [1], with the definition of F-FGBs, F-PGBs, and PSBs as potential positions for fatigue crack initiation. The basal energy is assigned to each position with different values, and it increases to a certain value after each cycle, which is calculated by Equations (1) and (2). When the basal energy increases to one, a crack is initiated. This approach assumes only one slip system in each ferrite grain, and stress and strain are not calculated. The crack initiation life and density calculated by this approach fit the experiment results very well. The simulation process is much faster than the finite element method because the equations solved in this model (Equations (1) and (2)) require a lower computational cost than equations in the finite element method. However, only fatigue crack behavior can be obtained from this approach; the stress and strain distributions are not available.

In the second approach, 2D geometric models similar to the first approach are firstly established, and a crystal plasticity constitutive model is developed to simulate the fatigue behavior of the first 3000 cycles. A new cyclic hardening/softening term is introduced in the model, which is able to shorten the simulation procedure. By adjusting parameter D, each cycle in CPFEM simulation is able to represent hundreds of cycles in the experiment. Crack density is calculated by three FIPs; the first two FIPs are adopted from previous works, and the third is developed in this paper. The results indicate that all three FIPs are capable of predicting the fatigue crack initiation positions, but the crack density calculated by $FI{P}_{ϵ}$ and $FI{P}_E$ is not in agreement with the experimental results; the reason is given in Section 3.2.4. The third FIP introduces the effective resolved shear stress (defined in Equation (19)) and neglects the slip rate under low resolved shear stress. Based on the results of crack density evolution, the third FIP is confirmed as valid in this work. Secondly, 3D models are generated containing 200 and 300 grains with coarse meshes to calculate the overall stress–strain response and cyclic hardening/softening effect using the constitutive model developed in this work. The results indicate that the constitutive model is suitable for both 2D and 3D geometric models.

As shown in Figure 9, the crack density evolution can be well predicted by the first approach and $FI{P}_{eff}$ in the second approach, and the crack initiation locations are also well predicted. Considering orientation’s influence, the critical value difference between the two orientation distributions is less than 5%; therefore, the influence of orientation on crack density is not significant. Comparing all crack density curves in Figure 9, $FI{P}_E$ and $FI{P}_{ϵ}$ are not able to predict the evolution of crack density, and the reason is discussed in the previous section. However, some drawbacks still exist in the current work. For example, even though $FI{P}_{eff}$ are able to reflect the crack density evolution, the critical value of $FI{P}_{eff}$ is determined according to the experimental results. For other materials, the critical value needs to be recalibrated by experiments, and the parameters in the CPFEM models and $FI{P}_{eff}$ also need to be recalibrated. However, once the parameters are determined for the material, this model is able to predict the crack density evolution.

Several conclusions are drawn:

• Both approaches are capable of predicting the fatigue crack initiation positions and initiation life in ferrite–pearlite steel; the crack density evolution is in good agreement with experiment results.
• The first approach determines whether a crack initiates by basal energy, which makes the simulation process much faster than the finite element method. However, stress and strain fields cannot be obtained from this approach.
• The second approach is based on CPFEM and more variables can be calculated, such as the stress and strain field, slip rate, resolved shear stress, etc., but a higher computational cost is required.
• Although the parameters are calibrated in 2D models, they can be applied to 3D models by recalibrating only one parameter.
• Both approaches are promising. For example, multiple slip systems can be introduced into the first approach, and more slip controlling parameters can be included in the second approach, such as temperature.