Effect of Thermal Load Caused by Tread Braking on Crack Propagation in Railway Wheels on Long Downhill Ramps

Abstract

To investigate the propagation behavior of thermal cracks on the wheel tread under the conditions of long downhill ramps, a three-dimensional finite element model of a 1/16 wheel, including an initial thermal crack, was developed using the finite element software ANSYS 17.0. The loading scenarios considered include mechanical wheel–rail loads, both with and without the superposition of thermal wheel–brake shoe friction loads. The virtual crack closure method (VCCM) is employed to analyze the variations in stress intensity factors (SIFs) for Modes I, II, and III (K_I, K_II, and K_III) at the 0°, mid, and 90° positions along the crack tip. The simulation results show that temperature is a critical factor for the propagation of thermal cracks. Among the SIFs, K_II (Mode II) is larger than K_I (Mode I) and K_III (Mode III). Specifically, the thermal load on the wheel tread during braking contributes up to 23.83% to K_II when the wheel tread reaches the martensitic phase transition temperature due to brake failure. These results are consistent with the observed radial propagation of thermal cracks in wheel treads under operational conditions.

1. Introduction

In regions characterized by mountainous landscapes and large differences in altitude, steep gradients are inevitably included in the design of railway routes. A prime example is the Baoji–Chengdu Railway through western China, which has a maximum gradient of 30% [1]. On long ramps, rolling contact fatigue (RCF) can easily occur on the wheel due to the rolling contact load between the wheel and the rail [2,3]. Rolling contact fatigue can cause cracks to initiate and propagate, ultimately leading to defects such as spalling [4], headchecks [5], and squats [6].

Rolling contact fatigue has been studied by many scholars. Jiang et al. investigated the rolling contact fatigue of trains at high speeds [7]. The crack propagation direction and crack propagation angles were calculated by finite element simulation. The results of the study showed that the crack propagation rate was faster at high train speeds. Zhang et al. investigated the effect of slip ratio and contact stress on rolling contact fatigue using a wheel/rail twin-disc machine [8]. The results showed that the slip ratio was the main reason for the change in crack angle. Zhang et al. [9] reviewed the prediction methods for the initiation of rolling contact fatigue cracks and proposed a research gap for the prediction of rolling contact fatigue crack initiation. Krishna et al. [10] compared four damage models for predicting the rolling contact fatigue of rails, including the wedge model [11], the KTH model [12], the surface fatigue index [13], and the whole-life rail model (WLRM) [10]. Rolling contact fatigue was quantified, and a method was proposed to measure the degree of agreement between the different models. Maglio et al. [14] investigated the effect of wheel tread damage on wheel impact load by creating a simulation model and validating it using on-site measured data. The results quantified the effects of wheel tread defects and provided an engineering reference for the fatigue damage behavior of the wheel tread. In addition, researchers have also investigated the influence of normal contact pressure [15,16,17,18], coefficient of friction [19], and curvature radius [20,21] on the initiation of fatigue cracks during rolling contact. However, temperature, an important factor, also significantly influences fatigue damage.

As higher speeds and higher loads are required for both freight and passenger trains, the kinetic energy of the vehicles inevitably increases. In the trend of high-speed and heavy-duty development, vehicle speed control is crucial, and the braking system plays a central role. Known for its simplicity and cost-effectiveness, the tread brake [22] is often used in braking systems. However, during tread braking, considerable frictional heat is generated between the wheels and the brake shoes, causing temperatures to rise rapidly, sometimes exceeding 700 °C [23], which is sufficient to promote austenite phase transformation. The forms of temperature transfer include conduction, convection, and radiation [24,25]. Srivastava et al. [26] conducted a comprehensive study of the evolution of the temperature field during train braking and the effects of thermomechanical coupling loads on thermal fatigue damage. The fatigue damage resulting from the interaction of thermal loads with other loads includes hot spots [27], shelling [28], spalling [29], corrugation [30], thermal cracks [31], etc. Therefore, the effect of friction-induced temperature rise from tread braking on the thermal fatigue damage of wheels is particularly important.

As the wheel is heated due to tread braking, the fatigue strength of the material decreases. The wheel tread generates residual tensile stresses during the cooling cycle, which can lead to thermal and mechanical cracks. Peng et al. [32] developed a three-dimensional finite element model of the wheel incorporating heat transfer and fracture mechanics theories. Using a semi-analytical method, they calculated the stress intensity factor (SIF) of thermal cracks on the wheel tread. The Generalized Frost–Dugdale method was used to predict the propagation behavior of thermal fatigue cracks on the wheel tread under braking conditions. The simulation results showed that the tensile stress in the circumferential direction resulting from thermal cycling promotes the propagation of fatigue cracks. Caprioli et al. [33] considered the temperature dependence of the wheel material and investigated the behavior of the material under simultaneous thermal and mechanical stress. The mechanical contact load was characterized as a moving Hertzian contact stress pattern coupled with a surface shear stress distribution that adapts to either full or partial slip conditions. Investigations found that partial slip scenarios induced larger plastic deformations within a thin layer adjacent to the contact interface. Esmaeili et al. [34] studied the material behavior of ER7-grade steel railway wheels, produced by SSAB in Luleå, Sweden, under the combination of tread brake thermal load and rolling contact load. They proposed a viscoplastic material model and validated the model based on experimental data. Fletcher [35] studied crack propagation behavior under coupled thermal and mechanical loads. The magnitude of the stress intensity factor was calculated at different temperatures using boundary element modeling. When the temperature reached 1000 °C, the calculated stress intensity factor was close to the fracture toughness values of typical rail steels. Duflot [36] considered the coupling of thermo-mechanical loads to simulate two- and three-dimensional stationary cracks using the extended finite element method (XFEM). A domain-formulated interaction integral was used to extract stress intensity factors from the extended finite element method solution without the need for direct integration across the crack surface. In a study by Bayat et al. [37], the effects of thermal nonlinearity and temperature dependence of materials on crack propagation paths were investigated using the extended finite element method. A fatigue crack growth algorithm was developed in MATLAB 2018 to perform fatigue propagation life prediction and crack shape analysis using Linear Elastic Fracture Mechanics (LEFM). Liu et al. [38] investigated the propagation behavior of fatigue cracks in rail welds using finite element analysis and calculated the stress intensity factor. For fatigue life estimation and crack shape analysis, a MATLAB-based algorithm was developed that utilizes the principles of linear elastic fracture mechanics to predict fatigue crack propagation. Nejad et al. [39] calculated the propagation behavior of fatigue cracks using the boundary element method and estimated the propagation rate of fatigue cracks based on the modified Paris model. Handa et al. [40] used a test bench to reproduce thermal cracks in the wheel tread under simultaneous loading of continuous rolling contact and cyclic frictional heat, finding that thermal cracks were caused by tensile stress due to plastic deformation and thermal stress near the wheel tread surface. Handa and Morimoto [41] stated that just 100 braking cycles under the influence of a normal contact force of 30 kN were enough to induce the formation of thermal cracks within the wheel tread.

In summary, researchers have studied the thermomechanically coupled fatigue behavior of wheel treads using both experimental and numerical methods, yielding significant results. However, previous studies have only focused on conditions to a limited extent, and the level of thermal load on the wheel tread was low. When a train runs on a long ramp, the thermal load on the wheel tread becomes exceptionally high, leading to increased susceptibility to fatigue damage. Figure 1a shows an example of a freight wheel that broke on a long downhill ramp in China due to the significant thermal loads caused by braking. Figure 1b shows the thermal cracks generated on the wheel tread due to thermal loads. Upon inspection by the rolling stock company, a martensite white etching layer was found on the wheel tread where the thermal crack started. This means that the area where the crack source was located was exposed to strong thermal impact. The thermal crack rapidly propagated in the radial direction under thermomechanical coupling loads. The operating conditions are as follows: the vehicle braking speed is 60 km/h, and the average slope and length of the ramp are 12‰ and 18 km, respectively.

This study uses these conditions as the background for the research. In this work, the influence of wheel tread braking on the SIF at the thermal fatigue crack tip during the martensitic phase transition temperature is investigated using the virtual crack closure method. The aim of this article is to analyze the effect of tread brake thermal loads on the propagation of thermal fatigue cracks in wheel treads under long ramp conditions. The study also aims to provide a technical reference for understanding thermal crack propagation under such conditions using the above operational route as an example.

2. Materials and Methods

2.1. Material Parameters

The material of the wheel is S-660 type, produced by China Northern Heavy Industries Group Co., Ltd. (NHI) in Shenyang, China, and the wheel material parameters used for the calculations in this work are listed in

Table 1. Wheel material parameters.

Temperature (°C)	Specific Heat Capacity (J/(kg·°C))	Thermal Conductivity (W/(m·°C))	Coefficient of Thermal Expansion (10−6/°C)	Density (kg/m3)	Poisson’s Ratio
0	434.02	48.30	1.0645	7833	0.3
100	473.21	46.43	1.1293
200	512.40	44.56	1.1944
300	551.59	42.69	1.2589
400	590.77	40.83	1.3237
500	629.96	38.96	1.3885

and

Table 2. Wheel Young’s modulus.

Temperature (°C)	Young’s Modulus (GPa)
25	200
100	190
200	180
300	170
400	160
500	150
600	140
700	130
800	120

. The material parameters of the brake shoe are shown in

Table 3. Brake shoe material parameters.

Temperature (°C)	Specific Heat Capacity (J/(kg·°C))	Thermal Conductivity (W/(m·°C))	Young’s Modulus (MPa)	Density (kg/m3)	Poisson’s Ratio
27	520	0.9970	400	2340	0.28
200		0.9606

. The material parameters for the wheel and brake shoe were obtained through experiments conducted by the rolling stock. Since the wheels are made of carbon steel, under the working conditions studied in this work, the friction between the wheel and the brake shoe generates a significant amount of heat that exceeds the martensitic phase transition temperature. As the temperature increases, the distance between the atoms increases, leading to a decrease in the Young’s modulus of the wheel material. For carbon steel, the Young’s modulus decreases by approximately 3% to 5% for every 100 °C increase in temperature. The variation of the Young’s modulus of the wheel material with temperature is shown in Table 2.

2.2. Calculation Process

Figure 2 shows the flowchart of the simulation analysis. To more accurately reflect the operating conditions during vehicle braking, it is necessary to obtain the actual route conditions and braking operation conditions during vehicle operation. Route conditions include the gradient and length of the slope. The braking operation conditions include the vehicle speed, the braking force, and the geometry and material parameters of the wheel and brake shoe. To study the propagation of thermal cracks on the wheel tread, it is essential to obtain the geometry (length and depth) of the thermal cracks and their position on the wheel tread.

Based on the above data and parameters, a finite element model of a 1/16 wheel with initial thermal cracks is created using the finite element software ANSYS 17.0, developed by Ansys, Inc. in Canonsburg, PA, USA. The parameters of the long ramp routes, the geometry and material parameters of the wheel and the brake shoe, as well as the location and geometric parameters of the thermal cracks, are input into the model. Taking into account the actual operating and route conditions, SIMPACK 2020, developed by Dassault Systèmes in Vélizy-Villacoublay, France, and MATLAB 2018, developed by MathWorks, Inc. in Natick, MA, USA, are used for the simulation analysis to determine the boundary conditions (thermal analysis boundary, static analysis boundary). MATLAB 2018 is used to determine the thermal boundary conditions for the wheel during vehicle braking operations, including heat flux boundary conditions and convective heat transfer boundary conditions. In addition, SIMPACK 2020 is employed to determine the wheel–rail contact positions when the vehicle runs along the track, facilitating the definition of boundary conditions for the static analysis.

These boundary conditions are then applied to the finite element model to calculate the temperature field of the wheel tread under friction braking conditions between the wheel and the brake shoe using the friction power method. The stress field of the wheel tread is then simulated and analyzed assuming a mechanical load acting between the wheel and the rail. To reduce the computational complexity of the finite element model, a moving load approach [42,43] is utilized to simulate the contact behavior between the wheel and the rail. The magnitude of the mechanical load between the wheel and the rail is calculated based on Hertzian contact theory and Coulomb’s law of friction. Another loading condition taken into account is the simultaneous application of the thermal load caused by the friction between the wheel and the brake shoe and the mechanical load between the wheel and the rail. Finally, based on the obtained stress and displacement data at the crack tip of the finite element model, the thermal crack propagation behavior is analyzed with MATLAB 2018 using the virtual crack closure method.

This study aims to fill this gap by simulating and analyzing the effects of elevated thermal loads on the wheel tread during the braking process, especially under conditions where martensitic phase transitions occur. To better understand the influence of such extreme conditions on crack propagation, a finite element model of the wheel is developed that takes into account the interactions between thermal and mechanical loads. The virtual crack closure method is used to investigate how thermal fatigue cracks develop under simultaneous thermal and mechanical loading. This provides insights that could improve the safety and durability of wheels in similar operating environments.

2.3. Thermomechanical Coupling Theory

Due to the long braking duration of the vehicle, the thermomechanical coupling behavior of the wheel tread during braking is solved using a decoupled approach between thermal and mechanical loads. The relationship between thermal and mechanical loads is considered weak coupling. The thermomechanical weak coupling theory is outlined as follows.

An idealized nonpolar continuum representing a solid deformable body occupies an open domain Ω in a three-dimensional Euclidean space characterized by a regular boundary ∂Ω at time t > 0. The motion of an object is described in the Cartesian coordinate system. The finite element method is used to solve the weak coupling form of the thermomechanical coupling problems. The weak coupling formulations of the thermomechanical equations referred to in [44] are expressed by Equations (1) and (2).

$${\displaystyle {\int }_{\Omega }\rho }\cdot \dot{\mathbf{v}}•\mathit{\delta }\mathbf{v}dv+{\displaystyle {\int }_{\Omega }\mathsf{\sigma }}:\mathit{\delta }\mathbf{d}dv={\displaystyle {\int }_{\Omega }\mathbf{b}}•\mathit{\delta }\mathbf{v}dv+{\displaystyle {\int }_{\partial \Omega }\mathbf{t}}•\mathit{\delta }\mathbf{v}ds,$$

(1)

$$\begin{array}{l}{\displaystyle {\int }_{\Omega }\mathit{\delta }T\cdot \rho }\cdot \dot{\mathbf{v}}•\mathbf{v}dv+{\displaystyle {\int }_{\Omega }\left(\nabla \mathit{\delta }T\right)•\left(\mathsf{\sigma }•\mathbf{v}\right)}dv+{\displaystyle {\int }_{\Omega }\mathit{\delta }T\cdot \dot{e}}dv-{\displaystyle {\int }_{\Omega }\left(\nabla \mathit{\delta }T\right)•\mathbf{q}}dv=\hfill \\ {\displaystyle {\int }_{\Omega }\mathit{\delta }T\cdot \mathbf{b}}•\mathbf{v}dv+{\displaystyle {\int }_{\partial \Omega }\mathit{\delta }T\cdot \mathbf{t}}•\mathbf{v}ds+{\displaystyle {\int }_{\partial \Omega }\mathit{\delta }T\cdot q_n}ds+{\displaystyle {\int }_{\Omega }\mathit{\delta }T\cdot r}dv\hfill \end{array},$$

(2)

where ρ is the material density; δv is the virtual velocity vector; σ is the Cauchy stress tensor; δd is the virtual strain rate tensor; b is the body force vector; t is the surface traction vector; δT is the virtual temperature; $\dot{e}$ is the time derivative of the internal energy per unit volume; q is the heat flux vector; r is the heat generation rate per unit volume; and q_n is the normal heat flux.

2.4. Virtual Crack Closure Method

Rybicki and Kanninen [45] pioneered the virtual crack closure method to solve the stress field strength factor at the crack tip. Originally the method was applied to solve a two-dimensional crack problem; Shivakumar et al. [46] later extended the method to solve a three-dimensional crack problem. Figure 3 shows a schematic diagram of the mesh in and around the oblique crack tip node i. In Figure 3, j and j* denote the two overlapping nodes in the crack plane closest to the tip node i; h₁, h₂, l₁, and l₂ denote the dimensions of the mesh near node i. The specific calculation involves two steps. The first step begins with determining the strain energy release rate G_I from Equation (3). This equation uses a local coordinate system o′x′y′z′, where the origin is defined at node i and its o′x′y′ plane coincides with the crack plane.

$$G_{\mathrm{I}}=-{\displaystyle \frac{F_{z^{\prime }i}\Delta u_{z^{\prime }(j,j^{\ast })}}{2\overline{A}}}{\displaystyle \frac{h_1}{h_2}},$$

(3)

where F_z′i is the nodal force at the crack tip node i in the z′ direction. Δu_{z′(j, j*)} is the displacement difference between nodes j and j* in the z′ direction. $\overline{A}$ is the correction related to the mesh inhomogeneity.

Based on the strain energy release rate, the Mode I SIF K_I is then derived from Equation (4).

$$K_{\mathrm{I}}=\sqrt{E^{\ast }{G}_{\mathrm{I}}},$$

(4)

where E is the Young’s modulus of the material, E* = E for the plane stress state, E* = E/(1 − ν²) for the plane strain state, and ν is the Poisson’s ratio.

In practice, the wheel is subjected to a multiaxial stress state resulting from the complicated wheel–rail contact loads. Crack propagation typically manifests as a combination of two or more propagation modes that form a composite crack. Under complex loading conditions, the propagation rate of a composite crack is determined by the equivalent SIF K_eff. This factor is calculated using Equation (5).

$$K_{\mathrm{eff}}=\sqrt{K_1^2+K_{\mathrm{II}}^2+(1-\nu )K_{\mathrm{III}}^2},$$

(5)

3. Finite Element Modeling

3.1. Physical Modeling of Thermal Cracks

In this study, a finite element model is constructed based on the actual location of thermal crack formation on the wheel. For thermal analysis, the element type of the finite element model is SOLID 70, and for structural analysis, the SOLID 70 element is converted to SOLID 185. In a specific operating scenario, the thermal cracks in the wheel tread are located 43 to 81 mm from the wheel flange side, with a length of 38 mm and a depth of 1.5 mm. To simplify the calculations while focusing on the stress concentration zone, a 1/16 wheel finite element model is created with the initial thermal crack. For static analysis, full constraints are applied on both sides of the 1/16 part of the wheel to simulate real operating conditions. The crack region of the finite element model is refined, as shown in Figure 4. The mesh size in the circumferential and axial directions is 0.5 mm and 0.3 mm, respectively. The finite element model consists of 239,868 elements and 258,207 nodes.

Figure 5 shows a schematic diagram of the initial thermal crack. The thermal crack image in Figure 5a is from an actual field wheel. The thermal crack has a length of 38 mm, a depth of 1.5 mm, and a semi-elliptical shape. Based on actual wheel fractures, it was observed that rapid fracture occurred in the left half of the thermal crack. Therefore, the points P1 (0° position), P2 (center position), and P3 (90° position) in Figure 5b are selected to study the propagation behavior of the thermal crack. Since water or contaminants may be present at the original crack location, the friction coefficient between the crack surfaces is set to 0 to simulate the working conditions when water is present. The contact stiffness of the cracked surfaces is set to 4 to prevent contact penetration.

3.2. Boundary Conditions of Thermal Analysis

3.2.1. Heat Flux Boundary

In the field of tribology and braking systems analysis, the friction power method is a well-established technique for estimating the heat flux generated during braking. This method is based on the principle that the heat flux (the heat transfer rate per unit area) can be derived from the friction work generated by the sliding contact between two surfaces. Specifically, the friction work during braking on the wheel tread is caused by the interaction between the wheel and the brake shoe, where mechanical energy is converted into thermal energy through friction. This friction work, which is crucial for determining the subsequent heat flux, is calculated using Equation (6). This equation takes into account various parameters such as the normal force acting between the wheel and the brake shoe, the coefficient of friction, and the relative sliding velocity. By solving Equation (6), the friction work generated during the braking process is obtained. The heat flux associated with the wheel and the brake shoe is then computed using Equations (7) and (8), respectively. These equations represent the mathematical framework that converts the friction work into heat flux values, allowing a quantitative assessment of the thermal load on each component. After frictional work is generated, the frictional energy is distributed between the wheels and the brake shoes. The heat partitioning factors are used to characterize the heat distribution between the wheels and the brake shoes. According to Ref. [47], the heat partitioning factors are defined by Equations (9)–(11). The heat flux action areas of the wheel and brake shoe are defined by Equations (12) and (13), as described in Refs. [48,49].

$$W_{f1/2}={\eta }_{1/2}\cdot F_{\tau }\cdot v={\eta }_{1/2}\cdot \mu \cdot F_N\cdot v,$$

(6)

$$q_1=\kappa \cdot {\displaystyle \frac{W_{f_1}}{A_{\mathrm{w}}}},$$

(7)

$$q_2=\kappa \cdot {\displaystyle \frac{W_{f_2}}{A_{\mathrm{s}}}},$$

(8)

$${\beta }_{\mathrm{w}/\mathrm{s}}=\sqrt{k_{\mathrm{w}/\mathrm{s}}\cdot {\rho }_{\mathrm{w}/\mathrm{s}}\cdot c_{\mathrm{w}/\mathrm{s}}},$$

(9)

$${\eta }_1={\displaystyle \frac{{\beta }_{\mathrm{w}}\cdot A_{\mathrm{w}}}{{\beta }_{\mathrm{w}}\cdot A_{\mathrm{w}}+{\beta }_{\mathrm{s}}\cdot A_{\mathrm{s}}}},$$

(10)

$${\eta }_2={\displaystyle \frac{{\beta }_{\mathrm{s}}\cdot A_{\mathrm{s}}}{{\beta }_{\mathrm{w}}\cdot A_{\mathrm{w}}+{\beta }_{\mathrm{s}}\cdot A_{\mathrm{s}}}},$$

(11)

$$A_{\mathrm{w}}=2\cdot \pi \cdot r\cdot \mathit{\delta },$$

(12)

$$A_{\mathrm{s}}=L\cdot \mathit{\delta },$$

(13)

where W_f1/2 is the friction work in the wheel/brake shoe, F_τ/(N) is the tangential force between the wheel and brake shoe, F_N/(N) is the normal force between the wheel and the brake shoe, μ is the coefficient of friction between the wheel and brake shoe, A_w/(m²) is the wheel heat flux area, A_s/(m²) is the brake shoe heat flux area, k_w/(W/(m·°C)) and k_s/(W/(m·°C)) are the thermal conductivity of the wheel and the brake shoe, ρ_w/(kg/m³) and ρ_s/(kg/m³) are the densities of the wheel and brake shoe; c_w/(J/(kg·°C)) and c_s/(J/(kg·°C)) are the specific heat capacities of the wheel and brake shoe, L/(m) is the length of the brake shoe, δ/(m) is the effective contact width of the wheel–brake shoe, r/(m) is the radius of the wheel, κ is the contact ratio of the wheel and the brake shoe, η₁ and η₂ are the heat partitioning factors of the wheel and the brake shoe, and q₁/(W/m²) and q₂/(W/m²) are the heat flux into the wheel and the brake shoe.

3.2.2. Convection Boundary

During vehicle braking, heat convection occurs between the wheel and the surrounding air, which dissipates the heat. According to Ref. [50], the convective heat transfer coefficient is defined by Equations (14)–(17). When a vehicle is braked at a constant speed, the rotational velocity of the wheels reaches a certain constant value. This wheel speed is crucial for the dynamic behavior of the fluid flow around the wheels, especially in the context of aerodynamic and hydrodynamic considerations. To gain insight into these fluid dynamics, the Reynolds number is used. The Reynolds number is a quantitative measure of the ratio of inertial forces to viscous forces in a fluid flow. It is derived from the calculation of the wheel speed during braking, and its computation is described in Equation (14). This equation captures the complicated relationship between wheel speed, fluid properties (such as kinematic viscosity), and geometric features of the wheel, allowing for a comprehensive analysis of the fluid dynamics around the wheels during braking.

$$Re={\displaystyle \frac{2R_{\mathrm{w}}{V}_h}{v_{\mathrm{a}}}},$$

(14)

where ν_a is the kinematic viscosity of the air, m²/s; V_h is the fluid velocity at a point of the wheel, m/s; R_w is the wheel radius, m.

The Prandtl constant is a pivotal parameter in aerodynamic theory, especially in the field of compressible flow. This constant is crucial for accurately predicting the aerodynamic performance of systems operating at high speeds, where compressibility effects become significant. The Prandtl constant is shown in Equation (15).

$$Pr={\displaystyle \frac{{\rho }_{\mathrm{a}}{v}_{\mathrm{a}}{c}_{\mathrm{a}}}{{\lambda }_{\mathrm{a}}}},$$

(15)

where ρ_a is the density of the air, kg/m³; c_a is the specific heat of the air, J/(kg·°C); λ_a is the thermal conductivity of the air, W/(m²·°C).

The Nusselt number, given in Equation (16), serves as a key performance indicator for evaluating the efficiency of convective heat transfer processes. According to Ref. [51], the Nusselt number is defined, while Equation (16) illustrates the relationship among the Nusselt number, Reynolds number, and Prandtl number. This equation, derived from the equations for convective heat transfer, illustrates the complicated interplay between fluid dynamics and the phenomena of heat transfer. Based on the above equations, the convective heat transfer coefficient is shown in Equation (17).

$$Nu=0.3+{\displaystyle \frac{0.62\sqrt[2]{Re}\sqrt[3]{Pr}}{\sqrt[4]{1+{\left(\frac{0.4}{Pr}\right)}^{\frac{2}{3}}}}}{[1+{({\displaystyle \frac{Re}{282000}})}^{{\displaystyle \frac{5}{8}}}]}^{{\displaystyle \frac{4}{5}}},$$

(16)

$$h_{\mathrm{a}}={\displaystyle \frac{Nu{\lambda }_{\mathrm{a}}}{2R_{\mathrm{w}}}},$$

(17)

where h_a is the convection coefficient [50], W/(m²·°C).

3.3. Boundary Conditions of Mechanical Analysis

The mechanical load between wheel and rail is determined based on Hertzian contact theory [52] and Coulomb’s law of friction. The wheel–rail contact follows Hertzian normal contact theory, whereby the contact patch has an elliptical shape. The normal pressure distribution in the contact patch is shown in Equation (18).

$$p_n\left(x,y\right)={\displaystyle \frac{3P}{2\pi ab}}\sqrt{1-{\displaystyle \frac{x^2}{a^2}}-{\displaystyle \frac{y^2}{b^2}}},$$

(18)

In Equation (18), x and y represent the longitudinal and lateral directions of the wheel–rail contact patch; P is the wheel load, N; a and b are the long half-axis and the short half-axis of the contact patch, respectively, in mm.

According to Coulomb’s law of friction, the tangential stress between the wheel and rail is shown in Equation (19).

$$p_{\tau }\left(x,y\right)={\mu }_{\mathrm{wr}}{p}_n\left(x,y\right),$$

(19)

In Equation (19), μ_wr represents the coefficient of friction between wheel and rail.

4. Numerical Results and Discussion

4.1. Calculated Working Conditions

4.1.1. Route Condition

The route profile is shown in Figure 6. In this study, the route conditions are derived from a real long ramp with a total length of 32.855 km, consisting of four sections with an average gradient of 12‰ and a ramp curve radius of 800 m. The vehicle brakes uniformly at a constant speed of 60 km/h, resulting in a total running time on the ramp of 32.2252 min. Based on the kinematic equation, the braking time is 28 min (1680 s), with the total brake release time being 4.2252 min (253.5 s). The uniform braking process of the vehicle is divided into four segments: T1 for 300 s, T2 for 720 s, T3 for 300 s, and T4 for 360 s. The vehicle’s brake release process is divided into three sections, each with an individual brake release time of 84.5 s.

4.1.2. Wheel–Rail Contact Position

The dynamics software SIMPACK 2020 is used to determine the wheel–rail contact position in real vehicle operation on the railway line. The axle load of the vehicle is 23 t, and the friction coefficient between wheel and rail is set to 0.5. The simulation is carried out to capture the state of wheel–rail interaction during operation in a curve with a radius of 800 m. The wheel is of the HEZD type with an LM tread profile, and the rail weight is 60 kg/m. Due to the “S” shape of the curve, the wheel can be positioned on either the outside or inside rail. From the dynamic results, two types of contact positions between wheel and rail are determined. Figure 7 shows a schematic diagram of the wheel–rail contact position and thermal crack location. The location of the crack is shown in the area of Figure 7, and the length is 38 mm. With the nominal rolling circle position to the right defined as the positive direction, the following two wheel–rail contact positions result:

(1) Condition 1: When the wheel–rail contact position in the transverse direction is −17.5 mm, the longitudinal creepage is −0.16%, and the lateral creepage is −437.04 × 10⁻⁶. At this contact point, the curvature radius of the wheel profile is 100 mm, and that of the rail profile is 80 mm. According to Hertz contact theory, half the width of the contact patch in the long and short axis is 6.8712 mm and 6.6503 mm, respectively.

(2) Condition 2: When the wheel–rail contact position in the transverse direction is −4.0625 mm, the longitudinal creepage is 0.117%, and the lateral creepage is 437.04 × 10⁻⁶. At this contact point, the curvature radius of the wheel profile is 100 mm, and that of the rail is 80 mm. According to Hertz contact theory, half the width of the contact patch in the long and short axis is 6.8712 mm and 6.6503 mm, respectively.

Figure 8a shows a schematic diagram of the moving load method used to simulate wheel–rail rolling contact. As shown in Figure 8b, the thermal crack is located at the position with coordinate 0° in the system shown. The wheel–rail rolling contact load rolls from position 1 to position 2. During this rolling process, the variation of the stress intensity factor at the tip of the thermal crack on the wheel tread is investigated.

4.2. Temperature Field of Wheel Tread

In this section, the boundary working condition for the generation of the martensite phase transition temperature is calculated. The heat flux between the wheel and the brake shoe is calculated using the friction power method described in Section 3.2.1, and a simulation of the frictional heat generation between the wheel and the brake shoe is carried out based on the uniform heat source method. The specific working conditions are set as follows: the ambient temperature is 25 °C, the contact force between the wheel and brake shoe is 7800 N, the friction coefficient is 0.25, the contact width between the wheel and brake shoe is 13 mm, the contact ratio is 15.3%, and the contact position is at a distance of 43–56 mm from the back of the wheel flange. A martensitic white layer forms on the tread surface of the wheel during actual running. At this point, the wheel tread temperature is above the martensitic phase transition temperature. Figure 9 shows the maximum temperature field of the wheel tread. Figure 10 shows the temperature–time curve of the highest temperature point on the wheel tread. The simulation results of the finite element model show that the maximum temperature of the wheel tread is 739.9 °C, exceeding the martensitic transformation temperature of 700 °C for a deformed (ferritic-pearlitic) microstructure [53].

4.3. Effect of Wheel–Rail Rolling Contact Load on Thermal Crack Propagation

To simplify the finite element model, the wheel–rail rolling contact is simulated with the translation of the normal pressure and the tangential traction across the contact surface in the FE mesh. The normal contact pressure is idealized as a Hertzian distribution, and full slip rolling is assumed, meaning that the tangential force is proportional to the normal pressure according to Coulomb’s law of friction.

The propagation of thermal cracks is investigated exclusively under the influence of the wheel–rail rolling contact load. For the wheel–rail contact scenarios represented by Condition 1 and Condition 2, the changes in SIFs (K_I, K_II, K_III) with the wheel rotation angle for the three types of crack tip positions (0°, center position, 90°) shown in Figure 5 are presented in Figure 11 and Figure 12. During the entire wheel rotation, the thermal cracks remain in a closed state, meaning that the SIF K_I at the crack tip position is always zero.

As shown in Figure 11 (for contact position Condition 1), the SIF K_III is highest at the 0° position P1 (see Figure 5b), and the thermal crack propagation is dominated by Mode III anti-plane shear. When only the mechanical wheel–rail load acts, the maximum values of the SIFs K_II, K_III, and K_eff are 1.01 MPa·m^1/2, 2.60 MPa·m^1/2, and 2.40 MPa·m^1/2, respectively. At the center position P2 (see Figure 5b), the SIF K_II is the highest and the thermal crack propagation is dominated by Mode II shear. The maximum values of the SIFs K_II, K_III, and K_eff are 29.59 MPa·m^1/2, 5.43 MPa·m^1/2, and 29.81 MPa·m^1/2, respectively. At the 90° position P3 (see Figure 5b), the SIF K_III is the highest, and the thermal crack propagation is dominated by Mode III anti-plane shear. The maximum values of the SIFs K_II, K_III, and K_eff are 1.80 MPa·m^1/2, 4.92 MPa·m^1/2, and 4.16 MPa·m^1/2, respectively.

As shown in Figure 12 (for contact position Condition 2), the SIF K_III is highest at the 0° position P1, and the thermal crack propagation is dominated by the anti-plane shear of Mode III. The maximum values of K_II, K_III, and K_eff are 0.15 MPa·m^1/2, 0.35 MPa·m^1/2, and 0.33 MPa·m^1/2, respectively. At the center position P2, the SIF K_III is the highest. The maximum values of K_II, K_III, and K_eff are 0.57 MPa·m^1/2, 4.29 MPa·m^1/2, and 3.64 MPa·m^1/2, respectively. At the 90° position P3, the SIF K_II is the highest, and the thermal crack propagation is dominated by Mode II shear. The maximum values of K_II, K_III, and K_eff are 20.78 MPa·m^1/2, 16.97 MPa·m^1/2, and 23.12 MPa·m^1/2, respectively.

4.4. Effect of Wheel–Rail Rolling Contact Load Superimposed on Wheel–Brake Shoe Frictional Heat Load on Thermal Crack Propagation

When the wheel–rail rolling contact load is combined with the frictional heat effect between the wheel and brake shoe, the propagation of the thermal crack is investigated. The changes in the SIFs at the tip of the thermal crack with the wheel rotation angle are shown in Figure 13 and Figure 14. During the entire rotation of the wheel, the thermal crack remains in a closed state. After being superimposed by the frictional heat between the wheel and brake shoe, the SIF K_I at the crack tip remains constant at zero due to thermal propagation and contraction.

As shown in Figure 13 (for contact position condition 1), at the 0° position P1, the maximum values of the SIFs K_II, K_III, and K_eff are 1.21 MPa·m^1/2, 3.09 MPa·m^1/2, and 2.85 MPa·m^1/2, respectively, when the wheel–rail contact load is superimposed on the thermal wheel–brake shoe frictional load. Compared to the results caused only by the wheel–rail rolling contact load (see Figure 11), the maximum values of SIFs K_II, K_III, and K_eff increase by 19.80%, 18.85%, and 18.75%, respectively, due to the wheel–brake shoe friction heat effect, with K_II making the largest contribution. At the center position P2, the maximum values of the SIFs K_II, K_III, and K_eff are 36.64 MPa·m^1/2, 6.14 MPa·m^1/2, and 36.86 MPa·m^1/2, respectively. Compared to the results caused only by the wheel–rail rolling contact load (see Figure 11), the maximum values of the SIFs K_II, K_III, and K_eff increase by 23.83%, 13.08%, and 23.65%, respectively, due to the wheel–brake shoe friction heat effect, with K_II making the largest contribution. At the 90° position P3, the maximum values of the SIFs K_II, K_III, and K_eff are 2.48 MPa·m^1/2, 7.01 MPa·m^1/2, and 5.91 MPa·m^1/2, respectively. Compared to the results caused only by the wheel–rail rolling contact load (see Figure 11), the maximum values of the SIFs K_II, K_III, and K_eff increase by 37.78%, 42.48%, and 42.07%, respectively, due to the wheel–brake shoe friction heat effect, with K_III making the largest contribution. The wheel–brake shoe frictional heat load contributes the most to the SIF K_III.

As shown in Figure 14 (for contact position condition 2), the maximum values of the SIFs K_II, K_III, and K_eff at the 0° position P1 are 0.16 MPa·m^1/2, 0.37 MPa·m^1/2, and 0.35 MPa·m^1/2, respectively, when the wheel–rail contact load is superimposed on the wheel–brake shoe frictional thermal load. Compared to the results caused only by the wheel–rail rolling contact load (see Figure 12), the maximum values of the SIFs K_II, K_III, and K_eff increase by 6.67%, 5.71%, and 6.06%, respectively, due to the wheel–brake shoe friction heat effect, with the frictional heat load between the wheel and brake shoe making the largest contribution to the SIF K_II. At the center position P2, the maximum values of the SIFs K_II, K_III, and K_eff are 0.66 MPa·m^1/2, 4.79 MPa·m^1/2, and 4.06 MPa·m^1/2, respectively. Compared with the results caused only by the wheel–rail rolling contact load (see Figure 12), the maximum values of the SIFs K_II, K_III, and K_eff increase by 15.79%, 11.66%, and 11.54%, respectively, due to the wheel–brake shoe friction heat effect, with the wheel–brake shoe friction heat load making the largest contribution to K_II.

At the 90° position P3, the maximum values of the SIFs K_II, K_III, and K_eff are 24.53 MPa·m^1/2, 19.41 MPa·m^1/2, and 25.39 MPa·m^1/2, respectively. Taking into account the effect of the wheel–brake shoe friction heat, the maximum values of the SIFs K_II, K_III and K_eff increase by 18.05%, 14.38%, and 9.82%, respectively, due to the wheel–brake shoe friction heat effect. The frictional heat load on the brake shoes makes the largest contribution to the SIF K_II.

To facilitate the analysis of crack propagation modes, the maximum SIFs and thermal crack propagation modes for two wheel–rail contact conditions are statistically summarized in Figure 15. For both conditions, the SIF K_I remains permanently closed at the 0° position when the wheel rolls.

For the wheel–rail contact position represented by Condition 1, the SIF is higher at the center position than at the 0° and 90° positions. In this scenario, Mode II crack propagation dominates, meaning that a shear type causes rapid radial depth propagation of the crack.

For the wheel–rail contact position represented by Condition 2, the SIF at the 90° position is higher than at the 0° and center positions. In this scenario, crack propagation in Mode II dominates near the 90° position, resulting in rapid propagation in the radial depth direction.

In summary, for the wheel–rail contact position represented by Condition 1, the SIF is largest at the center position and is dominated by crack propagation in Mode II, with a maximum equivalent SIF of 36.86 MPa·m^1/2. The frictional heat of the wheel–brake shoe contributes the most to the SIF K_II at 23.83% and causes the crack to propagate rapidly in the radial depth direction. For Condition 2, the SIF at the 90° position is largest, and the Mode II crack propagation shape dominates with a maximum equivalent SIF value of 25.39 MPa·m^1/2. The frictional heat effect of the wheel–brake shoe contributes the most to the SIF K_II at 18.05% and causes the crack to propagate rapidly in the radial depth direction.

In order to investigate the influence of the discretization of the system on solution accuracy, a mesh sensitivity analysis is performed. The finite element model is discretized with three different mesh densities: coarse, medium, and fine. Particular attention is paid to areas of stress concentration, such as the crack tip and the wheel–rail contact regions, where a refined mesh is used. The convergence of the numerical solution is evaluated by gradually refining the mesh until the differences in key results, such as the stress intensity factors (SIFs) and thermal stress distributions, become negligible. The study shows that the solution converges with increasing mesh density, with significant differences observed between the coarse and medium meshes. However, the differences between the medium and fine meshes are minimal, indicating that the medium mesh is sufficient to capture the main thermomechanical behaviors. This ensures that the computational cost is optimized without compromising the accuracy of the results. The final mesh configuration is selected based on this convergence study.

5. Conclusions

The study uses a real long downhill ramp line as a case study and evaluates the effect of thermal loading on the stress intensity factor at the tip of the thermal fatigue crack in the wheel tread when the wheel tread reaches the martensitic phase transition temperature due to the braking effect. Through a systematic analysis of the differences between these two methods, the main strengths and limitations of each approach are identified, which may be relevant for future research and applications in this field.

(1) When the superposition effect of the frictional heat of the wheel and brake shoe is taken into account, the two cases calculated in this paper show that Mode II SIF K_II is larger than Mode I SIF K_I and Mode III SIF K_III. The maximum contribution of the thermal load on the wheel tread during braking in K_II is 23.83%. The propagation mode of K_II is the slip type, in which the cracks propagate rapidly in the radial depth direction. This also explains the rapid fracture of the wheel in the radial direction in practice.

(2) For all types of stress intensity factors, the maximum increase in wheel tread brake thermal load is 42.48% (Mode III SIF K_III) when the superimposed effect of wheel–brake shoe friction heat reaching the martensitic phase transition temperature is taken into account. The thermal effect has a significant influence on the propagation of thermal fatigue cracks.

(3) The location of the thermal crack initiation and the position of the wheel–rail contact loads determine the mode and direction of crack propagation.

Although the preliminary simulations are promising, as the calculated outputs are in good agreement with the broken wheel from the actual site, the model still has the following limitations.

(1) While it accurately simulates the thermal and mechanical loading conditions, the use of a moving load approach to represent the wheel–rail contact behavior may oversimplify the complexities of real-world interactions. This simplification could lead to discrepancies between the simulated and actual performance under varying operational conditions.

(2) Although the virtual crack closure method is a valuable approach for crack propagation analysis, it cannot fully capture the dynamic behavior of thermal cracks, which limits the understanding of failure mechanisms in the wheel tread under different operating scenarios.

Future research should focus on addressing these limitations by incorporating more complex models that account for wheel–rail adhesion behavior and different operating conditions. Integrating real-time data acquisition and experimental validation into the modeling process could improve the accuracy of simulations and increase the understanding of thermal and mechanical responses during braking. Exploring alternative methods to analyze crack propagation, such as dynamic fracture mechanics, could also provide deeper insights into the failure mechanisms of wheel materials. Extending the study to a wider range of operating scenarios, such as different vehicle types and track conditions, would contribute to a more comprehensive understanding of thermal crack behavior and facilitate the development of more resilient wheel designs.