Gaps, challenges and possible solution for prediction of wheel–rail rolling contact fatigue crack initiation

The prediction of wheel/rail rolling contact fatigue (RCF) crack initiation during railway operations is an important task. Since RCF crack evolution is influenced by many factors, its prediction process is complex. This paper reviews the existing approaches to predict RCF crack initiation. The crack initiation region is predicted by the shakedown map. By combining the shakedown map with various initiation criteria and the critical plane method, the crack initiation life is calculated. The classification, methodologies, theories and applications of these approaches are included in this paper. The advantages and limitations of these methods are analyzed to provide recommendation for RCF crack initiation prediction. This review highlights that wheel/rail dynamic characteristic, complex working conditions, surface defects and wear all affect the RCF crack initiation. The optimal selection of criteria is essential in the crack initiation prediction. Based on the research gap regarding the challenging process of crack initiation prediction detailed in this review, a proposed prediction process of RCF crack initiation is proposed to achieve a more accurate result.


Introduction
Rolling contact between wheel and rail plays a vital role in many related wheel/rail damage modes, mainly including wear and rolling contact fatigue (RCF) [1][2][3]. With the growing demand of increased axle load and high speed, RCF, which is caused by repeated contact stress during the rolling motion, has become a severe and intractable problem in railway systems, leading to the reduction of wheel/rail service lifetime and the increase in maintenance budgets [4][5][6][7][8]. Furthermore, RCF is the incentive of surface and subsurface crack initiation, further propagation of the crack may result in fatigue phenomena such as spalling [9], squats [10] and headchecks [11]. For example, squats, which commonly occurred on the railhead as a typical type of RCF, could possibly lead to rail fracture failure. Cracks related to squats are long cracks with branches that some of might propagate upward to form spalling, and some of them might propagate downward into the depth of the rail [12]. Li et al. [13] mentioned that local plastic deformation was the basic characteristic of the squats, and the direct cause of squats was the excessive dynamic wheel-rail contact force. If the squats were left untreated, it might extend into the rail depth and break the rail, even leading to derailment [14]. Therefore, a better understanding of RCF crack initiation will help to prevent the risk of wheel/rail damage. There exist three main steps in the development of a surface initiated RCF crack (Fig. 1). During the evolution of the RCF process, a crack firstly initiates from one specific point at the contact surface as the consequence of plastic deformation accumulation (phase I). Then, a steady and continuous propagation will happen on the crack if there is no removal behavior utilized on the crack (phase II). Finally, branches occur in the tip of the crack aimed at causing significant fatigue damage, whereby some branches may propagate perpendicularly downwards into deeper areas, while others may go horizontally or even upward to the surface (phase III) [15,16]. Hence, the prevention of the initiation of RCF cracks is a key point to reduce the fatigue damage and prolong the service lifetime of the rail. Actually, the definition of crack initiation is indefinite. Some experiments or engineering projects are undertaken from the macroscopic aspect that defines the crack initiation as fracture or material separation of the specimen. However, most of the investigations mainly use a pre-determined initiated crack length, or a value of the reduction of the load in a strain-controlled test to determine the crack initiation [17]. It should be noted that many factors may contribute to RCF crack initiation.
According to previous studies on the initiation of RCF cracks, the influence factors can mainly be divided into three categories. The effect of working conditions such as normal contact pressure [18][19][20][21][22], wheel diameter [21], friction coefficient [18,19] and curvature radius [23,24] on RCF crack initiation is the emphasized part of the investigations. Among them, normal contact pressure is mostly being studied, not only for the magnitude of the load, but also for the form of contact between wheel and rail. For example, there are significant differences of contact condition between wheels running on a new rail and on a worn rail [25]. Meanwhile, internal properties are also important, such as material properties [20], material fatigue characteristics [26], material hardening [23] and residual stress [27]. All of these material parameters are used in the analysis of RCF crack initiation. McDowell et al. [26] reviewed the development of microstructure-sensitive computational methods for fatigue crack initiation in polycrystalline microstructures over the past few decades. And they pointed out that as a longer-term prospect in view of uncertainties in modeling mechanisms of cyclic slip, crack nucleation and growth, modeling can serve to support more quantitative predictions of fatigue lifetime as a function of microstructure. Furthermore, surface defects that occur during the rolling contact process may affect the initiation of RCF cracks. Defects such as squats [28,29] are being investigated regarding their influence on RCF crack initiation, also the size and location of defects [18] may affect the lifetime of wheels and rails. Based on the analysis of the influence of factors mentioned above, some prevention methods, appropriate combinations of input parameters and optimal material selection are proposed. For example, Ekberg et al. [30] put forward a predictive model of RCF crack initiation which established the relationship between the amount of surface fatigue and material fatigue resistance.
This paper provides an integrated overview of the prediction methods of RCF crack initiation, mainly focused on determining the location (the exact point on wheel and rail materials where it tended to initiate cracks) and the service period leading to RCF crack initiation. Firstly, the concept of RCF crack initiation prediction methods and the application of these methods on numerical simulation will be explained. Then, systematic approaches and an estimation methodology are described thoroughly, including the function of the shakedown map and critical plane methods for predicting RCF crack initiation. Development on shakedown map and review of fatigue calculation model are both mentioned. Furthermore, finite element analysis and dynamic simulation about wheel-rail modeling are the major subjects for the prediction models, and the comparison between their results with related experimental results will be illustrated. Finally, the research gap and a developed method for RCF crack initiation prediction are proposed.

Knowledge system in crack initiation prediction
RCF of railway wheel-rail materials are mainly affected by axle loads, wheel-rail profiles and friction coefficients, which can determine contact stresses, material responses and fatigue crack initiation in a great extent. A high axle load or friction coefficient relates to crack initiation due to low-cycle fatigue (LCF), and cracks initiated in railway wheel-rails in field observations and laboratory investigations all confirm that LCF is the domain fatigue life regime [19]. If the accumulation of plastic strain happens in a predominant direction and exceeds the critical level of material hardening and residual stresses, cracks will ultimately form at the higher fracture strain than a tensile test. This mode of fracture is called ratcheting [30]. In general, initiation of surface cracks is basically considered to be caused by LCF and ratcheting [31]. Also, LCF and ratcheting can be regarded as independent, competitive mechanisms of failure which need a deeper investigation. They also play a vital role in the prediction of RCF crack initiation. For example, if the predicted cycling number to failure by LCF is smaller than that by ratchetting, then the failure would be caused by LCF [32]. Significant challenges are faced in predicting RCF crack initiation, for the crack initiates under the complicated action of different factors (both internal and external). The critical task for the prediction is to determine a recognized method. The methods and analysis process of RCF crack initiation being reviewed in this paper mainly includes two parts: the location where crack initiate and the crack initiation life (Fig. 2). The shakedown map, which is mostly utilized in predicting the region of crack initiation as described in more detail in Sect. 3.2, is based on four material responses, namely elastic, elastic shakedown, plastic shakedown and ratcheting. Points in the plastic shakedown region and the ratcheting region, tend to initiate crack [33]. The predicted results showing plastic shakedown material response were used to evaluate the initiation criteria for LCF [19]. Note that the shakedown map is a template map directly used in the research, which originated from previous numerical calculations and empirical research activity. Thus, the accuracy of shakedown maps deserves further consideration.
However, as the shakedown map can only predict crack initiation in a qualitative way, the calculation equation of crack initiation life could be a quantitative way to predict crack initiation. The determination of the fatigue life equation mainly consists of three parts: the constitutive model, the critical plane (the plane where the FP max located) and the crack initiation criteria. It is noted that the criteria differ in the plastic deformation and ratcheting regions. Also, RCF life to crack initiation calculation criteria used in the research are various because researchers are focusing on different influencing factors (especially various kinds of stresses or strains). Although the criterion proposed by Jiang and Sehitoglu [34] is mostly used in recent research works, it is constantly being used in various updated ways to meet particular investigation requirements.
Given the complexity and diversity of crack initiation prediction, to propose a process of optimizing the prediction method according to the actual situation is important. And the authenticity and effectiveness of the prediction results need to be guaranteed. Each item in Fig. 2 is analytically reviewed in the following sections.

Prediction of the location where cracks initiate
The shakedown map comprehensively uses wheel/rail contact pressure and the tangential force to evaluate the bearing capacity of wheels and rails. It is widely used to qualitatively predict wheel/rail rolling contact fatigue, especially the crack initiation evaluation. Material response, which refers to the basic characteristic in establishing the shakedown map, is first reviewed in this section. Then, the concept, the input parameters and the application of the shakedown map is reviewed. Finally, due to the diversity of contact materials, contact conditions and effects of the external environment, it is necessary to update the shakedown map that has been established by the previous

Shakedown map
The shakedown map is commonly used by researchers to investigate the material response and the crack initiation. This part reviews the concept and the application of the shakedown map.

Concept of shakedown map
It is essential to establish a relationship between the friction coefficient and the maximum contact pressure among the four types of material response regions shown in Fig. 3 [24]. For this purpose, diagrams called "shakedown maps" that provide the material response showing the location for the greatest fatigue damage were first proposed by Johnson based on 2D numerical calculation ( Fig. 4) [35,36]. This shakedown map is aimed at the line contact situation, and the shakedown limit by the Tresca criterion is reached when P 0 /k = 1/μ. Then Ponter et al. [37] proposed the shakedown map for the more realistic case of point contact as shown in Fig. 5. This map is more representative of railway practice and is widely used in later research works. Surface or subsurface material damage is illustrated by the critical location of fatigue damage. The ordinate value, treated as a load factor, is the ratio of the maximum normal contact pressure p 0 to the ductility of material represented by the shear yield strength k. Meanwhile, the friction forces in the contact are characterized by traction coefficient μ on the abscissa of the map [15]. It should be mentioned that, although the shear yield strength k is commonly obtained using the Tresca criterion k = y ∕ √ 3 ( y is the material tensile yield strength) [38,39], some researchers pointed out that k can be estimated from the hardness measurement of the material. Jones et al. [40] considered the variation of the shear yield strength with depth after the material had been strain hardened, and put forward a relationship between k and the Vickers hardness HV: k = HV∕3 √ 3 . The crack tends to initiate when the investigated point falls into the plastic shakedown and ratcheting region. Because the shakedown map can be used to obtain a first prediction of the material response and the σ σ u σ ps σ es σ yl σ fl ε 1 ε i  [36] expected location of the most serious fatigue damage under a specific contact load and geometry, it is reliable to predict the location of the initiation of the RCF crack based on the results of various wheel/rail models. Some researchers have used the input parameter (normal contact pressure) and the output parameter (friction coefficient) from experimental studies to estimate the material response using the shakedown map [41,42]. However, as the field dynamic characteristics (traction force, track excitation, etc.) are difficult to determine in the experimental study, the finite element method (FEM) and multi-body system (MBS) analysis are the common and valuable tools used to analyze the wheel-rail contact problems in recent research studies.

Application of shakedown map
The 2D FE model is the basic and simplified wheel/rail simulation model established in ABAQUS as shown in Fig. 6 according to Ref. [19]. The model only includes a rail disc sample, driven by a moving load which represents the function of the wheel disc. The full-slip contact condition was considered because it represents more damaging to the material, and no changes in the contact patch caused by plasticity were introduced in the model. Different combinations of contact pressure and friction coefficient are performed by simulations to obtain elastic shakedown, plastic shakedown and ratcheting regions in the shakedown map (Fig. 4). Then, the crack initiation criteria (e.g., crack initiation angle) are evaluated in the different regions of the shakedown map [19]. Comparing with 2D FEM, the 3D model is closer to the real wheel and rail which can better express the stress state and the contact patch properties of the rolling contact. Figure 7a displays a 3D FEM of the wheel/rail interface proposed by Vo [33] in the finite element code ANASYS/LS-DYNA. Due to the consideration of the wheel vibration, wheelset-vehicle and the track-vehicle interaction in wheel/rail contact stress were analyzed since they might increase the magnitude of the stress. The wheel was simulated as a rigid body with a rotational speed under a vertical load exerted on it. Also, a constant coefficient of friction was used. Then the properties of wheel/rail materials and track components were applied into the FEM. It was necessary for the mesh size of the model to be fully taken into consideration as it would affect the accuracy of the simulation results [1]. The simulation results for the maximum contact load p 0 and the friction coefficient (0.45) were input together into the shakedown map ( Fig. 5) and plotted in Fig. 7b. This showed that all the investigated point fall into the ratcheting region (the area that contains the triangle in red). This meant that the damage tended to occur on the material, especially on the surface rather than in the sub-surface.
The main purpose of introducing the multi-body system (MBS) into the wheel-rail RCF analysis is that the vehicle dynamic characteristic can be taken into account in MBS, involving some parameters of vehicle components (e.g., mass, moment of inertia, suspension parameters, wheel profiles) and track parameters (curve radius, super elevation, rail profiles, track excitation). Meanwhile, the high efficiency and time-saving of MBS make it suitable to find a working condition which can lead to RCF initiation. The creep rate, creep force, normal force, contact patch size and maximum contact stress obtained by MBS provide the basic data for wheel/rail RCF prediction modeling. Lu [43] and Wang et al. [44] analyzed the location that is more likely to initiate cracks based on a vehicle system dynamic model using the dynamics software SIMPACK. Various 2D FE simulation model of wheel-rail contact (v is the moving direction of load distribution p, and p is the traction in the contact) [19] vehicle and track parameters were included in this model. The parameter of the abscissa on the shakedown map is the traction coefficient , calculated by Eq. (1). Maximum contact pressure P 0 , F x and F y are obtained from dynamic simulation while k is the constant value of the material.
where F x and F y are the vertical and lateral creep force in analyzed step respectively, and F n is the normal force between wheel and rail.
Anyway, the differences of the shakedown map input parameters (P 0 , k, ) in different research methods (experiment, 2D FE model, 3D FE model and MBS model) are summarized in Table 1. It also indicates the benefits and limitations of analyzing the location where cracks initiate by the shakedown map in different research methods.

Development of shakedown map
Although the shakedown maps proposed by Johnson and Ponter, respectively for 2D and 3D simulation models are commonly used in general wheel and rail RCF analysis, a further and deeper study of shakedown maps is needed. As a result of the nature of the shakedown phenomenon, both elastic and plastic shakedown limits are purely determined by elastic calculations, so that they are independent of materials [45]. Thus, the existing shakedown map is not suitable for some special situations such as different materials [38,46,47]. Some researchers have updated the shakedown maps to fit different situations using experimental  studies. Ringsberg et al. [46] modified the shakedown map by carrying out the twin-disc tests based on materials coated by 222 and 508 coatings. The traction coefficients obtained from the experiment when the eddy probe detected a crack under various input normal contact pressure were analyzed and the updated elastic and elastic shakedown limits for these coated materials were determined as shown in Fig. 8. Some researchers investigated the updating of shakedown maps using numerical calculation tools while others developed shakedown maps using simulation methods. Zhuang et al. [47] developed a shakedown solution for cohesivefriction materials (railway and pavement subgrade) under consideration of a general cohesive-friction yield criterion (as opposed to the purely cohesive criteria used in metal plasticity) based on the results from FE simulation and their updated shakedown map is shown in Fig. 9. Hasan [38] pointed out that the material characteristics of different steels (e.g., high-strength steel, heat-treated pearlitic steel, work-hardened steel, bainitic steel, etc.) are different from ordinary steels, so it is necessary to develop a shakedown map for a specific steel to help comparisons of steels to make a choice. The two commonly used shakedown maps are both calculated from numerical simulations, without any experimental parameters or empirical tips from the field. Therefore, in order to improve the reliability of the shakedown map, validation is needed using some experimental analysis. In recent studies, instead of using the previous application of Melan or Koiter theorems, Foletti and Desimone [48] reproduced the real stress-strain history of the material that had been subjected to rolling contact loads. Then, a new shakedown map ( Fig. 10) was established based on the previous one using the numerical simulation method. Comparing the new shakedown map with Johnson's, the agreement between the two shakedown solutions was complete for high friction coefficients (larger than =0.35). However, the new shakedown solution was between the linear kinematical hardening and elastic perfectly plastic solutions whereas both of Johnson's shakedown maps were in middle friction coefficients range (0 < <0.35). Therefore, in the future study on shakedown maps, this analysis regarding the validation of two widely used shakedown maps is significant. It means that the above experimental methods about establishing a new shakedown map can be used as references to determine whether points on the existing shakedown limit line are accurate or not during wheel/rail operation. Also, validation may be accompanied by some numerical simulation methods, for example, determining whether the stress-strain curve of the point is corresponding to the curve in Fig. 3 for different material response regions [18].
From this overview of shakedown maps, in addition to the function of predicting the location of crack initiation during wheel/rail rolling contact, two more applications  10 Comparison between the shakedown map calculated using the real rolling contact load with the traditional map using Melan theorem [48] of shakedown maps have been proposed. One is that it can qualitatively analyze the crack initiation combined with the existing damage function [43,44], the other is to provide the response region of the material to select the crack initiation life calculation method of the corresponding region [42].

Prediction of crack initiation life
According to the branch of the flowchart about the crack initiation life prediction shown in Fig. 2, four main determined items are included in the prediction processes. The accurate crack initiation calculation can be carried out by combining the determination of a cyclic constitutive model, the crack initiation criterion and the critical plane. Details of each item are reviewed in this section.

Cyclic constitutive model
The reliability and competency of the material model for simulating material responses (Fig. 3) are significant for fatigue analysis. For the difficult task of a material model to handle long-term ratcheting behaviors of a material, a proper plasticity model is needed to predict the decaying ratcheting effect [37]. Armstrong and Frederick [49] applied a kinematic hardening rule (A-F rule) which is a hardening rule that specifies the changes in the yield condition. Combining the A-F rule with the von Mises yield criterion, a realistic material model can be obtained within the framework of plasticity. However, later studies illustrated that the A-F rule is always lacking the ability to model cyclic material behavior. Thus, three elastic-plastic material models based on the A-F rule are currently used by researchers to model ratcheting material response, namely the Chaboche model [50], the Bower model [51] and the J-S model [52]. The comparison of these three models is shown in Table 2.
The Chaboche material model is mostly utilized in the research publications [17-19, 23, 24, 53, 54], all of which coupled the isotropic and kinematic hardening law. The coupling means that the yield surface is free to move in stress space and change shape or size [17]. The isotropic hardening law is aimed at simulating the decay in the ratcheting rate shown in Eq. (2): where t+Δt y is the updated yield stress at time t + Δt ; 0 y is the initial yield stress; Q and b are the material constants; e −p(t+Δt) governs the accumulated effective plastic strain at time t + Δt [23].
Meanwhile, Chaboche put forward a decomposed nonlinear kinematic hardening rule (Eq. 3) to improve the three critical segments that failed in the A-F rule. The updated hardening rule includes three decomposed hardening rules, namely the initial high modulus at the onset of yielding ( 1 ), the constant modulus segment at a higher strain range ( 2 ) and the transient nonlinear segment ( where C i and i are material constants; |d p | is the magnitude of the plastic strain increment tensor [23,57]. However, in RCF studies, the material model is used in the numerical simulation to obtain the stress-strain state: where X is the back stress; C i and i are material parameters, namely initial hardening of material and magnitude of the ratcheting rate, respectively [17,53]. Nevertheless, Ekh et al. [56] pointed out that the Chaboche model cannot synchronously simulate either the nonlinear behavior or the ratcheting rate. Also, Bari and Hassan [57] found that the Chaboche model can predict uniaxial ratchetting responses well, but it over predicts the biaxial and multiaxial ratchetting responses. This motivated many researchers to develop more advanced constitutive models [55]. The Bower model and J-S model use only a nonlinear kinematic law, without an isotropic rule. For the Bower model, the evolution processes for the A-F rule is the same as those of Chaboche model, except for an additional state variable Y in the development of the back-stress in Eq. (5). Also, the stress-strain curve of the Bower model is similar with the Chaboche model which failed to display the material nonlinear behavior and ratcheting rate.
where Y is obligated for the decreasing ratcheting rate which finally disappears with continued cycling; 2 determines the reduction in ratcheting rate with repeated cycles [17].
For the J-S model, the A-F rule has been modified by Jiang and Sehitoglu [52] to include several back-stresses X (i) : The total back-stress X is obtained by summation of all the back-stresses X (i) : where r (i) , c (i) and x (i) 0 (i = 1, 2,…, M) are material parameters obtained from experimental results [17].
Furthermore, the J-S model illustrates the best representation of both the shape and width of the stress-stain curve and the ratcheting rate. Therefore, Ekh et al. [56] mentioned that a more complex model such as the J-S model should be used in the RCF investigation in further studies. In addition, it should be noted that these material models had already been coded in the software ABAQUS. Other constitutive models can also be coded in the software using a Fortran subroutine if it is needed as Radim et al. [42] had already done.

Crack initiation criteria
In the wheel/rail rolling contact process, a specific region of material under the contact surface is subjected to nonproportional multi-axial load cycles. Factors such as the axle loads, wheel and rail profiles and friction coefficient can significantly influence the contact stress, material response and crack initiation during wheel/rail rolling and the fatigue life system is determined by the combination of these factors [15]. Additionally, the plastic flow of the material is concluded by the stress at the surface [19]. Therefore, a proper theoretical life prediction criterion is needed. As mentioned in the above section, the life prediction model differs in various material response regions.

Plastic shakedown
Various crack initiation criteria were used in the plastic shakedown region. The core of these criteria is to study how they perform when they are applied on RCF problems. The comparisons of criteria are illustrated in Table 3.
A simple form of SWT (Smith Watson Topper) damage parameter for predicting life under uni-axial loading was first proposed by Smith et al. [58] in 1976. However, the SWT parameter was used after a series of modifications in multiaxial proportional and non-proportional loading for materials with cracks. The most popular form of the SWT model is mentioned by Socie [34] with the consideration of mean stress as described in Eq. (8). The material plane with the largest value of the fatigue parameter, FP max , is recognized as the critical plane.
where max is the maximum normal stress in the plane with normal n (critical plane which is mentioned in Sect. 4.3); Δ 1 2 is the maximum principal strain amplitude (Fig. 3); � f is the fatigue strength coefficient; � f is the fatigue ductility coefficient; E is Young's modulus; b is the fatigue strength exponent; c is the fatigue ductility exponent; N f is the number of cycles to failure [59,60].
Therefore, a modification of max 1 was suggested by Beste in 1981. Only the stress amplitude Δ 1 2 is included and the compressive mean stress m 1 is neglected. This criterion was used as SWT/BESTE for predicting the crack initiation life [24].
Later in 1988, the additional hardening and mean stress effects were taken into account through the normal stress by Fatemi and Socie [61] based on axial fatigue data: where Δ max is the maximum shear strain range; max n is the maximum normal stress on the plane of Δ max ; y is the yield strength; k is a constant obtained from axial and torsional fatigue experimental data.
An adapted Coffin-Manson equation Eq. (10) was developed by Bannantine et al. [62] in 1990. In order to obtain the life predictions where material strain governs the life rather than stress, they calculated Δ 2 on the plane with the largest range of cyclic shear strain.
where Δ is the change in shear strain; � f is the shear fatigue strength coefficient; � f is the shear fatigue ductility coefficient; G is the shear elastic modulus.
In order to predict the crack initiation thoroughly, in 1999, Chen et al. [63] postulated that both normal and shear components of stress and strain on the critical plane is significant in the study of material damage [64]. Two criteria (CXH equations) about material characterized by normal fracture (mode I) and shear type failure (mode II) are described in Eq. (11).
where Δ max 1 is the maximum principal strain range; Δ 1 , Δ 1 and Δ 1 are the shear stress range that occur on the maximum principal strain range plane, respectively; Δ max is the maximum shear strain range; Δ , Δ n and Δ n are the shear stress range, normal strain range and normal stress range on the maximum shear strain plane, respectively [64].
Also in the year 1999, Jiang and Sehitoglu [34] proposed a crack initiation model which has a total strain energy density based formulation. It indicates that both normal and shear components of stress and strain contribute to the damage of the material on the critical plane. The model is expressed as where ⟨⟩ is the MacCauley bracket, and ⟨x⟩ = 0.5(�x� + x) ; J is a constant from tests and commonly use J = 0.2 [34].
Note that a proper selection of J plays a vital role in modeling a material that displaying mixed normal and shear cracking. For instance, J = 0 reduces the criterion to the SWT model for the normal cracking behavior. And the model depicts shear cracking fatigue when J is a large number. Meanwhile, Jiang and Sehitoglu [34] came up with the relationship between life N f and FP as where FP 0 , m, and C are material fatigue properties determined from baseline experiments [23,67].
Although some literature publications [23,29] predict crack initiation life using this original relationship, researchers prefer to combine the J-S model with SWT criterion [17,18,23,25,28,29,35,36,42,54,55,[70][71][72][73][74]. This updated criterion integrated the advantages in the above criterion. Also, it take into account of modeIand mode II fracture type according to CXH criterion. Therefore, the crack initiation life is more often calculated by where E and G are the elastic and shear modulus, respectively;

Ratcheting
The model for predicting the crack initiation life in the ratcheting region was postulated by Kapoor [32]. The number of cycles to crack initiation is calculated as where c is a material constant from experiments; Δ i and Δ i are the incremental normal and shear strain on the crack plane per load cycle. One of the material responses mentioned above occurs when the material is subjected to cyclic constant amplitude loading. The prediction of the RCF life to crack initiation is determined by one of these criteria related with the material response clarified [28].

Determination of critical plane
The critical plane method is proposed on the basis of physical observation of crack initiation and propagation, which has great advantages among other prediction methods. In the crack initiation stage, the crack surface is only affected by shear stress and shear strain, but part of the shear force is borne by the friction crack surface. Therefore, the shear force at the crack tip is reduced, and the crack growth progresses relatively slow. However, when the crack surface is affected by both shear stress and tensile stress, the crack surface will be separated by the tensile stress. In this case, the friction between the crack surfaces disappears, and all the shear forces are applied on the crack tip, then the cracks grow faster [70]. According to this mechanism, the critical plane method determines a physical quantity as the fatigue parameter (FP), and considers that the plane where the maximum fatigue parameter (FP max ) located is the critical plane. The demonstrations of critical angles generally can be displayed in different forms as shown in Fig. 11. The exact process of critical plane determination is displayed in Fig. 12.
All the stress and strain at the critical plane where FP max occurs are obtained from the simulation results of FEA. Except few literatures investigate the RCF using 2D FEM [17,18,75], most of researchers focus on 3D FEM which can better express the characteristics of real wheel and rail. Among those 3D FEA, the normal Hertz contact theory is utilized in most researches, while some researchers pointed out that the Hertz theory is not fit for the wheel/rail contact, especially in the curved sections [23,25,28]. Even if the wheel/rail contact surface is closely combined, Hertz theory may also fail to cause the non-conforming contact. Moreover, the large plastic deformation occurs as the result of severe friction on the contact patch, the highly concentrated stresses in the area near the contact patch and the non-proportional multi-axial loading subjected to the material lying under the contact region all make it difficult to determine the category and portion of stresses and strains which used in the evaluation of FP for crack initiation lifetime prediction [23].  Therefore, Akama tried to solve the full contact conditions between wheel and rail by applying constraints in the FEM simulation [28]. Also, Akama and his co-orkers [23] directly use the contact area and pressure distribution obtained from the FE simulation results of a certain step during rolling contact process, which are much different from that predicted by the traditional Hertz theory. Furthermore, Ringsberg et al. [25] used rail characteristics from the rail that have been in service for 16 months to study the non-Hertzian distribution condition, and compared the results with a newly manufactured rail that represented the Hertzian distribution condition.  [54] The comparison of normal contact pressure distribution for Hertzian and non-Hertzian condition is shown in Fig. 13.
Approaches of finding the critical plane where FP reaches the maximum value based on the stress and strain obtained from FEA is the key point of the life prediction. Most researchers preferred to determine the critical plane by surveying all the possible planes at a material point through a tensor rotation for the stress and strain [17, 18, 24, 25, 28, 29, 42, 53-55, 74, 75]. Figure 14 displays the schematic diagram of the tensor rotation plane. It indicates that critical plane approach incorporates the effects of all stress components in three dimensions. In tensor analysis, the spatial rotations of the stress matrices are accomplished by the application of Euler angles and mathematical transformations. Angle about the z-axis refers to the first rotation as shown in Fig. 14a, followed by the second rotation of angle happens about the new x-axis (which has already rotated automatically due to the first rotation), then the new-rotated coordinate system is generated after the last rotation of angle about the new y-axis. The transformation matrix of the stress tensor is related to the Euler angles , and as shown in Eqs. (16)- (19) [54]. However, some researchers [21,24,72] simplified the transformation matrix of the stress tensor, and the schematic diagram of the rotation plane displays in Fig. 15.
According to the above formula, there are numerous choices for the plane of the specific point on the material due to the variation of spatial angle. Among all the planes calculated by the criteria, the critical plane has the largest cos cos cos cos + sin cos sin sin sin − cos sin cos cos sin cos cos − sin sin sin sin sin − sin sin cos sin −sin cos cos cos FP value. Figure 16 shows a contour plot of the variation of FP under different material planes for the most critical point obtained by El-sayed [55]. Thus, the exact position of critical plane can be determined. However, Dollevoet et al. [54] pointed out that suggesting the obtained ranges of the critical plane angle seems to be more rational.

Crack initiation calculation
Due to the different life prediction model in various material response regions, this part reviews the crack initiation calculation in plastic shakedown region and ratcheting region.

Crack initiation in plastic shakedown region
In many publications [17-19, 23, 24, 28, 29, 53-55, 66, 72, 74, 75], after the determination of FP max , the crack initiation life is directly calculated using the prediction criteria mentioned above, and the value of stresses or strains needed in the criteria are consistent with the value in the critical plane. It is obvious that the maximum value of FP corresponds to the shortest life to crack initiation according to Ref. [29]. Otherwise, some researchers [25, 42, 71-73, 76, 77] put the fatigue damage effect into the consideration of the period to crack initiation prediction. Radim et al. [42] and Ringsberg and Josefson [25] assumed that fatigue damage is accumulated linearly, and the damage within one cycle as a result of RCF can be expressed as where FP 0 , K, and m are the material fatigue constants.
Furthermore, the evolution of the RCF cracks and the growth of wheel/rail wear reducing the material and changing the wheel/rail profiles. Thus, the coexistence of RCF and wear should be considered together. It should be mentioned that during simultaneous wear and RCF preexisting cracks will lead to stress concentrations which will influence the crack initiation parameter. The Archard model is a wear model that is commonly being used in research activities [78,79]: where V m is the volume of material worn away; D is the sliding distance; N is the wheel-rail normal force; H is the hardness of the material; K is a wear coefficient [28,80].
Bassevile and Gailletaud [65] investigated the evaluation of the competition between wear and crack initiation for alloy with Archard wear model based on SWT criterion. Zhou and his co-authors [71][72][73]76] proposed a calculation procedure of the period to crack initiation combined with fatigue damage effect as shown in Fig. 17. Firstly, the distribution of contact pressure and shear stresses were calculated using the software CONTACT based on the dynamic simulation results. And the damage analysis was according to the Archard damage theory. When wheels pass through, the equivalent wear of any point in the contact spot of the rail surface is calculated using Eq. (22): where Δz y is the abrasion amount of the longitudinal shadow cell of the contact spot when the abscissa is y; Δz(x, y) is the abrasion at cell (x, y); m and l are the number of vertical and horizontal cells of contact spot, respectively. Then, each rail profile and the distribution of the normal pressure and tangential traction are implemented on the 3D rail FE model. The FE model includes a long global rail and four wheels in two vehicle bogie. And the stress/strain calculated from FE model is used to estimate the crack initiation life. Later, once the accumulated wear at any point of the rail head reaches a certain value, the accumulated wear within the sliding area of the contact spot is superposed on the rail profile. And the profile is smoothed with cubic interpolation spline curve to obtain the wear profile [81]. Also, the accumulated number of wheel passing cycles n i in the phase of profile replacement is obtained. Then, the stage fatigue damage corresponding to wear stage n i is And according to the Miner theory [82], when the fatigue damage accumulation at any point of rail head in different wear stages reaches critical fatigue damage (Fig. 18), At this point, it is considered that the fatigue crack is initiated and the initiation life is where N f is the service life to crack initiation; and M is the replacement times of worn profile.
It should be noted that the Archard wear coefficient is being divided into four regions according to the sliding speed and contact pressure. And the specific value of each region needs to be determined by a large number of tests. Therefore, besides the Archard wear model, Bolton and Clayton [83] found that the wear rate of rail has a certain relationship with Tγ/A value (T tangential force at the wheel-rail interface; γ the creep coefficient; A the wheel-rail contact area) after rolling tests. Then, by analyzing a large number of wheel/rail material wear data from different test methods, Lewis et al. [84][85][86] established the relationship between wheel/rail material wear rate and Tγ/A value under different working conditions. It can be seen that the value of Tγ/A-wear rate curve changes continuously with the Tγ/A value, which could be more accurate than Archard wear model in simulation calculation [87]. Later, Meli and his co-authors [88][89][90] established an innovative model for the prediction of wheel/rail wear and RCF. Tγ/A was used and being improved to estimate the wheel/rail wear properly combined with the RCF model. The innovation model was used to update the wheel/rail profile, and provides a basis for the crack initiation life prediction when considering wear in future study.

Crack initiation in ratcheting region
Besides the plastic shakedown material response which caused by LCF as mentioned above, ratcheting effect is another main part that contribute to the fatigue damage. Generally, ratcheting deformation leads to the large shear deformation material microstructure, while the plastic shakedown produces less visible damage [15]. Therefore, Kapoor pointed out that LCF and ratcheting are the independent and competitive mechanisms for fatigue crack initiation. Equation (13) should be compared with a RCF crack initiation criterion described above [17]. It should be noted that Eq. (13) is aimed at the material that experiencing a constant ratcheting rate. However, for  Fig.18 Cumulative fatigue damage with rail profile evolution [76] most pearlitic rail steels, the ratcheting rate is typically not constant: when the cumulative deformation increases, the ratcheting rate decay can be observed. Then, a more general expression for the consideration of ratcheting-induced damage per cycle is proposed by Jiang [34] as where D r is the ratcheting damage; cri is a material constant about the ductility of the material; r is the orthogonal shear ratcheting strain; d r ∕dN represent the ratcheting rate.
Most of the literature publications [19,23,24,29] just predict the crack initiation life under ratcheting effect directly using the Kapoor criterion as Eq. (13) and compare the results with the LCF service life to crack initiation. The crack initiation criterion which is satisfied in the smaller number of cycles will make the material tend to failure. Ringsberg et al. [19] proposed a ratcheting model based on elastic analysis which is supplemented by the elastic-plastic material behavior based on empirical, and the plastic deformation in a single cycle is predicted: where ΔRd is the increment of plastic deformation in a single cycle; P r is the ratcheting load; p is the applied contact pressure; max is the maximum shear stress determined for each layer in the model; k ef f is the "current" shear yield stress of the material; C is a material constant which is independent of stress [91]. For each pass of the load (cycle), the increment of one-way strain ( ΔRd ) is added to the total cumulative plastic deformation Rd. Failure is defined as Rd exceeding the critical value of strain c determined by experiments. Furthermore, the results from this ratcheting model agree well with the results from the calculation using Kapoor criterion.
Other researchers [25,42,67] tended to take the ratcheting damage into account in the prediction of crack initiation life according to the assumption made by Jiang and Sehitoglo as shown in Eq. (26). If the material has low plasticity, LCF may occur in the form of crack initiation and propagation in a limited number of cycles. For ductile materials, extensive ratcheting strains can lead to the extrusion of thin surface slices. In both cases, the combination of ratcheting and fatigue can be used as a mean of failure prediction under these conditions. Jiang and Sehitoglo [34] proposed that the total damage is the sum of fatigue damage D RCF and ratcheting damage D r (Eq. (28)). When the total damage D reaches a unity, material fails.
However, Ringsberg and Josefson [25] mentioned that the damage of two mechanisms (LCF and ratcheting) should not be added in the same load cycle. It is accentuated if the damage caused by the two mechanisms is close in magnitude, which will lead to over conservative life prediction as Jiang and Sehitoglo did. Then a competitive approach is put forward in the summation of total damage. The largest damage component for the current load cycle is governed by the contribution of either RCF or ratcheting. The damage is assumed to accumulate linearly for every load cycle as Material failure occurs when damage D reaches unity. This prediction approach seems more reliable and comprehensive, and is often used in the service life to crack initiation prediction. Radim et al. [42] mentioned that the comparison of ratcheting failure life and LCF failure life is needed in a life study. Once the former is bigger than the latter, material would experience the LCF failure life. And they came to the conclusion that surface and subsurface cracks were mainly caused by ratcheting and LCF, respectively.
Hence, if the shakedown map is used in the first step of crack initiation prediction as mentioned in Sect. 3; then it can be clear that the material is under the ratcheting effect or plastic deformation, and the accurate criterion will be used in the corresponding region.

Validation of crack initiation life prediction
In order to confirm the validity of the numerical calculation method, most of the researchers chose to use the experimental results (the location or service life to crack initiation) to validate the reliability of RCF crack initiation prediction methods [15,25,28,42,69]. Radim et al. [42] pointed out that the maximum crack depth remained stable at about 43 μm after 30,000 cycles from the experimental analysis. This is similar with the numerical calculation results that indicated 30,858 cycles to crack initiation which occurred at the node placed 50 μm beneath the material surface. Zeng et al. [92] used a scaling factor arising from the relation between the actual wheel and twin-test sample dimensions to predict the actual wheel-rail RCF behavior based on the numerical simulation. And the morphology of RCF cracking of the defected rail material in the initiation stage from experimental analysis (Fig. 19) was compared with the evolution process of von Mises stress that occurred on the region around the defect in the simulation study (Fig. 20), where it can be seen that the crack that initiated near the middle edge of the defect was mainly caused by the maximum frictional shear stress. Also, some researchers have tried to compare the simulation results with the field observation [23,28,[53][54][55]. Xin et al. [53] plotted the stress distributions on the wing rail and whole crossing as shown in Fig. 21, it can be seen that the highest stress occurred at the position of 507 mm from the nose point in the simulation study, which is similar as in the field. Akama et al. [23] found that the predicted initiation cycles from numerical analysis was included in the range obtained from field observation. And the predicted crack angles were similar to those crack orientations of initiated head checks on the railhead in the field. However, Farley et al. [69] proposed that there exist differences between numerical calculation results and experimental results. The crack initiation life and location of the critical plane using the Coffin-Manson prediction criterion were approximately 500,000 cycles and 500 μm, respectively, while the experimental results were 100,000 cycles and 0.5-1 μm, respectively. The differing results pointed out that it is necessary to put forward a correction coefficient for   Fig. 21 Comparison of the plastic strain region in the simulation and experiment [53] the simulation results to fit the test results, and also to choose a proper life prediction criterion which is more similar with the experimental analysis. In order to select the most suitable calculation method from various criteria of RCF life to crack initiation, some researchers [19,60,64] compared the calculation results using different life criteria with the experimental results. Then the criterion closest to the experimental results was chosen. Ringsberg et al. [19] compared the initiation life calculated by the SWT/Beste approach, the Coffin-Manson relation, the Kapoor criterion and the ratcheting model with the results from the corresponding twin disc test. They pointed out the SWT/Beste approach was the most reliable one. Han et al. [64] examined six damage criteria for their ability to predict lives in different working conditions. The comparison between predicted and experimental lives for different criteria indicates that the FS criterion gives good life estimation. Also, Lykins et al. [60] pointed that both SWT and FS criteria predicted the failure location near the trailing edge of contact, which correlated well with the experimental results. Thus, it can also be concluded that the chosen crack criterion is different in each study, hence there is no unified standard.

Challenges and research gap in the prediction of RCF crack initiation
This review has described different numerical calculation methods of predicting the initiation of RCF cracks, mainly including the shakedown map and the critical plane theory with various criteria. The shakedown map is aimed at analyzing the region of crack initiation in the model and the material response, while the critical plane theory is for determining the service life to crack initiation. With these methods being utilized for wheel/rail system applications, these prediction approaches have been verified under normal wheel/rail contact conditions, or even for a new wheel/rail model. In order to obtain calculation results that are much more similar to the real wheel/rail contact scenarios, five main aspects should be taken into account: the wheel/rail dynamics boundary conditions, the complex working conditions, the surface defects, wear and the combination of shakedown map. Thus, further systematical and accurate studies about developing RCF crack initiation prediction approaches and characteristics of RCF cracks are needed. The multi-body system dynamics of wheel/rail contact used by researchers, which aimed at obtaining boundary conditions to reflect the actual wheel/rail rolling contact, can be divided into two types, namely vehicle system dynamics [23, 71-73, 76, 83] and train system dynamics [70]. An important feature that distinguishes vehicle system dynamics from train system dynamics is the consideration of dynamic in-train forces in the latter [93]. Thus, the differences of boundary conditions obtained from dynamics which are used as the input parameters in FE simulations should be mentioned. And the prediction of RCF crack initiation based on the contact stress distribution carried out by FEA will be influenced. Due to the various influence factors as mentioned before and the differences between the two types of dynamics, the research about considering dynamic characteristics is insufficient and non-systematic, which should be perfected in the further study.
Meanwhile, the complex working conditions are also a significant factor that affects the prediction of RCF crack initiation. As for shakedown map analysis, some researchers have already pointed out that the commonly used shakedown map is not fit for some special conditions such as coated material [46], cohesive material [47] or material under water lubrication conditions [41]. This indicates that there may be different shakedown maps corresponding to different working conditions, and the material response will be different from that using the common shakedown map. Therefore, it is important to validate the accuracy and update the 2D and 3D shakedown maps which are generally being used by researchers under specific working condition in the further study. In addition, the input parameters in the FE simulation are different for various working conditions. This leads to changes in the distribution of stresses and strains under the wheel/rail rolling contact process, and finally affects the determination of the crack plane and the calculation of the crack initiation life. Then, according to the literature publications mentioned in Sect. 4.5, in the future study, an appropriate life prediction criterion for a specific rolling contact process should be selected after comparing various life calculation results with the corresponding experimental results or observations in the field.
Furthermore, the occurrence of surface defects (e.g., squats [28,29,94]) on the interface between wheel and rail during the rolling contact process can significantly lead to RCF. Also, the size and location of defects [18] may affect the RCF crack initiation development of wheel/rail contact. Recent studies are mainly focused on the inner sub-surface defects [18,95,96] and the outer surface defects [92,[97][98][99][100]. These researchers analyzed the critical crack plane and the crack initiation life of the specific point around the defect. Further, the microstructure of the material around the defect in the cracking initiation stage in the experiment is compared with the simulation results as undertaken in [92,97]. However, the simulation analysis of material microstructure cannot clearly illustrate the crack initiation stage. Thus, it should focus on the evolution of crack initiation around the defects according to the microstructure in the cross-section of material samples from experiments undertaken in the investigation and improve the differences between test results and simulation results.
Also, wear plays an important role in the wheel/rail rolling contact process, which is competitive with RCF. Cantini et al. [101] mentioned that once the RCF cracks exist, they at the beginning tend to propagate faster than the wear. But then, the propagation speed of RCF cracks comes down while the crack gets deeper. If the wear rate is sufficient relative to the crack propagation speed, cracks will tend to disappear; if this does not occur, the crack will propagate deeper ending with a progressive detachment [101]. Therefore, the consideration of the wear rate in the study of RCF crack initiation prediction is essential. As the Sect. 4.4.1 mentioned, besides the traditional Archard wear model, the Tγ/A model could be a more appropriate wear model to reflect the wheel/ rail wear process.
Finally, because the calculation criteria of crack initiation life in the ratcheting effect region and the plastic deformation region are different, the crack initiation life can be calculated more efficiently and accurately by determining in advance which region the studied material belongs to. Then, this requirement can be realized by introducing a shakedown map to judge the region of the material before calculating the crack initiation life. Also, if the material is in the elastic or elastic deformation region in the shakedown map, then there is no need to calculate the crack initiation life.
The gap for this review is therefore to propose a prediction process for the crack initiation prediction method at the wheel/rail rolling contact. Thus, further study should focus on the following points: • Simulation models, both FEA and multi-body dynamics, should pay more attention to the dynamic characteristics and behavior of the field wheel/rail system. Such a consideration is limited in recent studies; • The optimal selection of a RCF crack initiation prediction approach in various working conditions needs further systematical study;

Fig. 22
Proposed prediction process of RCF crack initiation • It is necessary to determine the material response region based on the shakedown map before the undertaking the crack initiation life prediction.

A proposed prediction process of RCF crack initiation
For designing a reliable and effective process of RCF crack initiation prediction, a new and complete process is shown in Fig. 22. Details of each step are discussed in the following: Step 1 The simulation begins with the vehicle-track dynamic model with the standard unworn rail profile (i = 1) and the force, position, a-axis and b-axis of the wheel-rail contact patches are calculated. Also, the distribution of the normal pressure and tangential traction in the contact patches are obtained. It should be mentioned that the rail profile with irregularities or defects can also be used for the initial profile in this step to study some special working conditions; Step 2 Three-dimensional FE models of the rail in every profile (i = 1, 2, …,n) are established using the explicit mode in the Abaqus software, supplied with the distribution of the normal pressure and tangential traction. The stresses and strains at each point in the FE model are calculated; Step 3 Put the maximum normal contact stress combined with the shear yield strength and friction coefficient together into the shakedown map and determine which region of the shakedown map it will be located in. If the point falls in the crack initiation region (ratcheting and plastic deformation), the prediction process is linked to step 4. If it falls in the elastic and elastic shakedown region, the prediction process of RCF crack initiation is complete; Step 4 If the analyzed point falls in the plastic deformation region in the shakedown map, then this links to step 5. If the point falls in the ratcheting region, then it relates to step 7; Step 5 The stress/strain at all points in the rail head with every profile are applied to various criteria mentioned in Sect. 4.2 (C LCF1 , C LCF2 , …,C LCFn ) to identify the critical plane and the FP max ; Step 6 Each crack initiation life N ij corresponding to any point j of rail head at the ith profile according to various criteria is calculated. Then the fatigue damage D ij is calculated. Also, with the development of wear, when the rail profile develops from the first profile to the M profile, the cumulative fatigue damage D j of the rail head point j is obtained.
When the cumulative fatigue damage of point j reaches the critical fatigue damage (taken as 1 in this process), it is considered that the fatigue crack originates at point j, and the crack initiation life N f,RCF1 , N f, RCF2 , …, N f,RCFn corresponding to each criterion is calculated. Otherwise, combined with the Archard wear model and Tγ/A wear model, the wheel cycles passing number is updated and the rail profile is replaced by the worn profile (i = i + 1). Once the cumulative fatigue damage at each point of the rail head reaches the fatigue damage criterion, the crack is expected to initiate at the relevant point, and the crack initiation life N f,LRF1 , N f,LCF2 , …, N f,LCFn corresponding to each criterion is calculated. Mentioned that the cumulative fatigue damage should be calculated using various methods mentioned in Sect. 4.4.2, and compared the different calculation results with the corresponding experimental or field study and then determined an optimal method; Step 7 If the analyzed point falls in the ratcheting region, then the ratcheting life N f,R1 , N f,R2 , …, N f,Rn are calculated by the various ratcheting criteria (C R1 , C R2 , …,C Rn ) mentioned in Sect. 4.5. Also, the crack initiation life should be calculated using the plastic deformation criterion to compare which criterion was proper in different cases [25]. However, further study could focus on the crack initiation life prediction once the analyzed point is at the transition region between ratcheting and plastic deformation.
Step 8 After comparing the relevant experimental results (the working condition of wheel-rail in the experiment was the same with that in the simulation) with both the life (N f,LCF1 , N f,LCF2 , …, N f,LCFn ) and the ratcheting life (N f,R1 , N f,R2 , …,N f,Rn ), optimal criterion C RCFn or ratcheting criterion C Rn can be obtained.

Conclusion
A review of existing methods used for RCF crack initiation prediction is presented in this paper, particularly the crack initiation region detection and the crack initiation life calculation. The shakedown map is used to investigate the material response during the rolling contact process and the exact region of crack initiation. Meanwhile, four main steps are included in the crack initiation life prediction procedure: the cyclic constitutive model, the crack initiation criteria, the determination of the critical plane and the crack initiation calculation methods.
The shakedown map is commonly used to predict the location of crack initiation during wheel/rail rolling contact. Also, it can analyze the crack initiation combining with the damage function in a qualitative way. Four material response regions are included in the shakedown map, and the crack initiation life calculation method is different in these different regions.
As for the prediction of crack initiation life, life of material in the ratchetting effect region is calculated using the Kapoor criterion. When the material is in the plastic shakedown region, combined with the constitutive model and the distribution of stresses and strains obtained from simulation, the critical plane (FP max ) of material is calculated by various criteria in the plastic shakedown region. Among the various criteria, the criterion combining the J-S and SWT models is the most preferable criterion in research studies. Then the crack initiation life is calculated using the criteria on the critical plane. Besides, the method of considering the wear during rolling contact process is aimed at analyzing the coexistence for wear and RCF.
This review also proposed a developed prediction process of RCF crack initiation. The new prediction process considered the research gaps in the previous studies: the wheel/rail dynamics boundary conditions, the complex working conditions, the irregular profile (surface defects), and wear on wheel/rail operating surfaces. All of these factors would affect the crack initiation life in a great extent.